How is this expression for the energy of a free electron gas derived with the canonical ensemble (Ashcroft and Mermin)? Suppose we have some system of free electrons in a (say, 3D) box at a temperature $T$. From Ashcroft and Mermin equation (2.55), we can compute the energy of the system as
$$U=2\sum E(\mathbf{k})f(E(\mathbf{k}))$$
where the sum is over all allowed $\mathbf{k}$ states and where $f$ is the Fermi function. In particular, we are computing $U$, the average energy of the system, as the sum of the average energy of each possible electron level. Physically this makes sense, but I'm wondering how this is derived from the canonical ensemble considering the N electron system as a whole.
 A: I only sketch the answer here - for a more complete treatment I suggest consulting a statistical physics book.
The energy of a gas of free fermions is given by
$$E=\sum_\mathbf{k}\epsilon_\mathbf{k}n_\mathbf{k},$$
where $\epsilon_\mathbf{k}$ is the energy of a fermion in state $\mathbf{k}$, whereas $n_\mathbf{k}$ is the number of fermions in this state, which can be either zero or one (since these are fermions).
The partition function is then given by the sum of all possible combinations of the occupation numbers $n_\mathbf{k}$, which in the canonical ensemble would have to be restricted to summing only over the states for which $\sum_\mathbf{k}n_\mathbf{k}=N$.
$$
Z_C=\sum_{n_\mathbf{k_1}=0,1}...\sum_{n_\mathbf{k_M}=0,1} e^{-\beta
\sum_{j=1}^M\epsilon_{\mathbf{k}_j}n_{\mathbf{k}_j}}\delta_{N,\sum_j^M n_{\mathbf{k}_j}} 
$$
For practical purposes it is convenient to work with a finite number of momentum states, e.g., by imposing periodic boundary conditions, and taking the number of states to infinity only in the end of the calculation. Here $j$ enumerates the momentum states, whereas $\delta_{n,m}$ is the Kronecker symbol, restricting the sum to the states with $N$ fermions.
Calculations are a lot easier in the Grand canonical ensemble, where we need not constrain the sum:
$$
Z_{GC}=\sum_{n_\mathbf{k_1}=0,1}...\sum_{n_\mathbf{k_M}=0,1} e^{-\beta
\sum_{j=1}^M\epsilon_{\mathbf{k}_j}n_{\mathbf{k}_j} + \beta\mu N} =
\sum_{n_\mathbf{k_1}=0,1}...\sum_{n_\mathbf{k_M}=0,1} e^{-\beta
\sum_{j=1}^M\epsilon_{\mathbf{k}_j}n_{\mathbf{k}_j} + \beta\mu \sum_{j=1}^Mn_{\mathbf{k}_j}} =\\
\sum_{n_\mathbf{k_1}=0,1}...\sum_{n_\mathbf{k_M}=0,1} e^{-\beta
\sum_{j=1}^M(\epsilon_{\mathbf{k}_j}-\mu)n_{\mathbf{k}_j}}=
\prod_{j=1}^M\sum_{n_{\mathbf{k}_j}=0,1}e^{-\beta
(\epsilon_{\mathbf{k}_j}-\mu)n_{\mathbf{k}_j}}=
\prod_{j=1}^M\left[1+e^{-\beta
(\epsilon_{\mathbf{k}_j}-\mu)}\right]
$$
The energy of the system is by definition
$$
\langle E\rangle = \left\langle \sum_{j=1}^M\epsilon_{\mathbf{k}_j}n_{\mathbf{k}_j}\right\rangle =
\sum_{j=1}^M\epsilon_{\mathbf{k}_j}\langle n_{\mathbf{k}_j}\rangle,
$$
where
$$
\langle n_{\mathbf{k}_j}\rangle =Z_{GC}^{-1}
\sum_{n_\mathbf{k_1}=0,1}...\sum_{n_\mathbf{k_M}=0,1} n_{\mathbf{k}_j} e^{-\beta
\sum_{i=1}^M(\epsilon_{\mathbf{k}_i}-\mu)n_{\mathbf{k}_i}}=
Z_{GC}^{-1}\prod_{i=1,i\neq j}^M\sum_{n_{\mathbf{k}_i}=0,1}e^{-\beta
(\epsilon_{\mathbf{k}_i}-\mu)n_{\mathbf{k}_i}}\times 
\sum_{n_{\mathbf{k}_j}=0,1}n_{\mathbf{k}_j}e^{-\beta
(\epsilon_{\mathbf{k}_j}-\mu)n_{\mathbf{k}_j}}=
Z_{GC}^{-1} \prod_{i=1, i\neq j}^M\left[1+e^{-\beta
(\epsilon_{\mathbf{k}_i}-\mu)}\right]\times e^{-\beta
(\epsilon_{\mathbf{k}_j}-\mu)}=
\frac{e^{-\beta
(\epsilon_{\mathbf{k}_j}-\mu)}}{1+e^{-\beta
(\epsilon_{\mathbf{k}_j}-\mu)}}=\frac{1}{1+e^{\beta
(\epsilon_{\mathbf{k}_j}-\mu)}}=f_\mu(\epsilon_{\mathbf{k}_j})
$$
Therefore the total energy is
$$
\langle E\rangle = 
\sum_{j=1}^M\epsilon_{\mathbf{k}_j}f_\mu(\epsilon_{\mathbf{k}_j}).
$$
I would like to stress again that the result is obtained in any standard statistical physics textbook.
A: The simplest derivation of this expression from ensemble theory is to start from grand canonical ensemble to derive the Fermi distribution, and then apply this distribution as you have mentioned.
Your problem is equivalent to 'how to derive Fermi distribution from canonical ensemble'. You can refer to this Wikipedia page to get an answer.
A: Well, it doesn't follow from the canonical distribution. It follows from the grand canonical one. The two are not equivalent unless you take the thermodynamic limit $N\to\infty$. For any finite $N$ this would not be true in the canonincal ensemble.
To see how it follows in the grand canonical ensemble, see the answer by Roger Vadim, or open a statistical mechanics book.
