Density Matrix approach in Density Functional Theory - interpretation

Question

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?

Answer 1

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$.