Grünert Polynomials I was recently reading an article on a polynomial series defined by the following relation. 
$$\left(x\dfrac{\partial}{\partial x}\right)^ne^x=\phi_n(x)e^x$$
The same article also provided the following recurrence relation. 
$$\phi_{n+1}(x)=x\sum\limits_{k=0}^n\binom{n}{k}\phi_k(x)$$
The article says this can easily be proved inductively. Thus I request not only an inductive proof, but also some intuition or a means of understanding how this equation was arrived at in the first place, since induction often proves unsatisfactory in this regard. Does the generating function, 
$$e^{x(e^z-1)}=\sum\limits_{n=0}^\infty\phi_n(x)\dfrac{z^n}{n!}$$
make this relation apparent?
 A: For any smooth function $f(x)$, we have
$$\frac{d}{dx}(x f(x)) = x\frac{d}{dx} f(x) + f(x)$$
If we view differentiation by $x$, scalar multiplication by $x$ and constant $1$,
$$\mathcal{O}_{\frac{d}{dx}} : f(x) \mapsto \frac{d}{f(x)}{dx},\quad
\mathcal{O}_x : f(x) \mapsto xf(x)\quad\text{ and }\quad
\mathcal{O}_1 : f(x) \mapsto 1\cdot f(x) = f(x)$$ 
as linear maps on space of smooth functions, we have the identity
$$\mathcal{O}_{\frac{d}{dx}}\mathcal{O}_x = \mathcal{O}_x\mathcal{O}_{\frac{d}{dx}} + \mathcal{O}_1$$
In a typical computation, we will abuse notation and use the symbols $\frac{d}{dx}$, $x$ and $1$ to denote the corresponding operators. The above identity simplifies to
$$\frac{d}{dx} x = x \frac{d}{dx} + 1$$
Apply this to your expression, we have
$$\begin{align}\phi_{n+1}(x) e^x &= \left(x \frac{d}{dx}\right)^{n+1} e^x
= x \left(\frac{d}{dx} x\right)^n \frac{d}{dx}e^x
= x \left(x\frac{d}{dx} + 1\right)^n e^x\\
&\stackrel{*}{=} x \left[\sum_{k=0}^n \binom{n}{k} \left(x\frac{d}{dx}\right)^k \right] e^x
= x \left[\sum_{k=0}^n \binom{n}{k} \phi_{k}(x)\right] e^x\end{align}$$
The crucial step $(*)$ above is nothing but applying the binomial theorem to expand the product of operators${}^{\color{blue}{[1]}}$:
$$(\mathcal{O}_{\frac{d}{dx}} \mathcal{O}_x)^n = (\mathcal{O}_x \mathcal{O}_{\frac{d}{dx}} + \mathcal{O}_1)^n
= \sum_{k=0}^n \binom{n}{k}(\mathcal{O}_x\mathcal{O}_{\frac{d}{dx}})^k\mathcal{O}_1^{n-k}
= \sum_{k=0}^n \binom{n}{k}(\mathcal{O}_x\mathcal{O}_{\frac{d}{dx}})^k
$$
Notes


*

*$\color{blue}{[1]}$ - Expansion is legal because $\mathcal{O}_1$ is the identity operator which commutes with every operators. In general, if operators $A, B$ commutes with each other (i.e. $AB = BA$), we have:
$$(A+B)^n = \sum_{k=0}^n \binom{n}{k}A^kB^{n-k}$$

