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. 
 A: 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.
