Density matrices in a restricted formalism

Given a normalized reference state \(\require{physics} \ket{\Psi}\) \[\begin{equation} \ket{\Psi} = \sum_{\vb{k}} c_{\vb{k}} \ket{\vb{k}} % \thinspace , \end{equation}\] the expectation values of an arbitrary operator \(\hat{\Omega}\) (that is compatible with a restricted formalism) can be calculated by using density matrices: \[\begin{equation} \ev{\hat{\Omega}}{\Psi} = % \sum_{pq}^K D_{pq} \Omega_{pq} % + \frac{1}{2} \sum_{pqrs}^K d_{pqrs} \Omega_{pqrs} % + \Omega_0 % \thinspace , \end{equation}\] where \(\Omega_{pq}\) and \(\Omega_{pqrs}\) are the one- and two-electron integrals in the given orthonormal orbital basis and \(D_{pq}\) and \(d_{pqrs}\) are the 1- and 2-DMs. The total molecular energy can thus be written as \[\begin{equation} \label{eq:energy_DMs_spatial} E = \ev{ \hat{\mathcal{H}} }{\Psi} % = \sum_{pq}^K h_{pq} D_{pq} % + \frac{1}{2} \sum_{pqrs}^K g_{pqrs} d_{pqrs} % + h_{\text{nuc}} % \thinspace . \end{equation}\]

The one-electron density matrix, or 1-DM, \(\vb{D}\) is given by its elements \[\begin{equation} \label{eq:1-DM} D_{pq} = \ev{\hat{E}_{pq}}{\Psi} % \thinspace , \end{equation}\] whose diagonal elements are called the orbital occupation numbers: \[\begin{equation} \omega_p = D_{pp} % \thinspace , \end{equation}\] which give the occupation numbers for a single orbital. The trace is equal to the number of electrons in the reference state \(\ket{\Psi}\): \[\begin{equation} \tr \vb{D} = \sum_{p}^K D_{pp} = N % \thinspace . \end{equation}\] Since the 1-DM is Hermitian, it can be diagonalized, yielding natural orbitals and natural orbital occupation numbers: \[\begin{equation} \vb{D} = \vb{U} \thinspace \vb{D}^\text{NO} \thinspace \vb{U}^\dagger \thinspace . \end{equation}\] Sometimes, it is fruitful to calculate the separate spin contributions for the 1-DM, which we will define as \[\begin{align} D^{\alpha \alpha}_{pq} % &= \ev{ \hat{E}^\alpha_{pq} }{\Psi} \\ &= \ev{ \hat{a}^\dagger_{p \alpha} \hat{a}_{q \alpha} }{\Psi} % \thinspace , \end{align}\] and \[\begin{align} D^{\beta \beta}_{pq} % &= \ev{ \hat{E}^\beta_{pq} }{\Psi} \\ &= \ev{ \hat{a}^\dagger_{p \beta} \hat{a}_{q \beta} }{\Psi} % \thinspace , \end{align}\] such that their sum is the previously defined density matrix: \[\begin{align} D_{pq} % &= D^{\alpha \alpha}_{pq} + D^{\beta \beta}_{pq} \\ &= \sum_\sigma D^{\sigma \sigma}_{pq} % \thinspace , \end{align}\] in which: \[\begin{align} D^{\sigma \sigma}_{pq} % &= \ev{ \hat{E}^\sigma_{pq} }{\Psi} \\ &= \ev{ % \hat{a}^\dagger_{p \sigma} \hat{a}_{q \sigma} % }{\Psi} % \thinspace . \end{align}\]

The elements of the two-electron density matrix (2-DM) are \[\begin{equation} \label{eq:2-DM} d_{pqrs} = \ev{ \hat{e}_{pqrs} }{\Psi} % \thinspace , \end{equation}\] whose diagonal elements are \[\begin{equation} \omega_{pq} = d_{ppqq} % \thinspace , \end{equation}\] which are related to simultaneous pair occupations. The trace of the 2-DM becomes two times the number of pairs in the reference state: \[\begin{equation} \tr \vb{d} % = \sum_{pq}^K d_{ppqq} % = N(N-1) % \thinspace . \end{equation}\] There exists a relation between the 1- and 2-DM: \[\begin{equation} \label{eq:1-2-DM_spatial_orbitals} D_{pq} = % \frac{1}{N-1} \sum_r^K d_{pqrr} % \thinspace . \end{equation}\] Splitting up the previously defined (total) density matrix in its spin contributions, we have: \[\begin{align} d_{pqrs} % & = d^{\alpha \alpha \alpha \alpha}_{pqrs} % + d^{\alpha \alpha \beta \beta}_{pqrs} % + d^{\beta \beta \alpha \alpha}_{pqrs} % + d^{\beta \beta \beta \beta}_{pqrs} \\ & = \sum_{\sigma \tau} d^{\sigma \sigma \tau \tau}_{pqrs} % \thinspace , \end{align}\] in which \[\begin{equation} d^{\sigma \sigma \tau \tau}_{pqrs} % = \ev{ % \hat{a}^\dagger_{p \sigma} % hat{a}^\dagger_{r \tau} % \hat{a}_{s \tau} % \hat{a}_{q \sigma} % }{\Psi} % \thinspace . \end{equation}\] In this notation, we have chosen to let the superscript spin correspond to the subscript orbital index to resemble a kind of chemist’s notation. The elements of the mixed-spin DMs are related, because: \[\begin{equation} d^{\alpha \alpha \beta \beta}_{pqrs} % = d^{\beta \beta \alpha \alpha}_{rspq} % \thinspace . \end{equation}\] We also have: \[\begin{align} & d^{\alpha \alpha \beta \beta}_{pppp} \neq 0 \\ % & d^{\alpha \alpha \beta \beta}_{pppq} % = d^{\alpha \alpha \beta \beta *}_{ppqp} \\ % & d^{\alpha \alpha \beta \beta}_{pqpp} % = d^{\alpha \alpha \beta \beta *}_{qppp} % \thinspace . \end{align}\]

