Fermi-Dirac distribution derivation? I am trying to derive the Fermi-Dirac statistics using density matrix formalism.
I know that 
$$<A>= Tr \rho A.$$
So I started from
$$<n(\epsilon_i)>= Tr \rho n(\epsilon_i)=\frac {1}{Z} \sum e^{-\beta \epsilon_i n_i}n_i=\frac {1}{Z}  e^{-\beta \epsilon_i}. $$
In the last passage I used the pauli principle ($n_i=0,1$). Now to derive the correct Fermi-Dirac distribution I have to use for $Z=1 +e^{-\beta \epsilon_i}$.
Why I have not to use the general form of 
$$Z=\prod_i (1 +e^{-\beta \epsilon_i})~?$$
Can anybody give me a good explanation?
 A: 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.

