Eigenvalue of the first-order reduced density matrix (1-RDM) and condensation of bosons It is defined that a system of identical bosons is said to be condensed, if  the largest eigenvalue of the first-order reduced density matrix (1-RDM) is of the order of the number of particles in the system. Check Phys. Rev. 104, 576 (1956) and chapter 2 of Many-Body Schrödinger Dynamics of Bose-Einstein Condensates.
For a given wavefunction $\psi(r_{1},...,r_{N};t)$ of $\mathcal{N}$ identical spinless bosons with spatial coordinates $r_{i}$, the first-order reduced density matrix (1-RDM) is defined as
$$
\rho^{(1)}(r_{1},r^{'}_{1};t)=N\int\psi(r_{1},r_{2},...,r_{N};t)\psi^{*}(r^{'}_{1},r_{2},...,r_{N};t)dr_{2}dr_{3}...dr_{N}=\langle\psi(t)|\hat{\Psi}^{\dagger}(r{'}_{1})\hat{\Psi}(r_{1})|\psi(t)\rangle
$$
where $\hat{\Psi}(r)$ is the bosonic field operator and the wavefunctions are properly normalized $\langle\psi(t)|\psi(t)\rangle=1$. The first-order reduced density matrix can be regarded as the kernal of the operator
$$
\hat{\rho}^{(1)}=\mathcal{N}.{Tr}_{\mathcal{N}-1}[|\psi(t)\rangle\langle\psi(t)|]
$$
What exactly are the 1-RDM and its eigenvalues and how do I get a physical picture of it?
What is the physical meaning of the above statement based on the definition of the 1-RDM?
 A: The 1st order density matrix,
$$
\rho({\bf x}, {\bf x'}) = \langle \hat\psi^\dagger({\bf x'})\hat\psi({\bf x}) \rangle
$$
is a Hermitian operator since $\rho({\bf x}, {\bf x'})  = \rho^*({\bf x'}, {\bf x})$ and its diagonal entries,
$$
\rho({\bf x}, {\bf x}) = \langle \hat\psi^\dagger({\bf x})\hat\psi({\bf x}) \rangle = n({\bf x})
$$
give the number density of particles at location ${\bf x}$, such that 
$$
\int{d{\bf x} \;\rho({\bf x}, {\bf x})} \equiv \int{d{\bf x} \;n({\bf x})} = N
$$
is the (average) total number of particles in the system. Its orthonormal eigenfunctions, defined by
$$
\int{d{\bf x'} \rho({\bf x}, {\bf x'}) \phi_j({\bf x'})} = n_j \;\phi_j({\bf x})\\
\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'})
$$
represent natural orbitals. The eigenvalues 
$$
n_j = \int{\int{d{\bf x}d{\bf x'} \;\phi^*_j({\bf x})\; \rho({\bf x}, {\bf x'})\; \phi_j({\bf x'})}}
$$
satisfy
$$
\sum_k{n_k} = \int{d{\bf x} \;n({\bf x})} = N
$$
and represent the occupation numbers of the natural orbitals. 
In fact, if $a^\dagger_j$, $a_j$ are creation and annihilation operators of these orbitals, $[a_j, a^\dagger_k]=\delta_{jk}$, the field operators are expanded as
$$
\hat\psi({\bf x}) = \sum_j{\phi_j({\bf x})\;a_j}
$$
and we have
$$
\rho({\bf x}, {\bf x'}) = \sum_{jk}{\phi^*_j({\bf x'})\phi_k({\bf x}) \langle a^\dagger_j a_k \rangle}
$$
and 
$$
\rho_{jk} \equiv \langle a^\dagger_j a_k \rangle = \int{\int{d{\bf x}d{\bf x'} \phi^*_k({\bf x}) \rho({\bf x}, {\bf x'}) \phi_j({\bf x'})}} = \rho_{jj}\delta_{jk}\\
\rho_{jj} \equiv \langle a^\dagger_j a_j \rangle = n_j
$$
In the condensate state one orbital achieves a high occupation number $n_0 \sim N$, such that
$$
N = n_0 + \sum_{j\neq 0}{n_j}
$$
and the largest eigenvalue of $\rho({\bf x}, {\bf x'})$ becomes
$$
\rho_{00} = n_0
$$
See for example these lecture notes: "Dynamics of Bose-Einstein Condensates in Trapped Atomic Gases at Finite Temperature" 
