A General Formalism for Electromagnetic Response Based on Quasi-Energy Derivatives Within a Single-Determinant Framework

1. Introduction

Since the early development of quantum chemistry, one of its central goals has been the accurate prediction of molecular properties from first principles. Among these properties, those associated with the response of molecules to external electromagnetic fields—such as (hyper)polarizabilities and (hyper)magnetizabilities—play a critical role in understanding spectroscopic behavior, intermolecular interactions, and optical phenomena[1,2,3,4,5,6,7,8,9,10,11,12,13].

A straightforward way to evaluate such response properties is the finite-field approach[14,15], which computes changes in self-consistent field (SCF) energy under small variations of the applied field. In each field direction, this approach requires a few single-point calculations at slightly different field strengths to approximate a derivative of a given order. However, its accuracy is highly sensitive to the choice of field step size, and the accumulated numerical error becomes more pronounced for higher-order derivatives. This finite-field approach, moreover, is inherently limited to static field perturbations and cannot capture frequency-dependent response properties induced by time-dependent electromagnetic fields.

Another commonly used approach is the sum-over-states method[16,17,18,19], which expresses response properties as explicit sums over electronic excited states. This formulation is derived by taking derivatives of the expectation value of the dipole operator $\langle \widehat{\mathsf{\mu }}\rangle$. It originates from Rayleigh-Schrödinger perturbation theory, where the system’s response is expanded in terms of eigenstates of the unperturbed Hamiltonian. In principle, the expansion involves a complete set of many-electron excited states, rendering the method formally exact. In practice, however, the summation must be truncated to a manageable number of states in an infinite-dimensional Hilbert space, which introduces numerical approximations. While conceptually straightforward and potentially accurate when a sufficiently large and well-described excitation manifold is available, the sum-over-states method is computationally demanding and strongly dependent on the quality and completeness of the excited states. These limitations significantly hinder its routine applicability, especially for large systems and higher-order response properties.

A more widely applicable and systematically accurate alternative is the analytical derivative method[20,21,22,23,24,25], in which molecular response properties are obtained by expanding relevant physical quantities in a Taylor series with respect to external perturbation parameters. Within this formalism, the response of molecular orbitals to the applied field must be explicitly taken into account. If the ground-state orbitals are well-suited, it is often assumed that orbital responses of all orders can be represented as linear transformations of the unperturbed orbitals. Under this approximation, the resulting set of equations reduces to a linear system for the transformation matrix, which is often referred to as the coupled-perturbed (CP) Hartree-Fock (HF) equations. These equations, derived from the stationarity conditions of the underlying variational energy functionals, provide a rigorous route to compute response properties of arbitrary order. Because higher-order response properties are very sensitive to electron correlation effects[20,26], attempts to improve predictive accuracy have motivated the extension of analytical derivative theory beyond the HF level. Although post HF methods can offer systematic improvements, they are computationally demanding and impractical for large systems. The use of nonorthogonal localized molecular orbitals can reduce the cost in large-scale calculations[13,27,28]. In contrast, Kohn-Sham (KS) density functional theory (DFT) provides a more computationally affordable means of incorporating electron correlation, prompting several researchers[26,29,30,31,32,33] to extend the analytical derivative formalism to the KS framework. Compared with its HF counterpart, the KS-based formulation introduces additional complexity at higher orders—for instance, third-order response properties require explicit evaluation of the third derivative of the exchange-correlation (xc) functional with respect to the external perturbation.

In the historical development of the analytical derivative theory, there have emerged three different formulations, where the response properties are defined as the derivatives of (i) the expectation value of the dipole operator $\langle \widehat{\mathsf{\mu }}\rangle$, (ii) the energy $E=\langle \widehat{\mathcal{H}}\rangle$, or (iii) the so-called quasi-energy $W=\langle \widehat{\mathcal{H}}-\mathrm{i}{\partial }_t\rangle$.

Sekino and Bartlett[20] first introduced a systematic formulation at the HF level for evaluating static and frequency-dependent polarizabilities, based on derivatives of the expectation value of the dipole operator $\langle \widehat{\mathsf{\mu }}\rangle$. They thought that, whether the perturbation was static or dynamic, response properties derived from the dipole expectation value were equivalent to those obtained from energy derivatives. Building on this foundation, Rice et al. [21] and, independently, Karna and Dupuis[22] later proposed simplified formulations by leveraging the $(2n+1)$ rule, which streamlines the computation of higher-order responses. Their framework was subsequently adopted in numerous studies, including various extensions to KS-DFT[26,29,30,31], and has been widely implemented in modern electronic structure codes. In principle, response properties obtained from this framework are equivalent to those derived from the full, untruncated sum-over-states method. However, despite its popularity, the dipole-derivative-based framework suffers from conceptual and formal inconsistencies. First, the Hamiltonian does not commute with the dipole operator and hence cannot share its eigenstates, making the dipole expectation value, strictly speaking, not a well-defined observable. Moreover, the dipole operator is incompatible with periodic boundary conditions in crystalline systems—a well-known inconsistency that has motivated alternative formulations such as the KGB method, which introduces the vector potential to circumvent this issue[34]. As a consequence, response of $\langle \widehat{\mathsf{\mu }}\rangle$ to perturbations lacks rigorous physical meaning. Second, frequency-dependent response quantities, such as the first hyperpolarizability tensor ${\beta }^{abc}(-{\omega }_B-{\omega }_C;{\omega }_B,{\omega }_C)$, can violate the symmetry of mixed partial derivatives, and the frequency dependence associated with the first dimension is not treated in a physically meaningful way. Furthermore, the simplification of the second-order response density matrix necessitates the use of unphysical constructs involving the CP-HF/KS equations that do not exist under such electromagnetic perturbations. More fundamentally, for a time-dependent perturbation, the traditional Hellmann-Feynman theorem is no longer valid. There is no physical justification for the assumption that the electron density remains unaffected by the external field. In other words, the dipole moment, i.e., the first derivative of the energy E with respect to the external field, cannot be equated to $\langle \widehat{\mathsf{\mu }}\rangle$. The analytical forms of energy and dipole derivatives are no longer equivalent in the presence of dynamic fields.

Although, to the best of our knowledge, there have been no published formulations that rigorously express frequency-dependent response properties as energy derivatives under time-varying fields. Our work shows that, in the presence of time-varying fields, higher-order energy derivatives naturally contain additional terms with explicit frequency dependence, whose contributions become increasingly significant at high frequencies. Nevertheless, defining response properties via energy derivatives under time-dependent perturbations is not without its own theoretical challenges. For such systems, one must resort to the time-dependent Schrödinger equation, where energy ceases to be a well-defined quantity. The Hamiltonian expectation value only carries a time-averaged meaning and is not an observable in the conventional sense. In addition, for magnetic response properties, issues such as gauge-origin dependence become particularly pronounced. Our deeper analysis, however, reveals that the commonly perceived gauge problem stems from an incomplete treatment: while most existing methods account for the gauge transformation of the vector potential in the operators, they neglect the simultaneous transformation of both the basis functions and the scalar potential. In reality, electromagnetic fields exhibit gauge invariance, and the HF and KS equations should remain gauge covariant. Gauge transformations amount to translations of energy-momentum space and leave the physical system unchanged. Under such a transformation, the wavefunction merely acquires a phase factor ${\mathrm{e}}^{-\mathrm{i}g}$, with g being the gauge function associated with a chosen gauge origin. While this resolves the gauge-dependence of energy response in static systems, it also underscores the fact that energy is not a gauge-invariant quantity under time-dependent magnetic fields. As such, energy-derivative-based magnetic response properties remain inherently ill-defined in this context. This motivates a shift to quasi-energy-based formulations as a more consistent foundation for electromagnetic response theory.

The quasi-energy W is a well-defined, gauge-invariant observable in general, owing to the gauge covariance of the HF and KS equations. Notably, Rice and Handy[23], as well as Sasagane, Aiga, and Itoh[25], first suggested the idea of defining response properties as quasi-energy derivatives, though they did not provide a complete derivation. Later, Banerjee and Harbola[32], together with Aiga and coworkers[33], advocated using the Floquet-state[24] energy as a means to characterize response properties. However, the time-averaged Floquet energy is itself time-independent and therefore does not yield general and meaningful field derivatives. For these reasons, we adopt the quasi-energy derivative framework proposed by Sasagane et al., which offers a physically consistent foundation for treating time-dependent perturbations and recovers conventional response theories in the static limit ($\omega \to 0$).

In this work, we present a systematic derivation of analytical quasi-energy derivatives (up to third order) at the SCF level where the electronic wavefunction is described by a single Slater determinant. These analytical derivatives are based on a straightforward Taylor expansion in the external fields and can be applied to a wide range of electromagnetic response properties, including electric-field-induced polarizabilities and magnetic-field-induced magnetizabilities. The resulting formulation not only offers enhanced physical clarity for arbitrary-order derivatives but is also amenable to implementation. Moreover, we generalize the framework to arbitrary time-dependent electromagnetic perturbations, going beyond the conventional setting of static and monochromatic fields, and thereby extend the applicability of quasi-energy response theory to a broader class of physical phenomena such as second-harmonic generation (SHG). In addition, we address a commonly neglected yet essential aspect of analytical derivative theory: the explicit perturbation dependence of basis functions. In some cases, a basis set may be adapted to or parameterized by external fields, which introduces an additional layer of complexity. Under time-dependent perturbations, the operator $\mathrm{i}{\partial }_t$ also acts on such a basis. To account for this, we develop a generalized derivative theory that includes the explicit response of basis functions to external fields, ensuring that all sources of perturbation dependence are consistently treated.

2. General Formalism

2.1. Time-Harmonic Electromagnetic Fields

An electromagnetic field may be either static or dynamic, with the latter exhibiting time dependence. In most cases, a time-dependent field can be modeled as a time-harmonic oscillation. According to Maxwell’s equations, the time-dependent part of such a harmonic field typically consists of two linearly independent complex-conjugate solutions, ${\mathrm{e}}^{\mathrm{i}\omega t}$ and ${\mathrm{e}}^{-\mathrm{i}\omega t}$. Any linear combination of these components can be used to represent the time-harmonic field, i.e.,

$$\mathcal{F}(\mathit{r},t)={\mathcal{F}}^{+}\left(\mathit{r}\right){\mathrm{e}}^{\mathrm{i}\omega t}+{\mathcal{F}}^{-}\left(\mathit{r}\right){\mathrm{e}}^{-\mathrm{i}\omega t}.$$

(1)

Sometimes, multiple fields may be present simultaneously, possibly including a static component. The overall effect is then described by the vector sum of the individual fields,

$$\mathcal{F}(\mathit{r},t)=\sum _{\begin{array}{c}A\\ ({\omega }_A\ge 0)\end{array}}\left[{\mathcal{F}}_A^{+}\left(\mathit{r}\right){\mathrm{e}}^{\mathrm{i}{\omega }_{A}t}+{\mathcal{F}}_A^{-}\left(\mathit{r}\right){\mathrm{e}}^{-\mathrm{i}{\omega }_{A}t}\right]=\sum _A{\mathcal{F}}_A\left(\mathit{r}\right){\mathrm{e}}^{\mathrm{i}{\omega }_{A}t}.$$

(2)

In particular, the case of ${\omega }_A=0$ corresponds to a static field. Eq. (2) represents the general form of an electromagnetic field. Even in the case of an anharmonic time-varying field, expansion in a Fourier series allows it to be expressed as a linear combination of harmonic components, each with a specific frequency ${\omega }_A$.

We now turn to the question of how the quasi-energy of a system responds to these fields. For notational simplicity, the partial derivative of a physical quantity with respect to the a-th Cartesian component of an independent field ${\mathcal{F}}_A(\mathit{r},t;{\omega }_A)$, evaluated at zero field, is denoted by placing a superscript a and appending its characteristic frequency ${\omega }_A$ in parentheses, e.g.,

$$E^{Aa}\equiv E^a\left({\omega }_A\right)={\left.{\displaystyle \frac{\partial E}{\partial \left[{\mathcal{F}}_{Aa}\left(\mathit{r}\right){\mathrm{e}}^{\mathrm{i}{\omega }_{A}t}\right]}}\right|}_0.$$

(3)

It follows directly that

$$E^{*a}\left({\omega }_A\right)={\left.{\displaystyle \frac{\partial E^{*}}{\partial \left[{\mathcal{F}}_{Aa}\left(\mathit{r}\right){\mathrm{e}}^{\mathrm{i}{\omega }_{A}t}\right]}}\right|}_0={\left.{\displaystyle \frac{\partial E}{\partial \left[{\mathcal{F}}_{Aa}\left(\mathit{r}\right){\mathrm{e}}^{-\mathrm{i}{\omega }_{A}t}\right]}}\right|}_0^{*}=E^{a*}(-{\omega }_A),$$

(4)

where the asterisk (*) denotes complex conjugation and the spatial field distribution ${\mathcal{F}}_{Aa}\left(\mathit{r}\right)$ is assumed to be real.

2.2. Quasi-Energy Response

Under general time-dependent perturbations, the many-electron wavefunction is no longer an eigenstate of the Hamiltonian. Instead, its evolution is governed by the time-dependent Schrödinger equation (in atomic units)

$$\left[\widehat{\mathcal{H}}(\left\{\mathit{r}\right\},t)-\mathrm{i}{\partial }_t\right]\mathsf{\Psi }(\left\{\mathit{r}\right\},t)=0,$$

(5)

where $\left\{\mathit{r}\right\}$ represents the set of all electron coordinates. Following the treatment of Sasagane et al. [25], the wavefunction can be written as

$$\mathsf{\Psi }(\left\{\mathit{r}\right\},t)=\mathsf{\Phi }(\left\{\mathit{r}\right\},t)exp\left[-\mathrm{i}{\int }_{t_0}^t\mathrm{d}{t}^{\prime }W\left(t^{\prime }\right)\right],$$

(6)

where $\mathsf{\Phi }\left(\right\{\mathit{r}\},t)$ is represented as a single Slater determinant and W is the quasi-energy. Substitution of this ansatz into Eq. (5) transforms it into an eigenvalue equation

$$\left[\widehat{\mathcal{H}}(\left\{\mathit{r}\right\},t)-\mathrm{i}{\partial }_t\right]\mathsf{\Phi }(\left\{\mathit{r}\right\},t)=W\left(t\right)\mathsf{\Phi }(\left\{\mathit{r}\right\},t).$$

(7)

If the Hamiltonian carries no explicit time dependence, Eq. (7) reduces to the conventional time-independent Schrödinger equation. Assuming that $\mathsf{\Phi }\left(\right\{\mathit{r}\},t)$ remains normalized at all times, we can evaluate the quasi-energy $W\left(t\right)$ as

$$W\left(t\right)=\underset{E\left(t\right)}{\underbrace{\langle \mathsf{\Phi }\left|\widehat{\mathcal{H}}\right|\mathsf{\Phi }\rangle }}-\underset{T\left(t\right)}{\underbrace{\langle \mathsf{\Phi }\left|\mathrm{i}{\partial }_t\right|\mathsf{\Phi }\rangle }}.$$

(8)

The first term in Eq. (8) corresponds to the instantaneous energy of the system, which can be written in terms of the density matrix $D_{\nu \mu }$ as

$$E=D_{\nu \mu }{H}_{\mu \nu }+{\displaystyle \frac{1}{2}}{D}_{\nu \mu }{D}_{\lambda \kappa }\left(\mu \nu \right|\kappa \lambda )+E_{\mathrm{xc}},$$

(9)

where $H_{\mu \nu }$ denotes the core Hamiltonian matrix and $E_{\mathrm{xc}}$ is the xc energy. The Einstein summation convention is adopted throughout. The second term in Eq. (8) represents the dynamic phase arising from the time-dependence of the wavefunction, and evaluates to

$$T=D_{\nu \mu }{\mathcal{T}}_{\mu \nu }+\mathrm{i}{n}_i{C}_{\mu i}^{*}{S}_{\mu \nu }{\partial }_t{C}_{\nu i},$$

(10)

where $n_i$ is the occupation number of the i-th orbital, $C_{\mu i}$ denotes the matrix of orbital coefficients, $S_{\mu \nu }$ is the overlap matrix in the basis function representation, and ${\mathcal{T}}_{\mu \nu }$ is defined as

$${\mathcal{T}}_{\mu \nu }=\left(\mu \right|\mathrm{i}{\partial }_t\left|\nu \right).$$

(11)

Theoretically, the responses of ${\mathcal{T}}_{\mu \nu }$ to a static field should vanish at all orders, i.e.,

$$\begin{array}{cc}\hfill {\mathcal{T}}_{\mu \nu }^0& =0,\hfill \end{array}$$

(12a)

$$\begin{array}{cc}\hfill {\mathcal{T}}_{\mu \nu }^a\left(0\right)& =0,\hfill \end{array}$$

(12b)

$$\begin{array}{cc}\hfill {\mathcal{T}}_{\mu \nu }^{ab}(0,0)& =0,\hfill \\ \hfill & ⋮\hfill \end{array}$$

(12c)

In particular, if the basis functions are independent of the external perturbation, ${\mathcal{T}}_{\mu \nu }$ vanishes identically.

The energy response can be obtained directly by taking successive derivatives of Eq. (9). The zeroth-order energy reads

$$E^0=D_{\nu \mu }^0{H}_{\mu \nu }^0+{\displaystyle \frac{1}{2}}{D}_{\nu \mu }^0{D}_{\lambda \kappa }^0{\left(\mu \nu \right|\kappa \lambda )}^0+E_{\mathrm{xc}}^0,$$

(13)

which characterizes the unperturbed system. The first-order energy response is given by

$$E^{Aa}=D_{\nu \mu }^{Aa}{F}_{\mu \nu }^0+D_{\nu \mu }^0{H}_{\mu \nu }^{Aa}+{\displaystyle \frac{1}{2}}{D}_{\nu \mu }^0{D}_{\lambda \kappa }^0{\left(\mu \nu \right|\kappa \lambda )}^{Aa}+D_{\nu \mu }^0\left({\mu \nu }^{Aa}\right|v_{\mathrm{xc}}),$$

(14)

where $v_{\mathrm{xc}}$ denotes the xc potential, and $F_{\mu \nu }^0$ is the zeroth-order Fock matrix, defined as

$$F_{\mu \nu }^0=H_{\mu \nu }^0+D_{\lambda \kappa }^0{\left(\mu \nu \right|\kappa \lambda )}^0+\left({\mu \nu }^0\right|v_{\mathrm{xc}}).$$

(15)

Note that the derivative appearing in an xc integral acts on the basis function pair, e.g., $\left({\mu \nu }^{Aa}\right|v_{\mathrm{xc}})\equiv \left({\left(\mu \nu \right)}^{Aa}\right|v_{\mathrm{xc}})$. The second-order energy response takes the form

$$\begin{array}{cc}\hfill E^{AaBb}& =D_{\nu \mu }^{AaBb}{F}_{\mu \nu }^0+D_{\nu \mu }^{Aa}{F}_{\mu \nu }^{Bb}+D_{\nu \mu }^{Bb}{H}_{\mu \nu }^{Aa}+D_{\nu \mu }^0{H}_{\mu \nu }^{AaBb}\hfill \\ \hfill & +D_{\nu \mu }^{Bb}{D}_{\lambda \kappa }^0{\left(\mu \nu \right|\kappa \lambda )}^{Aa}+{\displaystyle \frac{1}{2}}{D}_{\nu \mu }^0{D}_{\lambda \kappa }^0{\left(\mu \nu \right|\kappa \lambda )}^{AaBb}\hfill \\ \hfill & +D_{\nu \mu }^{Bb}\left({\mu \nu }^{Aa}\right|v_{\mathrm{xc}})+D_{\nu \mu }^0\left({\mu \nu }^{AaBb}\right|v_{\mathrm{xc}})+D_{\nu \mu }^0({\mu \nu }^{Aa}|f_{\mathrm{xc}}\left|{\rho }^{Bb}\right),\hfill \end{array}$$

(16)

where $f_{\mathrm{xc}}$ is the xc kernel, i.e., the second functional derivative of $E_{\mathrm{xc}}$, and $\rho$ is not restricted to the electron density but depends on the variables entering the functional. The third-order energy response

$$\begin{array}{cc}\hfill E^{AaBbCc}& =D_{\nu \mu }^{AaBbCc}{F}_{\mu \nu }^0+D_{\nu \mu }^{AaBb}{F}_{\mu \nu }^{Cc}+D_{\nu \mu }^{AaCc}{F}_{\mu \nu }^{Bb}+D_{\nu \mu }^{Aa}{F}_{\mu \nu }^{BbCc}+D_{\nu \mu }^{BbCc}{H}_{\mu \nu }^{Aa}\hfill \\ \hfill & +D_{\nu \mu }^{Bb}{H}_{\mu \nu }^{AaCc}+D_{\nu \mu }^{Cc}{H}_{\mu \nu }^{AaBb}+D_{\nu \mu }^0{H}_{\mu \nu }^{AaBbCc}+D_{\nu \mu }^{BbCc}{D}_{\lambda \kappa }^0{\left(\mu \nu \right|\kappa \lambda )}^{Aa}\hfill \\ \hfill & +D_{\nu \mu }^{Bb}{D}_{\lambda \kappa }^{Cc}{\left(\mu \nu \right|\kappa \lambda )}^{Aa}+D_{\nu \mu }^{Bb}{D}_{\lambda \kappa }^0{\left(\mu \nu \right|\kappa \lambda )}^{AaCc}+D_{\nu \mu }^{Cc}{D}_{\lambda \kappa }^0{\left(\mu \nu \right|\kappa \lambda )}^{AaBb}\hfill \\ \hfill & +{\displaystyle \frac{1}{2}}{D}_{\nu \mu }^0{D}_{\lambda \kappa }^0{\left(\mu \nu \right|\kappa \lambda )}^{AaBbCc}+D_{\nu \mu }^{BbCc}\left({\mu \nu }^{Aa}\right|v_{\mathrm{xc}})+D_{\nu \mu }^{Bb}\left({\mu \nu }^{AaCc}\right|v_{\mathrm{xc}})\hfill \\ \hfill & +D_{\nu \mu }^{Bb}({\mu \nu }^{Aa}|f_{\mathrm{xc}}\left|{\rho }^{Cc}\right)+D_{\nu \mu }^{Cc}\left({\mu \nu }^{AaBb}\right|v_{\mathrm{xc}})+D_{\nu \mu }^0\left({\mu \nu }^{AaBbCc}\right|v_{\mathrm{xc}})\hfill \\ \hfill & +D_{\nu \mu }^0({\mu \nu }^{AaBb}|f_{\mathrm{xc}}|{\rho }^{Cc})+D_{\nu \mu }^{Cc}({\mu \nu }^{Aa}|f_{\mathrm{xc}}|{\rho }^{Bb})+D_{\nu \mu }^0({\mu \nu }^{AaCc}|f_{\mathrm{xc}}\left|{\rho }^{Bb}\right)\hfill \\ \hfill & +D_{\nu \mu }^0({\mu \nu }^{Aa}|f_{\mathrm{xc}}|{\rho }^{BbCc})+D_{\nu \mu }^0({\mu \nu }^{Aa}|k_{\mathrm{xc}}|{\rho }^{Bb}\left|{\rho }^{Cc}\right)\hfill \end{array}$$

(17)

involves $k_{\mathrm{xc}}$, the third functional derivative of $E_{\mathrm{xc}}$. All terms that involve $v_{\mathrm{xc}}$, $f_{\mathrm{xc}}$, and $k_{\mathrm{xc}}$ correspond to xc contributions of increasing order. Particularly, $k_{\mathrm{xc}}$ vanishes in the case of a pure HF exchange.

The response of T cannot be directly obtained by taking derivatives of Eq. (10), as the orbital coefficients $C_{\mu i}$ implicitly depend on the external perturbation. To resolve this, we expand both sides of Eq. (10) as Taylor series in powers of the external fields $\left\{\mathcal{F}\right\}$, which results in

$$\begin{array}{cc}\hfill \sum _{p=0}^{\infty }{\displaystyle \frac{T^{\left(p\right)}}{p!}}{\left\{\mathcal{F}\right\}}^p& =\sum _{p=0}^{\infty }\sum _{q=0}^{\infty }{\displaystyle \frac{D_{\nu \mu }^{\left(p\right)}{\mathcal{T}}_{\mu \nu }^{\left(q\right)}}{p!q!}}{\left\{\mathcal{F}\right\}}^p{\left\{\mathcal{F}\right\}}^q\hfill \\ \hfill & +\mathrm{i}{n}_i\sum _{p=0}^{\infty }\sum _{q=0}^{\infty }\sum _{r=0}^{\infty }{\displaystyle \frac{C_{\mu i}^{*\left(p\right)}{S}_{\mu \nu }^{\left(q\right)}{C}_{\nu i}^{\left(r\right)}}{p!q!r!}}{\left\{\mathcal{F}\right\}}^p{\left\{\mathcal{F}\right\}}^q{\partial }_t{\left\{\mathcal{F}\right\}}^r.\hfill \end{array}$$

(18)

Since each field in $\left\{\mathcal{F}\right\}$ is independent, we can collect terms of equal perturbation order on both sides. The expansion up to third order yields the following equations:

$$T^0=0,$$

(19a)

$$T^{Aa}=D_{\nu \mu }^0{\mathcal{T}}_{\mu \nu }^{Aa}-{\omega }_A{n}_i{C}_{\mu i}^{*0}{S}_{\mu \nu }^0{C}_{\nu i}^{Aa},$$

(19b)

$$\begin{array}{cc}\hfill T^{AaBb}& =D_{\nu \mu }^{Aa}{\mathcal{T}}_{\mu \nu }^{Bb}+D_{\nu \mu }^{Bb}{\mathcal{T}}_{\mu \nu }^{Aa}+D_{\nu \mu }^0{\mathcal{T}}_{\mu \nu }^{AaBb}\hfill \\ \hfill & -{\omega }_B{n}_i\left[C_{\mu i}^{*Aa}{S}_{\mu \nu }^0+C_{\mu i}^{*0}{S}_{\mu \nu }^{Aa}\right]C_{\nu i}^{Bb}\hfill \\ \hfill & -{\omega }_A{n}_i\left[C_{\mu i}^{*Bb}{S}_{\mu \nu }^0+C_{\mu i}^{*0}{S}_{\mu \nu }^{Bb}\right]C_{\nu i}^{Aa}\hfill \\ \hfill & -\left({\omega }_A+{\omega }_B\right)n_i{C}_{\mu i}^{*0}{S}_{\mu \nu }^0{C}_{\nu i}^{AaBb},\hfill \end{array}$$

(19c)

