Derivation of $P$ representation of the thermal density operator I'm trying to derive the P representation for the thermal state 
$$
\rho = \sum_{n=0}^\infty \frac{\mathrm{e}^{-\beta \omega n}}{Z} |n\rangle \langle n |
$$
where $\beta$ is the inverse temperature, $Z$ is the partition function, $\omega$ is the radiation frequency, $n$ is the quantum number.
The P function is defined as 
$$
\rho = \int \mathrm{d}^2 \alpha \,\, P(\alpha)|\alpha \rangle \langle \alpha |
$$
where $|\alpha \rangle$ is a coherent state.
I've looked through Glauber's original papers that work though some more general density operators, but not specifically thermal states. He gives a formula for calculating the P function by taking the inverse Fourier transform of the characteristic function.
$$
P(\alpha) = \frac{1}{\pi^2}\int \mathrm{d}^2 \eta \,\,\mathrm{e}^{\eta \alpha^{\ast}-\eta^{\ast}a}\chi(\eta)
$$
I have found this series for the characteristic function but can't express it in terms of elementary functions. The last sum looks annoyingly close to being a binomial series.
$$
\chi(\eta)= \mathrm{Tr}\{ \rho \, \mathrm{e}^{\eta a^{\dagger}}\mathrm{e}^{-\eta^{\ast} a}\} = \sum_{m=0}^\infty \frac{\mathrm{e}^{-\beta \omega m}}{Z} \sum_{k=0}^m (-|\eta|^2)^k \frac{m!}{k!k!(m-k)!}
$$
Can anyone help me here or suggest an alternative derivation?
 A: An alternative way is to go through the Husimi Q-function.
We can easily compute 
\begin{align}
Q(α)&=\left<α|ρ|α\right>\\
  &=\frac1Z \sum_{n=0}^∞ e^{-βωn} \left<α|n\right>\!\left<n|α\right>\\
  &=\frac1Z \sum_{n=0}^∞ e^{-βωn} e^{-|α|^2}\frac{|α|^{2n}}{n!}\\
  &=\frac1Z e^{-|α|^2} \exp \left(|α|^2e^{-βω}\right)\\
  &=\frac1Z\exp\left(-[1-e^{-β ω}]|α| ^2\right),
