Oscillator strengths and excited-state couplings for double excitations in time-dependent density functional theory

The challenge of describing excited states of double-excitation character has grown from being a relatively minor nuisance in certain regions of linear response spectra to being a real hindrance in our ability to simulate photo-excited dynamics in molecules. These states are loosely defined in the sense of having a significant proportion of doubly excited determinants with respect to a reference ground-state Slater determinant, where a doubly excited determinant promotes two occupied orbitals of the reference determinant to two virtual ones.¹ In ground-state spectra, these states have typically smaller transition strengths (“darker”) than predominantly singly excited states, but they can siphon off oscillator strength of neighboring single-excitations whose absorption peaks will be overestimated in simulations that do not account for double-excitations. Moreover, these states are readily accessible when a molecule is excited, offering photochemical/physical pathways that play a role in a range of applications, from photocatalytic design and photo-induced charge transport to laser-driven or polariton-enabled control. While the true correlated excited-state wavefunctions contain contributions from double excitations, these are either absent in approximate electronic structure methods or poorly approximated, resulting in states missing from the spectrum or with energies and transition strengths with typically larger errors than for states with predominantly single-excitation character especially when there is a large contribution from double-excitations,^2,3 wreaking havoc on predictions of excited-state dynamics.

The problem is particularly severe for single-reference methods in their usual modus operandi, such as time-dependent density functional theory (TDDFT) within the adiabatic approximation,^1,4–6 equation-of-motion coupled-cluster (EOM-CC) with singles and doubles only,^3,7 or the Bethe–Salpeter equation with the GW approximation (GW-BSE) with static screening.^8–10 While double-excitation character can be adequately included in EOM-CC if triple excitations are included in the cluster expansion and in GW-BSE with full frequency-dependence or equivalently in an expanded excitation space, these procedures increase the computational cost of these methods, both of which are already more expensive than TDDFT. Multiconfigurational methods may be better platforms for these states, but aside from their higher computational cost, other issues enter, such as gauge-sensitivity and the sizes of an adequate active space that can change during a dynamical process. For determination of oscillator strengths for these states (or any excited state), the difficulty increases significantly since increased accuracy requires cranking up the number of electronic configurations more so than for corresponding improvement in energies, which results in a large increase in the computational cost.^11,12

Fixing adiabatic TDDFT thus becomes an attractive option. Indeed, double-excitations provide a prime example where TDDFT in theory and in practice diverges: in theory, their excitation energies and oscillator strengths are captured exactly, but the approximate functionals used in practice cannot capture them. It is now well-known that the culprit is the adiabatic nature of the approximation:^1,4–6 the exact exchange–correlation (xc) kernel has a strong frequency-dependence (a pole) that effectively folds a double-excitation of the non-interacting Kohn–Sham (KS) system into TDDFT response equations, but adiabatic approximations use a frequency-independent kernel. The frequency-dependence implicit in orbital-dependence of a hybrid functional does not have the correct structure.¹³ Using first-principles arguments, the authors of Ref. 6 proposed a frequency-dependent kernel, sometimes called “dressed TDDFT,” that has been shown to successfully capture the energies of double-excitations in a range of molecules.^14–18 However, it is unable to yield information on the oscillator strengths because it is designed to operate within the framework of a Tamm–Dancoff-like approximation, which does not preserve oscillator strengths.

In this Communication, we derive a kernel designed to work within full TDDFT linear response, which yields both excitation energies and oscillator strengths of these states, preserving the oscillator strength sum-rule. We demonstrate its accuracy on a two-electron model system, the lowest ¹D excitations in the Be atom, and the LiH molecule over a range of interatomic separations. The kernel provides an approximation of how the transition density of a KS single excitation gets distributed into mixed single- and double-excitations of the true system. This enters into a formula to compute excited-to-excited state transition densities that has been recently derived¹⁹ to tame unphysical divergences in the adiabatic quadratic response function, giving an improvement in the transition dipole between two states in the LiH molecule as it begins to dissociate, with accuracy similar to that of the more ad hoc pseudo-wavefunction approximation.^20–24 

