How to decompose displaced Hermite-Gauss function into higher order HGs? The Hermite-Gauss functions appear commonly in physics.  These functions are formed from the product of a Hermite polynomial and a Gaussian:
$$ u_n(x) = \left(\frac{2}{\pi w_0^2}\right)^{1/4} \frac{1}{\sqrt{n! 2^n}} H_n\left(\frac{\sqrt{2}x}{w_0}\right)\exp\left\{-\left(\frac{x}{w_0}\right)^2\right\}$$
and are orthonormal:
$$ \int_{-\infty}^{\infty} u_n(x) u_m(x) dx = \delta_{n,m}$$
In a paper (2004 J. Opt. B 6 495) I found the following identity, which gives the decomposition of a displaced mode $u_0(x-a)$ in terms of a series over high-order Hermite-Gauss functions $u_n(x)$:
$$ \int_{-\infty}^{\infty} u_0(x - a) u_n(x) dx 
= \frac{a^n}{w_0^n \sqrt{n!}} \exp\left\{ -\frac{a^2}{2 w_0^2}\right\} $$
How is this derived? 
(In addition to deriving it by hand, I would like to know how to coax Mathematica into giving it.)
EDIT: Here is an animation showing the decomposition of a displaced Gaussian into higher-order Hermite-Gauss functions (modes):


For example via the generating function for $u_n$'s:
we have 
$$\exp(-(x-t)^2)=\sum_n H_n(x)\exp(-x^2) t^n/n!$$
Let $U_n(x)=H_n(x)\exp(-x^2)$ (i.e. $u_n$ up to rescaling $x$ and without normalization).
We thus have
$$\exp(-(x-t)^2+x^2/2)=\sum_n U_n(x) t^n/n!$$
Then we have 
$$\int_{-\infty}^{+\infty}\exp(-(x-a)^2/2)\exp(-(x-t)^2+x^2/2)dx=\int\exp(-(x-a/2-t)^2+at-a^2/4)dx=$$
$$=\exp(-a^2/4)\exp(at)\int\exp(-x^2)dx=\sqrt{2\pi}\exp(-a^2/4)\exp(at).$$
Taking the coefficient of $t^n/n!$ we get
$$\int_{-\infty}^{+\infty}\exp(-(x-a)^2/2)U_n(x)dx=\exp(-a^2/4)a^n\sqrt{2\pi}$$
which is your formula (modulo possible integration mistake I just produced :)
