# Expectation value of $a_i^\dagger a_i$ for thermal density matrix

Suppose we have some heat bath with Hamiltonian,

$$H=\sum_n \left(a^{\dagger}_na_n+\frac{1}{2}\right)\hbar\omega_n$$

and a density matrix $\rho=Z\exp(-\beta H)$ for some normalisation $Z$. Apparently we have the relation:

$$\langle a_k^\dagger a_k\rangle=\frac{1}{\exp(-\beta\hbar\omega_k)-1}$$

I however can't quite see how to derive this. My working so far (set $k=1$):

$$\langle a_1^\dagger a_1\rangle=Z\mathrm{Tr}\left( a_1^\dagger a_1 \exp\left(-\beta(a_1^\dagger a_1+1/2)\hbar\omega_1\right)\exp\left(-\beta(a_2^\dagger a_2+1/2)\hbar\omega_2\right)\cdots\right)$$

Let us denote $S_i=\sum_n\exp\left(-\beta(n+1/2)\hbar\omega_i\right)$, then $Z^{-1}=S_1S_2S_3\cdots$:

$$=\frac{\sum_nn\exp\left(-\beta(n+1/2)\hbar\omega_1\right)}{S_1}$$ This is however as far as I can go.

• Use the Fock basis, and that should be straight-forward. And $\omega_n$ should be replaced by $\omega_k$ in the last equation. – Adam Oct 8 '14 at 6:54
• I tried but wasn't able to reach the result. If you don't mind I added my working to the question, I would appreciate it if you could see if it is valid. I am a bit uneasy writing infinite products. – Ruvi Lecamwasam Oct 8 '14 at 9:30
• You're almost there. You just need to figure out what is the value of $\sum_n n e^{-a n}$. That's easy to do by computing $\sum_n e^{-a n}$ and then take a derivative with respect to $a$. – Adam Oct 8 '14 at 9:38

Let us first consider the trace of the density matrix: $$\text{Tr}\; \exp\left\{-\beta H\right\}=\sum_{n_1,n_2,\ldots}\langle n_1,n_2,\ldots|\exp\left\{-\beta\sum_k\left(a_k^\dagger a_k+\frac{1}{2}\right)\hbar\omega_n\right\}| n_1,n_2,\ldots\rangle\\ =\prod_k \sum_{n_1,n_2,\ldots}\langle n_1,n_2,\ldots|\exp\left\{-\beta\left(a_k^\dagger a_k+\frac{1}{2}\right)\hbar\omega_n\right\}| n_1,n_2,\ldots\rangle\\ =\prod_k\left(\sum_{n_k}\langle n_k|\exp\left\{-\beta\left(a_k^\dagger a_k+\frac{1}{2}\right)\hbar\omega_k\right\}|n_k\rangle\right)$$ Crucially, the different harmonic oscillators decouple and if you compute the expectation value of some operator $\mathcal{O}$ acting only on one of the oscillators all the other oscillator contributions will cancel in $\langle \mathcal{O} \rho\rangle/\langle \rho\rangle$. Let us ignore the vacuum energy for now and compute $$\sum_{n_k}\langle n_k|\exp\left\{-\beta a_k^\dagger a_k\hbar\omega_k\right\}|n_k\rangle\\ =\sum_{n_k} \exp\left\{-\beta n_k \hbar\omega_k\right\}\\ =\sum_{n_k} \left(\exp\left\{-\beta \hbar\omega_k\right\}\right)^{n_k}\\ =\frac{1}{1-e^{-\beta \hbar\omega_k}}.$$ Adding the vacuum energy you get $$\sum_{n_k}\langle n_k|\exp\left\{-\beta \left(a_k^\dagger a_k+\frac{1}{2}\right) \hbar\omega_k\right\}|n_k\rangle\\ =\frac{e^{-\frac{1}{2}\hbar\omega_k}}{1-e^{-\beta \hbar\omega_k}}.$$ Now, we tackle the actual expectation value you want to find: $$\langle a_1^\dagger a_1\rangle=\frac{\text{Tr}\;a_1^\dagger a_1 \rho}{\text{Tr}\;\rho}.$$
As stated above, only the contribution from the first oscillator will not cancel. We need $$\sum_{n_1}\langle n_1|a_1^\dagger a_1\exp\left\{-\beta\left(a_1^\dagger a_1+\frac{1}{2}\right)\hbar\omega_1\right\}|n_1\rangle\\ =\sum_{n_1} n_1 e^{-\beta\left(n_1+\frac{1}{2}\right)\hbar\omega_1}\\ =e^{-\frac{1}{2}\beta\hbar\omega_1} \left(-\frac{1}{\hbar\omega_1}\right)\frac{\partial}{\partial\beta}\sum_{n_1} e^{-\beta n_1 \hbar\omega_1}\\ =e^{-\frac{1}{2}\beta\hbar\omega_1}\left(-\frac{1}{\hbar\omega_1}\right)\frac{\partial}{\partial\beta}\frac{1}{1-e^{-\beta \hbar\omega_1}}\\ =e^{-\frac{1}{2}\beta\hbar\omega_1}\frac{e^{-\beta \hbar\omega_1}}{\left(1-e^{-\beta \hbar\omega_1}\right)^2}$$
Now we just divide by the contribution to the trace of the density matrix from the first oscillator, which is $$\frac{e^{-\frac{1}{2}\hbar\omega_k}}{1-e^{-\beta \hbar\omega_k}},$$ and we get $$\langle a_1^\dagger a_1\rangle = \frac{e^{-\beta \hbar\omega_1}}{1-e^{-\beta \hbar\omega_1}}=\frac{1}{e^{\beta \hbar\omega_1}-1}.$$ This differs from your result by a minus sign but I think mine should be correct since your result would give negative values for the expectation value of the number operator.
The crucial point was: Both the numerator and denominator in $$\langle\mathcal{O}\rangle=\frac{\text{Tr}\;\mathcal{O}\rho}{\text{Tr}\;\rho}$$ look somewhat like $$\prod_k\left(\sum_{n_k}\langle n_k|\exp\left\{-\beta\left(a_k^\dagger a_k+\frac{1}{2}\right)\hbar\omega_k\right\}|n_k\rangle\right)=\prod_k \frac{e^{-\frac{1}{2}\hbar\omega_k}}{1-e^{-\beta \hbar\omega_k}},$$ except the numerator gets changed in the factors where $\mathcal{O}$ is acting. All the other factors cancel.