Closed form expression for the density operator for the harmonic oscillator in thermal equilibrium I am looking for a closed form of the density operator of the quantum harmonic oscillator in thermal equilibrium, preferably in position representation. I am fairly sure it looks like a coherent state but I couldn't find it in any of my books or online sources that I skimmed. If my memory serves me well, then there  exists a closed form Wigner function for it too, so I would expect that there also exists a simple pure position representation. But please inform me if I am wrong, if there is no closed form.
A short derivation is also welcome.
EDIT
I have written down what i managed to do.
The density operator for a canonical ensemble with $\beta = \frac{1}{k_BT}$ should be $
\hat \rho = \exp(-\beta \hat H)$.Using the resolution of identity, we can write it as
$$\begin{aligned}
\hat \rho &=\exp(-\beta \hat H)\sum_{n=0}^\infty |n\rangle \langle n|\\
&= \sum_{n=0}^\infty \exp(-\beta  E_n)|n\rangle \langle n|\\
&=\sum_{n=0}^\infty \exp\left(-\beta  \left(\frac{1}{2}\hbar\omega(n+1) \right) \right)|n\rangle \langle n|\\
\end{aligned}$$
The eigenfunctions of the harmonic oscillator are
$$
\langle x | n \rangle =\phi_n(x) = \frac{1}{\sqrt{2^nn!}} \left(\frac{m\omega}{\pi\hbar}\right)^{\frac{1}{4}}\exp(-\frac{m\omega x^2}{2\hbar})H_n\left(\frac{m\omega}{\hbar} x \right)
$$
The position representation of the density operator is then
$$\begin{aligned}
\langle x|\hat \rho|x'\rangle &= \sum_{n=0}^\infty \exp\left(-\beta  \left(\frac{1}{2}\hbar\omega(n+1) \right) \right)\langle x|n\rangle \langle n|x'\rangle\\
&= \sum_{n=0}^\infty \exp\left(-\beta  \left(\frac{1}{2}\hbar\omega(n+1) \right) \right)\phi_n(x) \phi^*_n(x')\\
\end{aligned}$$
Can this sum be simplified ?
 A: If you want to write $\exp\{-\beta \hat H\}$ in closed form in the position representation, you can use Mehler's formula:
$$
\sum_{n=0}^\infty s^n \varphi_n(x)\varphi_n(y) =\\ \frac
1{\sqrt{\pi (1-s^2)}} \exp\left\{\frac{4xys
-(x^2+y^2)(1+s^2)}{2(1-s^2)}\right\}, \quad  0\le |s|<1.
$$
with
$$
s= e^{-\beta(n+1/2)}.
$$
Here
$$
\varphi_n(x)\equiv  \frac{1}{\sqrt{2^n n! \sqrt{\pi}}} H_n(x) e^{-x^2/2}
$$
is the normalized   harmonic oscillator wavefunction. I have set the frequency to be unity for convenience, but it is easy to generalise to arbitrary $\omega$.
A: For completeness sake and later reference i'll add the actual result here. The derivation requires besides the evaluation of the sum also a fair amount of hyperbolic trig function manipulation. Here it is,
$$\begin{aligned}
\rho(x,x') &\equiv \langle x | \hat \rho |x'\rangle \\[1.0em]
&= \frac{\langle x|\exp(-\beta \hat H)| x'\rangle }{Z}\\[1.0em]
&=\frac{1}{\sum^\infty_{n=0}\exp(-\beta E_n)}\sum^\infty_n \phi_n(x) \phi^*_n(x') \exp(-\beta E_n)\\[1.0em]
&= \sqrt{\frac{m\omega}{\hbar \pi} \tanh\left(\frac{1}{2}\beta \hbar \omega\right)  } \\
&\times \exp\left( -\frac{m\omega}{4\hbar } \left(   (x+x')^2\tanh\left(\frac{1}{2}\beta \hbar \omega\right) +(x-x')^2 \coth\left(\frac{1}{2}\beta \hbar \omega\right)  \right) \right)
\end{aligned}$$
Note that this $\hat \rho$ differs from the definition in my question and includes now the partition function which is missing in the original question.
A source with a fairly comprehensible derivation is found in
https://www.hep.phy.cam.ac.uk/theory/webber/tp2_06.pdf , Accessed 13.04.2021.
The final result is also given in an exercise in the book
Quantum Optics in Phase Space, Wolfgang P. Schleich, First Edition, page 64, Exercise 2.6.
