Operator expectation value for system of non-interacting particles (Fermions) I am reading the book "Electronic Structure" by Richard Martin which poses the following problem:
Show that the expectation value of an operator $\hat O$ in a system of identical, non-interacting Fermions (i.e. electrons with the independent particle approximation) has the following form:
$$ \left<\hat O\right> = \sum_{i,\sigma} f_i^\sigma \left<\psi_i^\sigma|\hat O|\psi_i^\sigma\right>$$
Where,
$f_i^\sigma = \frac{1}{e^{\beta(\epsilon_i^\sigma - \mu)} + 1}$ and $\psi_i^\sigma$ is the eigenvector of the ith single particle state with spin $\sigma$.
I have begun the derivation thusly (j is the jth eigenstate of the multibody system):
$$ \left<\hat O\right> = Tr\left(\hat \rho \hat O \right) = \sum_j \left<\Psi_j|\hat \rho \hat O |\Psi_j \right>$$
Since $\hat \rho$ is Hermitian:
$$ \sum_j \left<\Psi_j|\hat \rho \hat O |\Psi_j \right> = \sum_j \left<\hat \rho\Psi_j| \hat O |\Psi_j \right> $$
In the Grand Canoncial Ensemble: $\hat \rho = \frac{1}{Z}e^{-\beta(\hat H - \mu\hat N)} $
Further, considering that the jth eigenstate of the multibody system will be composed of N single particle wave functions, one can write:
$$ \sum_j \left<\hat \rho\Psi_j| \hat O |\Psi_j \right>  = \sum_{\{n_i,\sigma_i\}} \frac{1}{Z}e^{-\beta(\sum_i n_i\epsilon_i^{\sigma_i}-\mu \sum_i n_i)}\left<\Psi_j| \hat O |\Psi_j \right> $$
Based on my understanding of the other quantum mechanics reference I have been using, if $\hat O$ was simply the operator which returns the number of particles in the ith single particle state with spin $\sigma$, the sum reduces to the $f_i^\sigma$. So my guess as how to obtain the general result is to suppose that $\hat O$ is the combination of N single particle operators $\hat O = \sum_i \hat O_i$.
Hence,
$$\left<\Psi_j| \hat O |\Psi_j \right>  = \left<\Psi_j| \sum_i \hat O_i |\Psi_j \right>  = \sum_i \left<\psi_i^\sigma | \hat O_i | \psi_i^\sigma \right>$$
After this point, I am stuck and I do not see how this could reduce down to the given result. Further it would seem to me that the author should have placed a restriction on $\hat O$ to be only a single particle operator on the rhs of the target result.
 A: I think this is a situation where it is much easier to derive the desired result using second quantization. To start, let us consider a system of identical particles and let $O = \sum\limits_{k} o_k $ denote a generic one-body operator on the respective Fock space. In the language of second quantization, we can express this operator as
$$ O = \sum\limits_{ij}  \langle i|o|j\rangle \,a_i^\dagger a_j \quad . $$
The expectation value of an (one-body) operator in the state $\rho$ is defined by $\langle O \rangle_\rho \equiv \mathrm{Tr} \, \rho \, O$, where the trace is performed on the Fock space. By defining the elements of the one-body reduced density matrix in the state $\rho$ as
$$\gamma_{ij} \equiv \mathrm{Tr}\, \rho\, a_j^\dagger a_i  \quad ,$$
we see that we can write the expectation value of $O$ as
$$\langle O \rangle _\rho = \mathrm{tr}\, \gamma \,o \quad , $$
where now the trace is performed on the single-particle Hilbert space.

