Transition densities between excited states are key for nonlinear theoretical spectroscopy and multi-state non-adiabatic molecular dynamics (NAMD) simulations. In the framework of response theory, these transition densities are accessible from poles of the quadratic response function. It was shown recently that the thus obtained transition densities within time-dependent Hartree-Fock (TDHF) and adiabatic time-dependent density functional theory (TDDFT) exhibit unphysical divergences when the difference in excitation energy of the two states of interest matches another excitation energy. This unphysical behavior is a consequence of spurious poles in the quadratic response function. We show that the incorrect pole structure of the quadratic response is not limited to TDHF and adiabatic TDDFT, but is also present in many other approximate many-electron response functions, including those from coupled cluster and multiconfigurational self-consistent field response theory. The divergences appear in regions of the potential energy surface where the ground state is perfectly well behaved, and they are frequently encountered in NAMD simulations of photochemical reactions. The origin of the divergences is traced to an incorrect instantaneous time-dependence of the effective Hamiltonian. The implications for computations of frequency-dependent response properties are considerable and call into question the validity of conventional approximate many-electron response theories beyond linear response.

## I. INTRODUCTION

The prediction of electronic properties such as electrostatic moments, polarizabilities, molecular structures, as well as transition properties determining spectral intensities is a central goal of electronic structure theory. According to the principles of quantum mechanics,^{1} the expectation value of an observable property *O* associated with the self-adjoint operator $O\u02c6(t)$ is

where |Ψ(*t*)〉 denotes the time-dependent electronic wavefunction. In practice, (1) is less useful than one might expect because it assumes |Ψ(*t*)〉 to be the exact time-dependent wavefunction evolving according to the time-dependent Schrödinger equation

where $H\u02c6$ is the electronic Hamiltonian and atomic (Hartree) units are used throughout. For approximate wavefunctions, Eq. (1) yields inconsistent results; for example, dipole moments obtained from Eq. (1) disagree with energy derivatives in second-order Møller-Plesset (MP2) theory.^{2} In other settings such as density functional theory (DFT), the interacting many-electron wavefunction is entirely unavailable.

The limitations of Eq. (1) are overcome by response theory, which defines all properties by the response of the energy or action functional. According to Hellmann and Feynman,^{3,4} expectation values of variational wavefunctions can alternatively be obtained from energy derivatives. In the time-dependent case, the time-dependent wavefunction stationarizes the action functional^{5}

Time-dependent properties are defined by response of the action to a time-dependent perturbation $\lambda (t)O\u02c6(t)$ added to the Hamiltonian. By the time-dependent Hellmann-Feynman theorem^{6}

which recovers the expectation value of $O\u02c6(t)$. For example, the time-dependent dipole moment is the functional derivative of the action with respect to a spatially uniform time-dependent electric field at vanishing field strength.

The class of properties accessible by response theory is broad and goes beyond the electronic ground state. Eq. (4) includes the definition of the one-particle density matrix as the functional derivative of $A$ with respect to a one-particle perturbation. Second and higher order derivatives of the action with respect to multiple perturbations correspond to linear, quadratic, and higher order response properties. Excited state properties and transition properties are obtainable from the residues of the frequency-dependent linear and higher order response of the electronic ground state.^{7} For example, the frequency-dependent polarizability of a molecule has poles at electronic excitation energies, and the residues are related to transition moments.^{8–10} The quadratic response function provides access to excited state properties, transition moments between excited states, and two-photon transition moments;^{7,8,11,12} likewise, higher-order response functions contain multi-photon transition moments and excited state response properties.^{8,10,13,14}

As opposed to state-specific approaches, response theory does not require explicit knowledge of excited state wavefunctions. Nevertheless, excited state properties from response theory are often more accurate than those computed via Eq. (1) using approximate wavefunctions.^{15,16} Moreover, properties of many states can be obtained efficiently from response theory at once,^{17} and fundamental sum rules relating excited and ground state properties can be satisfied;^{18,19} achieving both simultaneously is difficult with state-specific approaches.

Due to its efficiency, consistency, and broad applicability, response theory has become the method of choice for on-the-fly non-adiabatic molecular dynamics (NAMD) simulations.^{20,21} Each time step of such simulations requires knowledge of ground and excited state energies and nuclear forces as well as derivative couplings between all electronic states. Separate calculations for each of the electronic states involved in a simulation would greatly increase the cost of a single time step. However, within response theory, the excitation energies and ground-to-excited state couplings are obtained simultaneously from the linear response function^{22} while couplings between all pairs of excited states (state-to-state) are obtained from the quadratic response function.^{23}

Recently, several groups independently reported^{23–25} a troubling aspect of state-to-state properties derived from quadratic response of time-dependent Hartree–Fock (TDHF) and adiabatic time-dependent DFT (TDDFT): transition properties between two excited states *M*, *N* with excitation energies Ω_{M} = *E _{M}* −

*E*