$$\begin{array}{cc}\hfill T^{AaBbCc}& =\left[D_{\nu \mu }^{AaBb}{\mathcal{T}}_{\mu \nu }^{Cc}+D_{\nu \mu }^{AaCc}{\mathcal{T}}_{\mu \nu }^{Bb}+D_{\nu \mu }^{BbCc}{\mathcal{T}}_{\mu \nu }^{Aa}+D_{\nu \mu }^{Aa}{\mathcal{T}}_{\mu \nu }^{BbCc}\right.\hfill \\ \hfill & +\left.D_{\nu \mu }^{Bb}{\mathcal{T}}_{\mu \nu }^{AaCc}+D_{\nu \mu }^{Cc}{\mathcal{T}}_{\mu \nu }^{AaBb}+D_{\nu \mu }^0{\mathcal{T}}_{\mu \nu }^{AaBbCc}\right]\hfill \\ \hfill & -{\omega }_C{n}_i\left[C_{\mu i}^{*AaBb}{S}_{\mu \nu }^0+C_{\mu i}^{*Aa}{S}_{\mu \nu }^{Bb}+C_{\mu i}^{*Bb}{S}_{\mu \nu }^{Aa}+C_{\mu i}^{*0}{S}_{\mu \nu }^{AaBb}\right]C_{\nu i}^{Cc}\hfill \\ \hfill & -{\omega }_B{n}_i\left[C_{\mu i}^{*AaCc}{S}_{\mu \nu }^0+C_{\mu i}^{*Aa}{S}_{\mu \nu }^{Cc}+C_{\mu i}^{*Cc}{S}_{\mu \nu }^{Aa}+C_{\mu i}^{*0}{S}_{\mu \nu }^{AaCc}\right]C_{\nu i}^{Bb}\hfill \\ \hfill & -{\omega }_A{n}_i\left[C_{\mu i}^{*BbCc}{S}_{\mu \nu }^0+C_{\mu i}^{*Bb}{S}_{\mu \nu }^{Cc}+C_{\mu i}^{*Cc}{S}_{\mu \nu }^{Bb}+C_{\mu i}^{*0}{S}_{\mu \nu }^{BbCc}\right]C_{\nu i}^{Aa}\hfill \\ \hfill & -\left({\omega }_B+{\omega }_C\right)n_i\left[C_{\mu i}^{*Aa}{S}_{\mu \nu }^0+C_{\mu i}^{*0}{S}_{\mu \nu }^{Aa}\right]C_{\nu i}^{BbCc}\hfill \\ \hfill & -\left({\omega }_A+{\omega }_C\right)n_i\left[C_{\mu i}^{*Bb}{S}_{\mu \nu }^0+C_{\mu i}^{*0}{S}_{\mu \nu }^{Bb}\right]C_{\nu i}^{AaCc}\hfill \\ \hfill & -\left({\omega }_A+{\omega }_B\right)n_i\left[C_{\mu i}^{*Cc}{S}_{\mu \nu }^0+C_{\mu i}^{*0}{S}_{\mu \nu }^{Cc}\right]C_{\nu i}^{AaBb}\hfill \\ \hfill & -\left({\omega }_A+{\omega }_B+{\omega }_C\right)n_i{C}_{\mu i}^{*0}{S}_{\mu \nu }^0{C}_{\nu i}^{AaBbCc}.\hfill \end{array}$$

(19d)

A detailed derivation is provided in Appendix A.

The quasi-energy responses up to third order are thus given by

$$\begin{array}{cc}\hfill W^0& =E^0,\hfill \end{array}$$

(20a)

$$\begin{array}{cc}\hfill W^{Aa}& =E^{Aa}-T^{Aa},\hfill \end{array}$$

(20b)

$$\begin{array}{cc}\hfill W^{AaBb}& =E^{AaBb}-T^{AaBb},\hfill \end{array}$$

(20c)

$$\begin{array}{cc}\hfill W^{AaBbCc}& =E^{AaBbCc}-T^{AaBbCc}.\hfill \end{array}$$

(20d)

2.3. Density Response

The density matrix $D_{\nu \mu }$ is defined as

$$D_{\nu \mu }=n_i{C}_{\nu i}{C}_{\mu i}^{*}.$$

(21)

Since the quasi-energy responses involve the density matrix at each order, it is necessary to evaluate how the density matrix (21) varies with respect to external perturbations. This can be directly achieved via differentiation. In the absence of perturbation, we recover the zeroth-order density matrix

$$D_{\nu \mu }^0=n_i{C}_{\nu i}^0{C}_{\mu i}^{0*},$$

(22)

which corresponds to the ground-state density. The first-order response of the density matrix is given by

$$\begin{array}{cc}\hfill D_{\nu \mu }^{Aa}& \equiv D_{\nu \mu }^a\left({\omega }_A\right)\hfill \\ \hfill & =n_i\left[C_{\nu i}^a\left({\omega }_A\right)C_{\mu i}^{0*}+C_{\nu i}^0{C}_{\mu i}^{a*}(-{\omega }_A)\right],\hfill \end{array}$$

(23)

the second-order response is

$$\begin{array}{cc}\hfill D_{\nu \mu }^{AaBb}& \equiv D_{\nu \mu }^{ab}({\omega }_A,{\omega }_B)\hfill \\ \hfill & =n_i\left[C_{\nu i}^{ab}({\omega }_A,{\omega }_B)C_{\mu i}^{0*}+C_{\nu i}^a\left({\omega }_A\right)C_{\mu i}^{b*}(-{\omega }_B)\right.\hfill \\ \hfill & +\left.C_{\nu i}^b\left({\omega }_B\right)C_{\mu i}^{a*}(-{\omega }_A)+C_{\nu i}^0{C}_{\mu i}^{ab*}(-{\omega }_A,-{\omega }_B)\right],\hfill \end{array}$$

(24)

and the third-order response reads

$$\begin{array}{cc}\hfill D_{\nu \mu }^{AaBbCc}& \equiv D_{\nu \mu }^{abc}({\omega }_A,{\omega }_B,{\omega }_C)\hfill \\ \hfill & =n_i\left[C_{\nu i}^{abc}({\omega }_A,{\omega }_B,{\omega }_C)C_{\mu i}^{0*}+C_{\nu i}^{ab}({\omega }_A,{\omega }_B)C_{\mu i}^{c*}(-{\omega }_C)\right.\hfill \\ \hfill & +C_{\nu i}^{ac}({\omega }_A,{\omega }_C)C_{\mu i}^{b*}(-{\omega }_B)+C_{\nu i}^a\left({\omega }_A\right)C_{\mu i}^{bc*}(-{\omega }_B,-{\omega }_C)\hfill \\ \hfill & +C_{\nu i}^{bc}({\omega }_B,{\omega }_C)C_{\mu i}^{a*}(-{\omega }_A)+C_{\nu i}^b\left({\omega }_B\right)C_{\mu i}^{ac*}(-{\omega }_A,-{\omega }_C)\hfill \\ \hfill & +\left.C_{\nu i}^c\left({\omega }_C\right)C_{\mu i}^{ab*}(-{\omega }_A,-{\omega }_B)+C_{\nu i}^0{C}_{\mu i}^{abc*}(-{\omega }_A,-{\omega }_B,-{\omega }_C)\right].\hfill \end{array}$$

(25)

For convenience in the subsequent derivations, we assume that the orbital coefficients at each order are related to the ground-state coefficients through a linear transformation. These transformation matrices $\mathit{U}$ are defined by

$$\begin{array}{cc}\hfill C_{\mu i}^a\left({\omega }_A\right)& =C_{\mu j}^0{U}_{ji}^a\left({\omega }_A\right),\hfill \end{array}$$

(26a)

$$\begin{array}{cc}\hfill C_{\mu i}^{ab}({\omega }_A,{\omega }_B)& =C_{\mu j}^0{U}_{ji}^{ab}({\omega }_A,{\omega }_B),\hfill \end{array}$$

(26b)

$$\begin{array}{cc}\hfill C_{\mu i}^{abc}({\omega }_A,{\omega }_B,{\omega }_C)& =C_{\mu j}^0{U}_{ji}^{abc}({\omega }_A,{\omega }_B,{\omega }_C),\hfill \\ \hfill & ⋮\hfill \end{array}$$

(26c)

Then, Eqs. (23)–(25) become

$$D_{\nu \mu }^a\left({\omega }_A\right)=n_i\left[C_{\nu j}^0{U}_{ji}^a\left({\omega }_A\right)C_{\mu i}^{0*}+C_{\nu i}^0{C}_{\mu j}^{0*}{U}_{ji}^{a*}(-{\omega }_A)\right],$$

(27)

$$\begin{array}{cc}\hfill D_{\nu \mu }^{ab}({\omega }_A,{\omega }_B)& =n_i\left[C_{\nu j}^0{U}_{ji}^{ab}({\omega }_A,{\omega }_B)C_{\mu i}^{0*}+C_{\nu j}^0{U}_{ji}^a\left({\omega }_A\right)C_{\mu k}^{0*}{U}_{ki}^{b*}(-{\omega }_B)\right.\hfill \\ \hfill & +\left.C_{\nu j}^0{U}_{ji}^b\left({\omega }_B\right)C_{\mu k}^{0*}{U}_{ki}^{a*}(-{\omega }_A)+C_{\nu i}^0{C}_{\mu j}^{0*}{U}_{ji}^{ab*}(-{\omega }_A,-{\omega }_B)\right],\hfill \end{array}$$

(28)

$$\begin{array}{cc}\hfill D_{\nu \mu }^{abc}({\omega }_A,{\omega }_B,{\omega }_C)& =n_i\left[C_{\nu j}^0{U}_{ji}^{abc}({\omega }_A,{\omega }_B,{\omega }_C)C_{\mu i}^{0*}+C_{\nu j}^0{U}_{ji}^{ab}({\omega }_A,{\omega }_B)C_{\mu k}^{0*}{U}_{ki}^{c*}(-{\omega }_C)\right.\hfill \\ \hfill & +C_{\nu j}^0{U}_{ji}^{ac}({\omega }_A,{\omega }_C)C_{\mu k}^{0*}{U}_{ki}^{b*}(-{\omega }_B)+C_{\nu j}^0{U}_{ji}^a\left({\omega }_A\right)C_{\mu k}^{0*}{U}_{ki}^{bc*}(-{\omega }_B,-{\omega }_C)\hfill \\ \hfill & +C_{\nu j}^0{U}_{ji}^{bc}({\omega }_B,{\omega }_C)C_{\mu k}^{0*}{U}_{ki}^{a*}(-{\omega }_A)+C_{\nu j}^0{U}_{ji}^b\left({\omega }_B\right)C_{\mu k}^{0*}{U}_{ki}^{ac*}(-{\omega }_A,-{\omega }_C)\hfill \\ \hfill & +\left.C_{\nu j}^0{U}_{ji}^c\left({\omega }_C\right)C_{\mu k}^{0*}{U}_{ki}^{ab*}(-{\omega }_A,-{\omega }_B)+C_{\nu i}^0{C}_{\mu j}^{0*}{U}_{ji}^{abc*}(-{\omega }_A,-{\omega }_B,-{\omega }_C)\right].\hfill \end{array}$$

(29)

It follows that

$$\begin{array}{cc}\hfill D_{\nu \mu }^a\left({\omega }_A\right)& =D_{\mu \nu }^{a*}(-{\omega }_A),\hfill \end{array}$$

(30a)

$$\begin{array}{cc}\hfill D_{\nu \mu }^{ab}({\omega }_A,{\omega }_B)& =D_{\mu \nu }^{ab*}(-{\omega }_A,-{\omega }_B),\hfill \end{array}$$

(30b)

$$\begin{array}{cc}\hfill D_{\nu \mu }^{abc}({\omega }_A,{\omega }_B,{\omega }_C)& =D_{\mu \nu }^{abc*}(-{\omega }_A,-{\omega }_B,-{\omega }_C),\hfill \\ \hfill & ⋮\hfill \end{array}$$

(30c)

2.4. Orthonormality

We assume that the molecular orbitals remain orthonormal at all times, i.e.,

$$S_{ij}=C_{\mu i}^{*}{S}_{\mu \nu }{C}_{\nu j}={\delta }_{ij}.$$

(31)

Note that matrices with subscripts i, j, k, etc. are represented in the molecular orbital basis. Since the ground-state orbitals are already orthonormal, we have

$$S_{ij}^0=C_{\mu i}^{0*}{S}_{\mu \nu }^0{C}_{\nu j}^0={\delta }_{ij}.$$

(32)

As a consequence, the higher-order orbital overlaps must vanish. For instance, the first-order condition reads

$$\begin{array}{cc}\hfill S_{ij}^{Aa}& =C_{\mu i}^{a*}(-{\omega }_A)S_{\mu \nu }^0{C}_{\nu j}^0+C_{\mu i}^{0*}{S}_{\mu \nu }^a\left({\omega }_A\right)C_{\nu j}^0+C_{\mu i}^{0*}{S}_{\mu \nu }^0{C}_{\nu j}^a\left({\omega }_A\right)\hfill \\ \hfill & =U_{ki}^{a*}(-{\omega }_A)S_{kj}^0+S_{ij}^a\left({\omega }_A\right)+S_{ik}^0{U}_{kj}^a\left({\omega }_A\right)\hfill \\ \hfill & =U_{ki}^{a*}(-{\omega }_A){\delta }_{kj}+S_{ij}^a\left({\omega }_A\right)+{\delta }_{ik}{U}_{kj}^a\left({\omega }_A\right)\hfill \\ \hfill & =U_{ji}^{a*}(-{\omega }_A)+S_{ij}^a\left({\omega }_A\right)+U_{ij}^a\left({\omega }_A\right)\hfill \\ \hfill & =0,\hfill \end{array}$$

(33)

where $S_{ij}^a\left({\omega }_A\right)=C_{\mu i}^{0*}{S}_{\mu \nu }^a\left({\omega }_A\right)C_{\nu j}^0$. In our notation, when perturbation indices appear in the superscript of a molecular orbital integral, such as $S_{ij}^{Aa}\equiv S_{ij}^a\left({\omega }_A\right)$, they denote a derivative acting solely on the matrix elements in the atomic orbital basis, not on the orbital coefficients. To indicate a total derivative, the perturbation indices are placed as a superscript on the entire matrix element, as in $S_{ij}^{Aa}$. It is worth noting that the derivatives of $S_{\mu \nu }$ with respect to different fields should be identical, such that $S_{\mu \nu }^{Aa}\equiv S_{\mu \nu }^a\left({\omega }_A\right)=S_{\mu \nu }^a$. Accordingly, we omit the frequency argument in parentheses for the derivatives of $S_{\mu \nu }$ throughout this work. Following a similar procedure as in the first-order case, one can derive the second- and third-order orthonormality conditions as

$$\begin{array}{cc}\hfill S_{ij}^{AaBb}& =U_{ji}^{ab*}(-{\omega }_A,-{\omega }_B)+U_{ki}^{a*}(-{\omega }_A)S_{kj}^b+U_{ki}^{a*}(-{\omega }_A)U_{kj}^b\left({\omega }_B\right)\hfill \\ \hfill & +U_{ki}^{b*}(-{\omega }_B)S_{kj}^a+S_{ij}^{ab}+S_{ik}^a{U}_{kj}^b\left({\omega }_B\right)\hfill \\ \hfill & +U_{ki}^{b*}(-{\omega }_B)U_{kj}^a\left({\omega }_A\right)+S_{ik}^b{U}_{kj}^a\left({\omega }_A\right)+U_{ij}^{ab}({\omega }_A,{\omega }_B)\hfill \\ \hfill & =0,\hfill \end{array}$$

(34)

$$\begin{array}{cc}\hfill S_{ij}^{AaBbCc}& =U_{ji}^{abc*}(-{\omega }_A,-{\omega }_B,-{\omega }_C)+U_{ki}^{ab*}(-{\omega }_A,-{\omega }_B)S_{kj}^c\hfill \\ \hfill & +U_{ki}^{ab*}(-{\omega }_A,-{\omega }_B)U_{kj}^c\left({\omega }_C\right)+U_{ki}^{ac*}(-{\omega }_A,-{\omega }_C)S_{kj}^b+U_{ki}^{a*}(-{\omega }_A)S_{kj}^{bc}\hfill \\ \hfill & +U_{ki}^{a*}(-{\omega }_A)S_{kl}^b{U}_{lj}^c\left({\omega }_C\right)+U_{ki}^{ac*}(-{\omega }_A,-{\omega }_C)U_{kj}^b\left({\omega }_B\right)\hfill \\ \hfill & +U_{ki}^{a*}(-{\omega }_A)S_{kl}^c{U}_{lj}^b\left({\omega }_B\right)+U_{ki}^{a*}(-{\omega }_A)U_{kj}^{bc}({\omega }_B,{\omega }_C)\hfill \\ \hfill & +U_{ki}^{bc*}(-{\omega }_B,-{\omega }_C)S_{kj}^a+U_{ki}^{b*}(-{\omega }_B)S_{kj}^{ac}+U_{ki}^{b*}(-{\omega }_B)S_{kl}^a{U}_{lj}^c\left({\omega }_C\right)\hfill \\ \hfill & +U_{ki}^{c*}(-{\omega }_C)S_{kj}^{ab}+S_{ij}^{abc}+S_{ik}^{ab}{U}_{kj}^c\left({\omega }_C\right)+U_{ki}^{c*}(-{\omega }_C)S_{kl}^a{U}_{lj}^b\left({\omega }_B\right)\hfill \\ \hfill & +S_{ik}^{ac}{U}_{kj}^b\left({\omega }_B\right)+S_{ik}^a{U}_{kj}^{bc}({\omega }_B,{\omega }_C)+U_{ki}^{bc*}(-{\omega }_B,-{\omega }_C)U_{kj}^a\left({\omega }_A\right)\hfill \\ \hfill & +U_{ki}^{b*}(-{\omega }_B)S_{kl}^c{U}_{lj}^a\left({\omega }_A\right)+U_{ki}^{b*}(-{\omega }_B)U_{kj}^{ac}({\omega }_A,{\omega }_C)\hfill \\ \hfill & +U_{ki}^{c*}(-{\omega }_C)S_{kl}^b{U}_{lj}^a\left({\omega }_A\right)+S_{ik}^{bc}{U}_{kj}^a\left({\omega }_A\right)+S_{ik}^b{U}_{kj}^{ac}({\omega }_A,{\omega }_C)\hfill \\ \hfill & +U_{ki}^{c*}(-{\omega }_C)U_{kj}^{ab}({\omega }_A,{\omega }_B)+S_{ik}^c{U}_{kj}^{ab}({\omega }_A,{\omega }_B)+U_{ij}^{abc}({\omega }_A,{\omega }_B,{\omega }_C)\hfill \\ \hfill & =0.\hfill \end{array}$$

(35)

These expressions impose constraints on the orbital transformation matrices $\mathit{U}$.

2.5. Determination of U: The CP-HF/KS Equations

To determine the transformation matrices $\mathit{U}$, we minimize the quasi-energy Eq. (8) with respect to orbital variations, subject to the constraint of orbital orthonormality. This leads to a Roothaan-like equation

$$\left[F_{\mu \nu }-{\mathcal{T}}_{\mu \nu }\right]C_{\nu i}-\mathrm{i}{S}_{\mu \nu }{\partial }_t{C}_{\nu i}=S_{\mu \nu }{C}_{\nu j}{\epsilon }_{ji}$$

(36)

in the atomic orbital basis. To facilitate the treatment of the time derivative ${\partial }_t$ acting on the orbital coefficients, we expand each quantity in the above expression as a Taylor series in the external fields:

$$\begin{array}{cc}\hfill & \sum _{p=0}^{\infty }\sum _{q=0}^{\infty }{\displaystyle \frac{{\left[F_{\mu \nu }-{\mathcal{T}}_{\mu \nu }\right]}^{\left(p\right)}{C}_{\nu i}^{\left(q\right)}}{p!q!}}{\left\{\mathcal{F}\right\}}^p{\left\{\mathcal{F}\right\}}^q-\mathrm{i}\sum _{p=0}^{\infty }\sum _{q=0}^{\infty }{\displaystyle \frac{S_{\mu \nu }^{\left(p\right)}{C}_{\nu i}^{\left(q\right)}}{p!q!}}{\left\{\mathcal{F}\right\}}^p{\partial }_t{\left\{\mathcal{F}\right\}}^q\hfill \\ \hfill & =\sum _{p=0}^{\infty }\sum _{q=0}^{\infty }\sum _{r=0}^{\infty }{\displaystyle \frac{S_{\mu \nu }^{\left(p\right)}{C}_{\nu j}^{\left(q\right)}{\epsilon }_{ji}^{\left(r\right)}}{p!q!r!}}{\left\{\mathcal{F}\right\}}^p{\left\{\mathcal{F}\right\}}^q{\left\{\mathcal{F}\right\}}^r.\hfill \end{array}$$

(37)

In analogy with the derivation of the responses of T, we obtain the original CP-HF/KS equations (up to second order)

$$F_{\mu \nu }^0{C}_{\nu i}^0=S_{\mu \nu }^0{C}_{\nu j}^0{\epsilon }_{ji}^0,$$

(38a)

$$\begin{array}{cc}\hfill & F_{\mu \nu }^a\left({\omega }_A\right)C_{\nu i}^0+F_{\mu \nu }^0{C}_{\nu i}^a\left({\omega }_A\right)-{\mathcal{T}}_{\mu \nu }^a\left({\omega }_A\right)C_{\nu i}^0+{\omega }_A{S}_{\mu \nu }^0{C}_{\nu i}^a\left({\omega }_A\right)\hfill \\ \hfill & =S_{\mu \nu }^a{C}_{\nu j}^0{\epsilon }_{ji}^0+S_{\mu \nu }^0{C}_{\nu j}^a\left({\omega }_A\right){\epsilon }_{ji}^0+S_{\mu \nu }^0{C}_{\nu j}^0{\epsilon }_{ji}^a\left({\omega }_A\right),\hfill \end{array}$$

(38b)

$$\begin{array}{cc}\hfill & F_{\mu \nu }^{ab}({\omega }_A,{\omega }_B)C_{\nu i}^0+F_{\mu \nu }^a\left({\omega }_A\right)C_{\nu i}^b\left({\omega }_B\right)+F_{\mu \nu }^b\left({\omega }_B\right)C_{\nu i}^a\left({\omega }_A\right)+F_{\mu \nu }^0{C}_{\nu i}^{ab}({\omega }_A,{\omega }_B)\hfill \\ \hfill & -{\mathcal{T}}_{\mu \nu }^{ab}({\omega }_A,{\omega }_B)C_{\nu i}^0-{\mathcal{T}}_{\mu \nu }^a\left({\omega }_A\right)C_{\nu i}^b\left({\omega }_B\right)-{\mathcal{T}}_{\mu \nu }^b\left({\omega }_B\right)C_{\nu i}^a\left({\omega }_A\right)\hfill \\ \hfill & +{\omega }_B{S}_{\mu \nu }^a{C}_{\nu i}^b\left({\omega }_B\right)+{\omega }_A{S}_{\mu \nu }^b{C}_{\nu i}^a\left({\omega }_A\right)+\left({\omega }_A+{\omega }_B\right)S_{\mu \nu }^0{C}_{\nu i}^{ab}({\omega }_A,{\omega }_B)\hfill \\ \hfill & =S_{\mu \nu }^{ab}{C}_{\nu j}^0{\epsilon }_{ji}^0+S_{\mu \nu }^a{C}_{\nu j}^b\left({\omega }_B\right){\epsilon }_{ji}^0+S_{\mu \nu }^a{C}_{\nu j}^0{\epsilon }_{ji}^b\left({\omega }_B\right)\hfill \\ \hfill & +S_{\mu \nu }^b{C}_{\nu j}^a\left({\omega }_A\right){\epsilon }_{ji}^0+S_{\mu \nu }^0{C}_{\nu j}^{ab}({\omega }_A,{\omega }_B){\epsilon }_{ji}^0+S_{\mu \nu }^0{C}_{\nu j}^a\left({\omega }_A\right){\epsilon }_{ji}^b\left({\omega }_B\right)\hfill \\ \hfill & +S_{\mu \nu }^b{C}_{\nu j}^0{\epsilon }_{ji}^a\left({\omega }_A\right)+S_{\mu \nu }^0{C}_{\nu j}^b\left({\omega }_B\right){\epsilon }_{ji}^a\left({\omega }_A\right)+S_{\mu \nu }^0{C}_{\nu j}^0{\epsilon }_{ji}^{ab}({\omega }_A,{\omega }_B),\hfill \end{array}$$

(38c)

where ${\epsilon }_{ji}$ is the Lagrange multiplier matrix, and the first- and second-order Fock matrices are, respectively,

$$F_{\mu \nu }^a\left({\omega }_A\right)=H_{\mu \nu }^a+D_{\lambda \kappa }^a\left({\omega }_A\right){\left(\mu \nu \right|\kappa \lambda )}^0+D_{\lambda \kappa }^0{\left(\mu \nu \right|\kappa \lambda )}^a+\left({\mu \nu }^a\right|v_{\mathrm{xc}})+({\mu \nu }^0|f_{\mathrm{xc}}\left|{\rho }^a\left({\omega }_A\right)\right),$$

(39)

and

$$\begin{array}{cc}\hfill F_{\mu \nu }^{ab}({\omega }_A,{\omega }_B)& =H_{\mu \nu }^{ab}+D_{\lambda \kappa }^{ab}({\omega }_A,{\omega }_B){\left(\mu \nu \right|\kappa \lambda )}^0+D_{\lambda \kappa }^a\left({\omega }_A\right){\left(\mu \nu \right|\kappa \lambda )}^b+D_{\lambda \kappa }^b\left({\omega }_B\right){\left(\mu \nu \right|\kappa \lambda )}^a\hfill \\ \hfill & +D_{\lambda \kappa }^0{\left(\mu \nu \right|\kappa \lambda )}^{ab}+\left({\mu \nu }^{ab}\right|v_{\mathrm{xc}})+({\mu \nu }^a|f_{\mathrm{xc}}|{\rho }^b\left({\omega }_B\right))+({\mu \nu }^b|f_{\mathrm{xc}}\left|{\rho }^a\left({\omega }_A\right)\right)\hfill \\ \hfill & +({\mu \nu }^0|f_{\mathrm{xc}}|{\rho }^{ab}({\omega }_A,{\omega }_B))+({\mu \nu }^0|k_{\mathrm{xc}}|{\rho }^a\left({\omega }_A\right)\left|{\rho }^b\left({\omega }_B\right)\right).\hfill \end{array}$$

(40)

The derivatives of $H_{\mu \nu }$ are typically independent of field frequencies, but may exhibit frequency dependence when the system is subject to a gauge transformation. This case is not considered here, but will be discussed in Subsec. 3.2.

2.5.1. The First-Order CP-HF/KS Equation

By substituting the definition (26a) into Eq. (38b) and contracting with $C_{\mu k}^{0*}$ over $\mu$ and $\nu$, Eq. (38b) can be transformed into the molecular orbital representation:

$$F_{ki}^a\left({\omega }_A\right)+{\epsilon }_{kj}^0{U}_{ji}^a\left({\omega }_A\right)-{\mathcal{T}}_{ki}^a\left({\omega }_A\right)+{\omega }_A{S}_{kj}^0{U}_{ji}^a\left({\omega }_A\right)=S_{kj}^a{\epsilon }_{ji}^0+S_{kl}^0{U}_{lj}^a\left({\omega }_A\right){\epsilon }_{ji}^0+S_{kj}^0{\epsilon }_{ji}^a\left({\omega }_A\right).$$

(41)

