Second-order response derivatives for variationally determined wave functions

In order to proceed to calculate the second-order derivatives for variationally determined wave functions, we can take the perturbational derivative of the first-order derivative and again apply the chain rule: \[\begin{align} \require{physics} \eval{ \frac{ \dd{^2}{\mathcal{E}(\boldsymbol{\eta})} }{\dd{\eta_m} \dd{\eta_n}} }_{\boldsymbol{\eta}_0} &= \eval{ \dv{\eta_m} \qty( \dv{ E(\boldsymbol{\eta}, \vb{p}^\star) }{\eta_n} ) }_{\boldsymbol{\eta}_0} \\ &= \eval{ \dv{\eta_m} \qty( \pdv{ E(\boldsymbol{\eta}, \vb{p}^\star) }{\eta_n} ) }_{\boldsymbol{\eta}_0} \\ &= \eval{ \pdv{ E (\boldsymbol{\eta}, \vb{p}^\star(\boldsymbol{\eta}_0) ) }{\eta_m}{\eta_n} }_{\boldsymbol{\eta}_0} + \sum_i^x \qty( \eval{ \pdv{ E(\boldsymbol{\eta}, \vb{p}) }{p_i}{\eta_n} }_{ \boldsymbol{\eta}_0, \vb{p}^\star(\boldsymbol{\eta}_0) } ) \qty( \eval{ \pdv{ p_i^\star(\boldsymbol{\eta}) }{\eta_m} }_{\boldsymbol{\eta}_0} ) \thinspace . \end{align}\] As can be seen, to calculate the second-order derivatives, we do have to calculate the first-order wave function response, even for variationally determined wave function parameters. In order to find the first-order wave function response, we have to solve the linear response equations.