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.
You can keep using the density matrix formalism but consider switching to Fock space and the grand canonical ensemble, where Fermi Dirac statistics are exactly derivable in about two lines.
– Nanite Mar 01 '14 at 14:29