Since ${\epsilon }_{ij}^0$ and $S_{ij}^0$ are diagonal, we can express the above equation as

$$F_{ji}^a\left({\omega }_A\right)+{\epsilon }_j^0{U}_{ji}^a\left({\omega }_A\right)-{\mathcal{T}}_{ji}^a\left({\omega }_A\right)+{\omega }_A{U}_{ji}^a\left({\omega }_A\right)=S_{ji}^a{\epsilon }_i^0+U_{ji}^a\left({\omega }_A\right){\epsilon }_i^0+{\epsilon }_{ji}^a\left({\omega }_A\right),$$

(42)

with no summation implied over the indices of the ground-state orbital energies ${\epsilon }_i^0$. If we further assume that ${\epsilon }_{ji}^a$ is block-diagonal, i.e.,

$${\epsilon }_{ji}^a\left({\omega }_A\right)=0$$

(43)

whenever $j\in \mathrm{vir}$ and $i\in \mathrm{occ}$, or $j\in \mathrm{occ}$ and $i\in \mathrm{vir}$, then we have

$$U_{ji}^a\left({\omega }_A\right)={\displaystyle \frac{F_{ji}^a\left({\omega }_A\right)-{\mathcal{T}}_{ji}^a\left({\omega }_A\right)-S_{ji}^a{\epsilon }_i^0}{{\epsilon }_i^0-{\epsilon }_j^0-{\omega }_A}}.$$

(44)

As for the occ-occ and vir-vir blocks of $U_{ji}^a\left({\omega }_A\right)$, they can be obtained from the orthonormality condition given in Eq. (33):

$$U_{ki}^{a*}(-{\omega }_A)+U_{ik}^a\left({\omega }_A\right)=-S_{ik}^a$$

(45)

for $i,k\in \mathrm{occ}$ or $i,k\in \mathrm{vir}$. The above equation also holds in the static limit of ${\omega }_A=0$, which implies that $S_{ik}^a$ must be Hermitian. For convenience, we further assume that

$$U_{ki}^{a*}(-{\omega }_A)=U_{ik}^a\left({\omega }_A\right),$$

(46)

so that the diagonal blocks of $U_{ik}^a\left({\omega }_A\right)$ can be directly obtained as

$$U_{ik}^a\left({\omega }_A\right)=-{\displaystyle \frac{1}{2}}{S}_{ik}^a.$$

(47)

It is important to note that the orthonormality condition (33) may no longer hold for the vir-occ and occ-vir blocks if the basis set depends on the perturbation. This can potentially violate the overall orthonormality of the orbitals, and care must be taken when constructing perturbation-dependent basis sets.

Eq. (44) allows for an iterative determination of the off-diagonal blocks of $U_{ji}^a\left({\omega }_A\right)$. Alternatively, a more direct approach is to substitute Eq. (39) into Eq. (44), and rearrange the terms such that those involving $U_{ji}^a\left({\omega }_A\right)$ appear on one side and all constant terms on the other. This yields

$$\begin{array}{cc}\hfill & \left[{\mathbb{A}}_{jikl}+{\omega }_A{\delta }_{ik}{\delta }_{jl}\right]U_{lk}^a\left({\omega }_A\right)+{\mathbb{B}}_{jilk}{U}_{lk}^{a*}(-{\omega }_A)\hfill \\ \hfill & =-H_{ji}^a-n_k{\left(ji\right|kk)}^a-\left({ji}^a\right|v_{\mathrm{xc}})-n_k({ji}^0|f_{\mathrm{xc}}\left|{kk}^a\right)+{\mathcal{T}}_{ji}^a\left({\omega }_A\right)+S_{ji}^a{\epsilon }_i^0,\hfill \end{array}$$

(48)

where

$${\mathbb{A}}_{jikl}=n_k{\left(ji\right|\left|kl\right)}^0+({\epsilon }_j^0-{\epsilon }_i^0){\delta }_{ik}{\delta }_{jl},$$

(49a)

$${\mathbb{B}}_{jilk}=n_k{\left(ji\right|\left|lk\right)}^0.$$

(49b)

The two-electron integral ${\left(ij\right|\left|kl\right)}^0$ is defined as

$${\left(ij\right|\left|kl\right)}^0=({ij}^0|r_{12}^{-1}+f_{\mathrm{xc}}\left|{kl}^0\right).$$

(50)

In Eq. (48), we may take $i,k\in \mathrm{occ}$ and $j\in \mathrm{vir}$, while the index l typically runs over the full orbital space. When the basis functions are independent of perturbations, the occ-occ block of $U_{lk}^a\left({\omega }_A\right)$ vanishes (see Eq. (47)), and thus only $l\in \mathrm{vir}$ contributes. However, if the basis set is perturbation-dependent, contributions from all orbitals may be nonzero. In this case, since $U_{lk}^a\left({\omega }_A\right)$ for $l\in \mathrm{occ}$ is already known from Eq. (47), we may split the summation over l into contributions from virtual and occupied orbitals, and move the known terms with $l\in \mathrm{occ}$ to the right-hand side as constants. This leads to

$$\left[{\mathbb{A}}_{jikl}+{\omega }_A{\delta }_{ik}{\delta }_{jl}\right]U_{lk}^a\left({\omega }_A\right)+{\mathbb{B}}_{jilk}{U}_{lk}^{a*}(-{\omega }_A)=R_{ji}^a\left({\omega }_A\right),$$

(51)

where $i,k\in \mathrm{occ}$, $j,l\in \mathrm{vir}$, and the right-hand-side constants $R_{ji}^a\left({\omega }_A\right)$ are defined as

$$\begin{array}{cc}\hfill R_{ji}^a\left({\omega }_A\right)& =-H_{ji}^a-n_k{\left(ji\right|kk)}^a-\left({ji}^a\right|v_{\mathrm{xc}})-n_k({ji}^0|f_{\mathrm{xc}}\left|{kk}^a\right)+{\mathcal{T}}_{ji}^a\left({\omega }_A\right)+S_{ji}^a{\epsilon }_i^0\hfill \\ \hfill & +{\displaystyle \frac{n_k}{2}}\left[{\left(ji\right|\left|km\right)}^0{S}_{mk}^a+{\left(ji\right|\left|mk\right)}^0{S}_{mk}^{a*}\right].\hfill \end{array}$$

(52)

The index m in the last term runs over the occupied orbitals. When the basis functions are field-independent, only $-H_{ji}^a$ survives on the right-hand side. By taking the complex conjugate of Eq. (51) and replacing ${\omega }_A$ with $-{\omega }_A$, we obtain another independent equation

$$\left[{\mathbb{A}}_{jikl}^{*}-{\omega }_A{\delta }_{ik}{\delta }_{jl}\right]U_{lk}^{a*}(-{\omega }_A)+{\mathbb{B}}_{jilk}^{*}{U}_{lk}^a\left({\omega }_A\right)=R_{ji}^{a*}(-{\omega }_A).$$

(53)

Together with Eq. (51), this forms a standard linear system:

$$\left[\begin{array}{cc}{\mathbb{A}}_{jikl}+{\omega }_A{\delta }_{ik}{\delta }_{jl}& {\mathbb{B}}_{jilk}\\ {\mathbb{B}}_{jilk}^{*}& {\mathbb{A}}_{jikl}^{*}-{\omega }_A{\delta }_{ik}{\delta }_{jl}\end{array}\right]\left[\begin{array}{c}{U}_{lk}^a\left({\omega }_A\right)\\ U_{lk}^{a*}(-{\omega }_A)\end{array}\right]=\left[\begin{array}{c}{R}_{ji}^a\left({\omega }_A\right)\\ R_{ji}^{a*}(-{\omega }_A)\end{array}\right],$$

(54)

whose solution directly determines the first-order response density. The above equation bears resemblance to the Casida equation[35], but with a key distinction: while Eq. (54) is solved for the response density, the Casida equation seeks the poles of this expression, i.e., the eigenvalues for which the operator matrix becomes singular, yielding the excitation energies.

2.5.2. The Second-Order CP-HF/KS Equation

In a similar manner, multiplication of both sides of Eq. (38c) by $C_{\mu k}^{0*}$ followed by summation over $\mu$ and $\nu$ yields

$$\begin{array}{cc}\hfill & F_{ki}^{ab}({\omega }_A,{\omega }_B)+F_{kj}^a\left({\omega }_A\right)U_{ji}^b\left({\omega }_B\right)+F_{kj}^b\left({\omega }_B\right)U_{ji}^a\left({\omega }_A\right)+{\epsilon }_{kj}^0{U}_{ji}^{ab}({\omega }_A,{\omega }_B)\hfill \\ \hfill & -{\mathcal{T}}_{ki}^{ab}({\omega }_A,{\omega }_B)-{\mathcal{T}}_{kj}^a\left({\omega }_A\right)U_{ji}^b\left({\omega }_B\right)-{\mathcal{T}}_{kj}^b\left({\omega }_B\right)U_{ji}^a\left({\omega }_A\right)\hfill \\ \hfill & +{\omega }_B{S}_{kj}^a{U}_{ji}^b\left({\omega }_B\right)+{\omega }_A{S}_{kj}^b{U}_{ji}^a\left({\omega }_A\right)+\left({\omega }_A+{\omega }_B\right)S_{kj}^0{U}_{ji}^{ab}({\omega }_A,{\omega }_B)\hfill \\ \hfill & =S_{kj}^{ab}{\epsilon }_{ji}^0+S_{kl}^a{U}_{lj}^b\left({\omega }_B\right){\epsilon }_{ji}^0+S_{kj}^a{\epsilon }_{ji}^b\left({\omega }_B\right)+S_{kl}^b{U}_{lj}^a\left({\omega }_A\right){\epsilon }_{ji}^0+S_{kl}^0{U}_{lj}^{ab}({\omega }_A,{\omega }_B){\epsilon }_{ji}^0\hfill \\ \hfill & +S_{kl}^0{U}_{lj}^a\left({\omega }_A\right){\epsilon }_{ji}^b\left({\omega }_B\right)+S_{kj}^b{\epsilon }_{ji}^a\left({\omega }_A\right)+S_{kl}^0{U}_{lj}^b\left({\omega }_B\right){\epsilon }_{ji}^a\left({\omega }_A\right)+S_{kj}^0{\epsilon }_{ji}^{ab}({\omega }_A,{\omega }_B).\hfill \end{array}$$

(55)

This is equivalent to the following expression, given that both ${\epsilon }_{ij}^0$ and $S_{ij}^0$ are diagonal matrices:

$$\begin{array}{cc}\hfill & F_{ji}^{ab}({\omega }_A,{\omega }_B)+F_{jl}^a\left({\omega }_A\right)U_{li}^b\left({\omega }_B\right)+F_{jl}^b\left({\omega }_B\right)U_{li}^a\left({\omega }_A\right)+{\epsilon }_j^0{U}_{ji}^{ab}({\omega }_A,{\omega }_B)\hfill \\ \hfill & -{\mathcal{T}}_{ji}^{ab}({\omega }_A,{\omega }_B)-{\mathcal{T}}_{jl}^a\left({\omega }_A\right)U_{li}^b\left({\omega }_B\right)-{\mathcal{T}}_{jl}^b\left({\omega }_B\right)U_{li}^a\left({\omega }_A\right)\hfill \\ \hfill & +{\omega }_B{S}_{jl}^a{U}_{li}^b\left({\omega }_B\right)+{\omega }_A{S}_{jl}^b{U}_{li}^a\left({\omega }_A\right)+\left({\omega }_A+{\omega }_B\right)U_{ji}^{ab}({\omega }_A,{\omega }_B)\hfill \\ \hfill & =S_{ji}^{ab}{\epsilon }_i^0+S_{kl}^a{U}_{li}^b\left({\omega }_B\right){\epsilon }_i^0+S_{jl}^a{\epsilon }_{li}^b\left({\omega }_B\right)+S_{kl}^b{U}_{li}^a\left({\omega }_A\right){\epsilon }_i^0+U_{ji}^{ab}({\omega }_A,{\omega }_B){\epsilon }_i^0\hfill \\ \hfill & +U_{jk}^a\left({\omega }_A\right){\epsilon }_{ki}^b\left({\omega }_B\right)+S_{jl}^b{\epsilon }_{li}^a\left({\omega }_A\right)+U_{jk}^b\left({\omega }_B\right){\epsilon }_{ki}^a\left({\omega }_A\right)+{\epsilon }_{ji}^{ab}({\omega }_A,{\omega }_B).\hfill \end{array}$$

(56)

We make a similar assumption as in the first-order case, namely that

$${\epsilon }_{ji}^{ab}({\omega }_A,{\omega }_B)=0$$

(57)

for the vir-occ and occ-vir blocks. Under this assumption, we obtain

$$U_{ji}^{ab}({\omega }_A,{\omega }_B)={\displaystyle \frac{F_{ji}^{ab}({\omega }_A,{\omega }_B)+G_{ji}^{ab}({\omega }_A,{\omega }_B)-{\mathcal{W}}_{ji}^{ab}({\omega }_A,{\omega }_B)}{{\epsilon }_i^0-{\epsilon }_j^0-{\omega }_A-{\omega }_B}},$$

(58)

where $G_{ji}^{ab}$ and ${\mathcal{W}}_{ji}^{ab}$ are defined as

$$G_{ji}^{ab}({\omega }_A,{\omega }_B)=F_{jl}^a\left({\omega }_A\right)U_{li}^b\left({\omega }_B\right)+F_{jl}^b\left({\omega }_B\right)U_{li}^a\left({\omega }_A\right)-U_{jk}^a\left({\omega }_A\right){\epsilon }_{ki}^b\left({\omega }_B\right)-U_{jk}^b\left({\omega }_B\right){\epsilon }_{ki}^a\left({\omega }_A\right),$$

(59a)

$$\begin{array}{cc}\hfill & {\mathcal{W}}_{ji}^{ab}({\omega }_A,{\omega }_B)=S_{ji}^{ab}{\epsilon }_i^0+S_{jl}^a{U}_{li}^b\left({\omega }_B\right){\epsilon }_i^0+S_{jk}^a{\epsilon }_{ki}^b\left({\omega }_B\right)+S_{jl}^b{U}_{li}^a\left({\omega }_A\right){\epsilon }_i^0+S_{jk}^b{\epsilon }_{ki}^a\left({\omega }_A\right)\hfill \\ \hfill & -{\omega }_B{S}_{jl}^a{U}_{li}^b\left({\omega }_B\right)-{\omega }_A{S}_{jl}^b{U}_{li}^a\left({\omega }_A\right)+{\mathcal{T}}_{ji}^{ab}({\omega }_A,{\omega }_B)+{\mathcal{T}}_{jl}^a\left({\omega }_A\right)U_{li}^b\left({\omega }_B\right)+{\mathcal{T}}_{jl}^b\left({\omega }_B\right)U_{li}^a\left({\omega }_A\right).\hfill \end{array}$$

(59b)

The diagonal blocks of $U_{ji}^{ab}({\omega }_A,{\omega }_B)$ can also be determined by exploiting the orthonormality condition (34):

$$\begin{array}{cc}\hfill & U_{ki}^{ab*}(-{\omega }_A,-{\omega }_B)+U_{ik}^{ab}({\omega }_A,{\omega }_B)=-U_{ji}^{a*}(-{\omega }_A)U_{jk}^b\left({\omega }_B\right)-U_{ji}^{b*}(-{\omega }_B)U_{jk}^a\left({\omega }_A\right)\hfill \\ \hfill & -S_{ik}^{ab}-U_{ji}^{a*}(-{\omega }_A)S_{jk}^b-U_{ji}^{b*}(-{\omega }_B)S_{jk}^a-S_{ji}^{a*}{U}_{jk}^b\left({\omega }_B\right)-S_{ji}^{b*}{U}_{jk}^a\left({\omega }_A\right).\hfill \end{array}$$

(60)

Choosing

$$U_{ki}^{ab*}(-{\omega }_A,-{\omega }_B)=U_{ik}^{ab}({\omega }_A,{\omega }_B)$$

(61)

leads to the following working expression:

$$\begin{array}{cc}\hfill U_{ik}^{ab}({\omega }_A,{\omega }_B)& =-{\displaystyle \frac{1}{2}}\left[U_{ji}^{a*}(-{\omega }_A)U_{jk}^b\left({\omega }_B\right)+U_{ji}^{b*}(-{\omega }_B)U_{jk}^a\left({\omega }_A\right)+S_{ik}^{ab}\right.\hfill \\ \hfill & +\left.U_{ji}^{a*}(-{\omega }_A)S_{jk}^b+U_{ji}^{b*}(-{\omega }_B)S_{jk}^a+S_{ji}^{a*}{U}_{jk}^b\left({\omega }_B\right)+S_{ji}^{b*}{U}_{jk}^a\left({\omega }_A\right)\right].\hfill \end{array}$$

(62)

To facilitate a more direct solution of the off-diagonal blocks of $U_{ji}^{ab}({\omega }_A,{\omega }_B)$, we insert Eq. (40) into Eq. (58), isolating the terms involving unknowns on one side to give

$$\begin{array}{cc}\hfill & \left[{\mathbb{A}}_{jikl}+{\omega }_A{\delta }_{ik}{\delta }_{jl}\right]U_{lk}^{ab}({\omega }_A,{\omega }_B)+{\mathbb{B}}_{jilk}{U}_{lk}^{ab*}(-{\omega }_A,-{\omega }_B)\hfill \\ \hfill & =-H_{ji}^{ab}-D_{\lambda \kappa }^a\left({\omega }_A\right){\left(ji\right|\kappa \lambda )}^b-D_{\lambda \kappa }^b\left({\omega }_B\right){\left(ji\right|\kappa \lambda )}^a-n_k{\left(ji\right|kk)}^{ab}\hfill \\ \hfill & -n_k{\left(ji\right|\left|nl\right)}^0\left[U_{lk}^a\left({\omega }_A\right)U_{nk}^{b*}(-{\omega }_B)+U_{lk}^b\left({\omega }_B\right)U_{nk}^{a*}(-{\omega }_A)\right]\hfill \\ \hfill & -\left({ji}^{ab}\right|v_{\mathrm{xc}})-({ji}^a|f_{\mathrm{xc}}|{\rho }^b\left({\omega }_B\right))-({ji}^b|f_{\mathrm{xc}}\left|{\rho }^a\left({\omega }_A\right)\right)\hfill \\ \hfill & -D_{\lambda \kappa }^a\left({\omega }_A\right)({ji}^0|f_{\mathrm{xc}}|{\kappa \lambda }^b)-D_{\lambda \kappa }^b\left({\omega }_B\right)({ji}^0|f_{\mathrm{xc}}|{\kappa \lambda }^a)-n_k({ji}^0|f_{\mathrm{xc}}\left|{kk}^{ab}\right)\hfill \\ \hfill & -({ji}^0|k_{\mathrm{xc}}|{\rho }^a\left({\omega }_A\right)\left|{\rho }^b\left({\omega }_B\right)\right)-G_{ji}^{ab}({\omega }_A,{\omega }_B)+{\mathcal{W}}_{ji}^{ab}({\omega }_A,{\omega }_B).\hfill \end{array}$$

(63)

The index l on the left-hand side is summed over all orbitals, as nonzero contributions arise regardless of whether the basis functions are perturbation-dependent. Following the same procedure, we partition the index l into virtual and occupied orbitals and move all constant terms to the right-hand side to obtain

$$\left[{\mathbb{A}}_{jikl}+\left({\omega }_A+{\omega }_B\right){\delta }_{ik}{\delta }_{jl}\right]U_{lk}^{ab}({\omega }_A,{\omega }_B)+{\mathbb{B}}_{jilk}{U}_{lk}^{ab*}(-{\omega }_A,-{\omega }_B)=R_{ji}^{ab}({\omega }_A,{\omega }_B),$$

(64)

where

$$\begin{array}{cc}\hfill R_{ji}^{ab}({\omega }_A,{\omega }_B)& =-H_{ji}^{ab}-D_{\lambda \kappa }^a\left({\omega }_A\right){\left(ji\right|\kappa \lambda )}^b-D_{\lambda \kappa }^b\left({\omega }_B\right){\left(ji\right|\kappa \lambda )}^a-n_k{\left(ji\right|kk)}^{ab}\hfill \\ \hfill & -n_k{\left(ji\right|\left|nl\right)}^0\left[U_{lk}^a\left({\omega }_A\right)U_{nk}^{b*}(-{\omega }_B)+U_{lk}^b\left({\omega }_B\right)U_{nk}^{a*}(-{\omega }_A)\right]\hfill \\ \hfill & -\left({ji}^{ab}\right|v_{\mathrm{xc}})-({ji}^a|f_{\mathrm{xc}}|{\rho }^b\left({\omega }_B\right))-({ji}^b|f_{\mathrm{xc}}\left|{\rho }^a\left({\omega }_A\right)\right)\hfill \\ \hfill & -D_{\lambda \kappa }^a\left({\omega }_A\right)({ji}^0|f_{\mathrm{xc}}|{\kappa \lambda }^b)-D_{\lambda \kappa }^b\left({\omega }_B\right)({ji}^0|f_{\mathrm{xc}}|{\kappa \lambda }^a)-n_k({ji}^0|f_{\mathrm{xc}}\left|{kk}^{ab}\right)\hfill \\ \hfill & -({ji}^0|k_{\mathrm{xc}}|{\rho }^a\left({\omega }_A\right)\left|{\rho }^b\left({\omega }_B\right)\right)-G_{ji}^{ab}({\omega }_A,{\omega }_B)+{\mathcal{W}}_{ji}^{ab}({\omega }_A,{\omega }_B)\hfill \\ \hfill & -n_k\left[{\left(ji\right|\left|km\right)}^0{U}_{mk}^{ab}({\omega }_A,{\omega }_B)+{\left(ji\right|\left|mk\right)}^0{U}_{mk}^{ab*}(-{\omega }_A,-{\omega }_B)\right].\hfill \end{array}$$

(65)

In the two equations above, we take $i,k,m\in \mathrm{occ}$ and $j,l\in \mathrm{vir}$. The index n may be restricted to virtual orbitals if the basis functions are field-independent; otherwise, it must run over all orbitals. Taking the complex conjugate of Eq. (64) and replacing $({\omega }_A,{\omega }_B)$ with $(-{\omega }_A,-{\omega }_B)$ yields

$$\left[{\mathbb{A}}_{jikl}^{*}-\left({\omega }_A+{\omega }_B\right){\delta }_{ik}{\delta }_{jl}\right]U_{lk}^{ab*}(-{\omega }_A,-{\omega }_B)+{\mathbb{B}}_{jilk}^{*}{U}_{lk}^{ab}({\omega }_A,{\omega }_B)=R_{ji}^{ab*}(-{\omega }_A,-{\omega }_B).$$

(66)

Combined with Eq. (64), this leads to a standard linear system:

$$\left[\begin{array}{cc}{\mathbb{A}}_{jikl}+\left({\omega }_A+{\omega }_B\right){\delta }_{ik}{\delta }_{jl}& {\mathbb{B}}_{jilk}\\ {\mathbb{B}}_{jilk}^{*}& {\mathbb{A}}_{jikl}^{*}-\left({\omega }_A+{\omega }_B\right){\delta }_{ik}{\delta }_{jl}\end{array}\right]\left[\begin{array}{c}{U}_{lk}^{ab}({\omega }_A,{\omega }_B)\\ U_{lk}^{ab*}(-{\omega }_A,-{\omega }_B)\end{array}\right]=\left[\begin{array}{c}{R}_{ji}^{ab}({\omega }_A,{\omega }_B)\\ R_{ji}^{ab*}(-{\omega }_A,-{\omega }_B)\end{array}\right],$$

(67)

which constitutes the working second-order CP-HF/KS equation. It depends on the solutions of the first-order equation (54), and together, the first- and second-order solutions determine the second-order response density.

2.6. Quasi-Energy Response: Working Form

Using the expressions derived in Subsecs. 2.3–2.5, the quasi-energy derivatives from Subsec. 2.2 can be simplified to yield the following explicit working formulas that are suitable for practical evaluation (see Appendix B for details).

The First Order:

$$\begin{array}{cc}\hfill W^a\left({\omega }_A\right)& =-n_i{S}_{ii}^a\left[{\epsilon }_i^0+{\displaystyle \frac{{\omega }_A}{2}}\right]+D_{\nu \mu }^0\left[H_{\mu \nu }^a-{\mathcal{T}}_{\mu \nu }^a\left({\omega }_A\right)\right]\hfill \\ \hfill & +{\displaystyle \frac{1}{2}}{D}_{\nu \mu }^0{D}_{\lambda \kappa }^0{\left(\mu \nu \right|\kappa \lambda )}^a+D_{\nu \mu }^0\left({\mu \nu }^a\right|v_{\mathrm{xc}}).\hfill \end{array}$$

(68)

The Second Order:

$$\begin{array}{cc}\hfill & W^{ab}({\omega }_A,{\omega }_B)\hfill \\ \hfill & =-{\displaystyle \frac{n_i}{2}}{U}_{ji}^a\left({\omega }_A\right)S_{ji}^{b*}\left({\epsilon }_i^0-{\omega }_A+{\omega }_B\right)-{\displaystyle \frac{n_i}{2}}{U}_{ji}^b\left({\omega }_B\right)S_{ji}^{a*}\left({\epsilon }_i^0+{\omega }_A-{\omega }_B\right)\hfill \\ \hfill & -{\displaystyle \frac{n_i}{2}}\left[U_{ji}^{a*}(-{\omega }_A)S_{ji}^b+U_{ji}^{b*}(-{\omega }_B)S_{ji}^a\right]\left({\epsilon }_i^0+{\omega }_A+{\omega }_B\right)\hfill \\ \hfill & -{\displaystyle \frac{n_i}{2}}{S}_{ii}^{ab}\left(2{\epsilon }_i^0+{\omega }_A+{\omega }_B\right)+{\displaystyle \frac{n_i}{4}}\left[S_{ki}^a{S}_{ki}^{b*}+S_{ki}^b{S}_{ki}^{a*}\right]\left({\epsilon }_k^0+3{\epsilon }_i^0\right)\hfill \\ \hfill & +{\displaystyle \frac{n_i}{8}}\left[\left({\omega }_A+3{\omega }_B\right)S_{ki}^a{S}_{ki}^{b*}+\left({\omega }_B+3{\omega }_A\right)S_{ki}^b{S}_{ki}^{a*}\right]\hfill \\ \hfill & -{\displaystyle \frac{n_i}{4}}\left\{S_{ki}^a\left[F_{ki}^{b*}(-{\omega }_B)-{\mathcal{T}}_{ki}^{b*}(-{\omega }_B)\right]+S_{ki}^{a*}\left[F_{ki}^b\left({\omega }_B\right)-{\mathcal{T}}_{ki}^b\left({\omega }_B\right)\right]\right.\hfill \\ \hfill & +\left.S_{ki}^b\left[F_{ki}^{a\phantom{b}*}(-{\omega }_A)-{\mathcal{T}}_{ki}^{a*}(-{\omega }_A)\right]+S_{ki}^{b*}\left[F_{ki}^{a\phantom{b}}\left({\omega }_A\right)-{\mathcal{T}}_{ki}^a\left({\omega }_A\right)\right]\right\}\hfill \\ \hfill & +{\displaystyle \frac{1}{2}}\left\{D_{\nu \mu }^a\left({\omega }_A\right)\left[H_{\mu \nu }^b-{\mathcal{T}}_{\mu \nu }^b\left({\omega }_B\right)\right]+D_{\nu \mu }^b\left({\omega }_B\right)\left[H_{\mu \nu }^a-{\mathcal{T}}_{\mu \nu }^a\left({\omega }_A\right)\right]\right.\hfill \\ \hfill & +D_{\nu \mu }^a\left({\omega }_A\right)D_{\lambda \kappa }^0{\left(\mu \nu \right|\kappa \lambda )}^b+D_{\nu \mu }^b\left({\omega }_B\right)D_{\lambda \kappa }^0{\left(\mu \nu \right|\kappa \lambda )}^a+D_{\nu \mu }^a\left({\omega }_A\right)\left({\mu \nu }^b\right|v_{\mathrm{xc}})\hfill \\ \hfill & +\left.D_{\nu \mu }^b\left({\omega }_B\right)\left({\mu \nu }^a\right|v_{\mathrm{xc}})+D_{\nu \mu }^0({\mu \nu }^b|f_{\mathrm{xc}}|{\rho }^a\left({\omega }_A\right))+D_{\nu \mu }^0({\mu \nu }^a|f_{\mathrm{xc}}\left|{\rho }^b\left({\omega }_B\right)\right)\right\}\hfill \\ \hfill & +D_{\nu \mu }^0\left[H_{\mu \nu }^{ab}-{\mathcal{T}}_{\mu \nu }^{ab}({\omega }_A,{\omega }_B)\right]+{\displaystyle \frac{1}{2}}{D}_{\nu \mu }^0{D}_{\lambda \kappa }^0{\left(\mu \nu \right|\kappa \lambda )}^{ab}+D_{\nu \mu }^0\left({\mu \nu }^{ab}\right|v_{\mathrm{xc}}),\hfill \end{array}$$

