Response density matrices
Suppose we have a normalized wave function \(\require{physics} \ket{\Psi(\vb{p}^\star)}\) that was variationally determined. If the usual energy function \[\begin{equation} \mathcal{E} ( \vb{h}, \vb{g} ) = E ( \vb{h}, \vb{g}, \vb{p}^\star(\vb{h}, \vb{g}) ) = \ev{ \hat{\mathcal{H}}(\vb{h}, \vb{g}) }{ \Psi( \vb{p}^\star(\vb{h}, \vb{g}) ) } \end{equation}\] is chosen, we define response 1- and 2-electron density matrices (T. Helgaker, Jørgensen, and Olsen 2000) as the response of the energy with respect to the 1- and 2-electron integrals. The “perturbation” we will examine is a change in one- and two-electron integrals. Since the energy is variationally determined, we can calculate the response 1-DM as: \[\begin{align} D^{\text{r}}_{PQ} &= \dv{ E(\vb{h}, \vb{g}, \vb{p}^\star(\vb{h}, \vb{g})) }{h_{PQ}} \\ &= \pdv{ E(\vb{h}, \vb{g}, \vb{p}^\star(\vb{h}, \vb{g})) }{h_{PQ}} \\ &= \ev{ \hat{E}_{PQ} }{\Psi(\vb{p}^\star(\vb{h}, \vb{g}))} = D_{PQ} \end{align}\] and analogously for the response 2-DM, we have: \[\begin{align} d^{\text{r}}_{PQRS} &= 2 \thinspace \dv{ E(\vb{h}, \vb{g}, \vb{p}^\star(\vb{h}, \vb{g})) }{g_{PQRS}} \\ &= 2 \thinspace \pdv{ E(\vb{h}, \vb{g}, \vb{p}^\star(\vb{h}, \vb{g})) }{g_{PQRS}} \\ &= \ev{ \hat{e}_{PQRS} }{\Psi(\vb{p}^\star(\vb{h}, \vb{g}))} = d_{PQRS} \thinspace , \end{align}\] which corresponds to usual definitions of the 1- and 2-DMs. Note that the extra prefactor for the 2-DM is chosen to obtain the same definition as in Helgaker (T. Helgaker, Jørgensen, and Olsen 2000). As synonyms for the response DMs, some authors use relaxed density matrices or variational density matrices.
By employing a Lagrangian formalism, we have succeeded in making non-variational wave functions variational/stationary. Therefore, we can also determine response DMs using the variational Lagrangian. For the response 1-DM, we have: \[\begin{align} D^{\text{r}}_{PQ} &= \dv{ \mathcal{L}( \vb{h}, \vb{g}, \vb{p}^\star(\vb{h}, \vb{g}), \boldsymbol{\lambda}^\star(\vb{h}, \vb{g}) ) }{h_{PQ}} \\ &= \pdv{ E( \vb{h}, \vb{g}, \vb{p}^\star(\vb{h}, \vb{g}), \boldsymbol{\lambda}^\star(\vb{h}, \vb{g}) ) }{h_{PQ}} + \sum_i^y \lambda^{\star}_i(\vb{h}, \vb{g}) \qty( \pdv{ f_i(\vb{h}, \vb{g}, \vb{p}^\star(\vb{h}, \vb{g})) }{h_{PQ}} ) \end{align}\] and for the 2-DM, we have: \[\begin{align} d^{\text{r}}_{PQRS} &= 2 \thinspace \dv{ E(\vb{h}, \vb{g}, \vb{p}^\star(\vb{h}, \vb{g})) }{g_{PQRS}} \\ &= 2 \thinspace \pdv{ E(\vb{h}, \vb{g}, \vb{p}^\star(\vb{h}, \vb{g})) }{g_{PQRS}} + \sum_i^y \lambda^{\star}_i(\vb{h}, \vb{g}) \qty( 2 \thinspace \eval{ \pdv{ % f_i(\vb{h}, \vb{g}, \vb{p}^\star(\vb{h}, \vb{g})) }{g_{PQRS}} }_{\vb{g}_0} ) \thinspace . \end{align}\] Note that when using a Lagrangian formalism to construct variationally determined wave functions, the resulting response density matrices are not necessarily \(N\)-representable.