_{0}, Ω

_{N}=

*E*−

_{N}*E*

_{0}diverge whenever the difference in energy between the two states, Ω

_{MN}= Ω

_{M}− Ω

_{N}, matches the excitation energy from the ground state to

*any other*excited state,

*K*. In other words, transition properties become strongly unphysical whenever

for any other state *K*. This matching condition does not correspond to any physical degeneracy or to any instability in the reference and can be encountered even when the reference and excited state energies are otherwise well-behaved.

These unphysical state-to-state transition properties can be traced to the incorrect pole structure of the approximate quadratic response function. Deficiencies in the pole structure of TDHF quadratic response functions were first noted in 1982 by Dalgaard.^{26} In particular, the TDHF quadratic response function contains terms proportional to a product of three poles [i.e., ∝1/(Ω_{N} − *ω*_{α})(Ω_{M} − *ω*_{β})(Ω_{K} − *ω*_{γ})] whereas all terms in the exact quadratic response function have products of at most two poles. Furthermore, the TDHF quadratic response function has terms proportional to, for example, 1/(Ω_{N} − *ω*_{α})(Ω_{M} − *ω*_{β}) that have no counterpart in exact theory. Since state-to-state properties are obtained as double residues of the response functions, the spurious third pole results in a single pole in these properties which in turn causes the unphysical divergences reported previously.

Here, we show that existence of the spurious pole is not unique to TDHF and adiabatic TDDFT, but it is also present in the response theories of many approximate electronic structure methods, including the response theories based on multi-configuration self-consistent field (MCSCF) and coupled-cluster (CC) references. We identify the artifacts that result from the incorrect pole structure of the approximate response functions and discuss the conditions under which these artifacts are most problematic. The electronic structure methods and their respective response theories, including sum-over-states expressions for all response functions and residues considered here, have been presented previously and extensively studied over the last three decades.^{7,26,27} However, a discussion of the impact of the spurious poles on properties of interest appears to be missing.

This paper is organized as follows: in Section II we review the response of a variational reference up through quadratic response in order to introduce response functions and residues of exact electronic eigenstates. We then use the response functions and residues from exact states to define transition moments and excited state properties in approximate variational methods and for CC. Section III provides numerical examples where the approximate response functions yield unphysical results. Section IV elaborates on the origins of the spurious poles. Section V discusses several avenues towards resolving the problem. Finally, we close in Section VI by discussing implications for the application of response theories.

## II. RESPONSE THEORY

In this section, we summarize the general working equations for response theory through second order and outline their derivation. For comprehensive derivations, the reader is referred to Refs. 28 and 29 in a general context, or to Refs. 7 and 26 for derivations in the TDHF and MCSCF contexts, respectively.

We start by partitioning the many-electron Hamiltonian into

where $H\u02c60$ is the time-independent field-free contribution; time-dependent perturbations are of the form

with *λ*_{α}, *ω*_{α}, and $v\u02c6(\alpha )$ the strength, frequency, and potential associated with the *α* component of the perturbation. In keeping with the Hermiticity of $v\u02c6\lambda $, the *λ*_{α} and *ω*_{α} are real and $(v\u02c6(\alpha ))\u2020=v\u02c6(\u2212\alpha )$. While the form of the field-free effective Hamiltonian will depend on the underlying electronic structure method, Eq. (7) is valid throughout this work.

### A. Response of a variational wavefunction

First, consider the response of a general variational wavefunction, which includes as special cases exact electronic states, TDHF, adiabatic TDDFT, and MCSCF response theory. The general working equations presented here are equally applicable to other variational methods, although the precise formulae provided for some of the matrix elements may be incomplete. For example, adiabatic TDDFT matrix elements should be augmented to include exchange-correlation functional derivatives.^{30} As the present focus is on the general features of approximate response theories, we refer the reader to the relevant literature for method-specific derivations.^{7,26,30}

The time-dependent perturbation induces timedependence in the reference which can be parametrized as a single-exponential

or a double-exponential

where $\kappa \u02c6(t)$ are anti-Hermitian operators that generate rotations between the static reference electronic state, |0〉, and other electronic states. The type of reference wavefunction used is determined by the underlying electronic structure method. In TDHF and TDDFT, |0〉 is a single Slater determinant whereas within the MCSCF theory, |0〉 is a linear combination of Slater determinants that is also an eigenstate of the configuration space. We write $\kappa \u02c6(t)$ in a general form as

where $\u2211\xb1\mu $ indicates a sum over both positive and negative indices of all operators (i.e., both *κ*_{1} and *κ*_{2}) and we have defined

and