(69)

where $i,k\in \mathrm{occ}$ and $j\in \mathrm{vir}$.

The Third Order:

$$\begin{array}{cc}\hfill & W^{abc}({\omega }_A,{\omega }_B,{\omega }_C)\hfill \\ \hfill & ={\displaystyle \frac{n_i}{2}}\left[U_{ji}^{ab}({\omega }_A,{\omega }_B)U_{ji}^{c*}(-{\omega }_C)-U_{ji}^{ab*}(-{\omega }_A,-{\omega }_B)U_{ji}^c\left({\omega }_C\right)\right.\hfill \\ \hfill & +U_{ji}^{ac}({\omega }_A,{\omega }_C)U_{ji}^{b*}(-{\omega }_B)-U_{ji}^{ac*}(-{\omega }_A,-{\omega }_C)U_{ji}^b\left({\omega }_B\right)\hfill \\ \hfill & +\left.U_{ji}^{bc}({\omega }_B,{\omega }_C)U_{ji}^{a*}(-{\omega }_A)-U_{ji}^{bc*}(-{\omega }_B,-{\omega }_C)U_{ji}^a\left({\omega }_A\right)\right]\left({\omega }_A+{\omega }_B+{\omega }_C\right)\hfill \\ \hfill & +n_i\left[F_{jl}^a\left({\omega }_A\right)-{\mathcal{T}}_{jl}^a\left({\omega }_A\right)\right]\left[U_{li}^b\left({\omega }_B\right)U_{ji}^{c*}(-{\omega }_C)+U_{li}^c\left({\omega }_C\right)U_{ji}^{b*}(-{\omega }_B)\right]\hfill \\ \hfill & +n_i\left[F_{jl}^b\left({\omega }_B\right)-{\mathcal{T}}_{jl}^b\left({\omega }_B\right)\right]\left[U_{li}^a\left({\omega }_A\right)U_{ji}^{c*}(-{\omega }_C)+U_{li}^c\left({\omega }_C\right)U_{ji}^{a*}(-{\omega }_A)\right]\hfill \\ \hfill & +n_i\left[F_{jl}^c\left({\omega }_C\right)-{\mathcal{T}}_{jl}^c\left({\omega }_C\right)\right]\left[U_{li}^a\left({\omega }_A\right)U_{ji}^{b*}(-{\omega }_B)+U_{li}^b\left({\omega }_B\right)U_{ji}^{a*}(-{\omega }_A)\right]\hfill \\ \hfill & -n_i\left[F_{ki}^{a\phantom{b}}\left({\omega }_A\right)-{\mathcal{T}}_{ki}^a\left({\omega }_A\right)\right]\left[U_{ji}^{b*}(-{\omega }_B)S_{jk}^c+U_{ji}^{b*}(-{\omega }_B)U_{jk}^c\left({\omega }_C\right)\right.\hfill \\ \hfill & +\left.U_{ji}^{c*}(-{\omega }_C)S_{jk}^b+S_{ik}^{bc}+S_{ji}^{b*}{U}_{jk}^c\left({\omega }_C\right)+U_{ji}^{c*}(-{\omega }_C)U_{jk}^b\left({\omega }_B\right)+S_{ji}^{c*}{U}_{jk}^b\left({\omega }_B\right)\right]\hfill \\ \hfill & -n_i\left[F_{ki}^b\left({\omega }_B\right)-{\mathcal{T}}_{ki}^b\left({\omega }_B\right)\right]\left[U_{ji}^{a*}(-{\omega }_A)S_{jk}^c+U_{ji}^{a*}(-{\omega }_A)U_{jk}^c\left({\omega }_C\right)\right.\hfill \\ \hfill & +\left.U_{ji}^{c*}(-{\omega }_C)S_{jk}^a+S_{ik}^{ac}+S_{ji}^{a*}{U}_{jk}^c\left({\omega }_C\right)+U_{ji}^{c*}(-{\omega }_C)U_{jk}^a\left({\omega }_A\right)+S_{ji}^{c*}{U}_{jk}^a\left({\omega }_A\right)\right]\hfill \\ \hfill & -n_i\left[F_{ki}^{c\phantom{b}}\left({\omega }_C\right)-{\mathcal{T}}_{ki}^c\left({\omega }_C\right)\right]\left[U_{ji}^{a*}(-{\omega }_A)S_{jk}^b+U_{ji}^{a*}(-{\omega }_A)U_{jk}^b\left({\omega }_B\right)\right.\hfill \\ \hfill & +\left.U_{ji}^{b*}(-{\omega }_B)S_{jk}^a+S_{ik}^{ab}+S_{ji}^{a*}{U}_{jk}^b\left({\omega }_B\right)+U_{ji}^{b*}(-{\omega }_B)U_{jk}^a\left({\omega }_A\right)+S_{ji}^{b*}{U}_{jk}^a\left({\omega }_A\right)\right]\hfill \\ \hfill & +{\mathcal{S}}^{abc}({\omega }_A,{\omega }_B,{\omega }_C)+{\mathcal{V}}^{abc}({\omega }_A,{\omega }_B,{\omega }_C),\hfill \end{array}$$

(70)

where $i,k\in \mathrm{occ}$, and $j\in \mathrm{vir}$ in the first bracketed term. The quantities ${\mathcal{S}}^{abc}({\omega }_A,{\omega }_B,{\omega }_C)$ and ${\mathcal{V}}^{abc}({\omega }_A,{\omega }_B,{\omega }_C)$ are defined as (with $i,k$ restricted to occupied orbitals)

$$\begin{array}{cc}\hfill & {\mathcal{S}}^{abc}({\omega }_A,{\omega }_B,{\omega }_C)\hfill \\ \hfill & =-{\displaystyle \frac{n_i}{2}}\left[U_{ki}^{ab}({\omega }_A,{\omega }_B)S_{ki}^{c*}\left({\epsilon }_k^0+{\epsilon }_i^0+{\omega }_A+{\omega }_B\right)+U_{ki}^{ab*}(-{\omega }_A,-{\omega }_B)S_{ki}^c\left({\epsilon }_k^0+{\epsilon }_i^0+{\omega }_C\right)\right.\hfill \\ \hfill & +U_{ki}^{ac}({\omega }_A,{\omega }_C)S_{ki}^{b*}\left({\epsilon }_k^0+{\epsilon }_i^0+{\omega }_A+{\omega }_C\right)+U_{ki}^{ac*}(-{\omega }_A,-{\omega }_C)S_{ki}^b\left({\epsilon }_k^0+{\epsilon }_i^0+{\omega }_B\right)\hfill \\ \hfill & +\left.U_{ki}^{bc}({\omega }_B,{\omega }_C)S_{ki}^{a*}\left({\epsilon }_k^0+{\epsilon }_i^0+{\omega }_B+{\omega }_C\right)+U_{ki}^{bc*}(-{\omega }_B,-{\omega }_C)S_{ki}^a\left({\epsilon }_k^0+{\epsilon }_i^0+{\omega }_A\right)\right]\hfill \\ \hfill & +{\displaystyle \frac{n_i}{4}}\left({\omega }_A+{\omega }_B+{\omega }_C\right)\left[U_{ki}^{ab}({\omega }_A,{\omega }_B)S_{ki}^{c*}+U_{ki}^{ac}({\omega }_A,{\omega }_C)S_{ki}^{b*}+U_{ki}^{bc}({\omega }_B,{\omega }_C)S_{ki}^{a*}\right.\hfill \\ \hfill & +\left.U_{ki}^{ab*}(-{\omega }_A,-{\omega }_B)S_{ki}^c+U_{ki}^{ac*}(-{\omega }_A,-{\omega }_C)S_{ki}^b+U_{ki}^{bc*}(-{\omega }_B,-{\omega }_C)S_{ki}^a\right]\hfill \\ \hfill & -{\displaystyle \frac{n_i}{2}}\left\{\left[U_{ji}^{b*}(-{\omega }_B)S_{jl}^c+U_{ji}^{c*}(-{\omega }_C)S_{jl}^b+S_{li}^{bc*}\right]U_{li}^a\left({\omega }_A\right)\left(2{\epsilon }_i^0-{\omega }_A+{\omega }_B+{\omega }_C\right)\right.\hfill \\ \hfill & +\left[U_{ji}^{a*}(-{\omega }_A)S_{jl}^c+U_{ji}^{c*}(-{\omega }_C)S_{jl}^a+S_{li}^{ac*}\right]U_{li}^b\left({\omega }_B\right)\left(2{\epsilon }_i^0+{\omega }_A-{\omega }_B+{\omega }_C\right)\hfill \\ \hfill & +\left.\left[U_{ji}^{a*}(-{\omega }_A)S_{jl}^b+U_{ji}^{b*}(-{\omega }_B)S_{jl}^a+S_{li}^{ab*}\right]U_{li}^c\left({\omega }_C\right)\left(2{\epsilon }_i^0+{\omega }_A+{\omega }_B-{\omega }_C\right)\right\}\hfill \\ \hfill & -{\displaystyle \frac{n_i}{2}}\left[\left({\omega }_C-{\omega }_A-{\omega }_B\right)U_{ji}^{ab}({\omega }_A,{\omega }_B)S_{ji}^{c*}+\left({\omega }_B-{\omega }_A-{\omega }_C\right)U_{ji}^{ac}({\omega }_A,{\omega }_C)S_{ji}^{b*}\right.\hfill \\ \hfill & +\left.\left({\omega }_A-{\omega }_B-{\omega }_C\right)U_{ji}^{bc}({\omega }_B,{\omega }_C)S_{ji}^{a*}\right]-{\displaystyle \frac{n_i}{2}}\left[U_{ji}^{ab*}(-{\omega }_A,-{\omega }_B)S_{ji}^c\right.\hfill \\ \hfill & +\left.U_{ji}^{ac*}(-{\omega }_A,-{\omega }_C)S_{ji}^b+U_{ji}^{bc*}(-{\omega }_B,-{\omega }_C)S_{ji}^a\right]\left({\omega }_A+{\omega }_B+{\omega }_C\right)\hfill \\ \hfill & -{\displaystyle \frac{n_i}{2}}\left[U_{ji}^{a*}(-{\omega }_A)S_{ji}^{bc}+U_{ji}^{b*}(-{\omega }_B)S_{ji}^{ac}+U_{ji}^{c*}(-{\omega }_C)S_{ji}^{ab}+S_{ii}^{abc}\right]\left(2{\epsilon }_i^0+{\omega }_A+{\omega }_B+{\omega }_C\right),\hfill \end{array}$$

(71)

and

$$\begin{array}{cc}\hfill & {\mathcal{V}}^{abc}({\omega }_A,{\omega }_B,{\omega }_C)\hfill \\ \hfill & =D_{\nu \mu }^a\left({\omega }_A\right)\left[H_{\mu \nu }^{bc}-{\mathcal{T}}_{\mu \nu }^{bc}({\omega }_B,{\omega }_C)+\left({\mu \nu }^{bc}\right|v_{\mathrm{xc}})+({\mu \nu }^b|f_{\mathrm{xc}}|{\rho }^c\left({\omega }_C\right))+({\mu \nu }^c|f_{\mathrm{xc}}\left|{\rho }^b\left({\omega }_B\right)\right)\right]\hfill \\ \hfill & +D_{\nu \mu }^b\left({\omega }_B\right)\left[H_{\mu \nu }^{ac}-{\mathcal{T}}_{\mu \nu }^{ac}({\omega }_A,{\omega }_C)+\left({\mu \nu }^{ac}\right|v_{\mathrm{xc}})+({\mu \nu }^a|f_{\mathrm{xc}}|{\rho }^c\left({\omega }_C\right))+({\mu \nu }^c|f_{\mathrm{xc}}\left|{\rho }^a\left({\omega }_A\right)\right)\right]\hfill \\ \hfill & +D_{\nu \mu }^c\left({\omega }_C\right)\left[H_{\mu \nu }^{ab}-{\mathcal{T}}_{\mu \nu }^{ab}({\omega }_A,{\omega }_B)+\left({\mu \nu }^{bc}\right|v_{\mathrm{xc}})+({\mu \nu }^a|f_{\mathrm{xc}}|{\rho }^b\left({\omega }_B\right))+({\mu \nu }^b|f_{\mathrm{xc}}\left|{\rho }^a\left({\omega }_A\right)\right)\right]\hfill \\ \hfill & +D_{\nu \mu }^0\left[H_{\mu \nu }^{abc}-{\mathcal{T}}_{\mu \nu }^{abc}({\omega }_A,{\omega }_B,{\omega }_C)+{\displaystyle \frac{1}{2}}{D}_{\lambda \kappa }^0{\left(\mu \nu \right|\kappa \lambda )}^{abc}+\left({\mu \nu }^{abc}\right|v_{\mathrm{xc}})\right.\hfill \\ \hfill & +D_{\lambda \kappa }^a\left({\omega }_A\right){\left(\mu \nu \right|\kappa \lambda )}^{bc}+D_{\lambda \kappa }^b\left({\omega }_B\right){\left(\mu \nu \right|\kappa \lambda )}^{ac}+D_{\lambda \kappa }^c\left({\omega }_C\right){\left(\mu \nu \right|\kappa \lambda )}^{ab}\hfill \\ \hfill & +\left.({\mu \nu }^{bc}|f_{\mathrm{xc}}|{\rho }^a\left({\omega }_A\right))+({\mu \nu }^{ac}|f_{\mathrm{xc}}|{\rho }^b\left({\omega }_B\right))+({\mu \nu }^{ab}|f_{\mathrm{xc}}\left|{\rho }^c\left({\omega }_C\right)\right)\right]\hfill \\ \hfill & +D_{\nu \mu }^a\left({\omega }_A\right)D_{\lambda \kappa }^b\left({\omega }_B\right){\left(\mu \nu \right|\kappa \lambda )}^c+D_{\nu \mu }^a\left({\omega }_A\right)D_{\lambda \kappa }^c\left({\omega }_C\right){\left(\mu \nu \right|\kappa \lambda )}^b\hfill \\ \hfill & +D_{\nu \mu }^b\left({\omega }_B\right)D_{\lambda \kappa }^c\left({\omega }_C\right){\left(\mu \nu \right|\kappa \lambda )}^a+({\rho }^a\left({\omega }_A\right)|k_{\mathrm{xc}}|{\rho }^b\left({\omega }_B\right)\left|{\rho }^c\left({\omega }_C\right)\right),\hfill \end{array}$$

(72)

respectively. The third-order derivative of the xc functional $\left({\rho }^a\left({\omega }_A\right)\right|k_{\mathrm{xc}}\left|{\rho }^b\left({\omega }_B\right)\right|{\rho }^c\left({\omega }_C\right))$ is computed numerically but in different forms for various types of functionals, as detailed in Appendix C.

If a perturbation-independent basis set is used, Eqs. (68)–(70) simplify to

$$W^a=D_{\nu \mu }^0{H}_{\mu \nu }^a,$$

(73a)

$$W^{ab}({\omega }_A,{\omega }_B)={\displaystyle \frac{1}{2}}\left[D_{\nu \mu }^a\left({\omega }_A\right)H_{\mu \nu }^b+D_{\nu \mu }^b\left({\omega }_B\right)H_{\mu \nu }^a\right]+D_{\nu \mu }^0{H}_{\mu \nu }^{ab},$$

(73b)

$$\begin{array}{cc}\hfill & W^{abc}({\omega }_A,{\omega }_B,{\omega }_C)\hfill \\ \hfill & ={\displaystyle \frac{n_i}{2}}\left[U_{ji}^{ab}({\omega }_A,{\omega }_B)U_{ji}^{c*}(-{\omega }_C)-U_{ji}^{ab*}(-{\omega }_A,-{\omega }_B)U_{ji}^c\left({\omega }_C\right)\right.\hfill \\ \hfill & +U_{ji}^{ac}({\omega }_A,{\omega }_C)U_{ji}^{b*}(-{\omega }_B)-U_{ji}^{ac*}(-{\omega }_A,-{\omega }_C)U_{ji}^b\left({\omega }_B\right)\hfill \\ \hfill & +\left.U_{ji}^{bc}({\omega }_B,{\omega }_C)U_{ji}^{a*}(-{\omega }_A)-U_{ji}^{bc*}(-{\omega }_B,-{\omega }_C)U_{ji}^a\left({\omega }_A\right)\right]\left({\omega }_A+{\omega }_B+{\omega }_C\right)\hfill \\ \hfill & +n_i{F}_{jl}^a\left({\omega }_A\right)\left[U_{li}^b\left({\omega }_B\right)U_{ji}^{c*}(-{\omega }_C)+U_{li}^c\left({\omega }_C\right)U_{ji}^{b*}(-{\omega }_B)\right]\hfill \\ \hfill & +n_i{F}_{jl}^b\left({\omega }_B\right)\left[U_{li}^a\left({\omega }_A\right)U_{ji}^{c*}(-{\omega }_C)+U_{li}^c\left({\omega }_C\right)U_{ji}^{a*}(-{\omega }_A)\right]\hfill \\ \hfill & +n_i{F}_{jl}^c\left({\omega }_C\right)\left[U_{li}^a\left({\omega }_A\right)U_{ji}^{b*}(-{\omega }_B)+U_{li}^b\left({\omega }_B\right)U_{ji}^{a*}(-{\omega }_A)\right]\hfill \\ \hfill & -n_i{F}_{ki}^{a\phantom{b}}\left({\omega }_A\right)\left[U_{ji}^{b*}(-{\omega }_B)U_{jk}^c\left({\omega }_C\right)+U_{ji}^{c*}(-{\omega }_C)U_{jk}^b\left({\omega }_B\right)\right]\hfill \\ \hfill & -n_i{F}_{ki}^b\left({\omega }_B\right)\left[U_{ji}^{a*}(-{\omega }_A)U_{jk}^c\left({\omega }_C\right)+U_{ji}^{c*}(-{\omega }_C)U_{jk}^a\left({\omega }_A\right)\right]\hfill \\ \hfill & -n_i{F}_{ki}^{c\phantom{b}}\left({\omega }_C\right)\left[U_{ji}^{a*}(-{\omega }_A)U_{jk}^b\left({\omega }_B\right)+U_{ji}^{b*}(-{\omega }_B)U_{jk}^a\left({\omega }_A\right)\right]\hfill \\ \hfill & +D_{\nu \mu }^a\left({\omega }_A\right)H_{\mu \nu }^{bc}+D_{\nu \mu }^b\left({\omega }_B\right)H_{\mu \nu }^{ac}+D_{\nu \mu }^c\left({\omega }_C\right)H_{\mu \nu }^{ab}+D_{\nu \mu }^0{H}_{\mu \nu }^{abc}\hfill \\ \hfill & +\left({\rho }^a\left({\omega }_A\right)\right|k_{\mathrm{xc}}\left|{\rho }^b\left({\omega }_B\right)\right|{\rho }^c\left({\omega }_C\right)).\hfill \end{array}$$

(73c)

3. Applications

The general formalism developed in the previous section can be directly applied to two commonly encountered types of perturbations: electric and magnetic fields. Electric fields give rise to (hyper)polarizabilities, while magnetic fields induce (hyper)magnetizabilities. In this section, we examine the explicit forms of these perturbations, which is essential for practical calculations and program implementation.

3.1. Polarizabilities

For a system with zero net charge, if the external electric field can be approximated as homogeneous, its interaction with the system may be described by the dipole approximation:

$${\mathcal{H}}^{\prime }=-{\widehat{\mathsf{\mu }}}_{\mathrm{e}}·\mathcal{E},$$

(74)

where $\mathcal{E}$ is the external electric field and ${\widehat{\mathsf{\mu }}}_{\mathrm{e}}$ is the electric dipole operator:

$${\widehat{\mathsf{\mu }}}_{\mathrm{e}}=\sum _I{Z}_I{\mathit{r}}_I-\sum _i{\mathit{r}}_i,$$

(75)

which includes contributions from both the nuclei and electrons. If the nuclear positions are considered fixed and unaffected by the external field, the nuclear contribution becomes a constant. As such, only the electronic part of the dipole operator needs to be taken into account in the SCF procedure.

To characterize the system’s response to such a perturbation, we define the electric dipole moment, polarizability, and first hyperpolarizability as follows:

$$\begin{array}{cc}\hfill {\mu }_{\mathrm{e}}^a\left({\omega }_A\right)& =-W^a\left({\omega }_A\right)\hfill \\ \hfill & =-D_{\nu \mu }^0{H}_{\mu \nu }^a,\hfill \end{array}$$

(76a)

$$\begin{array}{cc}\hfill {\alpha }^{ab}({\omega }_A,{\omega }_B)& =-W^{ab}({\omega }_A,{\omega }_B)\hfill \\ \hfill & =-{\displaystyle \frac{1}{2}}\left[D_{\nu \mu }^a\left({\omega }_A\right)H_{\mu \nu }^b+D_{\nu \mu }^b\left({\omega }_B\right)H_{\mu \nu }^a\right],\hfill \end{array}$$

(76b)

$$\begin{array}{cc}\hfill {\beta }^{abc}({\omega }_A,{\omega }_B,{\omega }_C)& =-W^{abc}({\omega }_A,{\omega }_B,{\omega }_C)\hfill \\ \hfill & =-{\displaystyle \frac{n_i}{2}}\left[U_{ji}^{ab}({\omega }_A,{\omega }_B)U_{ji}^{c*}(-{\omega }_C)-U_{ji}^{ab*}(-{\omega }_A,-{\omega }_B)U_{ji}^c\left({\omega }_C\right)\right.\hfill \\ \hfill & +U_{ji}^{ac}({\omega }_A,{\omega }_C)U_{ji}^{b*}(-{\omega }_B)-U_{ji}^{ac*}(-{\omega }_A,-{\omega }_C)U_{ji}^b\left({\omega }_B\right)\hfill \\ \hfill & +\left.U_{ji}^{bc}({\omega }_B,{\omega }_C)U_{ji}^{a*}(-{\omega }_A)-U_{ji}^{bc*}(-{\omega }_B,-{\omega }_C)U_{ji}^a\left({\omega }_A\right)\right]\hfill \\ \hfill & \times \left({\omega }_A+{\omega }_B+{\omega }_C\right)-({\rho }^a\left({\omega }_A\right)|k_{\mathrm{xc}}|{\rho }^b\left({\omega }_B\right)\left|{\rho }^c\left({\omega }_C\right)\right)\hfill \\ \hfill & -n_i{F}_{jl}^a\left({\omega }_A\right)\left[U_{li}^b\left({\omega }_B\right)U_{ji}^{c*}(-{\omega }_C)+U_{li}^c\left({\omega }_C\right)U_{ji}^{b*}(-{\omega }_B)\right]\hfill \\ \hfill & -n_i{F}_{jl}^b\left({\omega }_B\right)\left[U_{li}^a\left({\omega }_A\right)U_{ji}^{c*}(-{\omega }_C)+U_{li}^c\left({\omega }_C\right)U_{ji}^{a*}(-{\omega }_A)\right]\hfill \\ \hfill & -n_i{F}_{jl}^c\left({\omega }_C\right)\left[U_{li}^a\left({\omega }_A\right)U_{ji}^{b*}(-{\omega }_B)+U_{li}^b\left({\omega }_B\right)U_{ji}^{a*}(-{\omega }_A)\right]\hfill \\ \hfill & +n_i{F}_{ki}^{a\phantom{b}}\left({\omega }_A\right)\left[U_{ji}^{b*}(-{\omega }_B)U_{jk}^c\left({\omega }_C\right)+U_{ji}^{c*}(-{\omega }_C)U_{jk}^b\left({\omega }_B\right)\right]\hfill \\ \hfill & +n_i{F}_{ki}^b\left({\omega }_B\right)\left[U_{ji}^{a*}(-{\omega }_A)U_{jk}^c\left({\omega }_C\right)+U_{ji}^{c*}(-{\omega }_C)U_{jk}^a\left({\omega }_A\right)\right]\hfill \\ \hfill & +n_i{F}_{ki}^{c\phantom{b}}\left({\omega }_C\right)\left[U_{ji}^{a*}(-{\omega }_A)U_{jk}^b\left({\omega }_B\right)+U_{ji}^{b*}(-{\omega }_B)U_{jk}^a\left({\omega }_A\right)\right],\hfill \end{array}$$

(76c)

where $i,k\in \mathrm{occ}$ and $j,l\in \mathrm{vir}$. In practice, computation of such electric response properties can be performed with standard (field-independent) basis sets, for which the derivatives of the basis functions with respect to external fields vanish. The core Hamiltonian has nonzero field derivatives only at first order:

$$H_{\mu \nu }^a=\left(\mu \right|r^a\left|\nu \right).$$

(77)

