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

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)$$

• You're welcome! But I'd now still like to see a derivation through the characteristic function... Jan 17, 2018 at 12:49
• @FrédéricGrosshans Please have a look at my answer to see the derivation through characteristic function. Apr 27, 2018 at 13:40
• I should point out that this not a real Gaussian but a complex Gaussian en.wikipedia.org/wiki/…, so saying variance might be a little weird Mar 16 at 22:33

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}}$$

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.

• @glS I think the equation is fine. The extra factor comes from the change of variable.
• Happened to be going over this stuff again recently. All three answers were useful. This has got to be the speediest and most straightfoward. Just thought I'd note that the resolution of the identity for the coherent states here is missing a factor of $1/\pi$, which leaves the final answer missing that same factor. Oct 21, 2021 at 11:51
• The fact pi is included in my definition of the $d^2\alpha$! :)