A key point is that, whether with a frequency-dependent kernel or not, accurate oscillator strengths fulfilling the oscillator strength sum-rule can only be obtained when the full TDDFT matrix [Eq. (2)] is used, i.e., including both backward and forward transitions.^28,29 Backward transitions are neglected in the Tamm–Dancoff approximation where eigenvalues of the corresponding matrix are directly the frequencies ω_I (rather than $ω I 2$⁠). Although this approximation has been argued to sometimes give more accurate excitation energies with some approximate functionals,^29–31 the oscillator strengths are generally worse and the sum-rule is violated.

Equation (3), derived in Ref. 25, does not appear to be well known and, to our knowledge, used only once before³² in a detailed study of the linear response of the Hubbard dimer. An analogous redistribution of weights appears in many-body perturbation theory where dynamical kernels have also been explored.^33–35 In general, there has been much more attention paid to excitation energies than to oscillator strengths and the latter have only been studied with the adiabatic approximation^12,36–39 with eigenvectors normalized to 1. Most interesting for present purposes, Eq. (3) tells us how the oscillator strength of a KS single excitation that lies near a KS double excitation shatters into oscillator strengths of the mixed single and double excitations, when an appropriate frequency-dependent kernel is used.

We further note that since the matrix elements in Eq. (4) [Eq. (6)] involve a one-body operator, then in the limit that the true ground-state is dominated by a single Slater determinant (weak interaction limit), only the underlying single-excitation components of Ψ_I contribute to |G_I|². Since $∑ I G I 2 =1$⁠, this means that 1 − |G_I|² contains the double-excitation (and any higher-excitation) contribution, in this limit.

In Eq. (11), we also replaced E₀ in Eq. (9) by H₀₀ to better balance errors in the truncated matrix elements (as was done for the dressed SPA of Ref. 6).

Figure 1 shows the performance on an exactly solvable model system: two electrons in a one-dimensional harmonic plus linear potential $v ext ( x ) = 1 2 x 2 +γ|x|$⁠, where γ is a parameter in the range [0,1]; varying γ tunes the degree of mixing of the single and double-excitation character in the interacting excited states. The electrons interact via a soft-Coulomb interaction: $1 ( x 1 − x 2 ) 2 + 1$⁠. We observe in Fig. 1(a) that for small enough γ, the exact KS system has a double excitation lying very close to the second KS single excitation, with which it mixes strongly when the interaction is turned on as in the physical system, producing two states. As expected, the adiabatic approximation, here chosen to be exact exchange (AEXX), is blind to the double excitation and has only one excitation in this region of the spectrum, while all DSMAs correctly yield two excitations. (Note that the exact ground-state KS potential is used to find the bare KS orbitals and energies, while AEXX is used for $f H X C A$⁠.) Out of the DSMA’s outlined above, DSMA₀ performs the best, giving excitation frequencies very close to the exact. The excitation energies provided by each variant of DSMA are close to those given by the corresponding version of dressed SPA (DSPA) (see the supplementary material). However, the real and significant improvement of DSMA over DSPA lies in its capacity to accurately determine oscillator strengths. Figure 1(b) shows the fractions of oscillator strength shared by the two states. Again, we see that DSMA₀ is closest to the exact, but all give a significant improvement over the adiabatic result of 1 and 0 for any γ and also over the DSPA fractions shown in the supplementary material. All three flavors of DSMA respect the oscillator strength sum rule, Eq. (6) as shown in Fig. 1(b), in contrast to the DSPA (see the supplementary material for the analogous plot).

We next turn our attention to the beryllium atom. Here, the lowest two singlet ¹D states have a mixed single and double excitation character. The authors of Ref. 3 reported a roughly 30% single-excitation character to the predominantly doubly excited 1¹D state [1¹S(2s²) → 1¹D(2p²)]. As a result, adiabatic TDDFT fails to describe the states accurately. Figure 2 shows a plot of the two excitation frequencies of the 1¹D multiplet; the extreme right shows the exact reference from Ref. 5, with the lower one of predominantly 2p² nature, while the upper predominantly 2s3d (Ref. 3). Adiabatic TDDFT with the Perdew–Burke–Ernzerhof (PBE) functional gives only one state in this frequency region, which is closer in energy to the upper exact state, which would take up all the oscillator strength $( G q 2 = 1 )$ within the SMA. The three flavors of DSMA give two states, but overestimate their energy splitting. Whether and how the splitting improves with the inclusion of more diffuse basis functions is left to future study. Still, DSMA₀ and especially DSMA_S yield oscillator strengths (author request) which are close to that of Ref. 3 for the lower excitation, within the assumption that Be is weakly enough correlated that $G I 2$ meaningfully approximates the single-excitation component of the state I. On the other hand, the three flavors of DSPA also give two excitation frequencies but yield wrong oscillator strengths and violate the oscillator sum rule.