Derived from the quasi-energy formalism with a field-independent basis set, the first-order response (76a) coincides exactly with that obtained from the expectation value of the dipole operator, and the second-order form (76b) is also largely consistent. The only difference at second order lies in the quasi-energy formalism’s explicit incorporation of the symmetry with respect to the permutation of differentiation orders. However, discrepancies begin to appear at third order: the quasi-energy response (76c) includes additional terms involving the second-order response density and exhibits an explicit dependence on field frequencies. This indicates a breakdown of the so-called $(2n+1)$ rule for time-dependent perturbations. Nevertheless, it is worth noting that when the applied frequencies satisfy the equation ${\omega }_A+{\omega }_B+{\omega }_C=0$, the additional terms in Eq. (76c) cancel out. In other words, nonlinear optical phenomena such as SHG can still be correctly described. Eqs. (76a)–(76c) can thus reduce to the special cases presented in the works of Rice et al. [21] and of Karna and Dupuis[22]. However, their derivations are not sufficiently direct—for instance, interpreting SHG by invoking an artificial external field oscillating at $2\omega$. In contrast, our formulation is grounded in a more general framework that allows for arbitrary combinations of time-harmonic perturbations, leading to a clearer physical interpretation of the resulting response properties.

3.2. Magnetizabilities

When a magnetic field is applied, the system can also respond to it. If the source of the field is sufficiently distant, the magnetic field $\mathcal{B}$ can be approximated as uniform. Under this assumption, the nonrelativistic magnetic interaction takes the form

$$\begin{array}{cc}\hfill {\widehat{\mathcal{H}}}^{\prime }& =\sum _i\left[-\mathrm{i}\mathit{A}\left({\mathit{r}}_i\right)·{\nabla }_i+{\displaystyle \frac{1}{2}}{A}^2\left({\mathit{r}}_i\right)+{\displaystyle \frac{1}{2}}{\mathsf{\sigma }}_i·\mathcal{B}\right]\hfill \\ \hfill & =\sum _i\left[-{\displaystyle \frac{\mathrm{i}}{2}}\mathcal{B}\times {\mathit{r}}_i·{\nabla }_i+{\displaystyle \frac{1}{8}}{\left(\mathcal{B}\times {\mathit{r}}_i\right)}^2+{\displaystyle \frac{1}{2}}{\mathsf{\sigma }}_i·\mathcal{B}\right]\hfill \\ \hfill & =-{\widehat{\mathsf{\mu }}}_{\mathrm{m}}·\mathcal{B}+{\displaystyle \frac{1}{8}}\sum _i\left[{\mathcal{B}}^2{r}_i^2-{\left(\mathcal{B}·{\mathit{r}}_i\right)}^2\right],\hfill \end{array}$$

(78)

where $\mathit{A}$ is the vector potential, defined by $\nabla \times \mathit{A}=\mathcal{B}$, and ${\mathsf{\sigma }}_i\equiv \mathsf{\sigma }$ is the Pauli vector. The first term in Eq. (78) corresponds to the paramagnetic contribution of the electrons, which is precisely the classical magnetic dipole interaction, and the electronic magnetic moment operator ${\widehat{\mathsf{\mu }}}_{\mathrm{m}}$ arises from both orbital and spin angular momenta:

$${\widehat{\mathsf{\mu }}}_{\mathrm{m}}=-{\displaystyle \frac{1}{2}}\sum _i\left({\widehat{\mathit{l}}}_i+{\mathcal{g}}_{\mathrm{s}}{\widehat{\mathit{s}}}_i\right)\approx -{\displaystyle \frac{1}{2}}\sum _i\left(-\mathrm{i}{\mathit{r}}_i\times {\nabla }_i+{\mathsf{\sigma }}_i\right),$$

(79)

where the leading minus sign results from the negative charge of the electron, the factor of ${\displaystyle \frac{1}{2}}$ corresponds to the Bohr magneton, and ${\mathcal{g}}_{\mathrm{s}}\approx 2$ is the Landé g-factor for the electron spin, introduced to correct the gyromagnetic ratio. The second term in Eq. (78) represents the diamagnetic contribution of the electrons, which reflects the system’s tendency to oppose the applied magnetic field.

Strictly speaking, the magnetic interaction operator should also include paramagnetic and diamagnetic contributions from the nuclei. Even under the Born-Oppenheimer approximation, atomic nuclei possess spin magnetic moments and exhibit diamagnetic responses to an external magnetic field. However, these terms involve mass in the denominator, and since the nuclear mass is at least three orders of magnitude greater than the electronic mass, the nuclear contributions to the magnetic interaction are negligibly small compared to those of the electrons and are therefore commonly neglected.

We now define the (permanent) magnetic moment, magnetizability, and first hypermagnetizability of the system as

$$\begin{array}{cc}\hfill {\mu }_{\mathrm{m}}^a\left({\omega }_A\right)& =-W^a\left({\omega }_A\right)\hfill \\ \hfill & =-D_{\nu \mu }^0{H}_{\mu \nu }^a,\hfill \end{array}$$

(80a)

$$\begin{array}{cc}\hfill {\xi }^{ab}({\omega }_A,{\omega }_B)& =-W^{ab}({\omega }_A,{\omega }_B)\hfill \\ \hfill & =-{\displaystyle \frac{1}{2}}\left[D_{\nu \mu }^a\left({\omega }_A\right)H_{\mu \nu }^b+D_{\nu \mu }^b\left({\omega }_B\right)H_{\mu \nu }^a\right]-D_{\nu \mu }^0{H}_{\mu \nu }^{ab},\hfill \end{array}$$

(80b)

$$\begin{array}{cc}\hfill {\zeta }^{abc}({\omega }_A,{\omega }_B,{\omega }_C)& =-W^{abc}({\omega }_A,{\omega }_B,{\omega }_C)\hfill \\ \hfill & =-{\displaystyle \frac{n_i}{2}}\left[U_{ji}^{ab}({\omega }_A,{\omega }_B)U_{ji}^{c*}(-{\omega }_C)-U_{ji}^{ab*}(-{\omega }_A,-{\omega }_B)U_{ji}^c\left({\omega }_C\right)\right.\hfill \\ \hfill & +U_{ji}^{ac}({\omega }_A,{\omega }_C)U_{ji}^{b*}(-{\omega }_B)-U_{ji}^{ac*}(-{\omega }_A,-{\omega }_C)U_{ji}^b\left({\omega }_B\right)\hfill \\ \hfill & +\left.U_{ji}^{bc}({\omega }_B,{\omega }_C)U_{ji}^{a*}(-{\omega }_A)-U_{ji}^{bc*}(-{\omega }_B,-{\omega }_C)U_{ji}^a\left({\omega }_A\right)\right]\hfill \\ \hfill & \times \left({\omega }_A+{\omega }_B+{\omega }_C\right)-({\rho }^a\left({\omega }_A\right)|k_{\mathrm{xc}}|{\rho }^b\left({\omega }_B\right)\left|{\rho }^c\left({\omega }_C\right)\right)\hfill \\ \hfill & -n_i{F}_{jl}^a\left({\omega }_A\right)\left[U_{li}^b\left({\omega }_B\right)U_{ji}^{c*}(-{\omega }_C)+U_{li}^c\left({\omega }_C\right)U_{ji}^{b*}(-{\omega }_B)\right]\hfill \\ \hfill & -n_i{F}_{jl}^b\left({\omega }_B\right)\left[U_{li}^a\left({\omega }_A\right)U_{ji}^{c*}(-{\omega }_C)+U_{li}^c\left({\omega }_C\right)U_{ji}^{a*}(-{\omega }_A)\right]\hfill \\ \hfill & -n_i{F}_{jl}^c\left({\omega }_C\right)\left[U_{li}^a\left({\omega }_A\right)U_{ji}^{b*}(-{\omega }_B)+U_{li}^b\left({\omega }_B\right)U_{ji}^{a*}(-{\omega }_A)\right]\hfill \\ \hfill & +n_i{F}_{ki}^{a\phantom{b}}\left({\omega }_A\right)\left[U_{ji}^{b*}(-{\omega }_B)U_{jk}^c\left({\omega }_C\right)+U_{ji}^{c*}(-{\omega }_C)U_{jk}^b\left({\omega }_B\right)\right]\hfill \\ \hfill & +n_i{F}_{ki}^b\left({\omega }_B\right)\left[U_{ji}^{a*}(-{\omega }_A)U_{jk}^c\left({\omega }_C\right)+U_{ji}^{c*}(-{\omega }_C)U_{jk}^a\left({\omega }_A\right)\right]\hfill \\ \hfill & +n_i{F}_{ki}^{c\phantom{b}}\left({\omega }_C\right)\left[U_{ji}^{a*}(-{\omega }_A)U_{jk}^b\left({\omega }_B\right)+U_{ji}^{b*}(-{\omega }_B)U_{jk}^a\left({\omega }_A\right)\right]\hfill \\ \hfill & -D_{\nu \mu }^a\left({\omega }_A\right)H_{\mu \nu }^{bc}-D_{\nu \mu }^b\left({\omega }_B\right)H_{\mu \nu }^{ac}-D_{\nu \mu }^c\left({\omega }_C\right)H_{\mu \nu }^{ab},\hfill \end{array}$$

(80c)

where $i,k\in \mathrm{occ}$, $j,l\in \mathrm{vir}$, and a field-independent basis set is used. Due to the presence of nonlinear (quadratic) terms in the magnetic interaction operator, the one-electron core Hamiltonian has nonvanishing derivatives up to second order:

$$H_{\mu \nu }^a={\displaystyle \frac{1}{2}}\left(\mu \right|-\mathrm{i}{ϵ}_{\phantom{a}jk}^a{r}^j{\partial }^k+{\sigma }^a\left|\nu \right),$$

(81a)

$$H_{\mu \nu }^{ab}={\displaystyle \frac{1}{4}}\left(\mu \right|{\delta }^{ab}{r}^2-r^a{r}^b\left|\nu \right),$$

(81b)

where ${ϵ}_{\phantom{a}jk}^a$ is the Levi-Civita tensor. Eqs. (80b) and (80c) can be decomposed into paramagnetic and diamagnetic contributions, owing to the differing orders of the paramagnetic and diamagnetic interaction operators in the external magnetic field. Specifically, for magnetizability, we have

$${\xi }_{\mathrm{p}}^{ab}({\omega }_A,{\omega }_B)=-{\displaystyle \frac{1}{2}}\left[D_{\nu \mu }^a\left({\omega }_A\right)H_{\mu \nu }^b+D_{\nu \mu }^b\left({\omega }_B\right)H_{\mu \nu }^a\right],$$

(82a)

$${\xi }_{\mathrm{d}}^{ab}=-D_{\nu \mu }^0{H}_{\mu \nu }^{ab}.$$

(82b)

For the first hypermagnetizability, the paramagnetic and diamagnetic parts are given by

$$\begin{array}{cc}\hfill {\zeta }_{\mathrm{p}}^{abc}({\omega }_A,{\omega }_B,{\omega }_C)& =-{\displaystyle \frac{n_i}{2}}\left[U_{ji}^{ab}({\omega }_A,{\omega }_B)U_{ji}^{c*}(-{\omega }_C)-U_{ji}^{ab*}(-{\omega }_A,-{\omega }_B)U_{ji}^c\left({\omega }_C\right)\right.\hfill \\ \hfill & +U_{ji}^{ac}({\omega }_A,{\omega }_C)U_{ji}^{b*}(-{\omega }_B)-U_{ji}^{ac*}(-{\omega }_A,-{\omega }_C)U_{ji}^b\left({\omega }_B\right)\hfill \\ \hfill & +\left.U_{ji}^{bc}({\omega }_B,{\omega }_C)U_{ji}^{a*}(-{\omega }_A)-U_{ji}^{bc*}(-{\omega }_B,-{\omega }_C)U_{ji}^a\left({\omega }_A\right)\right]\hfill \\ \hfill & \times \left({\omega }_A+{\omega }_B+{\omega }_C\right)-({\rho }^a\left({\omega }_A\right)|k_{\mathrm{xc}}|{\rho }^b\left({\omega }_B\right)\left|{\rho }^c\left({\omega }_C\right)\right)\hfill \\ \hfill & -n_i{F}_{jl}^a\left({\omega }_A\right)\left[U_{li}^b\left({\omega }_B\right)U_{ji}^{c*}(-{\omega }_C)+U_{li}^c\left({\omega }_C\right)U_{ji}^{b*}(-{\omega }_B)\right]\hfill \\ \hfill & -n_i{F}_{jl}^b\left({\omega }_B\right)\left[U_{li}^a\left({\omega }_A\right)U_{ji}^{c*}(-{\omega }_C)+U_{li}^c\left({\omega }_C\right)U_{ji}^{a*}(-{\omega }_A)\right]\hfill \\ \hfill & -n_i{F}_{jl}^c\left({\omega }_C\right)\left[U_{li}^a\left({\omega }_A\right)U_{ji}^{b*}(-{\omega }_B)+U_{li}^b\left({\omega }_B\right)U_{ji}^{a*}(-{\omega }_A)\right]\hfill \\ \hfill & +n_i{F}_{ki}^{a\phantom{b}}\left({\omega }_A\right)\left[U_{ji}^{b*}(-{\omega }_B)U_{jk}^c\left({\omega }_C\right)+U_{ji}^{c*}(-{\omega }_C)U_{jk}^b\left({\omega }_B\right)\right]\hfill \\ \hfill & +n_i{F}_{ki}^b\left({\omega }_B\right)\left[U_{ji}^{a*}(-{\omega }_A)U_{jk}^c\left({\omega }_C\right)+U_{ji}^{c*}(-{\omega }_C)U_{jk}^a\left({\omega }_A\right)\right]\hfill \\ \hfill & +n_i{F}_{ki}^{c\phantom{b}}\left({\omega }_C\right)\left[U_{ji}^{a*}(-{\omega }_A)U_{jk}^b\left({\omega }_B\right)+U_{ji}^{b*}(-{\omega }_B)U_{jk}^a\left({\omega }_A\right)\right],\hfill \end{array}$$

(83a)

$${\zeta }_{\mathrm{d}}^{abc}({\omega }_A,{\omega }_B,{\omega }_C)=-D_{\nu \mu }^a\left({\omega }_A\right)H_{\mu \nu }^{bc}-D_{\nu \mu }^b\left({\omega }_B\right)H_{\mu \nu }^{ac}-D_{\nu \mu }^c\left({\omega }_C\right)H_{\mu \nu }^{ab}.$$

(83b)

Magnetic response properties bring in considerations related to gauge transformations, as the vector potential enters explicitly in the interaction Hamiltonian. According to Maxwell’s equations, a gauge transformation of the scalar and vector potentials,

$$\begin{array}{cc}\hfill \phi \to {\phi }^{\prime }& =\phi -{\partial }_{t}g,\hfill \end{array}$$

(84a)

$$\begin{array}{cc}\hfill \mathit{A}\to {\mathit{A}}^{\prime }& =\mathit{A}+\nabla g,\hfill \end{array}$$

(84b)

with g being the gauge function, leaves the physical electromagnetic fields unchanged. Similarly, the HF and KS equations must obey gauge covariance—i.e., they must retain their validity under such transformations. It would be unphysical for the HF/KS equations to break down simply because the gauge has changed. Under a gauge transformation, the one-electron wavefunction acquires a position- and time-dependent phase factor:

$${\psi }_i\to {\psi }_i^{\prime }={\psi }_i{\mathrm{e}}^{-\mathrm{i}g(\mathit{r},t)}.$$

(85)

As rigorously shown in Appendix D, The function g appearing in the exponent of the phase factor is exactly the gauge function.

When a gauge transformation is considered, both the transformed operators and wavefunctions (more specifically, the basis functions) must be taken into account. Due to the introduction of a phase factor, the basis functions become dependent on the external magnetic field (see Eq. (A63)). In this case, more general expressions for the quasi-energy derivatives are required to properly evaluate magnetic response properties. Fortunately, the added phase factor does not affect the overlap between basis functions, as it cancels out:

$$S_{\mu \nu }^{\prime }=\left({\mu }^{\prime }\right|{\nu }^{\prime })=\left(\mu \right|\nu )=S_{\mu \nu }.$$

(86)

It follows that all field derivatives of $S_{\mu \nu }$ vanish at any order. As a result, the magnetic response properties up to third order take the forms

$${\mu }_{\mathrm{m}}^a\left({\omega }_A\right)=D_{\nu \mu }^0\left[{\mathcal{T}}_{\mu \nu }^a\left({\omega }_A\right)-H_{\mu \nu }^a\left({\omega }_A\right)\right],$$

(87a)

$$\begin{array}{cc}\hfill {\xi }^{ab}({\omega }_A,{\omega }_B)& ={\displaystyle \frac{1}{2}}{D}_{\nu \mu }^a\left({\omega }_A\right)\left[{\mathcal{T}}_{\mu \nu }^b\left({\omega }_B\right)-H_{\mu \nu }^b\left({\omega }_B\right)\right]\hfill \\ \hfill & +{\displaystyle \frac{1}{2}}{D}_{\nu \mu }^b\left({\omega }_B\right)\left[{\mathcal{T}}_{\mu \nu }^a\left({\omega }_A\right)-H_{\mu \nu }^a\left({\omega }_A\right)\right]\hfill \\ \hfill & +D_{\nu \mu }^0\left[{\mathcal{T}}_{\mu \nu }^{ab}({\omega }_A,{\omega }_B)-H_{\mu \nu }^{ab}({\omega }_A,{\omega }_B)\right],\hfill \end{array}$$

(87b)

$$\begin{array}{cc}\hfill {\zeta }^{abc}({\omega }_A,{\omega }_B,{\omega }_C)& =-{\displaystyle \frac{n_i}{2}}\left[U_{ji}^{ab}({\omega }_A,{\omega }_B)U_{ji}^{c*}(-{\omega }_C)-U_{ji}^{ab*}(-{\omega }_A,-{\omega }_B)U_{ji}^c\left({\omega }_C\right)\right.\hfill \\ \hfill & +U_{ji}^{ac}({\omega }_A,{\omega }_C)U_{ji}^{b*}(-{\omega }_B)-U_{ji}^{ac*}(-{\omega }_A,-{\omega }_C)U_{ji}^b\left({\omega }_B\right)\hfill \\ \hfill & +\left.U_{ji}^{bc}({\omega }_B,{\omega }_C)U_{ji}^{a*}(-{\omega }_A)-U_{ji}^{bc*}(-{\omega }_B,-{\omega }_C)U_{ji}^a\left({\omega }_A\right)\right]\hfill \\ \hfill & \times \left({\omega }_A+{\omega }_B+{\omega }_C\right)-({\rho }^a\left({\omega }_A\right)|k_{\mathrm{xc}}|{\rho }^b\left({\omega }_B\right)\left|{\rho }^c\left({\omega }_C\right)\right)\hfill \\ \hfill & -n_i\left[F_{jl}^a\left({\omega }_A\right)-{\mathcal{T}}_{jl}^a\left({\omega }_A\right)\right]\left[U_{li}^b\left({\omega }_B\right)U_{ji}^{c*}(-{\omega }_C)+U_{li}^c\left({\omega }_C\right)U_{ji}^{b*}(-{\omega }_B)\right]\hfill \\ \hfill & -n_i\left[F_{jl}^b\left({\omega }_B\right)-{\mathcal{T}}_{jl}^b\left({\omega }_B\right)\right]\left[U_{li}^a\left({\omega }_A\right)U_{ji}^{c*}(-{\omega }_C)+U_{li}^c\left({\omega }_C\right)U_{ji}^{a*}(-{\omega }_A)\right]\hfill \\ \hfill & -n_i\left[F_{jl}^c\left({\omega }_C\right)-{\mathcal{T}}_{jl}^c\left({\omega }_C\right)\right]\left[U_{li}^a\left({\omega }_A\right)U_{ji}^{b*}(-{\omega }_B)+U_{li}^b\left({\omega }_B\right)U_{ji}^{a*}(-{\omega }_A)\right]\hfill \\ \hfill & +n_i\left[F_{ki}^{a\phantom{b}}\left({\omega }_A\right)-{\mathcal{T}}_{ki}^a\left({\omega }_A\right)\right]\left[U_{ji}^{b*}(-{\omega }_B)U_{jk}^c\left({\omega }_C\right)+U_{ji}^{c*}(-{\omega }_C)U_{jk}^b\left({\omega }_B\right)\right]\hfill \\ \hfill & +n_i\left[F_{ki}^b\left({\omega }_B\right)-{\mathcal{T}}_{ki}^b\left({\omega }_B\right)\right]\left[U_{ji}^{a*}(-{\omega }_A)U_{jk}^c\left({\omega }_C\right)+U_{ji}^{c*}(-{\omega }_C)U_{jk}^a\left({\omega }_A\right)\right]\hfill \\ \hfill & +n_i\left[F_{ki}^{c\phantom{b}}\left({\omega }_C\right)-{\mathcal{T}}_{ki}^c\left({\omega }_C\right)\right]\left[U_{ji}^{a*}(-{\omega }_A)U_{jk}^b\left({\omega }_B\right)+U_{ji}^{b*}(-{\omega }_B)U_{jk}^a\left({\omega }_A\right)\right]\hfill \\ \hfill & +D_{\nu \mu }^a\left({\omega }_A\right)\left[{\mathcal{T}}_{\mu \nu }^{bc}({\omega }_B,{\omega }_C)-H_{\mu \nu }^{bc}({\omega }_B,{\omega }_C)\right]\hfill \\ \hfill & +D_{\nu \mu }^b\left({\omega }_B\right)\left[{\mathcal{T}}_{\mu \nu }^{ac}({\omega }_A,{\omega }_C)-H_{\mu \nu }^{ac}({\omega }_A,{\omega }_C)\right]\hfill \\ \hfill & +D_{\nu \mu }^c\left({\omega }_C\right)\left[{\mathcal{T}}_{\mu \nu }^{ab}({\omega }_A,{\omega }_B)-H_{\mu \nu }^{ab}({\omega }_A,{\omega }_B)\right]\hfill \\ \hfill & +D_{\nu \mu }^0\left[{\mathcal{T}}_{\mu \nu }^{abc}({\omega }_A,{\omega }_B,{\omega }_C)-H_{\mu \nu }^{abc}({\omega }_A,{\omega }_B,{\omega }_C)\right].\hfill \end{array}$$

(87c)

Note that the gauge transformation of the scalar potential causes the core Hamiltonian matrix $H_{\mu \nu }$ to become explicitly dependent on the magnetic fields (or field frequencies), which must be properly accounted for in practical calculations. ${\mathcal{T}}_{\mu \nu }$ can now be explicitly expressed as

$$\begin{array}{cc}\hfill {\mathcal{T}}_{\mu \nu }& =\left({\mu }^{\prime }\right|\mathrm{i}{\partial }_t\left|{\nu }^{\prime }\right)\hfill \\ \hfill & =\left(\mu \right|\left({\partial }_{t}g\right)\left|\nu \right)\hfill \\ \hfill & =-{\displaystyle \frac{\mathrm{i}}{2}}\sum _A{\omega }_A{\mathcal{B}}_A{\mathrm{e}}^{\mathrm{i}{\omega }_{A}t}·\mathit{R}\times \left(\mu \right|\mathit{r}\left|\nu \right),\hfill \end{array}$$

(88)

which is, perhaps unexpectedly, identical to the result of the gauge-transformation operator for the scalar potential. Moreover, it involves only the first-order field derivative:

$${\mathcal{T}}_{\mu \nu }^a\left({\omega }_A\right)=-{\displaystyle \frac{\mathrm{i}{\omega }_A}{2}}{ϵ}_{\phantom{a}jk}^a{R}^j\left(\mu \right|r^k\left|\nu \right).$$

(89)

Eqs. (87a)–(87c) are exactly equivalent to Eqs. (80a)–(80c), thereby confirming that the response properties derived from quasi-energy derivatives are completely independent of the choice of gauge origin, as all gauge-dependent contributions cancel out. This is a direct consequence of the gauge covariance of the HF/KS equations.

3.2.1. Comparison with the GIAO Method

The magnetic response properties derived above are gauge-invariant, as the quasi-energy itself is a gauge-invariant quantity. In contrast, the energy and the dipole expectation are not conserved under gauge transformations. Consequently, frameworks based on their field derivatives inevitably lead to gauge-dependent magnetic properties. To address this so-called gauge-origin dependence, several gauge localization schemes[36,37,38,39,40,41] have been developed. Among them, the most widely adopted is the gauge-including/independent atomic orbitals (GIAOs)[39,40,41], originally referred to as London atomic orbitals[42]. In this method, the gauge in a gauge-transformed basis function ${\varphi }_{\nu }^{\prime }$ is modified as

$$g_{\nu }=-{\displaystyle \frac{1}{2}}\mathcal{B}·\left(\mathit{R}-{\mathit{R}}_{\nu }\right)\times \mathit{r},$$

(90)

where the gauge origin $\mathit{R}$ is fixed at the center ${\mathit{R}}_{\nu }$ of the $\nu$-th atomic orbital, i.e.,

$${\varphi }_{\nu }^{\prime }={\varphi }_{\nu }{\mathrm{e}}^{-\mathrm{i}{g}_{\nu }}={\varphi }_{\nu }{\mathrm{e}}^{-{\displaystyle \frac{\mathrm{i}}{2}}\mathcal{B}·\left({\mathit{R}}_{\nu }-\mathit{R}\right)\times \mathit{r}}.$$

(91)

However, the gauge function appearing in the field-dependent operators remains unchanged. This artificial treatment compromises the covariance of the HF/KS equations under gauge transformations.

In this case, the gauge-transformed overlap matrix $S_{\mu \nu }^{\prime }$ no longer vanishes:

$$S_{\mu \nu }^{\prime }=\left(\mu \right|{\mathrm{e}}^{\mathcal{B}·{\mathcal{U}}_{\mu \nu }}\left|\nu \right),$$

(92)

where we define

$${\mathcal{U}}_{\mu \nu }=-{\displaystyle \frac{\mathrm{i}}{2}}{\mathit{R}}_{\mu \nu }\times \mathit{r}=-{\displaystyle \frac{\mathrm{i}}{2}}\left({\mathit{R}}_{\nu }-{\mathit{R}}_{\mu }\right)\times \mathit{r}.$$

(93)

The matrix ${\mathcal{T}}_{\mu \nu }^{\prime }$ becomes

$${\mathcal{T}}_{\mu \nu }^{\prime }={\displaystyle \frac{\mathrm{i}}{2}}\sum _A{\omega }_A{\mathcal{B}}_A{\mathrm{e}}^{\mathrm{i}{\omega }_{A}t}·\left({\mathit{R}}_{\nu }-\mathit{R}\right)\times \left(\mu \right|\mathit{r}{\mathrm{e}}^{\mathcal{B}·{\mathcal{U}}_{\mu \nu }}\left|\nu \right).$$

(94)

Both Eq. (92) and Eq. (94) contribute nonzero higher-order field derivatives, which necessitates the use of the full formulation with field-dependent basis sets in magnetic response theory.

Moreover, one can show that gauge-origin dependent terms in Eqs. (87a)–(87c) such as

$$\begin{array}{cc}\hfill & H_{\mu \nu }^{\prime a}\left({\omega }_A\right)-{\mathcal{T}}_{\mu \nu }^{\prime a}\left({\omega }_A\right)\hfill \\ \hfill & =\left(\mu \left|{\mathcal{U}}_{\mu \nu }^a\left[-{\displaystyle \frac{1}{2}}{\nabla }^2-\sum _I{\displaystyle \frac{Z_I}{\left|\mathit{r}-{\mathit{r}}_I\right|}}\right]-{\displaystyle \frac{\mathrm{i}}{2}}{ϵ}_{\phantom{a}jk}^a{r}_{\nu }^j{\partial }^k+{\displaystyle \frac{1}{2}}{\sigma }^a-{\displaystyle \frac{\mathrm{i}}{2}}{\omega }_A{ϵ}_{\phantom{a}jk}^a{R}_{\nu }^j{r}^k\right|\nu \right)\hfill \end{array}$$

