# 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... Commented Jan 17, 2018 at 12:49
• @FrédéricGrosshans Please have a look at my answer to see the derivation through characteristic function. Commented 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 Commented Mar 16, 2022 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.
Commented Oct 24, 2020 at 6:52
• 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. Commented Oct 21, 2021 at 11:51
• The fact pi is included in my definition of the $d^2\alpha$! :)
There is in fact a rather direct method to get $$P,Q$$, and $$W$$ representations of thermal states, and more generally of displaced thermal states, by reasoning in terms of the $$T(\nu,s)$$ operators defined in the original Cahill and Glauber papers.
Definition of displaced thermal states — Define displaced thermal states as $$\rho(\alpha,x) \equiv D(\alpha)\rho(x) D(-\alpha),\\ \rho(x)\equiv (1-x)x^{a^\dagger a}\equiv(1-x)\sum_{n=0}^\infty x^n |n\rangle\!\langle n|,$$ where I'm parametrising temperatures via $$x\equiv e^{-\beta}$$ for notational convenience.
Closed formulas for quasiprobabilities — Generally speaking, the $$P$$ function of a state $$\rho$$ at a point $$\nu\in\mathbb{C}$$ can be written as $$P_\rho(\nu) = \frac1\pi \operatorname{tr}[T(\nu,1)\rho] \equiv \frac1\pi \lim_{s\to 1^-} \operatorname{tr}[T(\nu,s)\rho], \\ T(\nu,s) \equiv \int \frac{\mathrm d^2\beta}{\pi} e^{\alpha\bar\beta-\bar\alpha\beta} e^{\frac s2|\beta|^2} D(\beta) =\frac{2}{1-s} D(\nu)\left(\frac{s+1}{s-1}\right)^{a^\dagger a} D(-\nu),$$ for $$s\le1$$ (with the usual expressions with Dirac deltas arising for $$s=1$$). We are here using the expressions given in [CG1969] (see also Can we get quasiprobability distributions other than $P,Q,W$ from generalised characteristic functions?).
Main result — From $$T(\nu,s)=\rho(\nu,\frac{s+1}{s-1})$$ we get, using the properties of displacement operators and the shorthand $$t\equiv \frac{x+1}{x-1}$$, $$\operatorname{tr}[T(\nu,s) \rho(\alpha,x)] = \operatorname{tr}\left[T(\nu,s) T\left(\alpha,t\right)\right] = \int\frac{\mathrm d^2\gamma}{\pi} e^{\frac{s+t}{2}|\gamma|^2} e^{\gamma(\bar\nu-\bar\alpha)-\bar\gamma(\nu-\alpha)} \\ = \frac{-2}{s+t}\exp\left(\frac{2|\nu-\alpha|^2}{s+t}\right),$$ with the last step assuming $$\operatorname{Re}(s+t)<0$$. A similar calculation is also found in (6.39) of [CG1969]. We are interested in the case of $$\rho(\alpha,x)$$ being proper displaced thermal states, thus $$0\le x \le 1$$, i.e. $$t\le -1$$, which ensures $$s+t\le0$$ for all $$s\le 1$$.
The conclusion is now at hand. We can safely place $$s=1$$ in the derived expression, replace $$s+t=1+\frac{x+1}{x-1}=\frac{2x}{x-1}$$, and get $$P_{\rho(\alpha,x)}(\nu) \equiv \frac1\pi \lim_{s\to 1^-}\operatorname{tr}[T(\nu,s) \rho(\alpha,x)] = \frac{1-x}{\pi x} \exp\left(-\frac{1-x}{x}|\nu-\alpha|^2\right).$$ Note how for $$x\to0$$, corresponding to $$\beta\to\infty$$ or $$T\to0^+$$, we get $$\delta^2(\alpha-\nu)$$, consistently with $$\lim_{x\to0^+} \rho(\alpha,x)=|\alpha\rangle\!\langle\alpha|$$.
As a bonus, we also get for free the other quasiprobabilities. In particular, the Wigner reads $$W_{\rho(\alpha,x)}(\nu) = \frac1\pi \operatorname{tr}[T(\nu,0) \rho(\alpha,x)] = \frac{2(1-x)}{\pi(1+x)}\exp\left(-\frac{2(1-x)}{1+x}|\nu-\alpha|^2\right),$$ while the $$Q$$ reads $$Q_{\rho(\alpha,x)}(\nu) \equiv \frac1\pi \langle\nu|\rho(\alpha,x)|\nu\rangle = \frac1\pi \operatorname{tr}[T(\nu,-1) \rho(\alpha,x)] \\ = \frac1\pi(1-x)e^{-(1-x)|\alpha-\nu|^2}.$$ Note that we could have also directly gotten this formula using the general identity $$\lambda^{a^\dagger a}=N(e^{(\lambda-1)a^\dagger a})$$, discussed e.g. in Relationship between normal-ordered vacuum state and parity operator.