# Density Matrix approach in Density Functional Theory - interpretation

In a paper describing a Kohn-Sham Density Functional Theory implementation, the authors describe the use of the density matrix for e.g. the calculation of the electronic density and for efficiency improvements. With a basis expansion of the Kohn-Sham states $$|\Psi_n^\mathrm{KS}\rangle = \sum_i c_{in}\,|\phi_i\rangle$$ the density operator and matrix are \begin{align*} \hat{\rho} &= \sum_n f_n |\Psi_n\rangle\langle\Psi_n| =\sum_{ij} \rho_{ij}\, |\phi_i\rangle\langle\phi_j| \\ \rho_{ij} &= \sum_n f_n\, c^*_{in}\,c_{jn} \end{align*} with the occupation numbers $f_n$ - the electron density is then obtained via $$\langle\mathbf{r}|\hat{\rho}|\mathbf{r}\rangle = \rho = \sum_{ij} \phi_i\,\rho_{ij}\,\phi_j.$$ When comparing this with the standard definition of the density operator $$\hat{\rho} = \sum_n w_n |\Psi_n\rangle\langle\Psi_n|,$$ the weights $w_n$ have a different meaning, being the probabilities to find a subsystem in the state $|\Psi_n\rangle$ and fulfilling $\sum_n w_n = 1$. In my understanding, this has nothing to do with the density operator in above context where $f_n$ denote whether a single-particle orbital is occupied or not and where $\sum_n w_n \geq 1$ is possible.

Now to my question: is the density matrix in the context of DFT as explained above just a tool to ease technical/computational operations or can still any meaning in the 'classical' quantum-mechanical view be given to it?

While the Kohn-Sham states $$|n^{\rm KS}\rangle$$ are really only an auxiliary construct to compute the non-interacting kinetic energy, one often nevertheless goes ahead and interprets the $$|n^{\rm KS}\rangle$$ as single-particle states and the Kohn-Sham eigenvalues $$\epsilon_n$$ as quasi-particle energies (e.g. one uses the Kohn-Sham band structure as a first approximation to an experimentally observed quasi-particle band structure with such famous problems of density functional theory typically underestimating band gaps significantly).
If one does accept the interpretation of Kohn-Sham states as single-/quasi-particle states, the eigenvalues $$\lambda_n$$ of the density matrix $$\hat \rho = \sum_n \lambda_n |n^{\rm KS}\rangle\langle n^{\rm KS}|$$ are the probabilities to find a particle in state $$|n^{\rm KS}\rangle$$ (or simply the occupation of that state). $$\hat \rho$$ is therefore called one-particle reduced density matrix: $$N-1$$ particle degrees of freedom have been traced out. The eigenvectors of one-particle reduced density matrices are called natural orbitals (which in the case of Kohn-Sham density functional theory coincide with the Kohn-Sham states).
Since Kohn-Sham density functional theory only uses a single Slater determinant, at $$T=0$$K the $$\lambda_n$$ are either $$0$$ or $$1$$ (or $$0$$ or $$2$$ if spin degeneracy is factored into the occupation numbers); correlation beyond a single Slater determinant would lead to fractional occupation even at $$0$$K.
Besides this single-particle like interpretation of the above density matrix, a potential computational advantage comes about in the case of actually not decomposing $$\hat \rho$$ into its representation via the Kohn-Sham states, but expressing it in real space as $$\rho(\vec r, \vec r^\prime)$$. One can then take advantage of the fact that $$\rho(\vec r, \vec r^\prime)$$ typically decays 'relatively' quickly with $$|\vec r - \vec r^\prime|$$ and attempt to simply truncate $$\rho(\vec r, \vec r^\prime)$$ beyond a cutoff radius. This truncation can be used to construct linear-scaling approaches, avoiding the computational disadvantage of contructing orthogonal Kohn-Sham states which accounts for the general cubic scaling with system size of Kohn-Sham density functional theory (there are also basis set choice related techniques based on Kohn-Sham orbitals that allow for linear scaling computational complexity such as atom-centered atomic-orbital-like basis sets, where orbital overlap is truncated beyond a cutoff).
Finally a note on the trace of $$\hat \rho$$: unlike the trace of the pure density matrix of a single particle, which is $$1$$, the trace of the one-particle reduced density matrix $$\hat \rho$$ of an $$N$$-particle system is $$N$$; each particle is found with 100% probability somewhere, and there are $$N$$ particles. Hence this is not in contradiction to the eigenvalues of $$\rho$$ coinciding with single-particle occupation numbers: the sum of the single-particle occupation numbers must indeed be $$N$$.