(95)

and

$$\begin{array}{cc}\hfill & H_{\mu \nu }^{\prime ab}({\omega }_A,{\omega }_B)-{\mathcal{T}}_{\mu \nu }^{\prime ab}({\omega }_A,{\omega }_B)\hfill \\ \hfill & =\left(\mu \right|{\mathcal{U}}_{\mu \nu }^a{\mathcal{U}}_{\mu \nu }^b\left[-{\displaystyle \frac{1}{2}}{\nabla }^2-\sum _I{\displaystyle \frac{Z_I}{\left|\mathit{r}-{\mathit{r}}_I\right|}}\right]-{\displaystyle \frac{\mathrm{i}}{2}}\left[{\mathcal{U}}_{\mu \nu }^a{ϵ}_{\phantom{b}jk}^b+{\mathcal{U}}_{\mu \nu }^b{ϵ}_{\phantom{a}jk}^a\right]r_{\nu }^j{\partial }^k\hfill \\ \hfill & +{\displaystyle \frac{1}{2}}\left[{\mathcal{U}}_{\mu \nu }^a{\sigma }^b+{\mathcal{U}}_{\mu \nu }^b{\sigma }^a\right]-{\displaystyle \frac{\mathrm{i}}{2}}\left[{\omega }_B{\mathcal{U}}_{\mu \nu }^a{ϵ}_{\phantom{b}jk}^b+{\omega }_A{\mathcal{U}}_{\mu \nu }^b{ϵ}_{\phantom{a}jk}^a\right]R_{\nu }^j{r}^k+{\displaystyle \frac{1}{4}}\left[{\delta }^{ab}{r}_{\nu }^2-r_{\nu }^a{r}_{\nu }^b\right]\left|\nu \right)\hfill \end{array}$$

(96)

(with ${\mathit{r}}_{\nu }$ defined as $\mathit{r}-{\mathit{R}}_{\nu }$) merely redistribute the gauge dependence over the atoms. While the resulting expressions may appear free of explicit dependence on the gauge origin $\mathit{R}$, the computed magnetic properties remain strongly basis-set-dependent.

In this regard, the GIAO method does not fundamentally resolve the gauge problem but instead introduces additional complexity. If all the atomic centers ${\mathit{R}}_{\mu },{\mathit{R}}_{\nu },\dots$ in the GIAO expressions are set to zero, one immediately recovers the expressions under the conventional gauge transformation. The core issue is that GIAOs fail to recognize the inherent gauge covariance of the HF/KS equations. Since the quasi-energy is an observable, it should by construction be gauge-invariant—without resorting to artificial gauge localization schemes.

4. Results and Discussion

4.1. Numerical Tests

We have implemented the newly developed electromagnetic response theory based on the HF/KS quasi-energy framework within the PySCF platform[43,44,45,46] for computing general response properties. The core computational step is solving the CP-HF/KS equations (54) and (67), which are linear and can, in principle, be handled by standard linear solvers. However, such direct approaches are often inefficient for large-scale and sparse systems. To improve performance, our implementation includes interfaces to not only an exact solver but more efficient iterative techniques, including Krylov subspace methods[47,48] and Newton-Krylov (N-K) algorithms[49,50]. These methods enable faster convergence to a target accuracy in many cases. We benchmarked all implementations against the conventional solver as a reference, evaluating both their numerical accuracy and computational efficiency.

First-order response properties require no additional elaborate computation; once the unperturbed density is known, they can be directly evaluated. Therefore, these properties are not included in the following numerical tests. When the sum of the frequencies of the external fields acting on different components of the response tensor is zero, the second- and third-order electromagnetic response properties obtained from our formalism coincide with those given by existing frameworks. Accordingly, we computed the second- and third-order electromagnetic responses corresponding to several conventional optical properties of the H2O molecule (O-H bond length: 0.957 Å, H-O-H bond angle: 104.5°) at both the spin-restricted HF and KS-DFT levels, using the STO-3G basis set from the Basis Set Exchange (BSE) database[51]. The results obtained by different numerical solvers were compared.

For KS-DFT, which includes electron correlation effects to some extent, the numerical complexity of the $k_{\mathrm{xc}}$ term in Eq. (72) depends on the form of the xc functional (see Appendix C). For example, the local-density-approximation (LDA) functional depends only on the electron density $\rho$; the generalized-gradient-approximation (GGA) functional further includes its gradient $\nabla \rho$; and the meta-GGA (mGGA) functional introduces ${\nabla }^2\rho$ and the kinetic-energy density $\tau$ as additional variables. Non-collinear functionals may also depend on the spin magnetization vector $\mathit{m}$ and the kinetic spin-density vector $\mathit{u}$. Although our implementation is capable of numerically evaluating all these cases, the present tests were limited to collinear functionals, among which SVWN5, PBE, and TPSS were chosen as representatives of the LDA, GGA, and mGGA functionals, respectively. The radial and angular grids were set to 99 and 590 points for the numerical integration.

Table 1 and Table 2 summarize the computed results for the electric and magnetic response properties, respectively. The monochromatic field was assigned a frequency of 0.1 a.u. To make the data more compact, spherical averages of the response tensors were used—for example, the spherically averaged polarizability:

$$\overline{\alpha }={\displaystyle \frac{1}{3}}{\alpha }_{ii},$$

(97)

and the spherically averaged first hyperpolarizability:

$$\overline{\beta }=\sqrt{{\tilde{\beta }}_i{\tilde{\beta }}_i},\mathrm{where}{\tilde{\beta }}_i={\displaystyle \frac{1}{3}}\left({\beta }_{ijj}+{\beta }_{jij}+{\beta }_{jji}\right).$$

(98)

For the electric response properties, reference data from Gaussian 16 Rev. B.01[52] are provided for comparison, while for the magnetic responses, results from Dalton 2020.1[53] are included. There are only a few quantum chemistry packages that support the calculation of magnetizabilities. Even the Dalton program still lacks a complete library of mGGA functionals and is limited to the evaluation of static magnetizabilities for closed-shell systems. Since odd-order magnetic response properties of closed-shell systems are purely imaginary for one-component nonrelativistic theory, only the second-order results are presented here. In Table 2, the diamagnetic contribution depends solely on the unperturbed density and does not require solving the CP equations, whereas the paramagnetic term, analogous to the electric polarizability, involves the first-order CP equations. Accurate solutions of the CP-HF/KS equations rely on well-converged unperturbed molecular orbitals; therefore, the ground-state SCF energy was converged to ${10}^{-11}$ a.u. in all calculations. The Krylov solver used for the CP-HF/KS equations employed a convergence threshold of ${10}^{-9}$ a.u. within the subspace iterations.

As shown in Table 1 and Table 2, the results obtained from our implementation are largely reliable. In Table 1, the deviations from the Gaussian 16 results become slightly more pronounced with the increasing complexity of the xc functional, which is likely attributable to the different numerical integration grid pruning schemes employed in the two programs (NWChem-type grids in PySCF vs. SG-1 grids in Gaussian 16). In addition, the results produced by the various solvers implemented in our code are virtually identical. Considering the superior efficiency of the Krylov solver for large linear systems, it is therefore adopted as the default solver for the CP equations in subsequent calculations.

4.2. Non-Optical Properties

The electromagnetic response formalism developed in this work encompasses all conventional optical properties, which correspond to processes where the sum of the perturbation frequencies equals zero—i.e., the frequencies of the emitted fields match those of the incident ones. In practice, however, inelastic scattering processes can occur, where the frequency balance is broken due to energy exchange between the external field and the electronic system. Such situations give rise to non-optical response properties, including frequency-mixed or frequency-converted responses that are typically inaccessible in standard optical response theories. The present formalism enables direct prediction of these non-optical responses, which can exhibit pronounced dependence on the field frequencies. To illustrate this capability, we computed several frequency-dependent non-optical properties of the H2O molecule at $\omega =0.05,0.10,0.15,0.20$ a.u., and the results are summarized in Table 3.

To obtain nonzero third-order magnetic responses, the generalized HF/KS framework was employed. This approach allows each molecular orbital to mix $\alpha$- and $\beta$-spin components, thereby accounting for spin-magnetic coupling and non-collinear operator contributions that are essential for accurately describing magnetic correlation effects.

For non-optical properties, as indicated by Eqs. (76c) and (80c), the $(2n+1)$ rule breaks down, and the third-order responses show explicit frequency dependence. Consequently, they vary more strongly with the field frequency than the second-order responses under dynamic fields, as illustrated in Table 3.

4.3. Gauge-Invariant Magnetic Properties

Gauge transformations of the electromagnetic potentials introduce a gauge-dependent phase factor ${\mathrm{e}}^{-\mathrm{i}g}$ into the orbitals and simultaneously impose a unitary transformation on the Hamiltonian. Since a unitary transformation leaves all physical observables of the system unchanged, the quasi-energy and its derivatives are inherently gauge-invariant without any additional corrections.

In contrast, the magnetic properties derived from the GIAO method are based on energy derivatives. Under a dynamic electromagnetic field, the energy is not a gauge-conserved quantity, and the correction introduced by GIAO effectively adds a degree of arbitrariness by distributing the gauge origin among individual atomic centers. This leads to a stronger basis-set dependence of the computed properties. Nevertheless, as the basis set approaches completeness, the difference between GIAO and non-GIAO results diminishes, as can be seen from Table 4. Despite its widespread use, the GIAO approach is conceptually and computationally more involved, and offers limited physical insight.

5. Conclusions

Conventional electromagnetic response theories often suffer from internal inconsistencies and limited generality, particularly under time-dependent perturbations, where the expectation value of the Hamiltonian is no longer an observable. In such cases, frameworks based on energy derivatives become invalid, whereas the quasi-energy provides a more general and physically well-defined quantity. In this work, we have developed a general electromagnetic response theory based on quasi-energy derivatives, which enables the calculation of a wide range of electromagnetic response properties. First, the present theoretical framework provides a more general, consistent, and rigorous derivation, particularly for systems under time-dependent perturbations. Second, the formalism places no restriction on the type or combination of external fields, allowing the treatment of arbitrary dynamic or static perturbations. As a result, it not only encompasses all optical properties but also extends naturally to the prediction of non-optical responses. Additionally, the gauge dependence of magnetic properties is fully resolved within this framework, since the quasi-energy and its derivatives are themselves gauge-invariant quantities.

Appendix B.1. The First-Order W-Response

Let us first evaluate the first term in the first-order energy response (14):

$$\begin{array}{cc}\hfill D_{\nu \mu }^a\left({\omega }_A\right)F_{\mu \nu }^0& =n_i\left[C_{\nu j}^0{U}_{ji}^a\left({\omega }_A\right)C_{\mu i}^{0*}+C_{\nu i}^0{C}_{\mu j}^{0*}{U}_{ji}^{a*}(-{\omega }_A)\right]F_{\mu \nu }^0\hfill \\ \hfill & =n_i\left[U_{ji}^a\left({\omega }_A\right)F_{ij}^0+U_{ji}^{a*}(-{\omega }_A)F_{ji}^0\right].\hfill \end{array}$$

(A4)

Since the ground-state Fock matrix is diagonal in the molecular orbital basis, i.e., $F_{ij}^0={\delta }_{ij}{\epsilon }_j^0$, the summation over j yields

$$U_{ji}^a\left({\omega }_A\right)F_{ij}^0+U_{ji}^{a*}(-{\omega }_A)F_{ji}^0=\left[U_{ii}^a\left({\omega }_A\right)+U_{ii}^{a*}(-{\omega }_A)\right]{\epsilon }_i^0.$$

(A5)

By virtue of the orthonormality condition (33), we have

$$D_{\nu \mu }^a\left({\omega }_A\right)F_{\mu \nu }^0=-n_i{S}_{ii}^a{\epsilon }_i^0.$$

(A6)

Thus, the first-order response of the electronic energy becomes

$$E^a=-n_i{S}_{ii}^a{\epsilon }_i^0+D_{\nu \mu }^0{H}_{\mu \nu }^a+{\displaystyle \frac{1}{2}}{D}_{\nu \mu }^0{D}_{\lambda \kappa }^0{\left(\mu \nu \right|\kappa \lambda )}^a+D_{\nu \mu }^0\left({\mu \nu }^a\right|v_{\mathrm{xc}}).$$

(A7)

Now we turn to the first-order response of T, as given in Eq. (19b). The second term therein can be rewritten as

$$\begin{array}{cc}\hfill -{\omega }_A{n}_i{C}_{\mu i}^{0*}{S}_{\mu \nu }^0{C}_{\nu i}^a\left({\omega }_A\right)& =-{\omega }_A{n}_i{C}_{\mu i}^{0*}{S}_{\mu \nu }^0{C}_{\nu j}^0{U}_{ji}^a\left({\omega }_A\right)\hfill \\ \hfill & =-{\omega }_A{n}_i{\delta }_{ij}{U}_{ji}^a\left({\omega }_A\right)\hfill \\ \hfill & =-{\omega }_A{n}_i{U}_{ii}^a\left({\omega }_A\right)\hfill \\ \hfill & ={\displaystyle \frac{1}{2}}{\omega }_A{n}_i{S}_{ii}^a.\hfill \end{array}$$

(A8)

Then, Eq. (19b) is equivalent to

$$T^a\left({\omega }_A\right)=D_{\nu \mu }^0{\mathcal{T}}_{\mu \nu }^a\left({\omega }_A\right)+{\displaystyle \frac{1}{2}}{\omega }_A{n}_i{S}_{ii}^a.$$

(A9)

The total first-order quasi-energy response is therefore Eq. (68)

Appendix B.2. The Second-Order W-Response

Similarly, we begin by evaluating the first term in Eq. (16), using the definition of the second-order density matrix (28):

$$\begin{array}{cc}\hfill D_{\nu \mu }^{ab}({\omega }_A,{\omega }_B)F_{\mu \nu }^0& =n_i\left[U_{ji}^{ab}({\omega }_A,{\omega }_B)F_{ij}^0+U_{ji}^a\left({\omega }_A\right)U_{ki}^{b*}(-{\omega }_B)F_{kj}^0\right.\hfill \\ \hfill & +\left.U_{ji}^b\left({\omega }_B\right)U_{ki}^{a*}(-{\omega }_A)F_{kj}^0+U_{ji}^{ab*}(-{\omega }_A,-{\omega }_B)F_{ji}^0\right].\hfill \end{array}$$

(A10)

Owing to the diagonal structure of the zeroth-order Fock matrix in the molecular orbital basis, the summation over j can be carried out to produce

$$\begin{array}{cc}\hfill & U_{ji}^{ab}({\omega }_A,{\omega }_B)F_{ij}^0+U_{ji}^a\left({\omega }_A\right)U_{ki}^{b*}(-{\omega }_B)F_{kj}^0+U_{ji}^b\left({\omega }_B\right)U_{ki}^{a*}(-{\omega }_A)F_{kj}^0+U_{ji}^{ab*}(-{\omega }_A,-{\omega }_B)F_{ji}^0\hfill \\ \hfill & =\left[U_{ji}^a\left({\omega }_A\right)U_{ji}^{b*}(-{\omega }_B)+U_{ji}^b\left({\omega }_B\right)U_{ji}^{a*}(-{\omega }_A)\right]{\epsilon }_j^0+\left[U_{ii}^{ab}({\omega }_A,{\omega }_B)+U_{ii}^{ab*}(-{\omega }_A,-{\omega }_B)\right]{\epsilon }_i^0,\hfill \end{array}$$

(A11)

where the expression inside the second brackets can be reduced by virtue of the orthonormality condition (34), resulting in

$$\begin{array}{cc}\hfill & U_{ii}^{ab}({\omega }_A,{\omega }_B)+U_{ii}^{ab*}(-{\omega }_A,-{\omega }_B)\hfill \\ \hfill & =-U_{ji}^{a*}(-{\omega }_A)S_{ji}^b-U_{ji}^{a*}(-{\omega }_A)U_{ji}^b\left({\omega }_B\right)-U_{ji}^{b*}(-{\omega }_B)S_{ji}^a\hfill \\ \hfill & -S_{ii}^{ab}-S_{ji}^{a*}{U}_{ji}^b\left({\omega }_B\right)-U_{ji}^{b*}(-{\omega }_B)U_{ji}^a\left({\omega }_A\right)-S_{ji}^{b*}{U}_{ji}^a\left({\omega }_A\right).\hfill \end{array}$$

(A12)

Accordingly, we obtain

$$\begin{array}{cc}\hfill & D_{\nu \mu }^{ab}({\omega }_A,{\omega }_B)F_{\mu \nu }^0=n_i\left[U_{ji}^a\left({\omega }_A\right)U_{ji}^{b*}(-{\omega }_B)+U_{ji}^b\left({\omega }_B\right)U_{ji}^{a*}(-{\omega }_A)\right]\left({\epsilon }_j^0-{\epsilon }_i^0\right)\hfill \\ \hfill & -n_i\left[U_{ji}^{a*}(-{\omega }_A)S_{ji}^b+U_{ji}^{b*}(-{\omega }_B)S_{ji}^a+S_{ii}^{ab}+S_{ji}^{a*}{U}_{ji}^b\left({\omega }_B\right)+S_{ji}^{b*}{U}_{ji}^a\left({\omega }_A\right)\right]{\epsilon }_i^0.\hfill \end{array}$$

(A13)

For further simplification, we separate the summation over j in the above expression into contributions from virtual and occupied orbitals (the latter can be transformed using Eq. (47)), and write

$$\begin{array}{cc}\hfill & D_{\nu \mu }^{ab}({\omega }_A,{\omega }_B)F_{\mu \nu }^0=n_i\left[U_{ji}^a\left({\omega }_A\right)U_{ji}^{b*}(-{\omega }_B)+U_{ji}^b\left({\omega }_B\right)U_{ji}^{a*}(-{\omega }_A)\right]\left({\epsilon }_j^0-{\epsilon }_i^0\right)\hfill \\ \hfill & -n_i\left[U_{ji}^{a*}(-{\omega }_A)S_{ji}^b+U_{ji}^{b*}(-{\omega }_B)S_{ji}^a+S_{ii}^{ab}+S_{ji}^{a*}{U}_{ji}^b\left({\omega }_B\right)+S_{ji}^{b*}{U}_{ji}^a\left({\omega }_A\right)\right]{\epsilon }_i^0\hfill \\ \hfill & +{\displaystyle \frac{n_i}{4}}\left[S_{ki}^a{S}_{ki}^{b*}+S_{ki}^b{S}_{ki}^{a*}\right]\left({\epsilon }_k^0+3{\epsilon }_i^0\right),\hfill \end{array}$$

(A14)

where $i,k\in \mathrm{occ}$ and $j\in \mathrm{vir}$.

The second term in Eq. (16) can be recast into

$$D_{\nu \mu }^a\left({\omega }_A\right)F_{\mu \nu }^b\left({\omega }_B\right)=n_i\left[U_{ji}^a\left({\omega }_A\right)F_{ji}^{b*}(-{\omega }_B)+U_{ji}^{a*}(-{\omega }_A)F_{ji}^b\left({\omega }_B\right)\right].$$

(A15)

Similarly, by separating the summation over j into contributions from virtual and occupied orbitals, we obtain

$$\begin{array}{cc}\hfill D_{\nu \mu }^a\left({\omega }_A\right)F_{\mu \nu }^b\left({\omega }_B\right)& =n_i\left[U_{ji}^a\left({\omega }_A\right)F_{ji}^{b*}(-{\omega }_B)+U_{ji}^{a*}(-{\omega }_A)F_{ji}^b\left({\omega }_B\right)\right]\hfill \\ \hfill & -{\displaystyle \frac{n_i}{2}}\left[S_{ki}^a{F}_{ki}^{b*}(-{\omega }_B)+S_{ki}^{a*}{F}_{ki}^b\left({\omega }_B\right)\right]\hfill \\ \hfill & =n_i{U}_{ji}^a\left({\omega }_A\right)\left[U_{ji}^{b*}(-{\omega }_B)\left({\epsilon }_i^0-{\epsilon }_j^0+{\omega }_B\right)+{\mathcal{T}}_{ji}^{b*}(-{\omega }_B)+S_{ji}^{b*}{\epsilon }_i^0\right]\hfill \\ \hfill & +n_i{U}_{ji}^{a*}(-{\omega }_A)\left[U_{ji}^b\left({\omega }_B\right)\left({\epsilon }_i^0-{\epsilon }_j^0-{\omega }_B\right)+{\mathcal{T}}_{ji}^b\left({\omega }_B\right)+S_{ji}^b{\epsilon }_i^0\right]\hfill \\ \hfill & -{\displaystyle \frac{n_i}{2}}\left[S_{ki}^a{F}_{ki}^{b*}(-{\omega }_B)+S_{ki}^{a*}{F}_{ki}^b\left({\omega }_B\right)\right],\hfill \end{array}$$

(A16)

where $i,k\in \mathrm{occ}$ and $j\in \mathrm{vir}$, and Eq. (44) has been fully utilized.

It can be seen that the combination of Eqs. (A14) and (A16) leads to partial cancelation of terms:

$$\begin{array}{cc}\hfill & D_{\nu \mu }^{ab}({\omega }_A,{\omega }_B)F_{\mu \nu }^0+D_{\nu \mu }^a\left({\omega }_A\right)F_{\mu \nu }^b\left({\omega }_B\right)\hfill \\ \hfill & =n_i{U}_{ji}^a\left({\omega }_A\right)\left[{\mathcal{T}}_{ji}^{b*}(-{\omega }_B)+{\omega }_B{U}_{ji}^{b*}(-{\omega }_B)\right]\hfill \\ \hfill & +n_i{U}_{ji}^{a*}(-{\omega }_A)\left[{\mathcal{T}}_{ji}^b\left({\omega }_B\right)-{\omega }_B{U}_{ji}^b\left({\omega }_B\right)\right]\hfill \\ \hfill & -n_i\left[U_{ji}^{b*}(-{\omega }_B)S_{ji}^a+S_{ii}^{ab}+S_{ji}^{a*}{U}_{ji}^b\left({\omega }_B\right)\right]{\epsilon }_i^0\hfill \\ \hfill & +{\displaystyle \frac{n_i}{4}}\left[S_{ki}^a{S}_{ki}^{b*}+S_{ki}^b{S}_{ki}^{a*}\right]\left({\epsilon }_k^0+3{\epsilon }_i^0\right)\hfill \\ \hfill & -{\displaystyle \frac{n_i}{2}}\left[S_{ki}^a{F}_{ki}^{b*}(-{\omega }_B)+S_{ki}^{a*}{F}_{ki}^b\left({\omega }_B\right)\right].\hfill \end{array}$$

(A17)

Consequently, the second-order energy response takes the form

$$\begin{array}{cc}\hfill & E^{ab}({\omega }_A,{\omega }_B)\hfill \\ \hfill & =n_i{U}_{ji}^a\left({\omega }_A\right)\left[{\mathcal{T}}_{ji}^{b*}(-{\omega }_B)+{\omega }_B{U}_{ji}^{b*}(-{\omega }_B)\right]+n_i{U}_{ji}^{a*}(-{\omega }_A)\left[{\mathcal{T}}_{ji}^b\left({\omega }_B\right)-{\omega }_B{U}_{ji}^b\left({\omega }_B\right)\right]\hfill \\ \hfill & -n_i\left[U_{ji}^{b*}(-{\omega }_B)S_{ji}^a+S_{ii}^{ab}+S_{ji}^{a*}{U}_{ji}^b\left({\omega }_B\right)\right]{\epsilon }_i^0+{\displaystyle \frac{n_i}{4}}\left[S_{ki}^a{S}_{ki}^{b*}+S_{ki}^b{S}_{ki}^{a*}\right]\left({\epsilon }_k^0+3{\epsilon }_i^0\right)\hfill \\ \hfill & -{\displaystyle \frac{n_i}{2}}\left[S_{ki}^a{F}_{ki}^{b*}(-{\omega }_B)+S_{ki}^{a*}{F}_{ki}^b\left({\omega }_B\right)\right]+D_{\nu \mu }^b\left({\omega }_B\right)H_{\mu \nu }^a+D_{\nu \mu }^0{H}_{\mu \nu }^{ab}\hfill \\ \hfill & +D_{\nu \mu }^b\left({\omega }_B\right)D_{\lambda \kappa }^0{\left(\mu \nu \right|\kappa \lambda )}^a+{\displaystyle \frac{1}{2}}{D}_{\nu \mu }^0{D}_{\lambda \kappa }^0{\left(\mu \nu \right|\kappa \lambda )}^{ab}+D_{\nu \mu }^b\left({\omega }_B\right)\left({\mu \nu }^a\right|v_{\mathrm{xc}})\hfill \\ \hfill & +D_{\nu \mu }^0\left({\mu \nu }^{ab}\right|v_{\mathrm{xc}})+D_{\nu \mu }^0({\mu \nu }^a|f_{\mathrm{xc}}\left|{\rho }^b\left({\omega }_B\right)\right).\hfill \end{array}$$

(A18)

However, upon closer inspection, one can show that the first-order terms in Eq. (16) can be reorganized as

$$\begin{array}{cc}\hfill & D_{\nu \mu }^a\left({\omega }_A\right)F_{\mu \nu }^b\left({\omega }_B\right)+D_{\nu \mu }^b\left({\omega }_B\right)H_{\mu \nu }^a+D_{\nu \mu }^b\left({\omega }_B\right)D_{\lambda \kappa }^0{\left(\mu \nu \right|\kappa \lambda )}^a\hfill \\ \hfill & +D_{\nu \mu }^b\left({\omega }_B\right)\left({\mu \nu }^a\right|v_{\mathrm{xc}})+D_{\nu \mu }^0({\mu \nu }^a|f_{\mathrm{xc}}\left|{\rho }^b\left({\omega }_B\right)\right)\hfill \\ \hfill & =D_{\nu \mu }^b\left({\omega }_B\right)F_{\mu \nu }^a\left({\omega }_A\right)+D_{\nu \mu }^a\left({\omega }_A\right)H_{\mu \nu }^b+D_{\nu \mu }^a\left({\omega }_A\right)D_{\lambda \kappa }^0{\left(\mu \nu \right|\kappa \lambda )}^b\hfill \\ \hfill & +D_{\nu \mu }^a\left({\omega }_A\right)\left({\mu \nu }^b\right|v_{\mathrm{xc}})+D_{\nu \mu }^0({\mu \nu }^b|f_{\mathrm{xc}}\left|{\rho }^a\left({\omega }_A\right)\right)\hfill \end{array}$$

(A19)

according to Eq. (39). Proceeding in the same manner, we obtain an alternative expression for the second-order energy response:

