Derivation of $P$ representation of the thermal density operator

Question

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?

Answer 1

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

Answer 2

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

Answer 3

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.

Answer 4

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.