The linear response equations for variationally determined wave functions

In order to derive the linear response equations, we start by taking the (total) perturbational derivative of the variational conditions: \[\begin{equation} \require{physics} \eval{ \dv{\eta_m} \qty( \eval{ \pdv{ E(\boldsymbol{\eta}, \vb{p}) }{p_i} }_{\vb{p}^\star}) }_{\boldsymbol{\eta}_0} = 0 \thinspace . \end{equation}\] Applying the chain rule, we subsequentially find after some shuffling: \[\begin{equation} \sum_j^x \qty( \eval{ \pdv{ E(\boldsymbol{\eta}_0, \vb{p}) }{p_i}{p_j} }_{\vb{p}^\star(\boldsymbol{\eta}_0)} ) \qty( \eval{ \pdv{ p_j^\star(\boldsymbol{\eta}) }{\eta_m} }_{\boldsymbol{\eta}_0} ) = - \eval{ \pdv{ E(\boldsymbol{\eta}, \vb{p}) }{p_i}{\eta_m} }_{ \boldsymbol{\eta}_0, \vb{p}^\star(\boldsymbol{\eta}_0) } \thinspace , \end{equation}\] which are linear equations from which we can calculate the first-order wave function response and, accordingly, these equations are called the linear response equations. They can be cast into matrix-matrix multiplication form, reminescent of Hooke’s law: \[\begin{equation} \vb{k}_{\vb{p}} \thinspace \vb{x} = - \vb{F}_{\vb{p}} \thinspace , \end{equation}\] if we identify the second order parameter derivative of the energy function with a response constant: \[\begin{equation} \qty[ \vb{k}_{\vb{p}} ]_{ij} = \eval{ \pdv{ E(\boldsymbol{\eta}_0, \vb{p}) }{p_i}{p_j} }_{\vb{p}^\star(\boldsymbol{\eta}_0)} \end{equation}\] and the mixed parameter and perturbational derivative as the response force: \[\begin{equation} \qty[ \vb{F}_{\vb{p}} ]_{im} = \eval{ \pdv{ E(\boldsymbol{\eta}, \vb{p}) }{p_i}{\eta_m} }_{ \boldsymbol{\eta}_0, \vb{p}^\star(\boldsymbol{\eta}_0) } \thinspace . \end{equation}\] The linear wave function response then is the matrix \[\begin{equation} \qty[ \vb{x} ]_{im} = \eval{ \pdv{ p_i^\star(\boldsymbol{\eta}) }{\eta_m} }_{\boldsymbol{\eta}_0} \thinspace . \end{equation}\] We should note that, using this form of the matrix equations, both the wave function response \(\vb{x}\) and the response force \(\vb{F}_{\vb{p}}\) have as many columns as there are components of the external perturbation \(\boldsymbol{\eta}\). In the end, we can find the molecular Hessian by means of matrix-matrix multiplications as: \[\begin{equation} \eval{ \frac{ \dd{^2}{\mathcal{E}(\boldsymbol{\eta})} }{\dd{\eta_m} \dd{\eta_n}} }_{\boldsymbol{\eta}_0} = \eval{ \pdv{ E(\boldsymbol{\eta}, \vb{p}^\star) }{\eta_m}{\eta_n} }_{\boldsymbol{\eta}_0} + \qty( \vb{x}^\text{T} \thinspace \vb{F}_{\vb{p}} )_{mn} \thinspace . \end{equation}\]