$$\begin{array}{cc}\hfill & E^{ab}({\omega }_A,{\omega }_B)\hfill \\ \hfill & =n_i{U}_{ji}^b\left({\omega }_B\right)\left[{\mathcal{T}}_{ji}^{a*}(-{\omega }_A)+{\omega }_A{U}_{ji}^{a*}(-{\omega }_A)\right]+n_i{U}_{ji}^{b*}(-{\omega }_B)\left[{\mathcal{T}}_{ji}^a\left({\omega }_A\right)-{\omega }_A{U}_{ji}^a\left({\omega }_A\right)\right]\hfill \\ \hfill & -n_i\left[U_{ji}^{a*}(-{\omega }_A)S_{ji}^b+S_{ii}^{ab}+S_{ji}^{b*}{U}_{ji}^a\left({\omega }_A\right)\right]{\epsilon }_i^0+{\displaystyle \frac{n_i}{4}}\left[S_{ki}^a{S}_{ki}^{b*}+S_{ki}^b{S}_{ki}^{a*}\right]\left({\epsilon }_k^0+3{\epsilon }_i^0\right)\hfill \\ \hfill & -{\displaystyle \frac{n_i}{2}}\left[S_{ki}^b{F}_{ki}^{a*}(-{\omega }_A)+S_{ki}^{b*}{F}_{ki}^a\left({\omega }_A\right)\right]+D_{\nu \mu }^a\left({\omega }_A\right)H_{\mu \nu }^b+D_{\nu \mu }^0{H}_{\mu \nu }^{ab}\hfill \\ \hfill & +D_{\nu \mu }^a\left({\omega }_A\right)D_{\lambda \kappa }^0{\left(\mu \nu \right|\kappa \lambda )}^b+{\displaystyle \frac{1}{2}}{D}_{\nu \mu }^0{D}_{\lambda \kappa }^0{\left(\mu \nu \right|\kappa \lambda )}^{ab}+D_{\nu \mu }^a\left({\omega }_A\right)\left({\mu \nu }^b\right|v_{\mathrm{xc}})\hfill \\ \hfill & +D_{\nu \mu }^0\left({\mu \nu }^{ab}\right|v_{\mathrm{xc}})+D_{\nu \mu }^0({\mu \nu }^b|f_{\mathrm{xc}}\left|{\rho }^a\left({\omega }_A\right)\right).\hfill \end{array}$$

(A20)

This expression is mathematically equivalent to Eq. (A18), but to make the symmetry under the interchange of differentiation order explicit, we take the average of the two forms in the final result:

$$\begin{array}{cc}\hfill & E^{ab}({\omega }_A,{\omega }_B)\hfill \\ \hfill & ={\displaystyle \frac{n_i}{2}}{U}_{ji}^a\left({\omega }_A\right)\left[{\mathcal{T}}_{ji}^{b*}(-{\omega }_B)-S_{ji}^{b*}{\epsilon }_i^0\right]+{\displaystyle \frac{n_i}{2}}{U}_{ji}^{a*}(-{\omega }_A)\left[{\mathcal{T}}_{ji}^b\left({\omega }_B\right)-S_{ji}^b{\epsilon }_i^0\right]\hfill \\ \hfill & +{\displaystyle \frac{n_i}{2}}{U}_{ji}^b\left({\omega }_B\right)\left[{\mathcal{T}}_{ji}^{a*}(-{\omega }_A)-S_{ji}^{a*}{\epsilon }_i^0\right]+{\displaystyle \frac{n_i}{2}}{U}_{ji}^{b*}(-{\omega }_B)\left[{\mathcal{T}}_{ji}^a\left({\omega }_A\right)-S_{ji}^a{\epsilon }_i^0\right]\hfill \\ \hfill & +{\displaystyle \frac{n_i}{2}}\left({\omega }_B-{\omega }_A\right)\left[U_{ji}^a\left({\omega }_A\right)U_{ji}^{b*}(-{\omega }_B)-U_{ji}^b\left({\omega }_B\right)U_{ji}^{a*}(-{\omega }_A)\right]\hfill \\ \hfill & -n_i{S}_{ii}^{ab}{\epsilon }_i^0+{\displaystyle \frac{n_i}{4}}\left[S_{ki}^a{S}_{ki}^{b*}+S_{ki}^b{S}_{ki}^{a*}\right]\left({\epsilon }_k^0+3{\epsilon }_i^0\right)\hfill \\ \hfill & -{\displaystyle \frac{n_i}{4}}\left[S_{ki}^a{F}_{ki}^{b*}(-{\omega }_B)+S_{ki}^{a*}{F}_{ki}^b\left({\omega }_B\right)+S_{ki}^b{F}_{ki}^{a*}(-{\omega }_A)+S_{ki}^{b*}{F}_{ki}^a\left({\omega }_A\right)\right]\hfill \\ \hfill & +{\displaystyle \frac{1}{2}}\left[D_{\nu \mu }^a\left({\omega }_A\right)H_{\mu \nu }^b+D_{\nu \mu }^b\left({\omega }_B\right)H_{\mu \nu }^a+D_{\nu \mu }^a\left({\omega }_A\right)D_{\lambda \kappa }^0{\left(\mu \nu \right|\kappa \lambda )}^b\right.\hfill \\ \hfill & +D_{\nu \mu }^b\left({\omega }_B\right)D_{\lambda \kappa }^0{\left(\mu \nu \right|\kappa \lambda )}^a+D_{\nu \mu }^a\left({\omega }_A\right)\left({\mu \nu }^b\right|v_{\mathrm{xc}})+D_{\nu \mu }^b\left({\omega }_B\right)\left({\mu \nu }^a\right|v_{\mathrm{xc}})\hfill \\ \hfill & +\left.D_{\nu \mu }^0({\mu \nu }^b|f_{\mathrm{xc}}|{\rho }^a\left({\omega }_A\right))+D_{\nu \mu }^0({\mu \nu }^a|f_{\mathrm{xc}}\left|{\rho }^b\left({\omega }_B\right)\right)\right]\hfill \\ \hfill & +D_{\nu \mu }^0{H}_{\mu \nu }^{ab}+{\displaystyle \frac{1}{2}}{D}_{\nu \mu }^0{D}_{\lambda \kappa }^0{\left(\mu \nu \right|\kappa \lambda )}^{ab}+D_{\nu \mu }^0\left({\mu \nu }^{ab}\right|v_{\mathrm{xc}}),\hfill \end{array}$$

(A21)

where we choose $i,k\in \mathrm{occ}$ and $j\in \mathrm{vir}$.

The second derivative of T is also simplifiable. From Eq. (19c), we find that

$$\begin{array}{cc}\hfill & \left[C_{\mu i}^{a*}(-{\omega }_A)S_{\mu \nu }^0+C_{\mu i}^{0*}{S}_{\mu \nu }^a\right]C_{\nu i}^b\left({\omega }_B\right)\hfill \\ \hfill =& \left[C_{\mu k}^{0*}{U}_{ki}^{a*}(-{\omega }_A)S_{\mu \nu }^0+C_{\mu i}^{0*}{S}_{\mu \nu }^a\right]C_{\nu j}^0{U}_{ji}^b\left({\omega }_B\right)\hfill \\ \hfill =& \left[{\delta }_{kj}{U}_{ki}^{a*}(-{\omega }_A)+S_{ij}^a\right]U_{ji}^b\left({\omega }_B\right)\hfill \\ \hfill =& \left[U_{ji}^{a*}(-{\omega }_A)+S_{ji}^{a*}\right]U_{ji}^b\left({\omega }_B\right).\hfill \end{array}$$

(A22)

Similarly,

$$\begin{array}{cc}\hfill & \left[C_{\mu i}^{b*}(-{\omega }_B)S_{\mu \nu }^0+C_{\mu i}^{0*}{S}_{\mu \nu }^b\right]C_{\nu i}^a\left({\omega }_A\right)\hfill \\ \hfill =& \left[U_{ji}^{b*}(-{\omega }_B)+S_{ji}^{b*}\right]U_{ji}^a\left({\omega }_A\right).\hfill \end{array}$$

(A23)

For the last term in Eq. (19c), we have

$$\begin{array}{cc}\hfill & C_{\mu i}^{0*}{S}_{\mu \nu }^0{C}_{\nu i}^{ab}({\omega }_A,{\omega }_B)\hfill \\ \hfill & =S_{ij}^0{U}_{ji}^{ab}({\omega }_A,{\omega }_B)\hfill \\ \hfill & ={\delta }_{ij}{U}_{ji}^{ab}({\omega }_A,{\omega }_B)\hfill \\ \hfill & =U_{ii}^{ab}({\omega }_A,{\omega }_B)\hfill \\ \hfill & =-{\displaystyle \frac{1}{2}}\left[U_{ji}^{a*}(-{\omega }_A)S_{ji}^b+U_{ji}^{a*}(-{\omega }_A)U_{ji}^b\left({\omega }_B\right)\right.\hfill \\ \hfill & +U_{ji}^{b*}(-{\omega }_B)S_{ji}^a+S_{ii}^{ab}+S_{ji}^{a*}{U}_{ji}^b\left({\omega }_B\right)\hfill \\ \hfill & +\left.U_{ji}^{b*}(-{\omega }_B)U_{ji}^a\left({\omega }_A\right)+S_{ji}^{b*}{U}_{ji}^a\left({\omega }_A\right)\right].\hfill \end{array}$$

(A24)

Hence, $T^{ab}({\omega }_A,{\omega }_B)$ can be cast into

$$\begin{array}{cc}\hfill & T^{ab}({\omega }_A,{\omega }_B)\hfill \\ \hfill & ={\displaystyle \frac{n_i}{2}}\left[U_{ji}^a\left({\omega }_A\right){\mathcal{T}}_{ji}^{b*}(-{\omega }_B)+U_{ji}^{a*}(-{\omega }_A){\mathcal{T}}_{ji}^b\left({\omega }_B\right)\right.\hfill \\ \hfill & +\left.U_{ji}^b\left({\omega }_B\right){\mathcal{T}}_{ji}^{a*}(-{\omega }_A)+U_{ji}^{b*}(-{\omega }_B){\mathcal{T}}_{ji}^a\left({\omega }_A\right)\right]\hfill \\ \hfill & +{\displaystyle \frac{1}{2}}\left[D_{\nu \mu }^a\left({\omega }_A\right){\mathcal{T}}_{\mu \nu }^b\left({\omega }_B\right)+D_{\nu \mu }^b\left({\omega }_B\right){\mathcal{T}}_{\mu \nu }^a\left({\omega }_A\right)\right]+D_{\nu \mu }^0{\mathcal{T}}_{\mu \nu }^{ab}({\omega }_A,{\omega }_B)\hfill \\ \hfill & +{\displaystyle \frac{n_i}{2}}\left({\omega }_A-{\omega }_B\right)\left\{\left[U_{ji}^{a*}(-{\omega }_A)+S_{ji}^{a*}\right]U_{ji}^b\left({\omega }_B\right)-\left[U_{ji}^{b*}(-{\omega }_B)+S_{ji}^{b*}\right]U_{ji}^a\left({\omega }_A\right)\right\}\hfill \\ \hfill & +{\displaystyle \frac{n_i}{2}}\left({\omega }_A+{\omega }_B\right)\left[U_{ji}^{a*}(-{\omega }_A)S_{ji}^b+S_{ii}^{ab}+U_{ji}^{b*}(-{\omega }_B)S_{ji}^a\right].\hfill \end{array}$$

(A25)

Splitting the summation over j into virtual and occupied orbitals gives

$$\begin{array}{cc}\hfill & T^{ab}({\omega }_A,{\omega }_B)\hfill \\ \hfill & ={\displaystyle \frac{n_i}{2}}\left[U_{ji}^a\left({\omega }_A\right){\mathcal{T}}_{ji}^{b*}(-{\omega }_B)+U_{ji}^{a*}(-{\omega }_A){\mathcal{T}}_{ji}^b\left({\omega }_B\right)\right.\hfill \\ \hfill & +\left.U_{ji}^b\left({\omega }_B\right){\mathcal{T}}_{ji}^{a*}(-{\omega }_A)+U_{ji}^{b*}(-{\omega }_B){\mathcal{T}}_{ji}^a\left({\omega }_A\right)\right]\hfill \\ \hfill & -{\displaystyle \frac{n_i}{4}}\left[S_{ki}^a{\mathcal{T}}_{ki}^{b*}(-{\omega }_B)+S_{ki}^{a*}{\mathcal{T}}_{ki}^b\left({\omega }_B\right)+S_{ki}^b{\mathcal{T}}_{ki}^{a*}(-{\omega }_A)+S_{ki}^b{\mathcal{T}}_{ki}^a\left({\omega }_A\right)\right]\hfill \\ \hfill & +{\displaystyle \frac{1}{2}}\left[D_{\nu \mu }^a\left({\omega }_A\right){\mathcal{T}}_{\mu \nu }^b\left({\omega }_B\right)+D_{\nu \mu }^b\left({\omega }_B\right){\mathcal{T}}_{\mu \nu }^a\left({\omega }_A\right)\right]+D_{\nu \mu }^0{\mathcal{T}}_{\mu \nu }^{ab}({\omega }_A,{\omega }_B)\hfill \\ \hfill & +{\displaystyle \frac{n_i}{2}}\left({\omega }_A-{\omega }_B\right)\left\{\left[U_{ji}^{a*}(-{\omega }_A)+S_{ji}^{a*}\right]U_{ji}^b\left({\omega }_B\right)-\left[U_{ji}^{b*}(-{\omega }_B)+S_{ji}^{b*}\right]U_{ji}^a\left({\omega }_A\right)\right\}\hfill \\ \hfill & +{\displaystyle \frac{n_i}{2}}\left({\omega }_A+{\omega }_B\right)\left[U_{ji}^{a*}(-{\omega }_A)S_{ji}^b+S_{ii}^{ab}+U_{ji}^{b*}(-{\omega }_B)S_{ji}^a\right]\hfill \\ \hfill & -{\displaystyle \frac{n_i}{8}}\left[\left(3{\omega }_A+{\omega }_B\right)S_{ki}^{a*}{S}_{ki}^b+\left(3{\omega }_B+{\omega }_A\right)S_{ki}^{b*}{S}_{ki}^a\right],\hfill \end{array}$$

(A26)

where $i,k\in \mathrm{occ}$ and $j\in \mathrm{vir}$.

Consequently, the second-order quasi-energy response can be written as Eq. (69).

Appendix B.3. The Third-Order W-Response

By analogy with previous derivations, we make use of the canonical form of the ground-state molecular orbitals and the orthonormality condition (35), which allows us to simplify the first term in the third-order energy response (17) as

$$\begin{array}{cc}\hfill & D_{\nu \mu }^{abc}({\omega }_A,{\omega }_B,{\omega }_C)F_{\mu \nu }^0\hfill \\ \hfill & =n_i\left[U_{ji}^{ab}({\omega }_A,{\omega }_B)U_{ji}^{c*}(-{\omega }_C)+U_{ji}^c\left({\omega }_C\right)U_{ji}^{ab*}(-{\omega }_A,-{\omega }_B)\right.\hfill \\ \hfill & +U_{ji}^{ac}({\omega }_A,{\omega }_C)U_{ji}^{b*}(-{\omega }_B)+U_{ji}^b\left({\omega }_B\right)U_{ji}^{ac*}(-{\omega }_A,-{\omega }_C)\hfill \\ \hfill & +\left.U_{ji}^{bc}({\omega }_B,{\omega }_C)U_{ji}^{a*}(-{\omega }_A)+U_{ji}^a\left({\omega }_A\right)U_{ji}^{bc*}(-{\omega }_B,-{\omega }_C)\right]\left({\epsilon }_j^0-{\epsilon }_i^0\right)\hfill \\ \hfill & -n_i\left[U_{ji}^{ab*}(-{\omega }_A,-{\omega }_B)S_{ji}^c+U_{ji}^{ac*}(-{\omega }_A,-{\omega }_C)S_{ji}^b+U_{ji}^{a*}(-{\omega }_A)S_{ji}^{bc}\right.\hfill \\ \hfill & +U_{ji}^{a*}(-{\omega }_A)S_{jk}^b{U}_{ki}^c\left({\omega }_C\right)+U_{ji}^{a*}(-{\omega }_A)S_{jk}^c{U}_{ki}^b\left({\omega }_B\right)+U_{ji}^{bc*}(-{\omega }_B,-{\omega }_C)S_{ji}^a\hfill \\ \hfill & +U_{ji}^{b*}(-{\omega }_B)S_{ji}^{ac}+U_{ji}^{b*}(-{\omega }_B)S_{jk}^a{U}_{ki}^c\left({\omega }_C\right)+U_{ji}^{c*}(-{\omega }_C)S_{ji}^{ab}\hfill \\ \hfill & +S_{ii}^{abc}+S_{ji}^{ab*}{U}_{ji}^c\left({\omega }_C\right)+U_{ji}^{c*}(-{\omega }_C)S_{jk}^a{U}_{ki}^b\left({\omega }_B\right)+S_{ji}^{ac*}{U}_{ji}^b\left({\omega }_B\right)\hfill \\ \hfill & +S_{ji}^{a*}{U}_{ji}^{bc}({\omega }_B,{\omega }_C)+U_{ji}^{b*}(-{\omega }_B)S_{jk}^c{U}_{ki}^a\left({\omega }_A\right)+U_{ji}^{c*}(-{\omega }_C)S_{jk}^b{U}_{ki}^a\left({\omega }_A\right)\hfill \\ \hfill & +\left.S_{ji}^{bc*}{U}_{ji}^a\left({\omega }_A\right)+S_{ji}^{b*}{U}_{ji}^{ac}({\omega }_A,{\omega }_C)+S_{ji}^{c*}{U}_{ji}^{ab}({\omega }_A,{\omega }_B)\right]{\epsilon }_i^0.\hfill \end{array}$$

(A27)

According to the definition of the second-order Fock matrix (40), we find in Eq. (17) that

$$\begin{array}{cc}\hfill & D_{\nu \mu }^a\left({\omega }_A\right)F_{\mu \nu }^{bc}({\omega }_B,{\omega }_C)+D_{\nu \mu }^{bc}({\omega }_B,{\omega }_C)\left[H_{\mu \nu }^a+D_{\lambda \kappa }^0{\left(\mu \nu \right|\kappa \lambda )}^a+\left({\mu \nu }^a\right|v_{\mathrm{xc}})\right]\hfill \\ \hfill & +D_{\nu \mu }^0\left[({\mu \nu }^a|f_{\mathrm{xc}}|{\rho }^{bc}({\omega }_B,{\omega }_C))+({\mu \nu }^a|k_{\mathrm{xc}}|{\rho }^b\left({\omega }_B\right)\left|{\rho }^c\left({\omega }_C\right)\right)\right]\hfill \\ \hfill & =D_{\nu \mu }^{bc}({\omega }_B,{\omega }_C)F_{\mu \nu }^a\left({\omega }_A\right)+D_{\nu \mu }^a\left({\omega }_A\right)\left[H_{\mu \nu }^{bc}+D_{\lambda \kappa }^b\left({\omega }_B\right){\left(\mu \nu \right|\kappa \lambda )}^c+D_{\lambda \kappa }^c\left({\omega }_C\right){\left(\mu \nu \right|\kappa \lambda )}^b\right.\hfill \\ \hfill & +\left.D_{\lambda \kappa }^0{\left(\mu \nu \right|\kappa \lambda )}^{bc}+\left({\mu \nu }^{bc}\right|v_{\mathrm{xc}})+({\mu \nu }^b|f_{\mathrm{xc}}|{\rho }^c\left({\omega }_C\right))+({\mu \nu }^c|f_{\mathrm{xc}}\left|{\rho }^b\left({\omega }_B\right)\right)\right]\hfill \\ \hfill & +D_{\nu \mu }^b\left({\omega }_B\right)({\mu \nu }^c|f_{\mathrm{xc}}|{\rho }^a\left({\omega }_A\right))+D_{\nu \mu }^c\left({\omega }_C\right)({\mu \nu }^b|f_{\mathrm{xc}}\left|{\rho }^a\left({\omega }_A\right)\right)\hfill \\ \hfill & +D_{\nu \mu }^0({\mu \nu }^{bc}|f_{\mathrm{xc}}|{\rho }^a\left({\omega }_A\right))+({\rho }^a\left({\omega }_A\right)|k_{\mathrm{xc}}|{\rho }^b\left({\omega }_B\right)\left|{\rho }^c\left({\omega }_C\right)\right),\hfill \end{array}$$

(A28)

where

$$\begin{array}{cc}\hfill & D_{\nu \mu }^{bc}({\omega }_B,{\omega }_C)F_{\mu \nu }^a\left({\omega }_A\right)\hfill \\ \hfill & =n_i\left[U_{ji}^{bc}({\omega }_B,{\omega }_C)F_{ji}^{a*}(-{\omega }_A)+U_{ji}^{bc*}(-{\omega }_B,-{\omega }_C)F_{ji}^a\left({\omega }_A\right)\right]\hfill \\ \hfill & +n_i{F}_{jl}^a\left({\omega }_A\right)\left[U_{li}^b\left({\omega }_B\right)U_{ji}^{c*}(-{\omega }_C)+U_{li}^c\left({\omega }_C\right)U_{ji}^{b*}(-{\omega }_B)\right].\hfill \end{array}$$

(A29)

The summation over j in the first brackets above can be partitioned into contributions from the occupied and virtual subspaces. For the occupied part, using the orthonormality condition (34), we have

$$\begin{array}{cc}\hfill & n_i\left[U_{ki}^{bc}({\omega }_B,{\omega }_C)F_{ki}^{a*}(-{\omega }_A)+U_{ki}^{bc*}(-{\omega }_B,-{\omega }_C)F_{ki}^a\left({\omega }_A\right)\right]\hfill \\ \hfill & =n_i{F}_{ki}^a\left({\omega }_A\right)\left[U_{ik}^{bc}({\omega }_B,{\omega }_C)+U_{ki}^{bc*}(-{\omega }_B,-{\omega }_C)\right]\hfill \\ \hfill & =-n_i{F}_{ki}^a\left({\omega }_A\right)\left[U_{ji}^{b*}(-{\omega }_B)S_{jk}^c+U_{ji}^{b*}(-{\omega }_B)U_{jk}^c\left({\omega }_C\right)+U_{ji}^{c*}(-{\omega }_C)S_{jk}^b\right.\hfill \\ \hfill & +\left.S_{ik}^{bc}+S_{ji}^{b*}{U}_{jk}^c\left({\omega }_C\right)+U_{ji}^{c*}(-{\omega }_C)U_{jk}^b\left({\omega }_B\right)+S_{ji}^{c*}{U}_{jk}^b\left({\omega }_B\right)\right],\hfill \end{array}$$

(A30)

where $i,k\in \mathrm{occ}$, and j runs over the full orbital space. For the virtual part, by virtue of the first-order CP-HF/KS equation (44), we obtain

$$\begin{array}{cc}\hfill & n_i\left[U_{ji}^{bc}({\omega }_B,{\omega }_C)F_{ji}^{a*}(-{\omega }_A)+U_{ji}^{bc*}(-{\omega }_B,-{\omega }_C)F_{ji}^a\left({\omega }_A\right)\right]\hfill \\ \hfill & =n_i{U}_{ji}^{bc}({\omega }_B,{\omega }_C)\left[U_{ji}^{a*}(-{\omega }_A)\left({\epsilon }_i^0-{\epsilon }_j^0+{\omega }_A\right)+{\mathcal{T}}_{ji}^{a*}(-{\omega }_A)+S_{ji}^{a*}{\epsilon }_i^0\right]\hfill \\ \hfill & +n_i{U}_{ji}^{bc*}(-{\omega }_B,-{\omega }_C)\left[U_{ji}^a\left({\omega }_A\right)\left({\epsilon }_i^0-{\epsilon }_j^0-{\omega }_A\right)+{\mathcal{T}}_{ji}^a\left({\omega }_A\right)+S_{ji}^a{\epsilon }_i^0\right],\hfill \end{array}$$

(A31)

where $i,k\in \mathrm{occ}$ and $j\in \mathrm{vir}$. Consequently, Eq. (A29) becomes

$$\begin{array}{cc}\hfill & D_{\nu \mu }^{bc}({\omega }_B,{\omega }_C)F_{\mu \nu }^a\left({\omega }_A\right)\hfill \\ \hfill & =n_i{U}_{ji}^{bc}({\omega }_B,{\omega }_C)\left[U_{ji}^{a*}(-{\omega }_A)\left({\epsilon }_i^0-{\epsilon }_j^0+{\omega }_A\right)+{\mathcal{T}}_{ji}^{a*}(-{\omega }_A)+S_{ji}^{a*}{\epsilon }_i^0\right]\hfill \\ \hfill & +n_i{U}_{ji}^{bc*}(-{\omega }_B,-{\omega }_C)\left[U_{ji}^a\left({\omega }_A\right)\left({\epsilon }_i^0-{\epsilon }_j^0-{\omega }_A\right)+{\mathcal{T}}_{ji}^a\left({\omega }_A\right)+S_{ji}^a{\epsilon }_i^0\right]\hfill \\ \hfill & -n_i{F}_{ki}^a\left({\omega }_A\right)\left[U_{ji}^{b*}(-{\omega }_B)S_{jk}^c+U_{ji}^{b*}(-{\omega }_B)U_{jk}^c\left({\omega }_C\right)+U_{ji}^{c*}(-{\omega }_C)S_{jk}^b\right.\hfill \\ \hfill & +\left.S_{ik}^{bc}+S_{ji}^{b*}{U}_{jk}^c\left({\omega }_C\right)+U_{ji}^{c*}(-{\omega }_C)U_{jk}^b\left({\omega }_B\right)+S_{ji}^{c*}{U}_{jk}^b\left({\omega }_B\right)\right]\hfill \\ \hfill & +n_i{F}_{jl}^a\left({\omega }_A\right)\left[U_{li}^b\left({\omega }_B\right)U_{ji}^{c*}(-{\omega }_C)+U_{li}^c\left({\omega }_C\right)U_{ji}^{b*}(-{\omega }_B)\right],\hfill \end{array}$$

(A32)

where $i,k\in \mathrm{occ}$, and the indice j is restricted to the virtual orbitals in the first two terms, but otherwise runs over all orbitals. Similarly, we have

$$\begin{array}{cc}\hfill & D_{\nu \mu }^{ac}({\omega }_A,{\omega }_C)F_{\mu \nu }^b\left({\omega }_B\right)\hfill \\ \hfill & =n_i{U}_{ji}^{ac}({\omega }_A,{\omega }_C)\left[U_{ji}^{b*}(-{\omega }_B)\left({\epsilon }_i^0-{\epsilon }_j^0+{\omega }_B\right)+{\mathcal{T}}_{ji}^{b*}(-{\omega }_B)+S_{ji}^{b*}{\epsilon }_i^0\right]\hfill \\ \hfill & +n_i{U}_{ji}^{ac*}(-{\omega }_A,-{\omega }_C)\left[U_{ji}^b\left({\omega }_B\right)\left({\epsilon }_i^0-{\epsilon }_j^0-{\omega }_B\right)+{\mathcal{T}}_{ji}^b\left({\omega }_B\right)+S_{ji}^b{\epsilon }_i^0\right]\hfill \\ \hfill & -n_i{F}_{ki}^b\left({\omega }_B\right)\left[U_{ji}^{a*}(-{\omega }_A)S_{jk}^c+U_{ji}^{a*}(-{\omega }_A)U_{jk}^c\left({\omega }_C\right)+U_{ji}^{c*}(-{\omega }_C)S_{jk}^a\right.\hfill \\ \hfill & +\left.S_{ik}^{ac}+S_{ji}^{a*}{U}_{jk}^c\left({\omega }_C\right)+U_{ji}^{c*}(-{\omega }_C)U_{jk}^a\left({\omega }_A\right)+S_{ji}^{c*}{U}_{jk}^a\left({\omega }_A\right)\right]\hfill \\ \hfill & +n_i{F}_{jl}^b\left({\omega }_B\right)\left[U_{li}^a\left({\omega }_A\right)U_{ji}^{c*}(-{\omega }_C)+U_{li}^c\left({\omega }_C\right)U_{ji}^{a*}(-{\omega }_A)\right],\hfill \end{array}$$

