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.

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...
• Thank you very much for the detailed answer. Lately, I've been trying to solve the integral with the generating function of the Laguerre polynomials as well, as you suggested earlier, but I was not able to arrive at the correct factor of $1/2\pi\hslash$ in the end. Otherwise it works really well though, with just a few lines of computation. I am going to compare both solutions and maybe post my solution here as well. Oct 23 '17 at 17:54
• Here, by virtue of the re-dimensionalization, $8\pi \to 2\pi \hbar$. The basic Weierstrass transform of the Gaussian is given in the WP article. Oct 23 '17 at 18:05
• One more question: The most crucial steps seems to be what you did after you wrote "consequently". When I instead use the original mapping from $\exp(ax^2)$ to $\exp(ax^2/(1-4a))/\sqrt{1-4a}$, with $a= \frac{1}{2}\frac{s-1}{1+s}$, I don't arrive at the correct result. Why do we need that second formulation of the mapping with $\exp(2\xi^2a)$, and how do you actually obtain it? You don't simply do a substitution like $a \rightarrow a^2$, but square both sides. Oct 24 '17 at 8:54
• Before "consequently" I copy the WP formula with the WP normalizations of the convolving Gaussian filter, which you may easily derive by completing Gaussian squares. I apply it independently on the x and y integrations of our ξ^2 expressions, adjusting the normalizations so as to get the proper a w.r.t. the filter gaussians, and, of course, the overall normalization sort squares since we have the 2d product of two such, the x and the y one! Oct 24 '17 at 10:38