By using a restricted formalism, we can show that the second-quantized electron density operator becomes \[\begin{equation} \label{eq:electron_density_operator_spatial} \hat{\rho}(\vb{r}) % = \sum_{pq}^K % \phi^*_p(\vb{r}) \phi_q(\vb{r}) % \thinspace \hat{E}_{pq} % \thinspace , \end{equation}\] whose expectation value is the electron density at the point \(\vb{r}\): \[\begin{align} \rho(\vb{r}) % &= \ev{\hat{\rho}(\vb{r})}{\Psi} \\ % &= \sum_{pq}^K % D_{pq} \thinspace % \phi^*_p(\vb{r}) \phi_q(\vb{r}) % \thinspace , \end{align}\] in which \(D_{pq}\) is the 1-DM expressed in spatial orbital basis. Using the field operators \[\begin{align} & \hat{\psi}_\sigma(\vb{r}) = % \sum_p^K % \phi_p(\vb{r}) \thinspace % \hat{a}_{p \sigma} \\ % & \hat{\psi}_\sigma^\dagger(\vb{r}) = % \sum_p^K % \phi^*_p(\vb{r}) \thinspace % \hat{a}^\dagger_{p \sigma} % \thinspace , \end{align}\] we can alternatively write: \[\begin{equation} \hat{\rho}(\vb{r}) = % \sum_\sigma % \hat{\psi}_\sigma^\dagger(\vb{r}) % \hat{\psi}_\sigma(\vb{r}) % \thinspace . \end{equation}\]

Let us, by a natural extension of the previous formula, also introduce the first-order reduced density matrix operator as \[\begin{align} \hat{\rho}(\vb{r}_1 ; \vb{r}_1') % &= \sum_\sigma % \hat{\psi}_\sigma^\dagger(\vb{r}') % \hat{\psi}_\sigma(\vb{r}) \\ &= \sum_{pq}^K % \phi^*_p(\vb{r}_1') \phi_q(\vb{r}_1) % \thinspace \hat{E}_{pq} % \label{eq:1-RDM_operator_spatial} % \thinspace . \end{align}\] The first-order reduced density matrix (1-DM) is then given by \[\begin{equation} \rho(\vb{r}_1 ; \vb{r}_1') % = \sum_{pq}^K % D_{pq} \thinspace % \phi^*_p(\vb{r}_1) \phi_q(\vb{r}_1) % \thinspace , \end{equation}\] whose diagonal is of course the electron density \(\rho(\vb{r}_1)\). Furthermore, we will define the second-order reduced density matrix operator as \[\begin{align} \hat{\rho}(\vb{r}_1, \vb{r}_2 ; \vb{r}_1', \vb{r}_2') % &= \frac{1}{2} \sum_{\sigma \tau} % \hat{\psi}^\dagger_\sigma(\vb{r}_1') % \hat{\psi}^\dagger_\tau(\vb{r}_2') % \hat{\psi}_\tau(\vb{r}_2) % \hat{\psi}_\sigma(\vb{r}_1) \\ % &= \frac{1}{2} \sum_{pqrs}^K % \phi^*_p(\vb{r}_1') % \phi_r^*(\vb{r}_2') % \phi_s(\vb{r}_2) % \phi_q(\vb{r}_1) % \thinspace \hat{e}_{pqrs} % \label{eq:2-RDM_operator_spatial} % \thinspace , \end{align}\] such that the second-order reduced density matrix (2-DM) is given by the expectation value: \[\begin{align} \rho(\vb{r}_1, \vb{r}_2 ; \vb{r}_1', \vb{r}_2') % &= \ev{ % \hat{\rho}(\vb{r}_1, \vb{r}_2 ; \vb{r}_1', \vb{r}_2') % }{\Psi} \\ &= \frac{1}{2} \sum_{pqrs}^K % d_{pqrs} % \phi^*_p(\vb{r}_1') % \phi_r^*(\vb{r}_2') % \phi_s(\vb{r}_2) % \phi_q(\vb{r}_1) % \thinspace . \end{align}\] The pair density (operator) is then the diagonal of the 2-DM (operator): \[\begin{align} & \hat{\rho}(\vb{r}_1, \vb{r}_2) % = \hat{\rho}(\vb{r}_1, \vb{r}_2 ; \vb{r}_1, \vb{r}_2) \\ % & \rho(\vb{r}_1, \vb{r}_2) % = \rho(\vb{r}_1, \vb{r}_2 ; \vb{r}_1, \vb{r}_2) % \thinspace . \end{align}\]