thermodynamics trace calculation I'm trying to calculate a trace to get the average energy.
The Hamiltonoperator is $H = \sum\limits_k \varepsilon_k a_k^\dagger a_k$ and $N = \sum\limits_k a_k^\dagger a_k$. The trace to calculate is
$$
E = \operatorname{Tr}(H \exp( -\beta( H - \mu N ) )
$$
I usually split the trace to a product of the elements of each Hilbert space and use the fact that
$\operatorname{Tr}(A \otimes B) = \operatorname{Tr}(A) \cdot \operatorname{Tr}(B)$.
But that doesn't work here because of the $H$ that's multiplied with the exponential function. I'm aware of the fact that it's not always possible to separate the Hilbert spaces and I think that's the case here.
Any hints on how to solve this?
 A: In the following, I assume that $a_k^\dagger, a_k$ are bosonic operators. The argument is similar in the fermionic case.
First, use that
$$
E=\operatorname{Tr}(H\exp(-\beta(H-\mu N))=\sum_k\epsilon_k\operatorname{Tr}(a_k^\dagger a_k\exp(-\beta(H-\mu N)).
$$
You can then perform the trace by tracing over the Hilbert space of each mode $k$ separately:
$$
\operatorname{Tr}(a_p^\dagger a_p\exp(-\beta(H-\mu N))\\
=-\frac{1}{\beta}\operatorname{Tr}(\frac{\partial}{\partial \epsilon_p}\exp(-\beta(H-\mu N))\\
=-\frac{1}{\beta}\frac{\partial}{\partial \epsilon_p}\prod_k \sum_{n_k=0}^{\infty}\langle n_k|\exp(-\beta(\epsilon_k a_k^\dagger a_k-\mu a_k^\dagger a_k))|n_k\rangle\\
=-\frac{1}{\beta}\frac{\partial}{\partial \epsilon_p}\prod_k \sum_{n_k=0}^{\infty}\exp(-\beta(\epsilon_k- \mu) n_k)\\
=-\frac{1}{\beta}\frac{\partial}{\partial \epsilon_p}\prod_k\frac{1}{1-e^{-\beta(\epsilon_k- \mu)}},
$$
where I used that $N=\sum_k a_k^\dagger a_k$ and $\sum_{n=0}^\infty q^n=\frac{1}{1-q}$.
Now, take the derivative with respect to $\epsilon_p$:
$$
-\frac{1}{\beta}\frac{\partial}{\partial \epsilon_p}\prod_k\frac{1}{1-e^{-\beta(\epsilon_k- \mu)}}=-\frac{1}{\beta}\prod_k\frac{1}{1-e^{-\beta(\epsilon_k- \mu)}}\frac{-\beta e^{-\beta(\epsilon_p-\mu)}}{1-e^{-\beta(\epsilon_p-\mu)}}\\
=\frac{1}{e^{\beta(\epsilon_p-\mu)}-1}\prod_k\frac{1}{1-e^{-\beta(\epsilon_k- \mu)}}
$$
Usually, one normalizes
$$
\operatorname{Tr}(\exp(-\beta(H-\mu N))=\prod_k\frac{1}{1-e^{-\beta(\epsilon_k- \mu)}}=1.
$$
Consequently,
$$
E=\sum_k\epsilon_k\operatorname{Tr}(a_k^\dagger a_k\exp(-\beta(H-\mu N))\\
=\sum_k \frac{\epsilon_k}{e^{\beta(\epsilon_k-\mu)}-1},
$$
which is exactly the sum of the energies of the modes weighted with the occupation numbers of the Bose-Einstein distribution.
