# 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!

## 1 Answer

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.

It looks like this page

http://functions.wolfram.com/Polynomials/HermiteH/23/02/

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