We now use DSMA to compute the proper weight G_c associated with the first and second terms in Eq. (15) for the transition dipole moment between the first and fourth singlet excited states in LiH, $μ 14 = ⟨ Ψ 1 | x ̂ | Ψ 4 ⟩ =∫x n c a App ( r ) dr$⁠, where x is along the internuclear axis. In Fig. 3, all calculations are done in the aug-cc-pVDZ basis set using PySCF,⁴⁰ while the pseudo-wavefunction calculation is performed with the Turbomole package.^41,42 As demonstrated in the supplementary material, both the TDDFT and reference Full Configuration Interaction (FCI) results show significant basis-set dependence; the transition dipole moments more so than the excitation energies, but convergence appears to be reached with the aug-cc-pVDZ basis. The top panel shows the intersection of the energy of the KS double-excitation $2 ν 1 APBE0$ with the fourth KS excitation frequency, justifying the need for the DSMA approach. The black lines demonstrate the “resonance condition” for FCI, Ω₄ = 2Ω₁, while near the intersection of the corresponding curves for APBE0 (blue lines), $μ 14 APBE0$ diverges (see Refs. 19 and 24 and the supplementary material). The middle panel shows the fourth and fifth singlet excitation energies of the exact and the one APBE0 excitation in this frequency-range, consistent with its blindness to the KS double-excitation in the top panel. Applying our DSMA to this fourth state, we retrieve two states, $Ω 4 ( 1 )$ and $Ω 4 ( 2 )$⁠, but with some error compared to the exact; both DSMA_S and DSMA_a overshoot the splitting between the two states, with the upper level of the latter being off the scale of the graph. The lowest panel compares μ₁₄ computed from Eq. (15) $( Ω c = Ω 4 ( 1 ) )$ using G_c obtained from DSMA₀ and DSMA_S, against the reference FCI and the pseudo-wavefunction approximation.^20–24,43 The latter is often used to tame the divergence of the raw adiabatic transition dipole, but is somewhat ad hoc in that its underlying kernel structure is not yet known. In the indicated region of strong mixing between the single and double excitation $R∼ [ 2.6 − 3.4 A ̊ ]$⁠, although they overestimate the dipole, DSMA_S and especially DSMA₀ both have a better trend as a function of R and are closer to the FCI reference than the pseudo-wavefunction approach. The dotted curves indicate the importance of the double-excitation contribution.

Note that the analysis here applies not just to transition dipole moments between excited states but also to any one-body multiplicative coupling matrix element since the TDDFT formalism is based on the one-body density. In particular, it also applies to non-adiabatic couplings between excited states that arise due to coupled electron–nuclear motion, $d I J , ν = ⟨ Ψ I | ∇ ν Ψ J ⟩ = ⟨ Ψ I | ∇ ν V e n | Ψ J ⟩ E J − E I$ (inner product is over the electronic Hilbert space only), where ∇_ν is the gradient with respect to the R_νth nuclear coordinate and V_en is the electron–nuclear Coulomb interaction.

By providing oscillator strengths and excitation energies of states of double-excitation character, the new kernel presented here improves the reliability of TDDFT for linear response spectra, as well as for quadratic response properties where it provides an otherwise missing contribution to the excited-to-excited state transition densities. Future work involves extensive tests on a range of systems including to cases when more than one single-excitation couples strongly with a double-excitation. In addition, ramifications of Eq. (3) outside the problem of double excitations will be explored. For example, calculating oscillator strengths using an orbital-dependent functional (exact exchange, hybrids, meta-GGAs) within pure DFT implies a frequency-dependence, and an analogous formula for the generalized KS framework may be needed.