Unphysical divergences in response theory

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,¹ the expectation value of an observable property O associated with the self-adjoint operator $Oˆ(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ˆ$ 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.² 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⁵ 

Time-dependent properties are defined by response of the action to a time-dependent perturbation $λ(t)Oˆ(t)$ added to the Hamiltonian. By the time-dependent Hellmann-Feynman theorem⁶ 

which recovers the expectation value of $Oˆ(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.⁷ 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,¹⁷ 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²² while couplings between all pairs of excited states (state-to-state) are obtained from the quadratic response function.²³ 

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₀, Ω_N = E_N − E₀ 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.²⁶ 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.

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ˆ0$ is the time-independent field-free contribution; time-dependent perturbations are of the form

with λ_α, ω_α, and $vˆ(α)$ the strength, frequency, and potential associated with the α component of the perturbation. In keeping with the Hermiticity of $vˆλ$⁠, the λ_α and ω_α are real and $(vˆ(α))†=vˆ(−α)$⁠. While the form of the field-free effective Hamiltonian will depend on the underlying electronic structure method, Eq. (7) is valid throughout this work.

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.³⁰ 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 $κˆ(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 $κˆ(t)$ in a general form as

where $∑±μ$ indicates a sum over both positive and negative indices of all operators (i.e., both κ₁ and κ₂) and we have defined

and

Depending on the reference electronic structure method, several choices for the $Tˆμ$ 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ˆn†≡|n〉〈0|$⁠.⁷ 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ˆpq=aˆp†aˆq$⁠.^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.⁷ 

The response functions in terms of the Fourier components of derivatives of $κˆ$ 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λα≡f(α)$⁠, $d2fdλαdλβ≡f(αβ)$⁠, 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²⁹ 

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

where $Aμν=E−μν[2]$⁠, $Bμν=E−μ−ν[2]$⁠, $Σμν=S−μν[2]$⁠, and $Ξμν=S−μ−ν[2]$⁠. The linear response equations can then be compactly written as²⁹ 

where

and

In addition, we define the combined column vector $|x,y)≡xy$ 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²⁹ 

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 $〈0|[[[Hˆ0,Oˆμ],Oˆν],Oˆζ]|0〉$⁠. 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 2n + 1 rule that is obtained from applying the inverse linear response operator contained in $X(αβ),Y(αβ)$ 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.

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ˆn†$ is a linear combination of many-body excitations (up to N_e-body excitations where N_e is the number of electrons) that transforms the ground-state wavefunction |0〉 into eigenstate |n〉. Evaluating Eqs. (15)-(17) and (20), we find that A_nm = Ω_nδ_nm, Σ_nm = δ_nm, B_nm = Ξ_nm = 0, and $Pn(α)=−vn0(α)$⁠, 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̄nm=vnm−δ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.

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.

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²⁶ 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.

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, $〈R|τ̄μi†=〈μi|$⁠,

is the cluster operator, and t_μi is the cluster amplitude for the μ-th component of the i-fold excitation operator $τˆμi$⁠. The excitation operators, $τˆμi$⁠, and de-excitation operators, $τ̄μi†$⁠, are chosen to satisfy $〈R|τ̄μi†τˆνj|R〉=δμνδij$ and, $τ̄μi†$ and $τˆμ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.²⁹ 

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^N are the 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,³⁵ but has been shown to be removable by orbital optimization³⁶ (which has other deficiencies, however³⁷). Interestingly, only $v0N(α)$ 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.³⁴ Furthermore, due caution is necessary whenever symmetries that are respected in exact theory, $v0N(α)−vN0(α)∗=0$⁠, are catastrophically broken such as is the case when $v0N(α)−vN0(α)∗→∞$⁠. The relevant observable in the case of CC response theory is the symmetrized transition strength $S0Nαβ=12(v0N(β)vN0(α)+v0N(α)∗vN0(β)∗)$ 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