Depending on the reference electronic structure method, several choices for the $T\u02c6\mu $ are possible. For example, the response of exact electronic states is conveniently expressed using a single exponential in which state-transfer operators project from the reference eigenstate, |0〉, to electronic state *n*, $R\u02c6n\u2020\u2261|n\u3009\u30080|$.^{7} On the other hand, in the case of the time-dependent Slater determinant in TDHF and TDDFT, single-electron excitation operators are used, i.e., $E\u02c6pq=a\u02c6p\u2020a\u02c6q$.^{26,30} MCSCF response theory employs a double-exponential in which state-transfer operators are responsible for rotations within the configuration space, and single-particle excitation operators are responsible for orbital rotations.^{7}

The response functions in terms of the Fourier components of derivatives of $\kappa \u02c6$ are found by differentiating the action, Eq. (3), with respect to the perturbation strengths, e.g.,

In the previous equation, we use a superscript with Greek indices to indicate differentiation such that $dfd\lambda \alpha \u2261f(\alpha )$, $d2fd\lambda \alpha d\lambda \beta \u2261f(\alpha \beta )$, and so on. Next, response equations are found at a given order by enforcing the stationarity of the action. For example, the linear response equations are determined by requiring

The linear response operator is written as **E**^{[2]} − *ω***S**^{[2]} using^{29}

Before writing the linear response equations, we first reorder **E**^{[2]} and **S**^{[2]} into the form

where $A\mu \nu =E\u2212\mu \nu [2]$, $B\mu \nu =E\u2212\mu \u2212\nu [2]$, $\Sigma \mu \nu =S\u2212\mu \nu [2]$, and $\Xi \mu \nu =S\u2212\mu \u2212\nu [2]$. The linear response equations can then be compactly written as^{29}

where

and

In addition, we define the combined column vector $|x,y)\u2261xy$ and the combined row vector (*x*, *y*| ≡ |*x*, *y*)^{†} for vectors *x* and *y*. The linear response equation obtained by evaluating Eqs. (15) and (16) within the TDHF or TDDFT context (single Slater determinant reference, single exponential parametrized by single excitations) differs superficially from the linear response equation encountered in the literature:^{30–32} the two formulations differ by a unitary transformation and therefore lead to identical properties and eigenvalues. The linear response function is then given as

The second-order equations can be expressed in a similar form as^{29}

where *X*^{(αβ)} and *Y*^{(αβ)} are defined in analogy to *X*^{(α)} and *Y*^{(α)}. The right-hand-side is given by

The response quantities **F**^{(α)} and **G** are response tensors of rank 2 and 3, respectively, defined as in Ref. 29 and repeated in the supplementary material. The **G** tensor is comprised of third derivatives of the action with respect to the rotation parameters

which is equivalent to a sum of triple commutators such as $\u30080|[[[H\u02c60,O\u02c6\mu ],O\u02c6\nu ],O\u02c6\zeta ]|0\u3009$. Elements of the **G** tensor play a central role in the spurious poles, as shown below.

The quadratic response function for the general case is then

Eq. (28) is a realization of the 2*n* + 1 rule that is obtained from applying the inverse linear response operator contained in $X(\alpha \beta ),Y(\alpha \beta )$ to the left. Although the two formulae are equivalent, we include both as certain residues are more easily formulated proceeding from one form or the other. In addition, the two forms imply different computational strategies. For example, computing the hyperpolarizability from Eq. (28) is advantageous computationally since no second-order parameters, *X*^{(αβ)} and *Y*^{(αβ)}, are necessary and therefore all of the first-order parameters can be computed simultaneously. On the other hand, Eq. (27) is more convenient for computing transition moments between excited states.

### B. Response of an exact state

Now we turn to evaluate several important residues of the linear and quadratic response functions using state-transfer operators between exact electronic eigenstates, i.e., states that satisfy

In this case, the state transfer operator $R\u02c6n\u2020$ is a linear combination of many-body excitations (up to *N _{e}*-body excitations where

*N*is the number of electrons) that transforms the ground-state wavefunction |0〉 into eigenstate |

_{e}*n*〉. Evaluating Eqs. (15)-(17) and (20), we find that

*A*= Ω

_{nm}_{n}

*δ*, Σ

_{nm}_{nm}=

*δ*,

_{nm}*B*= Ξ

_{nm}_{nm}= 0, and $Pn(\alpha )=\u2212vn0(\alpha )$, such that the linear response function becomes

We thus recognize that the poles of the linear response function correspond to excitation energies, and the associated residues,

provide transition moments. Furthermore, both the excitation energies and transition moments are found by solving the eigenvalue problem

subject to the normalization condition

with

To continue the analysis to second order, we first consider the elements of **G**. An important feature of *G*_{μνζ} is that it vanishes when *μ*, *ν*, and *ζ* all label state-transfer operators (both excitation and de-excitation operators), but is in general non-zero for many-body excitation and de-excitation operators. Combining this result with Eq. (27) we arrive at the usual sum-over-states expression,

with $v\u0304nm=vnm\u2212\delta nmv00$ and *ω*_{α} + *ω*_{β} + *ω*_{γ} = 0.

Two residues of the exact quadratic response function are especially important. First, the single residue