\end{align}
The Q-function $Q(α)$ is therefore a Gaussian of variance $σ_Q^2=\frac{2}{1-e^{-βω}}$ centred in $α=0$. 
Now we can move on to the original question, computing the Glauber-Sudarshan P function. The Q-function can be computed as the $Q$ function convolved with a Gaussian of variance 2:
$$Q(β)=∫P(α)e^{-|β-α|^2} d^2α.$$
In general, inverting the convolution to compute $P$ from $Q$ can be a bit heavy, and it is rarely the best way to compute $P$, I guess. But here, $Q$ is a Gaussian of variance $σ_Q\frac{2}{1-e^{-βω}}≥2$, and it helps a lot: $P(α)$ is a Gaussian of variance $σ_P^2=σ_Q^2-2=\frac{2}{e^{βω}-1}$:
$$P(α)∝\exp
\left(-[e^{βω}-1]|α|^2\right)$$
A: A brief sketch of obtaining Glauber-Sudarshan $P$-function using normal ordered moment generating function for a signle thermalized bosonic mode.
$P[\alpha_{}^{},\alpha_{}^{*}]$ is defined as the double Fourier transform of normal ordered moment generating function ($G[\eta_{}^{},\eta_{}^{*}]$) as,
$$\mathbf{P[\alpha_{}^{},\alpha_{}^{*}]=\int_{\eta \in\mathbb{C}}^{}\frac{d_{}^{2}\eta}{\pi_{}^{2}}e_{}^{-i(\alpha_{}^{*} \eta_{}^{} + \alpha_{}^{} \eta_{}^{*})}G[\eta_{}^{},\eta_{}^{*}], \tag{1}}$$
with normal ordered moment generating function defined as,
$$\mathbf{G[\eta_{}^{},\eta_{}^{*}]=\mathbf{Tr}\left[e_{}^{i\eta_{}^{} \hat{a}_{}^{\dagger}} e_{}^{i\eta_{}^{*} \hat{a}_{}^{}} \hat{\rho}_{}^{}\right]. \tag{2}}$$
Use Baker-Campbell-Hausdorff formula twice, to put $G[\eta_{}^{},\eta_{}^{*}]$ in a convenient form as,
$$\mathbf{G[\eta_{}^{},\eta_{}^{*}]=e_{}^{\eta_{}^{*}\eta_{}^{}} \underbrace{\mathbf{Tr}\left[e_{}^{i\eta_{}^{} \hat{a}_{}^{\dagger}}  \hat{\rho}_{}^{} e_{}^{i\eta_{}^{*} \hat{a}_{}^{}}\right]}_{\text{anti-normal ordered MGF}}^{}. \tag{3}}$$
Perform trace in the coherent state representation as,
$$\mathbf{G[\eta_{}^{},\eta_{}^{*}]=e_{}^{\eta_{}^{*}\eta_{}^{}} \int_{\gamma \in \mathbb{C}}^{} \frac{d_{}^{2}\gamma}{\pi}\langle \gamma|e_{}^{i\eta_{}^{} \hat{a}_{}^{\dagger}}  \hat{\rho}_{}^{} e_{}^{i\eta_{}^{*} \hat{a}_{}^{}}|\gamma\rangle=e_{}^{\eta_{}^{*}\eta_{}^{}} \int_{\gamma \in \mathbb{C}}^{} \frac{d_{}^{2}\gamma}{\pi}e_{}^{i(\eta_{}^{*} \gamma_{}^{} +  \eta_{}^{} \gamma_{}^{*})} \langle \gamma|  \hat{\rho} |\gamma\rangle. \tag{4}}$$
For the case OP is interested in, $\hat{\rho}_{}^{}$ is the density matrix of the thermal state given by
$$\mathbf{\hat{\rho}_{}^{}=\frac{e_{}^{-\beta_{}^{}\epsilon_{}^{}\hat{a}_{}^{\dagger}\hat{a}_{}^{}}}{\mathbf{Tr}[e_{}^{-\beta_{}^{}\epsilon_{}^{}\hat{a}_{}^{\dagger}\hat{a}_{}^{}}]}.}$$
Now $\langle \gamma|  \hat{\rho}_{}^{} |\gamma\rangle$ can be evaluated 
by  deriving a differential equation satisfied by $\langle \gamma|  \hat{\rho}_{}^{}(k)_{}^{} |\gamma\rangle$ with 
$$\mathbf{\hat{\rho}_{}^{}(k)_{}^{}=\frac{e_{}^{-k\beta_{}^{}\epsilon_{}^{}\hat{a}_{}^{\dagger}\hat{a}_{}^{}}}{\mathbf{Tr}[e_{}^{-\beta_{}^{}\epsilon_{}^{}\hat{a}_{}^{\dagger}\hat{a}_{}^{}}]},}$$ as follows :
Take the derivative of $\langle \gamma|  \hat{\rho}_{}^{} |\gamma\rangle$ with respect to $k$ to get
$$\mathbf{\frac{\partial}{\partial k}\langle \gamma|  \hat{\rho}_{}^{}(k)_{}^{} |\gamma\rangle=-\beta_{}^{}\epsilon_{}^{}\langle \gamma| \hat{a}_{}^{\dagger}\hat{a}_{}^{} \hat{\rho}_{}^{}(k)_{}^{} |\gamma\rangle \underbrace{=}_{\hat{a}_{}^{}e_{}^{-k\beta_{}^{}\epsilon_{}^{}\hat{a}_{}^{\dagger}\hat{a}_{}^{}}=e_{}^{-k\beta_{}^{}\epsilon_{}^{}}e_{}^{-\beta_{}^{}\epsilon_{}^{}\hat{a}_{}^{\dagger}\hat{a}_{}^{}}\hat{a}_{}^{}}^{} -\beta_{}^{}\epsilon_{}^{}e_{}^{-k\beta_{}^{}\epsilon_{}^{}}\langle \gamma| \hat{a}_{}^{\dagger} \hat{\rho}_{}^{}(k)_{}^{} \hat{a}_{}^{} |\gamma\rangle = -\beta_{}^{}\epsilon_{}^{}e_{}^{-k\beta_{}^{}\epsilon_{}^{}} \gamma_{}^{*}\gamma_{}^{}\langle\gamma| \hat{\rho}_{}^{}(k)_{}^{} |\gamma\rangle.}$$
The solution of the resultant differential equation satisfied by $\langle \gamma|  \hat{\rho}_{}^{} |\gamma\rangle$
$$\mathbf{\frac{\partial}{\partial k}\langle \gamma|  \hat{\rho}_{}^{}(k)_{}^{} |\gamma\rangle=-\beta_{}^{}\epsilon_{}^{}e_{}^{-k\beta_{}^{}\epsilon_{}^{}} \gamma_{}^{*}\gamma_{}^{}\langle\gamma| \hat{\rho}_{}^{}(k)_{}^{} |\gamma\rangle, \tag{5}}$$
along with the obvious boundary condition $\langle \gamma|  \hat{\rho}_{}^{}(k)_{}^{} |\gamma\rangle|_{k=0}^{}=\frac{1}{\mathbf{Tr}[e_{}^{-\beta\epsilon\hat{a}_{}^{\dagger}\hat{a}_{}^{}}]}=(1-e_{}^{-\beta\epsilon})$ is given as,
$$\mathbf{\langle \gamma|  \hat{\rho}_{}^{}(k)_{}^{} |\gamma\rangle =  (1-e_{}^{-\beta\epsilon})e_{}^{(e_{}^{-k\beta_{}^{}\epsilon_{}^{}}-1)\gamma_{}^{*}\gamma_{}^{}}.\tag{6}}$$
Use $\langle \gamma|  \hat{\rho}_{}^{} |\gamma\rangle = \langle \gamma|  \hat{\rho}_{}^{}(k)_{}^{} |\gamma\rangle|_{k=1}^{}$ in Eq.$(4)$ to get 
$$\mathbf{G[\eta_{}^{},\eta_{}^{*}]=(1-e_{}^{-\beta\epsilon})e_{}^{\eta_{}^{*}\eta_{}^{}} \int_{\gamma \in \mathbb{C}}^{} \frac{d_{}^{2}\gamma}{\pi}e_{}^{i(\eta_{}^{*} \gamma_{}^{} +  \eta_{}^{} \gamma_{}^{*})}e_{}^{(e_{}^{-\beta_{}^{}\epsilon_{}^{}}-1)\gamma_{}^{*}\gamma_{}^{}}.}$$
Now perform $\gamma$ integrals (clearly the integrand is integrable for $\beta_{}^{} \epsilon_{}^{} > 0$) to get,
$$\mathbf{G[\eta_{}^{},\eta_{}^{*}]=e_{}^{-\frac{1}{e_{}^{\beta_{}^{}\epsilon_{}^{}}-1}\eta_{}^{*}\eta_{}^{}} \tag{7}.}$$
Finally use Eq.$(7)$ in Eq.$(1)$ and perform $\eta_{}^{}$ integrals (integrand is integrable for $\beta_{}^{}\epsilon_{}^{} > 0$) to get the sought after Sudarshan-Glauber distribution function aka $P$-function for single thermalized bosonic mode as
$$\mathbf{P[\alpha_{}^{},\alpha_{}^{*}]= \frac{1}{\pi}\frac{e_{}^{-\frac{\alpha_{}^{*}\alpha_{}^{}}{f_{be}^{}(\epsilon_{}^{})}}}{f_{be}^{}(\epsilon_{}^{})}, \tag{8}}$$
with the Bose-Einstein distribution function defined as,
$$\mathbf{f_{be}^{}(\epsilon_{}^{})=\frac{1}{e_{}^{\beta_{}^{}\epsilon_{}^{}}-1}. \tag{9}}$$
A: Another derivation, just using the properties of the coherent states.
Writing $\rho=\frac{e^{-\beta\omega\, a^\dagger a}}Z$, with $Z^{-1}=1-e^{-\beta\omega}$, and using the (over)-completeness of the coherent states $1=\int d^2\alpha |\alpha\rangle\langle\alpha|$, with $d^2\alpha\equiv\frac{d Re(\alpha)\,dIm(\alpha)}{\pi}$, we have
$$
\rho=\frac1Z e^{-\frac{\beta\omega}2 a^\dagger a}e^{-\frac{\beta\omega}2 a^\dagger a}
=\frac1Z\int d^2\alpha\, e^{-\frac{\beta\omega}2 a^\dagger a}|\alpha\rangle\langle\alpha|e^{-\frac{\beta\omega}2 a^\dagger a}.
$$
Using the definition of the coherent states, we have
$$
e^{-\frac{\beta\omega}2 a^\dagger a}|\alpha\rangle=e^{-\frac{|\alpha|^2}2(1-e^{-\beta\omega})}|e^{-\frac{\beta\omega}2}\alpha\rangle.
$$
Using the change of variable $\alpha'=e^{-\frac{\beta\omega}2}\alpha$ and relabeling $\alpha'$ into $\alpha$, we obtain
$$
\rho=\frac{e^{\beta\omega}}Z \int d^2\alpha\,e^{-(e^{\beta\omega}-1)|\alpha|^2}|\alpha\rangle\langle\alpha|=\int d^2\alpha\,\frac{e^{-\frac{|\alpha|^2}{n_B(\omega)}}}{n_B(\omega)}|\alpha\rangle\langle\alpha|,
$$
with $n_B(\omega)=\frac1{e^{\beta\omega}-1}$ the Bose function.
