Reduced density matrices
In a natural extension to density operators in second quantization, we introduce the first-order reduced density matrix operator as: \[\begin{align} \require{physics} \hat{\rho}(\vb{r}_1 ; \vb{r}_1') &= \hat{\Psi}^\dagger(\vb{r}_1') \hat{\Psi}(\vb{r}_1) \\ &= \sum_{PQ} \phi_P^\dagger(\vb{r}_1') \phi_Q(\vb{r}_1) \thinspace \hat{E}_{PQ} \thinspace , \end{align}\] whose expectation value is called the first-order reduced density matrix (1-DM): \[\begin{align} \rho(\vb{r}_1 ; \vb{r}_1') &= \ev{\hat{\rho}(\vb{r}_1 ; \vb{r}_1')}{\Psi} \\ &= \sum_{PQ}^M D_{PQ} \thinspace \phi^\dagger_P(\vb{r}_1') \phi_Q(\vb{r}_1) \thinspace . \end{align}\] We should note that the one-electron density matrix appears as expansion coefficients in the Cartesian product basis \(\qty{\phi_P(\vb{r}')} \times \qty{\phi_P(\vb{r})}\). We can show that we can use the 1-DM to calculate the expectation value of any one-electron operator as: \[\begin{equation} \ev{\hat{f}}{\Psi} = \int \dd{\vb{r}_1} \qty[f^c(\vb{r}_1) \rho(\vb{r}_1 ; \vb{r}_1')]_{\vb{r}_1' = \vb{r}_1} \thinspace . \end{equation}\]
We can see that the electron density (operator) is the diagonal element (i.e. \(\vb{r}' = \vb{r}\)) of the first-order reduced density matrix (operator): \[\begin{align} & \hat{\rho}(\vb{r}) = \hat{\rho}(\vb{r} ; \vb{r}) \\ & \rho(\vb{r}) = \rho(\vb{r} ; \vb{r}) \thinspace . \end{align}\]
What about a second-order reduced density matrix? If we define the corresponding operator as \[\begin{align} \hat{\rho}(\vb{r}_1, \vb{r}_2 ; \vb{r}_1', \vb{r}_2') &= \frac{1}{2} \thinspace \hat{\Psi}^\dagger(\vb{r}_1') \hat{\Psi}^\dagger(\vb{r}_2') \hat{\Psi}(\vb{r}_2) \hat{\Psi}(\vb{r}_1) \\ &= \frac{1}{2} \sum_{PQRS}^M \phi_P^\dagger(\vb{r}_1') \phi_R^\dagger(\vb{r}_2') \phi_S(\vb{r}_2) \phi_Q(\vb{r}_1) \thinspace \hat{e}_{PQRS} \thinspace , \end{align}\] we can calculate the 2-DM as the expectation value of that operator: \[\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}^M d_{PQRS} \thinspace \phi_P^\dagger(\vb{r}_1') \phi_R^\dagger(\vb{r}_2') \phi_S(\vb{r}_2) \phi_Q(\vb{r}_1) \thinspace . \end{align}\]
As in the case of the 1-DM, the 2-DM appears as some sort of expansion coefficients, and the expectation value of any two-electron operator can be calculated by means of the 2-DM: \[\begin{equation} \ev{\hat{g}}{\Psi} = \int \int \dd{\vb{r}_1} \dd{\vb{r}_2} \qty[ g^c(\vb{r}_1, \vb{r}_2) \rho(\vb{r}_1, \vb{r}_2 ; \vb{r}_1', \vb{r}_2') ]_{ \vb{r}_1' = \vb{r}_1, \vb{r}_2' = \vb{r}_2 } \thinspace . \end{equation}\] The diagonal element of the 2-DM (operator) (i.e. \(\vb{r}_1' = \vb{r}_1\) and \(\vb{r}_2' = \vb{r}_2\)) is equal to the electron pair density (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}\]