yields the two-photon absorption amplitude

The second important residue is the double residue

which yields expectation values of excited states and transition properties between excited states. The derivations of both of the quadratic response residues implicitly assume that Ω_{mn}≠Ω_{k} for any *k*.

### C. Spectral representation

Before discussing the excited state properties and residues obtained from approximate methods, it is instructive to introduce the spectral representation of the linear response operator

which is written in terms of the excitation energies and vectors from Eq. (33). In the previous equation and for the remainder of this paper, we label solutions of the response eigenvalue problem with uppercase *N*, *M*, *K* and reserve lowercase *n*, *m*, *k* for results from exact eigenstates. The spectral representation is most helpful in determining residues. For instance, it is easy to see from the spectral representation and the linear response equations that

and

With the aim of simplifying the proceeding section, we introduce the following notation:

We further introduce diagonal elements of response operators as

where **W** is a rank *k* operator and the *N _{k}* label eigenvectors of the response eigenvalue problem.

### D. Residues of approximate variational states

Following Refs. 7, 28, and 29, the results of Sec. II C may be used to *define* excited state properties for approximate electronic structure methods, especially those without clearly identifiable excited states such as TDHF and adiabatic TDDFT. For example, excitation energies and transition moments between the ground state and an excited state are defined according to Eqs. (33) and (35).

By inserting the spectral representation of the linear response operator into the definition of the quadratic response function in Eq. (27), the quadratic response function is expressed as

The previous equation shows that in the general case, the approximate quadratic response function contains terms with products of three poles, in contrast to the exact case in Eq. (36). Furthermore, while the residue at *ω*_{α}, *ω*_{β} → Ω_{n}, Ω_{m} vanishes in the exact case, it can be seen to be non-zero for the approximate response function. This generalizes Dalgaard’s observations^{26} from TDHF to a general variational response method.

The residues of the quadratic response function may be used to define state-to-state one-photon transition properties, excited-state expectation values, and the ground-to-excited state two-photon transition moments. Consider, for instance, the two-photon absorption cross section which we first rewrite using the permutation and time-reversal symmetry of the quadratic response function as

From the previous equation and Eq. (28) it follows that

However, this yields a sum-over-states expression

which exhibits a striking difference with the sum-over-states expression using exact electronic states: the approximate two-photon absorption amplitudes contain terms with *two* poles whereas the exact amplitudes have only terms with a *single* pole.

Starting from Eq. (27), we can conveniently write the state-to-state transition moments and excited state expectation values as

with

and

Inserting the spectral representation of the inverse linear response operator

we again see that the approximate transition moments contain an extra pole relative to the exact transition moments. The consequence of this extra pole is clear: state-to-state transition properties diverge whenever |Ω_{MN}| approaches any other excitation energy Ω_{K}. This is a key result of this paper.

For both the two-photon absorption and the state-to-state transition moments, the residues of the spurious poles are proportional to elements of **G**. As a consequence, the unphysical divergences will only manifest in methods for which **G** is non-zero. Since triple commutators of state-transfer operators always identically vanish, methods that employ only state-transfer operators—full configuration interaction (CI) or orbital unrelaxed CI—recover the correct pole structure. On the other-hand, methods that parametrize the time-dependent state using (at least in part) an incomplete set of many-body operators—TDHF, adiabatic TDDFT, and MCSCF response—will diverge erroneously.

### E. CC response theory

In this section, we analyze the behavior of coupled cluster response theory.^{29,33,34} For a thorough derivation and for definitions of all relevant terms, the reader is referred to Ref. 29. The general forms of the response equations (and thus the primary results) in this section are applicable to both canonical truncated coupled-cluster (CCSD, CCSDT,...) and iterative models such as CC2 and CC3, although the specific representations of the terms provided in Ref. 29 and in the supplementary material are valid only for truncated CC.

The CC Lagrangian is written as

where

is the CC wavefunction and

is the set of Lagrange multipliers. |R〉 is the reference (usually Hartree–Fock) determinant, $\u3008R|\tau \u0304\mu i\u2020=\u3008\mu i|$,

is the cluster operator, and *t*_{μi} is the cluster amplitude for the *μ*-th component of the *i*-fold excitation operator $\tau \u02c6\mu i$. The excitation operators, $\tau \u02c6\mu i$, and de-excitation operators, $\tau \u0304\mu i\u2020$, are chosen to satisfy $\u3008R|\tau \u0304\mu i\u2020\tau \u02c6\nu j|R\u3009=\delta \mu \nu \delta ij$ and, $\tau \u0304\mu i\u2020$ and $\tau \u02c6\mu i$ are not necessarily adjoints of each other.

The linear response equations are

from which both the linear response function and the quadratic response function are obtained, where **I** is the identity matrix, and *ξ*^{(α)}, *η*^{(α)} are rank 1 tensors, and **A** (the CC Jacobian) and **F** are rank 2 tensors defined as in Ref. 29 and repeated in the supplementary material. The CC quadratic response function has previously been shown to include terms with products of four poles and to have non-zero residues at *ω*_{α}, *ω*_{β} → Ω_{N}, Ω_{M}.^{29}