In the following, $i,j$ denote elements of the basis in which the single-particle Hamiltonian is diagonal.
For a system of non-interacting fermions and in equilibrium in the grand canonical ensemble, it follows that
$$ \gamma_{ij} = \delta_{ij} \,\langle n_i \rangle_\rho \quad . $$
This can be derived by e.g. applying Wick's theorem.
Here,
$$\langle n_i \rangle_\rho = \frac{1}{e^{(\epsilon_i - \mu)/{k_{\mathrm{B}}T}}+1} $$
is the well-known expression for the average occupation number of the single-particle state $i$.
Finally, this shows that indeed
$$ \langle O \rangle_\rho  = \sum\limits_i \frac{1}{e^{(\epsilon_i - \mu)/{k_{\mathrm{B}}T}}+1}\, \langle i|o|i\rangle \quad \quad ,$$
for a system of non-interacting fermions in grand canonical equilibrium.
As a last point, note that an analogous result holds for the case of non-interacting identical bosons.
A: So I believe it can also be derived using 1st quantization as follows:
$$\langle \hat O \rangle = \sum_{\left\{ n \right\}_j \in N^{\infty}} \langle \Psi_j|\hat \rho \hat O |\Psi_j\rangle = \sum_{\left\{ n \right\}_j \in N^{\infty}} \langle \hat \rho \Psi_j|\hat O |\Psi_j\rangle =\sum_{\left\{ n \right\}_j \in N^{\infty}} \tilde \rho_j \langle \Psi_j| \hat O |\Psi_j\rangle$$
Where $\tilde \rho_j$ is the eigenvalue of $\hat \rho$ associated with state j.
Presuming that $\hat O$ can be written as the sum of one body operators:
$$\hat O = \sum_k \hat o_k$$
Where k is the index of particle k. One can now write:
$$\langle \hat O \rangle = \sum_{\left\{ n \right\}_j \in N^{\infty}} \tilde \rho_j \langle \Psi_j| \sum_k \hat o_k |\Psi_j\rangle = \sum_{\left\{ n \right\}_j \in N^{\infty}} \tilde \rho_j \sum_i n_i\langle \psi_i| \hat o |\psi_i\rangle$$
Where I have used the fact that the expectation value of the single body operator $\hat o_k$ acting on the the multibody wave function $\Psi_j$ is:
$$ \langle \Psi_j | \hat o_k | \Psi_j \rangle = \sum_i n_i \frac{(N-1)!}{N!} \langle\psi_i|\hat o |\psi_i \rangle$$
Where i sums over all the single-body eigenstates present in $\Psi_j$, the $N!$ in the denominator comes from the fact that $\Psi_j$ is represented by a Slater determinant of single-body eigenstates. The $(N-1)!$ in the numerator represents the fact that each particle index is repeated $N!/N$ times. Finally, $n_i$ represents the number of times that eigenstate i is repeated.
Thus when this sum is performed N times, the factorials are eliminated.
Returning to $\langle \hat O \rangle$, one has:
$$\langle \hat O \rangle = \sum_{\left\{ n \right\}_j \in N^{\infty}} \tilde \rho_j \sum_i n_i\langle \psi_i| \hat o |\psi_i\rangle = \sum_{\left\{ n \right\}_j \in N^{\infty}}  \sum_i \tilde \rho_j n_i\langle \psi_i| \hat o |\psi_i\rangle = \sum_{i} \sum_{\left\{ n \right\}_j \in N^{\infty}} \tilde \rho_j n_i\langle \psi_i| \hat o |\psi_i\rangle $$
Which can be achieved by re-indexing the sum appropriately (I think?).
Continuing:
$$\langle \hat O \rangle = \sum_i  \left( \sum_{\left\{ n \right\}_j \in N^{\infty}} \tilde \rho_j n_i \right) \langle \psi_i| \hat o |\psi_i\rangle = \sum_i  \langle n_i \rangle \langle \psi_i| \hat o |\psi_i\rangle$$
Using the fact that:
$$ \langle n_i \rangle = \frac{1}{e^{\beta(\epsilon_i-\mu)} \pm 1} = f_i$$
One has:
$$\langle \hat O \rangle = \sum_i  f_i \langle \psi_i| \hat o |\psi_i\rangle $$
As originally desired.
