# Thermal wave function of the harmonic oscillator - proving that it's a gaussian?

I'm a bit stumped trying to prove this. I've computed the probability density for a thermal density matrix for the quantum harmonic oscillator, namely

$$\rho(x) = \frac{\sum_n^\infty e^{-\frac{\hbar\omega}{2kT}(2n+1)}\frac{1}{2^nn!}\left(\frac{m\omega}{\pi\hbar}\right)^{1/2}e^{-\frac{m\omega}{\hbar}x^2}H_n^2\left(\sqrt{\frac{m\omega}{\hbar}}x\right)}{\sum_n^\infty e^{-\frac{\hbar\omega}{2kT}(2n+1)}}$$

Now, I can compute the expectation value of $\langle x^2 \rangle$ for this distribution making use of the properties of Hermite polynomials. It turns out to be $\langle x^2 \rangle = \frac{\hbar}{m\omega}\frac{1+\xi^2}{1-\xi^2}$ with $\xi = e^{-\frac{\hbar\omega}{2kT}}$. I have the strong impression the overall function is really just a Gaussian with the corresponding variance. I tried calculating it numerically for a number of temperatures, and it always fits to a very high precision. However I can't prove it theoretically. I've tried multiple lines of attack - trying to prove that all momenta are equivalent to the normal distribution's by expanding the $x^l$ term in Hermite polynomials and making use of the triple product symbol, trying to express the polynomials as a Taylor series, Fourier transforms... the problem remains too hard to bring back to a simple analytical form. Basically the core of it is proving that:

$$e^{-z^2}\sum_n^{\infty}\xi^{2n+1}\frac{1}{2^nn!}H_n^2(z)$$

is still Gaussian in $z$, albeit with different width. Any ideas? Thanks!

According to Feynman's Statistical Mechanics, equation 2.84, it is indeed a Gaussian:

$$\rho(x)=\sqrt{\frac{m\omega}{2\pi\hbar\sinh{2f}}}\exp(-\frac{m\omega}{\hbar}x^2\tanh{f})$$

where

$$f=\frac{\hbar\omega}{2kT}$$

However, Feynman derives this by solving a differential equation, not by doing the sum you're trying to do.

of formulas for infinite summations of Hermite polynomials has what you are looking for as the 10th formula if you set $z=z_1$.
• Hmm... I don't see how the normalization of Feynman's formula is correct. The integral over $x$ doesn't give 1. When I apply the Wolfram formula to your sum I get a different normalization for $\rho(x)$ which does integrate to 1. – G. Smith Jul 2 '18 at 21:04
• I get $$\rho(x)=\sqrt{\frac{\alpha}{\pi}}\exp(-a x^2)$$ where $$\alpha=\frac{m\omega}{\hbar}\tanh{\frac{\hbar\omega}{2kT}}$$ – G. Smith Jul 2 '18 at 21:20
• By the way, I think your calculation of $\langle x^2 \rangle$ is off by a factor of 2. According to my calculation and Feynman's eqn 2.85, it is $$\langle x^2\rangle=\frac{1}{2\alpha}=\frac{\hbar}{2m\omega}\coth{\frac{\hbar\omega}{2kT}}$$ Also, Feynman's 2.89 makes clear that by his definition of $\rho(x)$ it doesn't integrate to 1. – G. Smith Jul 2 '18 at 21:42
• Thanks! It's all very useful. I can redo the calculations (I'm sure what I have for $\langle x^2 \rangle$ works since I've been doing the numerical checks so I guess I might have made a mistake copying here), but even just knowing that it's a Gaussian is very helpful to me! – Okarin Jul 2 '18 at 22:36