To examine the pole structure of the CC linear and quadratic response functions, one can diagonalize the non-Hermitian CC Jacobian

where the rows of **L** and the columns of **R** contain the left and right eigenvectors, respectively. This leads to a spectral representation for the CC inverse response operator of

where *R ^{N}* and

*L*are the

^{N}*N*-th column and row of

**R**and

**L**, respectively.

The non-Hermitian nature of the CC ansatz complicates the identification of transition moments from response functions. However, by defining

the residue of the linear response function may be expressed as

which matches the expectation from exact response theory. The response properties in the CC theory have an incorrect pole structure already for linear response: the linear response function diverges near a conical intersection involving the ground state (Ω_{N} → 0). This well-known aberration in the linear response function is also found in Brueckner coupled cluster,^{35} but has been shown to be removable by orbital optimization^{36} (which has other deficiencies, however^{37}). Interestingly, only $v0N(\alpha )$ exhibits this divergence. Although one may be thus tempted to prefer the left transition moment over the right transition moment, this is not physically justifiable as only observables are well-defined in response theory and individual transition moments are not observables.^{34} Furthermore, due caution is necessary whenever symmetries that are respected in exact theory, $v0N(\alpha )\u2212vN0(\alpha )\u2217=0$, are catastrophically broken such as is the case when $v0N(\alpha )\u2212vN0(\alpha )\u2217\u2192\u221e$. The relevant observable in the case of CC response theory is the symmetrized transition strength $S0N\alpha \beta =12(v0N(\beta )vN0(\alpha )+v0N(\alpha )\u2217vN0(\beta )\u2217)$ which diverges when either the left or the right transition moment diverges.

Just as one-photon transition moments were defined such that the CC linear response function would resemble the exact linear response function, the two-photon absorption amplitudes may be defined as

where *ω*_{β} = Ω_{K} − *ω*_{α} and $P(\alpha ),(\beta )$ symmetrizes with respect to (*α*) and (*β*). The response tensors **A**^{[α]}, **F**^{[α]}, **B**, and **G** are defined in Ref. 29 and repeated in the supplementary material. Similarly to the case of the linear response transition moments, the symmetrized transition strength

diverges when either $v0N(\alpha \beta )(\u2212\omega )$ or $vN0(\alpha \beta )(\omega )$ diverges.

Finally, state-to-state transition properties are expressed as^{29}

which has a similar structure as in the variational theories with one major difference: assuming Ω_{N}≠Ω_{M} and all Ω_{K} > 0, only one of $vNM(\alpha )$ or $vMN(\alpha )$ diverges since Ω_{MN} = − Ω_{NM}. The symmetrized transition strength, $SNM\alpha \beta =12(vNM(\beta )vMN(\alpha )+vNM(\alpha )\u2217vMN(\beta )\u2217)$, diverges when either transition moment diverges. This is reminiscent of the situation for linear transition moments (see above). Finally, the state-to-state properties derived from Brueckner CC response^{35} exhibit the same spurious poles.

. | Variational . | CC . | ||
---|---|---|---|---|

Residue . | Condition . | Order . | Condition . | Order . |

Linear response | ||||

$v0N(\alpha )$ | 0 | −Ω_{N} | 1 | |

Quadratic response | ||||

$v0N(\alpha \beta )(\omega )$ | ω, Ω_{N} − ω | 1 | ω, Ω_{N} − ω, − Ω_{N} | 2 |

$vNM(\alpha )$ | |Ω_{MN}| | 1 | Ω_{MN} | 1 |

. | Variational . | CC . | ||
---|---|---|---|---|

Residue . | Condition . | Order . | Condition . | Order . |

Linear response | ||||

$v0N(\alpha )$ | 0 | −Ω_{N} | 1 | |

Quadratic response | ||||

$v0N(\alpha \beta )(\omega )$ | ω, Ω_{N} − ω | 1 | ω, Ω_{N} − ω, − Ω_{N} | 2 |

$vNM(\alpha )$ | |Ω_{MN}| | 1 | Ω_{MN} | 1 |

## III. NUMERICAL EXAMPLES

In this section, we illustrate the consequences of the incorrect pole structure of the CC linear response function and the quadratic response function for several approximate response methods (Table I). For quadratic response, we examine exclusively the state-to-state transition properties for which the erroneous behavior is most clear and striking. We expect the spurious poles found in the two-photon absorption amplitudes to most strongly affect the resonant two-photon absorption and lead to amplitudes that are “overly resonant.” In particular, we focus on the transition dipole moments, but the same behavior is expected for any transition property such as higher-order electric moments, magnetic moments, or nonadiabatic couplings.

### A. PSB3

