Expression for $ \left( \frac{d^2}{dx^2} + \frac{c}{x} \frac{d}{dx} \right)^k \phi(x)$ Consider the following differential operator
\begin{align}
D_c= \frac{d^2}{dx^2} + \frac{c}{x}  \frac{d}{dx}
\end{align}
defined on $x>0$ for some given non-zero constant $c$.
Let $\phi(x)=\exp(-x^2/2)$ be the Gaussian function. For every positive integer $k$, we are interested in characterizing
\begin{align}
D^k_c \phi(x)
\end{align}
i.e., repeated application of this transform.
What I did
First, note that
\begin{align}
 \frac{d^k}{dx^k} \phi(x)= (-1)^k \operatorname{He}_k(x) \phi(x)
\end{align}
where $\operatorname{He}_k(x)$ is the Hermite polynomial.
Here is the expression for the first few expressions
\begin{align}
D_c \phi(x)&=  \operatorname{He}_2(x) \phi(x)- \frac{c}{x}  \operatorname{He}_1(x) \phi(x)\\
&= \left(\operatorname{He}_2(x) - c \right) \phi(x)
\end{align}
and (if correct)
\begin{align}
D_c^2 \phi(x)= \left( \operatorname{He}_2(x)^2-c \operatorname{He}_2(x)-4  \operatorname{He}_1^2(x)+2 + c(2-( \operatorname{He}_2(x)-c))  \right) \phi(x)
\end{align}
So the expression appears to be of the form
\begin{align}
D_c^k \phi(x)= \text{Poly}(2k) \phi(x)
\end{align}
where $\text{Poly}(2k)$ is a polynomial of degree $2k$.  But I couldn't really find the expression for the polynomial.
One more thing this operator looks like Bessel operator but is slightly different.
 A: $\newcommand{\ddx}{\frac{d}{dx}}$
You can show that $D_c^n e^{\frac{-x^2}{2}} = p_n(x)e^{\frac{-x^2}{2}}$ where $p_n$ is a polynomial of degree $2n$ with only even degree terms.
This is certainly true for $n=0$, where $p_0(x)=1$, and by doing some simple computations you can show that
$$
p_n(x) = \left(x^2-(1+c)\right)p_{n-1}(x) + \left(\frac{c}{x} - 2x\right)p_{n-1}'(x) + p_{n-1}''(x),
$$
so by induction it is true for all $n$.
I don't know of a closed form for $p_n(x)$, and I'm not sure if there is one.
Let $a_{n,k}$ be the k'th coefficient of $p_n(x) = \sum a_{n,k}x^k$.
It's not hard to show that $a_{n,2n} = 1$ and, maybe more interestingly,
$$
a_{n,2n-2} = -n(2n-1 + c).
$$
In particular, it follows from the above formula that
$$
a_{n,k} = (k+2)(k+1+c)a_{n-1,k+2} - \left( 2k + 1 + c\right)a_{n-1,k}  +   a_{n-1,k-2}.
$$
I also feel that
$$
a_{n,0} = (-1)^n\prod_{i=0}^{n-1} (c+2i+1),
$$
but I haven't tried to prove it. (EDIT: I doubt that it is true)

Let's rename the coefficients to $b_{n,i} = a_{n,2(n-i)}$, so $p_n(x) = \sum_{i=0}^n b_{n,n-i}x^{2i} = \sum_{i=0}^n b_{n,i}x^{2n-2i}$.
In other words $b_{n,i}$ is the $i$'th largest coefficient, starting from the coefficient of $x^{2n}$ and going down in steps of two. We also set $b_{n,i} = 0$ for $i < 0$ and $i > n$.
We have shown that $b_{n,0} = 1$ and that $b_{n,1} = -n(2n-1+c)$.
The inductive formula for coefficients $i \leq n$ translates to
$$
b_{n,i} = (2n-2i-2)(2n-2i+1+c)b_{n-1,i-2} - (4n-4i + 1+c) b_{n-1, i-1} + b_{n-1,i}.
$$
By solving these equations we can show taht

*

*$b_{n,0} = 1$  ,

*$b_{n,1} = -n(2n-1-c)$,

*$b_{n,2} = 2n^4 + O(n^3)$ (it gets rather complicated).

In particular, for $i=O(1)$ we can show with Faulhaber's formula that
$$
b_{n,i} = \frac{(-2)^i}{i!}n^{2i} + O(n^{2i-1}),
$$
it is true for $i=0$ and follows for the rest from the inductive formula. It doesn't follow for all $i$, since the formula is recursive and the recursion depth depends on $i$, but it works for $i=O(1)$ (and I believe even for $i=o(n)$).
You can also go forward and show, for example that
$$
b_{n,i} = \frac{(-2)^i}{i!}n^{2i} + \frac{(-2)^i}{12}\left(  \frac{4}{(i-2)!} -\frac{3(c+1)4^i i!}{(2i)!}  \right)n^{2i-1} + O(n^{2i-2}),
$$
and I guess you can go forward as much as you wont, but it will only get more complicated.
Again, I see no closed form, what you could try is expand $\left(\left(\ddx + \frac{c}{x}\right) \ddx \right)^n$, and then see how it applies to $e^{\frac{-x^2}{2}}$, but I doubt it will get any simpler.
