The derivation of the Fermi-Dirac distribution using the density matrix formalism proceeds as follows:
The setup.
We assume that the single-particle hamiltonian has a discrete spectrum, so the single-particle energy eigenstates are labeled by an index $i$ which runs over some finite or countably infinite index set $I$. A basis for the Hilbert space of the system is the occupation number basis
\begin{align}
|\mathbf n\rangle = |n_0, n_1, \dots\rangle
\end{align}
where $n_i$ denotes the number of particles occupying the single-particle energy eigenstate $i$. For a system of non-interacting identical fermions, the set $\mathscr N_-$ of admissible occupation sequences $\mathbf n$ consists of those sequences with each $n_i$ equal to either $0$ or $1$. Let $H$ be the hamiltonian for such a system, and let $N$ be the number operator, then we have
\begin{align}
H|\mathbf n\rangle = \left(\sum_{i\in I}n_i\epsilon_i\right)|\mathbf n\rangle, \qquad N|\mathbf n\rangle = \left(\sum_{i\in I} n_i\right) |\mathbf n\rangle
\end{align}
where $\epsilon_i$ is the energy of eigenstate $i$. We can also define an observable $N_i$ which tells us the occupation number of the $i^\mathrm{th}$ single-particle energy state;
\begin{align}
N_i|\mathbf n\rangle = n_i|\mathbf n\rangle
\end{align}
Note that we are attempting to determine the ensemble average occupation number of the $j^\mathrm{th}$ energy eigenstate. In the density matrix formalism, this is given by
\begin{align}
\langle n_j\rangle =\mathrm{tr}(\rho N_i)
\end{align}
where
\begin{align}
\rho = \frac{e^{-\beta(H-\mu N)}}{Z}, \qquad Z = \mathrm {tr}\big(e^{-\beta(H-\mu N)}\big)
\end{align}
The proof.
- Show that
\begin{align}
Z = \sum_{\mathbf n\in \mathscr N_-}\prod_{i\in I}x_i^{n_i}
\end{align}
where $x_j = e^{-\beta(\epsilon_j-\mu)}$, the sum is over admissible sequences $\mathbf n$ of occupation numbers of single-particle energy states, and the product is over indices $i$ labeling an orthonormal basis of single particle energy eigenstates.
- Show that the ensemble average occupation number of the $j^\mathrm{th}$ state can be computed as follows:
\begin{align}
\langle n_j\rangle = x_j\frac{\partial}{\partial x_j}\ln Z
\end{align}
- Show that the product and the sum in the partition function can be "exchanged" to give
\begin{align}
Z = \prod_{i\in I}\sum_{n=0}^1 x_i^n
\end{align}
where the product is now over single-particle energy eigenstates, and the sum is over admissible occupation numbers of a single-particle state.
- Combine the results of steps 2 and 3 to show that
\begin{align}
\langle n_j\rangle = \frac{1}{e^{\beta(\epsilon_j-\mu)}+1}
\end{align}
which is the desired result.
