Fourier Transform of Hermite polynomial $H_n(x)e^{-\frac{x^2}{2}}$ I'm trying to solve the 2D paraxial equation $2i\partial_zu=-\partial_x^2u$, for the initial condition $u(x,z=0)=H_n(x)e^{-x^2/2}$, with $x$ and $z$ both real and $n\geq0$. 
For $n=0$, I used the Fourier transform—defined as  $\mathscr{F}_x\big\{f(x)\big\}(k)=\int_{-\infty}^{+\infty}dxf(x)e^{-ikx}$ — to get 
$$\tilde{u}(k,z)=\sqrt{2\pi}e^{-\frac{k^2}{2}(1+iz)}$$
$${u}(x,z)=\dfrac{1}{\sqrt{1+iz}}e^{-\frac{x^2}{2(1+iz)}}$$
in $k$ and $x$ domains respectively, i.e., I correctly get the propagation of a normalized Gaussian beam.
For the general case of $n\geq1$


*

*How can I prove that $\mathscr{F}\big\{H_n(x)e^{-x^2/2}\big\}=\sqrt{2\pi}(-i)^nH_n(k)e^{-k^2/2}$?

*What would be the inverse Fourier transform $\mathscr{F}^{-1}\big\{\sqrt{2\pi}(-i)^nH_n(k)e^{-\frac{k^2}{2}(1+iz)}\big\}$.
 A: We can prove that:
$$I(k,a)=\int_{-\infty}^{\infty}H_n(ax)e^{-\frac{x^2}{2}}e^{-ikx}dx=\sqrt{2\pi}(-i)^n\Big(2a^2-1\Big)^{n/2}H_n{\Big(\frac{ak}{\sqrt{2a^2-1}}\Big)}e^{-\frac{k^2}{2}}$$
This shows that when $a=1$ we recover the desired result, that the Hermite polynomials are eigenfunctions of the Fourier transform. This also allows to compute the second Fourier transform in question since:
$$\int_{-\infty}^{\infty}H_n(k)e^{-\frac{k^2}{2b^2}}e^{ikx}\frac{dk}{2\pi}=\int_{-\infty}^{\infty}H_n(x)e^{-\frac{x^2}{2b^2}}e^{ikx}\frac{dx}{2\pi}=\frac{b}{2\pi}\int_{-\infty}^{\infty}H_n(b x)e^{-\frac{x^2}{2}}e^{ikb x}dx=(-1)^n\frac{b}{2\pi}I(kb,b)$$
and therefore we obtain the slightly more general result, valid for complex $b$ in general for which the integral converges:
$$I(x,a,b)=\int_{-\infty}^{\infty}H_n(ak)e^{-\frac{k^2}{2b^2}}e^{ikx}dk=i^nb\sqrt{2\pi}(2a^2b^2-1)^{n/2}H_n\Big(\frac{a~b^2 ~x}{\sqrt{2a^2b^2-1}}\Big)e^{-b^2x^2/2}$$
from which we finally obtain that:
$$\mathcal{F^{-1}}(\sqrt{2\pi}(-i)^nH_n(k)e^{-(1+iz)k^2/2})=\frac{(-1)^n}{\sqrt{1+iz}}\Big(\frac{1-iz}{1+iz}\Big)^{n/2}H_n\Big(\frac{x}{\sqrt{1+z^2}}\Big)e^{-\frac{x^2}{2(1+iz)}}$$
Note for calculation of $I$ to be added soon.
$\textbf{EDIT:} ~~\small\text{Calculation of $I(k,a,b)$ valid for all $b\in\mathbb{C}$}$
First, write $H_n(x)=\sum_{l}c_{nl}x^l$. Substitute this in and perform the integrals:
$$I(k,a,b)=\sum_{l}c_{nl}a^l\int_{-\infty}^{\infty}x^le^{-x^2/2b^2}e^{-ikx}dx\\=b\sqrt{2\pi}\sum_{l}c_{nl}a^l\Big(i\frac{d}{dk}\Big)^le^{-k^2b^2/2}\\=b\sqrt{2\pi}\sum_{l}c_{nl}(iab)^le^{-k^2b^2/2}\Big[e^{k^2b^2/2}\Big(\frac{d}{d(kb)}\Big)^le^{-k^2b^2/2}\Big]\\=b\sqrt{2\pi}\sum_{l}c_{nl}(-iab)^le^{-k^2b^2/2}He_n(kb)$$
where $He_n(x)$ are the probabilists Hermite's polynomials as defined on the wikipedia page.
Now utilize the represenations in terms of differential operators$$He_n(ax)=a^ne^{-D^2/2a^2}x^n, H_n(ax)=(2a)^ne^{-D^2/4a^2}x^n, D\equiv\frac{d}{dx}$$  repeatedly  to rewrite again as:
$$\begin{align}I(k,a,b)&=b\sqrt{2\pi}e^{-k^2b^2/2}e^{-D^2/2b^2}\sum_{l}c_{nl}(-iab^2k)^l\\&=b\sqrt{2\pi}e^{-k^2b^2/2}e^{-D^2/2b^2}H_n(-iab^2k)\\&=b\sqrt{2\pi}e^{-k^2b^2/2}e^{-D^2/2b^2}(2iab^2)^n e^{D^2/4a^2b^4}k^n\\&=(2iab^2)^n e^{-k^2b^2/2}b\sqrt{2\pi}e^{-k^2b^2/2}e^{-D^2(1/2b^2-1/4a^2b^4)}k^n\\&=b\sqrt{2\pi}e^{-k^2b^2/2}\Big(-iab\sqrt{2-\frac{1}{a^2b^2}}\Big)^nH_n\Big(\frac{kb}{\sqrt{2-\frac{1}{a^2b^2}}}\Big)\end{align}$$
which upon some basic algebra yields the quoted result:
$$I(k,a,b)=b\sqrt{2\pi}(-i)^n\Big(2a^2b^2-1\Big)^{n/2}H_n\Big(\frac{kab^2}{\sqrt{2a^2b^2-1}}\Big)e^{-k^2b^2/2}$$
A: Perhaps this is not what you wanted, but it is clearest to me: some details about the quantum harmonic oscillator include that all eigenvectors are generated from $u_1(x)=e^{-x^2/2}$ by the "raising operator" $R=i{\partial\over \partial x}ix$, and that the $n+1$ eigenspace is exactly scalar multiples of $R^n u_1$. Since $u_1$ is annihilated by the lowering operator $L=i{\partial\over \partial x}+ix$, the interaction of Fourier transform with such operators shows that $\widehat{u_1}=u_1$ (at least with suitable normalization of Fourier transform.
Then an induction on $n$ shows that Fourier transform on $R^n u_1$ is multiplication by $i^{-n}$.
