# 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](#variationally-determined-wave-functions):
\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}
