Calculation of hyperpolarizabilities with auxiliary density perturbation theory

A new approach for the calculation of static first hyperpolarizabilities within the auxiliary density perturbation theory is derived. It takes advantage of Wigner's 2n+1 rule. As a byproduct quadratic response equations are obtained. The final expressions are similar to those found in other approaches but present subtle differences related to the use of the auxiliary density. Validation of our new approach with respect to finite field and literature results demonstrates its numerical stability and accuracy.