Deriving the Husimi Function of Harmonic Oscillator Eigenstates by Convolution In phase space formulation of quantum mechanics, the so-called Husimi function can be defined as the convolution of the Wigner function by an appropriate Gaussian. There are apparently alternative definitions and more elegant treatments, but I was unable to find an explicit derivation of the Husimi function of the n-th Hermite function by convolution. My ansatz reads$$ Q_n(x,p) = (W_n\ast W_0)(x,p),$$
where $Q_n$ and $W_n$ denote the Husimi and Wigner functions of the n-th eigenstate, respectively. Starting with (see p.106 of Wolfgang Schleich's Quantum Optics in Phase Space)
$$W_n(x,p) = \frac{(-1)^n}{\pi\hslash} \text{exp}
\left\lbrace - \left[ \left( \frac{p}{\hslash \kappa} \right)^2 + \left(\kappa x \right)^2 \right]\right\rbrace
L_n \left\lbrace 2 \left[ \left( \frac{p}{\hslash \kappa} \right)^2 + \left(\kappa x \right)^2 \right]\right\rbrace, \enspace \kappa = \sqrt{m\omega/\hslash},$$
one can define a radial variable $\xi^2 = 2 \left[ \left(\frac{p}{\hslash \kappa} \right)^2 + \left(\kappa x \right)^2 \right]$ and obtain
\begin{align}
(W_n \ast W_0) (\xi) &= \int d \eta \frac{(-1)^n}{\pi\hslash} \text{exp}\left(- \frac{1}{2}\eta^2 \right) L_n(\eta^2) \frac{1}{\pi \hslash}\exp\left(-\frac{1}{2}(\xi - \eta)^2\right) \\
&= \frac{(-1)^n}{(\pi\hslash)^2} \exp\left(-\frac{1}{4}\xi^2 \right) \int d\eta \exp \left[-\left(\eta -\frac{1}{2}\xi \right)^2 \right] L_n(\eta^2) .
\end{align}
Of course, the Husimi function is also known and reads (taken from a publication)
$$Q_n(\xi)=\frac{1}{2\pi\hslash n!} \left(\frac{1}{2}\xi \right)^{2n} \exp\left(-\frac{1}{4}\xi^2\right).
$$
It seems that I already got the exponential correctly. Then I can solve the integral by defining $\beta = \eta -\frac{1}{2}\xi$, applying the definition of Laguerre polynomials and then using 
$$\left(\beta+ \frac{1}{2}\xi\right)^{2k} = \sum_{j=0}^{k} \binom{2k}{2j} \left(\frac{1}{2}\xi \right)^{2k-2j} \beta^{2j},$$
instead of the standard binomial theorem since the integral vanishes if $\beta^j$ is odd. Eventually I get 
$$(W_n \ast W_0) (\xi) = \frac{(-1)^n}{(\pi\hslash)^2} \exp\left(-\frac{1}{4}\xi^2 \right)\sum_{k=0}^n \frac{(-1)^k}{k!}\binom{n}{k}\sum_{j=0}^{k} \binom{2k}{2j} \left(\frac{1}{2}\xi \right)^{2k-2j} \\
\times \int d\beta \exp (-\beta^2) \beta^{2j},$$
in which $\int d\beta \exp (-\beta^2) \beta^{2j} = \sqrt{\pi} \frac{(2j-1)!!}{2^j}.$ The difficulty is in simplifying this expression so that is equals $Q_n(\xi)$, which I haven't been able to do.  I would also appreciate if someone could come up with an alternative way of solving the integral involving the Laguerre polynomial.
 A: It might pay to avoid getting lost in the Laguerre trees, to appreciate the magnificent forest. I'll first simplify the units (whose de-nondimensionalized superfluous constants, $1/\kappa^2=\hbar=4$ , you may reinstate uniquely in the end), so that 
$$
\xi^2\equiv \frac{x^2+p^2}{2}~.
$$ 
As reviewed in our book, A Concise Treatise of Quantum Mechanics in Phase Space, eqn (122) et seq, star-multiplication with $W_0$ amounts to just a 2d Weierstrass transform, a diffusive convolution with a Gaussian, $e^{-x^2/4}/\sqrt{4\pi}$, which thus maps $e^{ax^2}$ to $\exp (\frac{ax^2}{1-4a})/\sqrt{1-4a}~$.  
Consequently,
$$
e^{2\xi^2 a} \qquad \mapsto \qquad \frac{e^{2\xi^2 a/(1-4a)}}{1-4a} ~,
$$ 
as the products of the independent x and y transforms combine through summation of their respective exponents.
The generating functional of the oscillator Wigner stargen-functions is
$$
\mathfrak{W} (\xi, s)\equiv \sum_n s^n W_n(\xi)= \sum_n s^n \frac{(-)^n}{4\pi} e^{-\xi^2/2} L_n(\xi^2) \\
=\frac{e^{-\xi^2/2}}{4\pi } \frac{1}{s+1} e^{s\xi^2/(1+s)}= \frac{1}{4\pi(1+s)} e^{\frac {\xi^2}{2}\frac{s-1}{1+s}}, 
$$
where the generating function of the Laguerres is 
$$\sum_n s^n L_n(z)= \exp(-sz/(1-s)) ~ /~(1-s) .$$
The technique is illustrated in the book linked. 
The Weierstrass transform of $\mathfrak{W}$ then is $\mathfrak{Q}$, with 
$$
a=\frac{1}{4}\frac{s-1}{1+s} ~~~\Longrightarrow 1-4a=\frac{2}{1+s} ~,
$$
so that 
$$
\mathfrak{Q}=\sum_n s^n Q_n=\frac{1+s}{4\pi (1+s) 2 }~ \exp \left ( \xi^2 \frac{1}{2}\frac{s-1}{1+s}\frac{1+s}{2} \right ) = \frac{1}{8\pi} e^{-\xi^2(1-s)/4} \\ =
\frac{e^{-\xi^2/4}}{8\pi}\sum_n s^n \frac{(\xi^2/4)^n}{n!} ~.
$$
It is then concluded that 
$$
Q_n= \frac{e^{-\xi^2/4}}{8\pi n!} (\xi^2 /a)^n ~,
$$
the standard expression you strive for, into which you may readily reinstate the re-dimensionalized units uniquely, so 4$\to \hbar$ in the normalization and the original form in $\xi$.


*

*Appreciate how this generating functional immediately shows the s dependence collapse at the origin of phase space, $\xi=0$, so that all Husimi stargen-functions beyond the Gaussian ground state must vanish at the origin---the location of the most negative part of the Wigner stargenstates has been low-pass filtered to zero. But at what a terrible price... (I assume you are familiar with the pitfalls of the Husimi.) So each $Q_n$ has  minimum at the origin and just one maximum; all of them... ugh... I assume you have experienced the nastiness of utilizing those to take matrix elements with a non-ignorable star product needed inside the phase space integrals...