First, we examine the spurious poles of the CC linear response function near a ground state degeneracy by computing transition dipole moments between the reference state and the lowest energy response state of the penta-2,4-dieniminium cation (PSB3) using RI-CC2^{38,39} as implemented in Turbomole 7.0.^{40} PSB3, shown in Fig. 1, is an extensively studied computational model of the retinal protonated Schiff base chromophore (rPSB) of visual pigments.^{41,42} Here, we scan along the bond length alternation path reported by Olivucci and co-workers.^{43} The path is defined by interpolating between two structures denoted as TS_{CT} (for charge-transfer transition state) and TS_{DIR} (for diradical transition state).

A conical intersection between the ground state and the first excited state can be found with RI-CC2 by extrapolating beyond the TS_{DIR} structure.^{44} From Eqs. (63) and (64), we expect the right moment (*μ*_{01}) to diverge as Ω_{1} → 0 but not the left moment (*μ*_{10}). Fig. 1 shows the total energies and transition moments computed along the (extrapolated) bond length alternation coordinate using RI-CC2 with the def2-SVP^{45} basis set. Both the structure interpolation and extrapolation were performed using Linear Synchronous Transit^{46} (LST). Even though the energies all vary smoothly across the whole range, *both* transition moments become unreliable when the ground state becomes nearly degenerate: the right transition moment diverges as expected while the left moment is affected by the instability in the CC2 equations themselves as evidenced by the fact that the calculations failed to converge close to the degeneracy.

### B. LiH

Next, we consider the transition moment between two excited states of LiH as a function of the bond length. We compare the results of adiabatic TDDFT (PBE0 functional^{47}), complete active space self-consistent field (CASSCF) (4 electrons in the 6 core and valence orbitals), and RI-CC2, all using the def2-SVP basis set. The RI-CC2 calculations were performed using Turbomole 7.0 and the adiabatic TDDFT calculations were performed using a development version of Turbomole. The CASSCF calculations were performed using the Dalton program package.^{48} In all methods, the transition moments between the first and the fourth excited states (first and second excited states in *A*_{1} irreducible representation of *C*_{2v} point group symmetry) were computed. Fig. 2 compares the results to full CI calculations with the same basis set. Importantly, the excitation energies computed from all methods are fairly accurate for the entire range considered. This is especially apparent for the CASSCF energies which are virtually indistinguishable from the full CI energies. Furthermore, the transition dipole moments are essentially correct near the equilibrium bond length. Still, all approximate methods fail catastrophically when Ω_{10} = Ω_{41}. Note that even though *μ*_{41} for CC is well behaved, the large discrepancy between *μ*_{41} and *μ*_{14} should nonetheless be considered an unphysical result.

### C. 1,3-cyclohexadiene

We now turn to 1,3-cyclohexadiene (CHD), an organic molecule that undergoes conrotatory ring-opening to form hexatriene in accordance with the Woodward–Hoffmann rules.^{49} To define a path, we first optimized the ground state (S_{0}) geometry and the first excited state (S_{1}) geometry using PBE0 with the def2-SVP basis set. The optimized geometries are included in the supplementary material. The path is obtained by interpolating between the S_{0} and S_{1} geometries using LST. Compared to the ground state geometry, the carbon-carbon single bond in the relaxed excited state geometry is elongated (1.53 Å to 1.77 Å), and the dihedral angle defined by the three consecutive single bonds is flattened (14.4° to 4.6°). Fig. 3 shows the energies and transition moments between the first two excited states along this path using PBE0 and RI-CC2. Results from CASSCF response theory are not shown because the matching condition (Ω_{21} = Ω_{1}) is never satisfied, an artifact of the lack of dynamical correlation in CASSCF.^{50} We again notice that the energies are well behaved throughout the path and qualitatively similar with both methods. However, near the S_{1} geometry the transition dipole moments from both methods diverge. It is especially problematic that the divergences observed with both methods occur so close to the excited state relaxed geometry. We find the crossings to approximately occur at interpolation fractions of 0.82 and 1.02 for PBE0 and RI-CC2, respectively.

## IV. DIAGNOSIS

To understand the origin of the spurious poles in the state-to-state transition properties, we examine the conditions under which the spurious poles vanish, focusing on the case of a variational method to facilitate comparison with the full CI limit. The response tensor **G** is the third derivative of the action with respect to the response parameters [Eq. (26)],

where $P\mu \nu \zeta \alpha \beta \gamma $ symmetrizes with permutations of the pairs ($O\u02c6\mu ,\omega \alpha $), ($O\u02c6\nu ,\omega \beta $), and ($O\u02c6\zeta ,\omega \gamma $).^{29} For static response functions, this reduces to the straightforward third derivative of the energy as a function of the rotation parameters in ** κ**. The extra poles in the quadratic response function (i.e., the terms containing products of three poles) vanish if the diagonal elements of

**G**take the form

