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

## 3 Answers

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... – Frédéric Grosshans Jan 17 '18 at 12:49
• @FrédéricGrosshans Please have a look at my answer to see the derivation through characteristic function. – Sunyam Apr 27 '18 at 13:40

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.