Calculating thermal average of an observable in quantum mechanics Single-band Hamiltonian:
If we have a Hamiltonian given as
$$
\hat H_s = \sum_k \mathcal{E}_k c_k^\dagger c_k 
$$
then the thermal average of operator $c_k^\dagger c_k$ is
$$
\langle c_k^\dagger c_k\rangle = \frac{\text{Tr}  [\hat\rho c_k^\dagger c_k]  }{\text{Tr} [\hat\rho] }
$$
here $\hat\rho=e^{-\beta \hat H_s}$ is density matrix, and $(c_k,c_k^\dagger)$ are Fermi operators. If we use bases states $|k\rangle $ then we can prove that
$$
\langle c_k^\dagger c_k\rangle = f(\mathcal{E}_k)
$$
where $f$ is Fermi distribution function.
My question is about the multi-band Hamiltonian.
Multi-band Hamiltonian: Now, we have Hamiltonian
$$
\hat H_m = \sum_k \sum_{\alpha\beta}\mathcal{E}_k^{\alpha\beta} c_{k,\alpha}^\dagger c_{k,\beta} 
$$
here $(\alpha,\beta)$ goes from $1$ to $N$, and $\mathcal{E}_k^{\alpha\beta}$ is a $N\times N$ matrix (Hamiltonian is not diagonalized). I want to calculate the thermal average of operator $c_{k,\alpha}^\dagger c_{k,\beta}$, i.e.
$$
\langle c_{k,\alpha}^\dagger c_{k,\beta} \rangle = ??
$$
 A: Let's move to the basis which diagonalizes the hamiltonian, whose fermionic operators shall be denoted by $\gamma_{k,\alpha}$ so that
\begin{equation}
H_m = \sum_{k,\alpha} \epsilon_{k,\alpha} \gamma^\dagger_{k,\alpha}\gamma_{k,\alpha}
\end{equation}
where
\begin{equation}
\gamma_{k,\alpha} = \sum_\beta u_{k,\beta} c_{k,\beta}
\end{equation}
for some unitary matrix whose elements are $u_{k,\beta}$. The problem of finding $\langle c^\dagger_{k,\alpha}c_{k,\beta} \rangle$ therefore boils down to calculating $\langle{\gamma^\dagger_{k,\alpha}\gamma_{k,\beta}}\rangle$. The best we can do is
\begin{align}
\langle \gamma^\dagger_{k,\alpha}\gamma^\dagger_{k,\beta}\rangle &= \frac{1}{Z}\text{Tr}[e^{-\beta H} \gamma^\dagger_{k,\alpha}\gamma_{k,\beta}]\\
&=\frac{1}{Z}\sum_{m,n} e^{-\beta E_m} \langle E_m|\gamma^\dagger_{k,\alpha}|E_n\rangle  \langle E_n|\gamma_{k,\beta}|E_m\rangle 
\end{align}
It may seem at first that we haven't really achieved much, a very similar result would have followed had we used a non-diagonal basis $c_{k,\alpha}$. However note that in the diagonal basis, $H_m$ conserves particle number $N_{k,\alpha}$ if we define the number operator as $N_{k,\alpha}=\gamma^\dagger_{k,\alpha} \gamma_{k,\alpha}$. This means that the eigenstates of our system have a definite particle number $|{E_m,N_m}\rangle$, and thus the calculation of the matrix elements of $\gamma$ are much easier since $\langle E_m, N_m|\gamma^\dagger_{k,\alpha}|E_n, N_n\rangle  \langle E_n, N_n|\gamma_{k,\beta}|E_m, N_m\rangle $ is only non-zero if $N_m=N_n+1$.