which can be verified by substituting the previous equation into Eqs. (49) and (54). If one assumes that the diagonal excitation operators satisfy $[H\u02c60,\kappa \u02c6N]=\Omega N\kappa \u02c6N$, a condition that is equivalent to assuming the full CI limit, then

and thus Eq. (71) is satisfied: there are no spurious poles in either the quadratic response function or its residues. On the other hand, in the full CI limit, the excitation operator, $\kappa \u02c6N$, is equivalent to the state-transfer operator, |*n*〉〈0| and all of the commutators in the previous equation vanish. Short of the full CI limit—i.e., employing a truncated set of many-body excitation operators for ${T\u02c6\mu}$—the commutation condition is not satisfied, and the remainder causes the divergences seen in Sec. III.

An equivalent interpretation is possible in terms of the non-linearity of the reference. Since the linear response operator is the second derivative of the action with respect to rotation parameters

**G** can be seen as the change in the linear response operator stemming from a change in the reference wavefunction. In full CI, the reference is the vacuum state and thus the linear response operator is independent of the reference. For all other approximate methods, however, the response operator depends parametrically on the reference and thus **G** is non-zero.

The spurious poles in the quadratic response function of adiabatic TDDFT are closely related to the peak-shifting phenomenon of real-time (RT)-TDDFT, as both arise from the non-linearity of the reference.^{51,52} However, peak-shifting is only observed in the non-perturbative regime, whereas the spurious poles occur already in second-order response theory. An intriguing possibility may be to use Eq. (71) to inform exact conditions on frequency-dependent exchange-correlation kernels, *f ^{xc}*, since

**G**contains contributions from both

*f*and the exchange–correlation hyperkernel

^{xc}where the last line represents the adiabatic approximation to the exchange–correlation kernel.

## V. REMEDIES

We can envision several routes towards circumventing the unphysical behavior of the state-to-state transition properties. The first possible route is to set **G** = 0. We refer to this as the “unrelaxed approximation” since it ignores the relaxation of the wavefunction parameters in response to the perturbation. This route has the advantage of reproducing the correct pole structure while eliminating the need to solve for the second-order response. It is, however, not fully consistent in the sense that the resulting response properties in general do not agree with the derivatives of the action. For example, the excited state difference dipole moment, $\mu \u0304NN=\mu NN\u2212\mu 00$, is computed either directly from the residue of the quadratic response function or from the derivative of the excitation energy with respect to the external field, $d\Omega Nd\epsilon $. Within the unrelaxed approximation, the residue of the quadratic response function will not match $d\Omega Nd\epsilon $ computed from fully relaxed finite difference calculations.^{53}

The second option is to solve the second-order response equation at zero frequency, i.e., set Ω_{MN} = 0 in Eq. (51). In the context of TDHF and adiabatic TDDFT, this is equivalent to the “pseudo-wavefunction approximation” proposed by Subotnik and co-workers.^{54,55} Similar to the unrelaxed approximation, the pseudo-wavefunction approximation avoids the spurious pole at |Ω_{MN}| = Ω_{K}. Meanwhile, in contrast to the unrelaxed approximation, the pseudo-wavefunction approximation is consistent for certain excited state properties, for which Ω_{NN} = 0. However, rather than *removing* the spurious pole at |Ω_{MN}| = Ω_{K}, it is simply shifted to Ω_{K} = 0; thus any excited state property *and* any state-to-state property become unphysical near a degenerate ground state. For transition properties between excited states or excited state dynamical polarizabilities, however, the pseudo-wavefunction approximation is also inconsistent with action derivatives.

Another possibility is to obviate the quadratic response function entirely and obtain state-to-state properties from the linear response function *from a different reference*. For example, *μ _{NM}* could be just as well computed from the linear response function using |

*N*〉 as the stationary reference. This is equivalent to a redefinition of “state-to-state” property and “ground-to-excited state” property or “reference-to-excited state” property. Such an approach is similar in spirit to spin-flip methods in which the ground state and excited states are expressed as “response states” that result from spin-flipping excitations from a high-spin reference.

^{56–58}To illustrate this possibility, we again compute the transition dipole moment between the first and fourth excited states of LiH with CASSCF, this time by optimizing the orbitals for the first excited state and then computing excitations from the first excited state with linear response theory. Fig. 4 compares the results to those discussed previously using quadratic response from the ground state reference. Using the linear response function from the excited state reference, the transition moment no longer diverges. However, obtaining

*all possible*state-to-state transition moments in this manner would require separate non-linear optimizations for each excited state which may be difficult and costly. Furthermore, such a procedure would introduce ambiguities since transition moments and energies computed from different references will not agree, in general. As an extension of this idea, one could use an ensemble reference and exploit the fact that all possible state-to-state properties are contained in the ensemble linear response function

where *w _{k}* is the weight associated with the

*k*-th component of the ensemble. Such a strategy resembles state-averaged CASSCF and exhibits a similar ambiguity arising from the weight dependence.

