Condensate fraction and single-particle density matrix In Bose–Einstein condensation (BEC), how to prove the largest eigenvalue
of the single-particle density matrix $$\rho_{ij}=\frac{\langle\Psi|a_i^{\dagger}a_j|\Psi\rangle}{N}$$ is $$\frac{1}{N}\sum_{i}{|\langle \Psi|a_i|\Psi\rangle|^2},$$ which is the condensate fraction f of the total particle number N in the thermodynamic limit?
 A: If operators $a_j$, $a^\dagger_j$ correspond to the system's orthonormal natural orbitals $\phi_j({\bf x})$, such that
$$ 
\int{d{\bf x}\; \phi^*_j({\bf x})\phi_k({\bf x})} = \delta_{jk}, \;\;\;\sum_j{\phi^*_j({\bf x})\phi_j({\bf x'})} = \delta({\bf x} - {\bf x'})\\
\hat\psi({\bf x}) = \sum_j{\phi_j({\bf x})\;a_j},\;\;\; a_j = \int{d{\bf x} \;\phi^*_j({\bf x})\;\hat\psi({\bf x})}
$$
then the 1st order density matrix is already diagonal (see this related answer),
$$
\rho_{ij} = \rho_{jj}\delta_{ij}
$$
and, with the given normalization, the eigenvalues are simply the occupation fractions of orbitals $\phi_j$,
$$
\rho_{jj} \equiv \frac{n_j}{N} = \frac{\langle a^\dagger_j a_j \rangle}{N}, \;\;\;
\sum_j{\frac{n_j}{N}} = \sum_j{\frac{\langle a^\dagger_j a_j \rangle}{N}} = 1
$$
The condensate phase is defined as the macroscopic occupation of a dominant orbital $\phi_0$, such that 
$$
N \rho_{00} = n_0 \sim N \;\; >> N \rho_{jj} = n_j, \;\; j\neq 0
$$ 
Formally this amounts to a macroscopic component of the field operators, such that
$$
\hat\psi({\bf x}) = \sqrt{n_0}\phi_0({\bf x}) + \bar \psi({\bf x})\\
\langle \bar\psi({\bf x}) \rangle = 0, \;\;\;\langle \hat\psi({\bf x}) \rangle = \sqrt{n_0}\phi_0({\bf x}) 
$$
Equivalently, for orbital operators we have
$$
a_j = \int{d{\bf x} \phi^*_j({\bf x})\hat\psi({\bf x})} = \int{d{\bf x} \phi^*_j \left( \sqrt{n_0}\phi_0({\bf x}) + \bar\psi({\bf x}) \right)}
= \sqrt{n_0} \delta_{0j} + \bar a_j\\
\langle \bar a_j \rangle = 0, \;\;\; \langle a_j \rangle = \sqrt{n_0} \delta_{0j},\;\;\; \langle a^\dagger_j a_j \rangle = n_j
$$
and so
$$
\frac{1}{N} \sum_j{|\langle a_j \rangle |^2} = \frac{1}{N} \sum_j{| \sqrt{n_0} \delta_{0j}  |^2} = \frac{n_0}{N} = \frac{\langle a^\dagger_0 a_0 \rangle}{N} \equiv \rho_{00}
$$