(A33)

and

$$\begin{array}{cc}\hfill & D_{\nu \mu }^{ab}({\omega }_A,{\omega }_B)F_{\mu \nu }^c\left({\omega }_C\right)\hfill \\ \hfill & =n_i{U}_{ji}^{ab}({\omega }_A,{\omega }_B)\left[U_{ji}^{c*}(-{\omega }_C)\left({\epsilon }_i^0-{\epsilon }_j^0+{\omega }_C\right)+{\mathcal{T}}_{ji}^{c*}(-{\omega }_C)+S_{ji}^{c*}{\epsilon }_i^0\right]\hfill \\ \hfill & +n_i{U}_{ji}^{ab*}(-{\omega }_A,-{\omega }_B)\left[U_{ji}^c\left({\omega }_C\right)\left({\epsilon }_i^0-{\epsilon }_j^0-{\omega }_C\right)+{\mathcal{T}}_{ji}^c\left({\omega }_C\right)+S_{ji}^c{\epsilon }_i^0\right]\hfill \\ \hfill & -n_i{F}_{ki}^c\left({\omega }_C\right)\left[U_{ji}^{a*}(-{\omega }_A)S_{jk}^b+U_{ji}^{a*}(-{\omega }_A)U_{jk}^b\left({\omega }_B\right)+U_{ji}^{b*}(-{\omega }_B)S_{jk}^a\right.\hfill \\ \hfill & +\left.S_{ik}^{ab}+S_{ji}^{a*}{U}_{jk}^b\left({\omega }_B\right)+U_{ji}^{b*}(-{\omega }_B)U_{jk}^a\left({\omega }_A\right)+S_{ji}^{b*}{U}_{jk}^a\left({\omega }_A\right)\right]\hfill \\ \hfill & +n_i{F}_{jl}^c\left({\omega }_C\right)\left[U_{li}^a\left({\omega }_A\right)U_{ji}^{b*}(-{\omega }_B)+U_{li}^b\left({\omega }_B\right)U_{ji}^{a*}(-{\omega }_A)\right].\hfill \end{array}$$

(A34)

By separating the summation over j in the second-order $\mathit{U}$-dependent terms of Eq. (A27) into contributions from virtual and occupied orbitals, we find that Eqs. (A27) and (A32)–(A34) can be combined in such a way that certain terms cancel out:

$$\begin{array}{cc}\hfill & D_{\nu \mu }^{abc}({\omega }_A,{\omega }_B,{\omega }_C)F_{\mu \nu }^0+D_{\nu \mu }^{ab}({\omega }_A,{\omega }_B)F_{\mu \nu }^c\left({\omega }_C\right)\hfill \\ \hfill & +D_{\nu \mu }^{ac}({\omega }_A,{\omega }_C)F_{\mu \nu }^b\left({\omega }_B\right)+D_{\nu \mu }^{bc}({\omega }_B,{\omega }_C)F_{\mu \nu }^a\left({\omega }_A\right)\hfill \\ \hfill & =n_i\left\{U_{ji}^{ab}({\omega }_A,{\omega }_B)\left[{\mathcal{T}}_{ji}^{c*}(-{\omega }_C)+{\omega }_C{U}_{ji}^{c*}(-{\omega }_C)\right]\right.\hfill \\ \hfill & +U_{ji}^{ab*}(-{\omega }_A,-{\omega }_B)\left[{\mathcal{T}}_{ji}^c\left({\omega }_C\right)-{\omega }_C{U}_{ji}^c\left({\omega }_C\right)\right]\hfill \\ \hfill & +U_{ji}^{ac}({\omega }_A,{\omega }_C)\left[{\mathcal{T}}_{ji}^{b*}(-{\omega }_B)+{\omega }_B{U}_{ji}^{b*}(-{\omega }_B)\right]\hfill \\ \hfill & +U_{ji}^{ac*}(-{\omega }_A,-{\omega }_C)\left[{\mathcal{T}}_{ji}^b\left({\omega }_B\right)-{\omega }_B{U}_{ji}^b\left({\omega }_B\right)\right]\hfill \\ \hfill & +U_{ji}^{bc}({\omega }_B,{\omega }_C)\left[{\mathcal{T}}_{ji}^{a*}(-{\omega }_A)+{\omega }_A{U}_{ji}^{a*}(-{\omega }_A)\right]\hfill \\ \hfill & +\left.U_{ji}^{bc*}(-{\omega }_B,-{\omega }_C)\left[{\mathcal{T}}_{ji}^a\left({\omega }_A\right)-{\omega }_A{U}_{ji}^a\left({\omega }_A\right)\right]\right\}\hfill \\ \hfill & +n_i{F}_{jl}^a\left({\omega }_A\right)\left[U_{li}^b\left({\omega }_B\right)U_{ji}^{c*}(-{\omega }_C)+U_{li}^c\left({\omega }_C\right)U_{ji}^{b*}(-{\omega }_B)\right]\hfill \\ \hfill & +n_i{F}_{jl}^b\left({\omega }_B\right)\left[U_{li}^a\left({\omega }_A\right)U_{ji}^{c*}(-{\omega }_C)+U_{li}^c\left({\omega }_C\right)U_{ji}^{a*}(-{\omega }_A)\right]\hfill \\ \hfill & +n_i{F}_{jl}^c\left({\omega }_C\right)\left[U_{li}^a\left({\omega }_A\right)U_{ji}^{b*}(-{\omega }_B)+U_{li}^b\left({\omega }_B\right)U_{ji}^{a*}(-{\omega }_A)\right]\hfill \\ \hfill & -n_i{F}_{ki}^a\left({\omega }_A\right)\left[U_{ji}^{b*}(-{\omega }_B)S_{jk}^c+U_{ji}^{b*}(-{\omega }_B)U_{jk}^c\left({\omega }_C\right)+U_{ji}^{c*}(-{\omega }_C)S_{jk}^b\right.\hfill \\ \hfill & +\left.S_{ik}^{bc}+S_{ji}^{b*}{U}_{jk}^c\left({\omega }_C\right)+U_{ji}^{c*}(-{\omega }_C)U_{jk}^b\left({\omega }_B\right)+S_{ji}^{c*}{U}_{jk}^b\left({\omega }_B\right)\right]\hfill \\ \hfill & -n_i{F}_{ki}^b\left({\omega }_B\right)\left[U_{ji}^{a*}(-{\omega }_A)S_{jk}^c+U_{ji}^{a*}(-{\omega }_A)U_{jk}^c\left({\omega }_C\right)+U_{ji}^{c*}(-{\omega }_C)S_{jk}^a\right.\hfill \\ \hfill & +\left.S_{ik}^{ac}+S_{ji}^{a*}{U}_{jk}^c\left({\omega }_C\right)+U_{ji}^{c*}(-{\omega }_C)U_{jk}^a\left({\omega }_A\right)+S_{ji}^{c*}{U}_{jk}^a\left({\omega }_A\right)\right]\hfill \\ \hfill & -n_i{F}_{ki}^c\left({\omega }_C\right)\left[U_{ji}^{a*}(-{\omega }_A)S_{jk}^b+U_{ji}^{a*}(-{\omega }_A)U_{jk}^b\left({\omega }_B\right)+U_{ji}^{b*}(-{\omega }_B)S_{jk}^a\right.\hfill \\ \hfill & +\left.S_{ik}^{ab}+S_{ji}^{a*}{U}_{jk}^b\left({\omega }_B\right)+U_{ji}^{b*}(-{\omega }_B)U_{jk}^a\left({\omega }_A\right)+S_{ji}^{b*}{U}_{jk}^a\left({\omega }_A\right)\right]\hfill \\ \hfill & -{\displaystyle \frac{n_i}{2}}\left[U_{ki}^{ab}({\omega }_A,{\omega }_B)S_{ki}^{c*}+U_{ki}^{ab*}(-{\omega }_A,-{\omega }_B)S_{ki}^c+U_{ki}^{ac}({\omega }_A,{\omega }_C)S_{ki}^{b*}\right.\hfill \\ \hfill & +\left.U_{ki}^{ac*}(-{\omega }_A,-{\omega }_C)S_{ki}^b+U_{ki}^{bc}({\omega }_B,{\omega }_C)S_{ki}^{a*}+U_{ki}^{bc*}(-{\omega }_B,-{\omega }_C)S_{ki}^a\right]\left({\epsilon }_k^0+{\epsilon }_i^0\right)\hfill \\ \hfill & -n_i\left[U_{ji}^{a*}(-{\omega }_A)S_{ji}^{bc}+U_{ji}^{a*}(-{\omega }_A)S_{jl}^b{U}_{li}^c\left({\omega }_C\right)+U_{ji}^{a*}(-{\omega }_A)S_{jl}^c{U}_{li}^b\left({\omega }_B\right)\right.\hfill \\ \hfill & +U_{ji}^{b*}(-{\omega }_B)S_{ji}^{ac}+U_{ji}^{b*}(-{\omega }_B)S_{jl}^a{U}_{li}^c\left({\omega }_C\right)+U_{ji}^{c*}(-{\omega }_C)S_{ji}^{ab}+S_{ii}^{abc}\hfill \\ \hfill & +S_{ji}^{ab*}{U}_{ji}^c\left({\omega }_C\right)+U_{ji}^{c*}(-{\omega }_C)S_{jl}^a{U}_{li}^b\left({\omega }_B\right)+S_{ji}^{ac*}{U}_{ji}^b\left({\omega }_B\right)\hfill \\ \hfill & +\left.U_{ji}^{b*}(-{\omega }_B)S_{jl}^c{U}_{li}^a\left({\omega }_A\right)+U_{ji}^{c*}(-{\omega }_C)S_{jl}^b{U}_{li}^a\left({\omega }_A\right)+S_{ji}^{bc*}{U}_{ji}^a\left({\omega }_A\right)\right]{\epsilon }_i^0,\hfill \end{array}$$

(A35)

with $i,k$ restricted to occupied orbitals, and $j,l$ running over all orbitals except in the first term in braces, where the index j refers only to virtual orbitals. Therefore, the third-order energy response reads

$$\begin{array}{cc}\hfill & E^{abc}({\omega }_A,{\omega }_B,{\omega }_C)\hfill \\ \hfill & =D_{\nu \mu }^{abc}({\omega }_A,{\omega }_B,{\omega }_C)F_{\mu \nu }^0+D_{\nu \mu }^{ab}({\omega }_A,{\omega }_B)F_{\mu \nu }^c\left({\omega }_C\right)+D_{\nu \mu }^{ac}({\omega }_A,{\omega }_C)F_{\mu \nu }^b\left({\omega }_B\right)\hfill \\ \hfill & +D_{\nu \mu }^{bc}({\omega }_B,{\omega }_C)F_{\mu \nu }^a\left({\omega }_A\right)+D_{\nu \mu }^a\left({\omega }_A\right)D_{\lambda \kappa }^b\left({\omega }_B\right){\left(\mu \nu \right|\kappa \lambda )}^c\hfill \\ \hfill & +D_{\nu \mu }^a\left({\omega }_A\right)D_{\lambda \kappa }^c\left({\omega }_C\right){\left(\mu \nu \right|\kappa \lambda )}^b+D_{\nu \mu }^b\left({\omega }_B\right)D_{\lambda \kappa }^c\left({\omega }_C\right){\left(\mu \nu \right|\kappa \lambda )}^a\hfill \\ \hfill & +D_{\nu \mu }^a\left({\omega }_A\right)\left[H_{\mu \nu }^{bc}+\left({\mu \nu }^{bc}\right|v_{\mathrm{xc}})+({\mu \nu }^b|f_{\mathrm{xc}}|{\rho }^c\left({\omega }_C\right))+({\mu \nu }^c|f_{\mathrm{xc}}\left|{\rho }^b\left({\omega }_B\right)\right)\right]\hfill \\ \hfill & +D_{\nu \mu }^b\left({\omega }_B\right)\left[H_{\mu \nu }^{ac}+\left({\mu \nu }^{ac}\right|v_{\mathrm{xc}})+({\mu \nu }^a|f_{\mathrm{xc}}|{\rho }^c\left({\omega }_C\right))+({\mu \nu }^c|f_{\mathrm{xc}}\left|{\rho }^a\left({\omega }_A\right)\right)\right]\hfill \\ \hfill & +D_{\nu \mu }^c\left({\omega }_C\right)\left[H_{\mu \nu }^{ab}+\left({\mu \nu }^{bc}\right|v_{\mathrm{xc}})+({\mu \nu }^a|f_{\mathrm{xc}}|{\rho }^b\left({\omega }_B\right))+({\mu \nu }^b|f_{\mathrm{xc}}\left|{\rho }^a\left({\omega }_A\right)\right)\right]\hfill \\ \hfill & +D_{\lambda \kappa }^0\left[D_{\nu \mu }^a\left({\omega }_A\right){\left(\mu \nu \right|\kappa \lambda )}^{bc}+D_{\nu \mu }^b\left({\omega }_B\right){\left(\mu \nu \right|\kappa \lambda )}^{ac}+D_{\nu \mu }^c\left({\omega }_C\right){\left(\mu \nu \right|\kappa \lambda )}^{ab}\right]\hfill \\ \hfill & +D_{\nu \mu }^0\left[({\mu \nu }^{bc}|f_{\mathrm{xc}}|{\rho }^a\left({\omega }_A\right))+({\mu \nu }^{ac}|f_{\mathrm{xc}}|{\rho }^b\left({\omega }_B\right))+({\mu \nu }^{ab}|f_{\mathrm{xc}}\left|{\rho }^c\left({\omega }_C\right)\right)\right]\hfill \\ \hfill & +D_{\nu \mu }^0\left[H_{\mu \nu }^{abc}+{\displaystyle \frac{1}{2}}{D}_{\lambda \kappa }^0{\left(\mu \nu \right|\kappa \lambda )}^{abc}+\left({\mu \nu }^{abc}\right|v_{\mathrm{xc}})\right]+({\rho }^a\left({\omega }_A\right)|k_{\mathrm{xc}}|{\rho }^b\left({\omega }_B\right)\left|{\rho }^c\left({\omega }_C\right)\right).\hfill \end{array}$$

(A36)

For $T^{abc}({\omega }_A,{\omega }_B,{\omega }_C)$ as given in Eq. (19d), one can observe that

$$\begin{array}{cc}\hfill & \left[C_{\mu i}^{bc*}(-{\omega }_B,-{\omega }_C)S_{\mu \nu }^0+C_{\mu i}^{b*}(-{\omega }_B)S_{\mu \nu }^c+C_{\mu i}^{c*}(-{\omega }_C)S_{\mu \nu }^b+C_{\mu i}^{0*}{S}_{\mu \nu }^{bc}\right]C_{\nu i}^a\left({\omega }_A\right)\hfill \\ \hfill =& \left[U_{ki}^{bc*}(-{\omega }_B,-{\omega }_C)+U_{ji}^{b*}(-{\omega }_B)S_{jk}^c+U_{ji}^{c*}(-{\omega }_C)S_{jk}^b+S_{ki}^{bc*}\right]U_{ki}^a\left({\omega }_A\right),\hfill \end{array}$$

(A37a)

$$\begin{array}{cc}\hfill & \left[C_{\mu i}^{ac*}(-{\omega }_A,-{\omega }_C)S_{\mu \nu }^0+C_{\mu i}^{a*}(-{\omega }_A)S_{\mu \nu }^c+C_{\mu i}^{c*}(-{\omega }_C)S_{\mu \nu }^a+C_{\mu i}^{0*}{S}_{\mu \nu }^{ac}\right]C_{\nu i}^b\left({\omega }_B\right)\hfill \\ \hfill =& \left[U_{ki}^{ac*}(-{\omega }_A,-{\omega }_C)+U_{ji}^{a*}(-{\omega }_A)S_{jk}^c+U_{ji}^{c*}(-{\omega }_C)S_{jk}^a+S_{ki}^{ac*}\right]U_{ki}^b\left({\omega }_B\right),\hfill \end{array}$$

(A37b)

$$\begin{array}{cc}\hfill & \left[C_{\mu i}^{ab*}(-{\omega }_A,-{\omega }_B)S_{\mu \nu }^0+C_{\mu i}^{a*}(-{\omega }_A)S_{\mu \nu }^b+C_{\mu i}^{b*}(-{\omega }_B)S_{\mu \nu }^a+C_{\mu i}^{0*}{S}_{\mu \nu }^{ab}\right]C_{\nu i}^c\left({\omega }_C\right)\hfill \\ \hfill =& \left[U_{ki}^{ab*}(-{\omega }_A,-{\omega }_B)+U_{ji}^{a*}(-{\omega }_A)S_{jk}^b+U_{ji}^{b*}(-{\omega }_B)S_{jk}^a+S_{ki}^{ab*}\right]U_{ki}^c\left({\omega }_C\right).\hfill \end{array}$$

(A37c)

Similarly, we can obtain

$$\begin{array}{cc}\hfill & \left[C_{\mu i}^{c*}(-{\omega }_C)S_{\mu \nu }^0+C_{\mu i}^{0*}{S}_{\mu \nu }^c\right]C_{\nu i}^{ab}({\omega }_A,{\omega }_B)\hfill \\ \hfill =& \left[U_{ji}^{c*}(-{\omega }_C)+S_{ji}^{c*}\right]U_{ji}^{ab}({\omega }_A,{\omega }_B),\hfill \end{array}$$

(A38a)

$$\begin{array}{cc}\hfill & \left[C_{\mu i}^{b*}(-{\omega }_B)S_{\mu \nu }^0+C_{\mu i}^{0*}{S}_{\mu \nu }^b\right]C_{\nu i}^{ac}({\omega }_A,{\omega }_C)\hfill \\ \hfill =& \left[U_{ji}^{b*}(-{\omega }_B)+S_{ji}^{b*}\right]U_{ji}^{ac}({\omega }_A,{\omega }_C),\hfill \end{array}$$

(A38b)

$$\begin{array}{cc}\hfill & \left[C_{\mu i}^{a*}(-{\omega }_A)S_{\mu \nu }^0+C_{\mu i}^{0*}{S}_{\mu \nu }^a\right]C_{\nu i}^{bc}({\omega }_B,{\omega }_C)\hfill \\ \hfill =& \left[U_{ji}^{a*}(-{\omega }_A)+S_{ji}^{a*}\right]U_{ji}^{bc}({\omega }_B,{\omega }_C).\hfill \end{array}$$

(A38c)

The third-order coefficients $C_{\nu i}^{abc}({\omega }_A,{\omega }_B,{\omega }_C)$ appearing in $C_{\mu i}^{0*}{S}_{\mu \nu }^0{C}_{\nu i}^{abc}({\omega }_A,{\omega }_B,{\omega }_C)$ deserve particular attention, as this term can be recast as

$$\begin{array}{cc}\hfill & C_{\mu i}^{0*}{S}_{\mu \nu }^0{C}_{\nu i}^{abc}({\omega }_A,{\omega }_B,{\omega }_C)\hfill \\ \hfill & =S_{ij}^0{U}_{ji}^{abc}({\omega }_A,{\omega }_B,{\omega }_C)\hfill \\ \hfill & ={\delta }_{ij}{U}_{ji}^{abc}({\omega }_A,{\omega }_B,{\omega }_C)\hfill \\ \hfill & =U_{ii}^{abc}({\omega }_A,{\omega }_B,{\omega }_C),\hfill \end{array}$$

(A39)

to which we need to apply a further assumption. Imposing the relation $U_{ki}^{abc}(-{\omega }_A,-{\omega }_B,-{\omega }_C)=U_{ik}^{abc}({\omega }_A,{\omega }_B,{\omega }_C)$ for the occ-occ block gives

$$\begin{array}{cc}\hfill & U_{ii}^{abc}({\omega }_A,{\omega }_B,{\omega }_C)\hfill \\ \hfill & =-{\displaystyle \frac{1}{2}}\left[U_{ji}^{ab*}(-{\omega }_A,-{\omega }_B)S_{ji}^c+U_{ji}^{ab*}(-{\omega }_A,-{\omega }_B)U_{ji}^c\left({\omega }_C\right)+U_{ji}^{ac*}(-{\omega }_A,-{\omega }_C)S_{ji}^b\right.\hfill \\ \hfill & +U_{ji}^{a*}(-{\omega }_A)S_{ji}^{bc}+U_{ji}^{a*}(-{\omega }_A)S_{jk}^b{U}_{ki}^c\left({\omega }_C\right)+U_{ji}^{ac*}(-{\omega }_A,-{\omega }_C)U_{ji}^b\left({\omega }_B\right)\hfill \\ \hfill & +U_{ji}^{a*}(-{\omega }_A)S_{jk}^c{U}_{ki}^b\left({\omega }_B\right)+U_{ji}^{a*}(-{\omega }_A)U_{ji}^{bc}({\omega }_B,{\omega }_C)+U_{ji}^{bc*}(-{\omega }_B,-{\omega }_C)S_{ji}^a\hfill \\ \hfill & +U_{ji}^{b*}(-{\omega }_B)S_{ji}^{ac}+U_{ji}^{b*}(-{\omega }_B)S_{jk}^a{U}_{ki}^c\left({\omega }_C\right)+U_{ji}^{c*}(-{\omega }_C)S_{ji}^{ab}\hfill \\ \hfill & +S_{ii}^{abc}+S_{ji}^{ab*}{U}_{ji}^c\left({\omega }_C\right)+U_{ji}^{c*}(-{\omega }_C)S_{jk}^a{U}_{ki}^b\left({\omega }_B\right)+S_{ji}^{ac*}{U}_{ji}^b\left({\omega }_B\right)\hfill \\ \hfill & +S_{ji}^{a*}{U}_{ji}^{bc}({\omega }_B,{\omega }_C)+U_{ji}^{bc*}(-{\omega }_B,-{\omega }_C)U_{ji}^a\left({\omega }_A\right)+U_{ji}^{b*}(-{\omega }_B)S_{jk}^c{U}_{ki}^a\left({\omega }_A\right)\hfill \\ \hfill & +U_{ji}^{b*}(-{\omega }_B)U_{ji}^{ac}({\omega }_A,{\omega }_C)+U_{ji}^{c*}(-{\omega }_C)S_{jk}^b{U}_{ki}^a\left({\omega }_A\right)+S_{ji}^{bc*}{U}_{ji}^a\left({\omega }_A\right)\hfill \\ \hfill & +\left.S_{ji}^{b*}{U}_{ji}^{ac}({\omega }_A,{\omega }_C)+U_{ji}^{c*}(-{\omega }_C)U_{ji}^{ab}({\omega }_A,{\omega }_B)+S_{ji}^{c*}{U}_{ji}^{ab}({\omega }_A,{\omega }_B)\right].\hfill \end{array}$$

(A40)

Therefore, $T^{abc}({\omega }_A,{\omega }_B,{\omega }_C)$ finally takes the form

$$\begin{array}{cc}\hfill & T^{abc}({\omega }_A,{\omega }_B,{\omega }_C)\hfill \\ \hfill & =\left[D_{\nu \mu }^{ab}({\omega }_A,{\omega }_B){\mathcal{T}}_{\mu \nu }^c\left({\omega }_C\right)+D_{\nu \mu }^{ac}({\omega }_A,{\omega }_C){\mathcal{T}}_{\mu \nu }^b\left({\omega }_B\right)+D_{\nu \mu }^{bc}({\omega }_B,{\omega }_C){\mathcal{T}}_{\mu \nu }^a\left({\omega }_A\right)\right.\hfill \\ \hfill & +D_{\nu \mu }^a\left({\omega }_A\right){\mathcal{T}}_{\mu \nu }^{bc}({\omega }_B,{\omega }_C)+D_{\nu \mu }^b\left({\omega }_B\right){\mathcal{T}}_{\mu \nu }^{ac}({\omega }_A,{\omega }_C)+D_{\nu \mu }^c\left({\omega }_C\right){\mathcal{T}}_{\mu \nu }^{ab}({\omega }_A,{\omega }_B)\hfill \\ \hfill & +\left.D_{\nu \mu }^0{\mathcal{T}}_{\mu \nu }^{abc}({\omega }_A,{\omega }_B,{\omega }_C)\right]-{\omega }_A{n}_i\left[U_{ki}^{bc*}(-{\omega }_B,-{\omega }_C)+U_{ji}^{b*}(-{\omega }_B)S_{jk}^c\right.\hfill \\ \hfill & +\left.U_{ji}^{c*}(-{\omega }_C)S_{jk}^b+S_{ki}^{bc*}\right]U_{ki}^a\left({\omega }_A\right)-{\omega }_B{n}_i\left[U_{ki}^{ac*}(-{\omega }_A,-{\omega }_C)+U_{ji}^{a*}(-{\omega }_A)S_{jk}^c\right.\hfill \\ \hfill & +\left.U_{ji}^{c*}(-{\omega }_C)S_{jk}^a+S_{ki}^{ac*}\right]U_{ki}^b\left({\omega }_B\right)-{\omega }_C{n}_i\left[U_{ki}^{ab*}(-{\omega }_A,-{\omega }_B)+U_{ji}^{a*}(-{\omega }_A)S_{jk}^b\right.\hfill \\ \hfill & +\left.U_{ji}^{b*}(-{\omega }_B)S_{jk}^a+S_{ki}^{ab*}\right]U_{ki}^c\left({\omega }_C\right)-\left({\omega }_A+{\omega }_B\right)n_i\left[U_{ji}^{c*}(-{\omega }_C)+S_{ji}^{c*}\right]U_{ji}^{ab}({\omega }_A,{\omega }_B)\hfill \\ \hfill & -\left({\omega }_A+{\omega }_C\right)n_i\left[U_{ji}^{b*}(-{\omega }_B)+S_{ji}^{b*}\right]U_{ji}^{ac}({\omega }_A,{\omega }_C)-\left({\omega }_B+{\omega }_C\right)n_i\left[U_{ji}^{a*}(-{\omega }_A)\right.\hfill \\ \hfill & +\left.S_{ji}^{a*}\right]U_{ji}^{bc}({\omega }_B,{\omega }_C)-\left({\omega }_A+{\omega }_B+{\omega }_C\right)n_i{U}_{ii}^{abc}({\omega }_A,{\omega }_B,{\omega }_C).\hfill \end{array}$$

(A41)

Subtracting Eq. (A41) from Eq. (A36), one arrives at Eq. (70).