## VI. CONCLUSIONS

The quadratic response function computed from TDHF and adiabatic TDDFT exhibits spurious poles that result in “overly resonant” two-photon absorption amplitudes and unphysical divergences in the transition properties between two excited states. We have shown that these spurious poles are deficiencies of approximate response theory that plague not only TDHF and adiabatic TDDFT but also response theories based on MCSCF and CC references. Furthermore, we showed that the CC linear response function suffers from similar deficiencies. A major result of the present paper is that despite the successes of response theories in terms of providing accurate excitation energies and linear and non-linear polarization properties, none of the approximate many-electron response theories discussed here are free of unphysical defects that render them inappropriate for the description of non-adiabatic dynamics.

The consequences of the incorrect pole structure are profound. Considering the transition dipole vector, *μ*_{nm}, and applying the Cauchy–Schwarz inequality

it follows that the transition dipole moment for square-integrable electronic wavefunctions must be bounded. This bound is violated in the state-to-state transition moments obtained from all approximate methods discussed here. Furthermore, our numerical results demonstrate that state-to-state transition moments from approximate quadratic response functions can fail dramatically *even when the energies and methods are otherwise well-defined and well-behaved*. The transition moments between the ground state and an excited state obtained from the CC linear response function also violate this bound when the optical gap approaches zero. This problem, however, is less concerning since single-reference CC is unreliable^{59,60} near a degenerate ground state. Nonetheless, this suggests one must be cautious when computing linear transition properties using CC with a nearly degenerate ground state since even well converged and smoothly varying energies do not guarantee physically meaningful transition properties. This is in addition to the incorrect dimensionality of conical intersections derived from CC.^{61}

The assumption that the pole structure of approximate response functions is qualitatively correct underlies much of modern response theory.^{7} The present work shows that this assumption is not justified for the most common electronic structure methods beyond linear response.

The implications of the spurious poles for NAMD simulations are also considerable. In NAMD simulations, the various electronic states are typically coupled through one-electron properties (e.g., derivative couplings or electrostatic potentials) obtained through linear and quadratic response. Unphysical divergences in the coupling matrix elements will result in incorrect electronic wavefunctions and indirectly in fictitious electronic transitions. The occurrence of a single such fictitious transition would cast doubt on the validity of the whole trajectory. Furthermore, we expect such accidental crossings to be ubiquitous during NAMD simulations. Consider, for example, a trajectory in which the initial configuration satisfies 0 < Ω_{1} < Ω_{2} < 2Ω_{1} and, in addition, the first two excited states are far from degenerate. If the trajectory proceeds through a configuration that exhibits a conical intersection between the ground state and the first excited state (i.e., if Ω_{1} → 0), then there must exist an intermediate configuration where Ω_{2} − Ω_{1} = Ω_{1}. These conditions on the initial configuration are common; they are consistent with the observation underpinning Kasha’s rule^{62} that the S_{1}-S_{0} energy gap tends to be larger than the S_{n}-S_{n−1} gap (*n* ≥ 2) and were satisfied in the equilibrium configurations of all of the organic chromophores we tested. The PSB3 and CHD examples in Sec. III demonstrated that the unphysical couplings are likely to occur in regions of interest for NAMD simulations. Since the matching condition does not correspond to any true physical degeneracy, predicting or avoiding their locations in advance is not in general feasible.

The present conclusions apply only to excited states and couplings derived from response theory. Excited states and couplings derived from, for example, state-averaged CASSCF or its multi-state perturbative extensions do not diverge because the electronic states are explicitly parametrized.

The spurious poles are a consequence of the (i) nonlinear and (ii) instantaneous dependence of the effective Hamiltonian on variational parameters in approximate time-dependent electronic structure theories. These parts of the effective Hamiltonian become unphysically resonant when the difference of two excitation energies equals another excitation energy of the system, causing spurious poles in quadratic and higher-order response theory. In exact response theory, the action is a quadratic functional of the variational parameters, and the Hamiltonian is time-dependent only through external fields; thus there are only physical poles.

Retardation due to “memory effects” also underlies the lack of double and higher excitations and other spurious behavior in adiabatic TDDFT.^{63–67} Physically, this retardation may be related to dissipation, i.e., damping, due to the electron interaction.^{68} Our results suggest that retardation and the resulting frequency dependence of response kernels are crucial features of time-dependent effective many-electron theories and need to be present in a qualitatively correct nonlinear response theory.

## SUPPLEMENTARY MATERIAL

See the supplementary material for explicit expressions of all response quantities mentioned as well as ground and excited state geometries of 1,3-cyclohexadiene.

## Acknowledgments

This material is based on work supported by the U.S. Department of Energy under Award No. DE-SC0008694.

Principal investigator Filipp Furche has an equity interest in Turbomole GmbH. The terms of this arrangement have been reviewed and approved by the University of California, Irvine, in accordance with its conflict of interest policies.