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 Maitra (2022).

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 neighbouring 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, 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 larger errors than for states with predominantly single-excitation character Loos et al. (2019), 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 Jamorski, Casida, and Salahub (1996); Tozer and Handy (2000); Maitra et al. (2004); Maitra (2022), equation-of-motion coupled-cluster (EOM-CC) with singles and doubles only Hirata, Nooijen, and Bartlett (2000); Loos et al. (2019), or the Bethe-Salpeter equation with the GW approximation (GW-BSE) with static screening Romaniello et al. (2009); Sangalli et al. (2011); Bintrim and Berkelbach (2022). 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 which 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 costDamour et al. (2023); Sarkar et al. (2021).

Fixing adiabatic TDDFT thus becomes an attractive option. Indeed, double-excitations provide a prime example where TDDFT in theory and in practise diverge: in theory, their excitation energies and oscillator strengths are captured exactly, but the approximate functionals used in practise cannot capture them. It is now well-known that the culprit is the adiabatic nature of the approximation Jamorski, Casida, and Salahub (1996); Tozer and Handy (2000); Maitra et al. (2004); Maitra (2022): 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 the 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, Ref. Maitra et al. (2004) 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 Cave et al. (2004); Elliott et al. (2011); Huix-Rotllant et al. (2011); Mazur and Włodarczyk (2009); Mazur et al. (2011). 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 ${}^{1}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 was recently derived Dar, Roy, and Maitra (2023) 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 pseudowavefunction approximation Li and Liu (2014); Zhang and Herbert (2015); Ou, Alguire, and Subotnik (2015); Ou et al. (2015); Parker, Roy, and Furche (2016).

In TDDFT, excitation spectra are extracted from the linear response of the density, usually cast in the form of a matrix equation in the space of single KS excitations, $q=i\to a$, with $i(a)$ labelling occupied (unoccupied) orbitals Casida (1995); Petersilka, Gossmann, and Gross (1996); Grabo, Petersilka, and Gross (2000); Gross and Maitra (2012):

where

Here $f_{{\scriptscriptstyle\rm HXC}\,qq^{\prime}}=\int d{\bf r}d{\bf r}^{\prime}\Phi_{q}({\bf r})f_{\scriptscriptstyle\rm HXC}({\bf r},{\bf r}^{\prime},\omega)\Phi_{q^{\prime}}({\bf r}^{\prime})$ is the matrix element of the Hartree-exchange-correlation kernel, $f_{\scriptscriptstyle\rm HXC}({\bf r},{\bf r}^{\prime},\omega)=\frac{1}{|{\bf r}-{\bf r}^{\prime}|}+f_{\scriptscriptstyle\rm XC}({\bf r},{\bf r}^{\prime},\omega)$ with $\Phi_{q}=\phi_{i}^{*}\phi_{a}$, being the product of occupied and unoccupied orbitals. The KS frequencies $\nu_{q}=\epsilon_{a}-\epsilon_{i}$ are corrected to the true ones $\omega_{I}$ through solving Eq. (1). While the square-root of the (pseudo-)eigenvalues, $\omega_{I}$, give the excitation frequencies of the physical system, the oscillator strength of the transition can be obtained from ${\bf G}_{I}$ when normalized according to Casida (1995)

Oscillator strengths $f_{I}$ (or any one-body transition properties) can be obtained from these eigenvectors, through

where on the right ${\bf r}=(x,y,z)$ with $x$ being a vector in single-excitation space with matrix elements $x_{q}=\langle\phi_{i}|x|\phi_{a}\rangle$ and $S$ is a diagonal matrix $S_{qq}=\frac{1}{\nu_{q}}=\frac{1}{\epsilon_{a}-\epsilon_{i}}$. The oscillator strength sum-rule, $\sum_{I}f_{I}=N$ is satisfied by both the true and KS systems and Eq. (4) gives the contributions of the KS oscillator strengths to a given true $f_{I}$ through the normalized eigenvector ${\bf G}_{I}$. In the adiabatic approximation $\Omega(\omega)=\Omega$ has no frequency-dependence so the (pseudo-)eigenvectors are just normalized to $1$, but a frequency-dependent kernel causes a redistribution of oscillator strengths through Eq. (3).

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 Gross and Maitra (2012); Hirata and Head-Gordon (1999). Backward transitions are neglected in the Tamm-Dancoff approximation where eigenvalues of the corresponding matrix are directly the frequencies $\omega_{I}$ (rather than $\omega_{I}^{2}$). Although this approximation has been argued to sometimes give more accurate excitation energies with some approximate functionals Liang et al. (2022); Casida et al. (2000); Hirata and Head-Gordon (1999), the oscillator strengths are generally worse, and the sum-rule violated.

Eq. (3), derived in Ref. Casida (1995), does not appear to be well known and, to our knowledge, used only once before Carrascal et al. (2018) in a detailed study of the linear response of the Hubbard dimer. 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 Sarkar et al. (2021); Chrayteh et al. (2021); Robinson (2018); Caricato et al. (2011); Laurent and Jacquemin (2013) 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.

Consider the limit that a KS excitation is isolated enough in energy from neighbouring excitations, such that the diagonal corrections in the TDDFT matrix dominates, and the full TDDFT matrix Eq. (2) locally reduces to the $1\times 1$ “small-matrix approximation (SMA)",

Setting this equal to $\omega^{2}$, we observe that frequency-dependence of the kernel $f_{\scriptscriptstyle\rm XC}({\bf r},{\bf r}^{\prime},\omega)$ creates more than one TDDFT excitation from this single KS excitation. The eigenvectors $G_{I}$ reduce to one number for each eigenvalue, and indicate the fraction of the KS oscillator strength that goes into the corresponding true transition, with the second equality in Eq. (4) reducing to

In this case, the oscillator strength associated with this isolated KS transition will be preserved over the true excitations in this subspace provided that the approximation for $f_{\scriptscriptstyle\rm HXC}(\omega)$ that will go into $\Omega(\omega)$ gives eigenvectors with the property $\sum_{I}G_{I}^{2}=1$:

where $I$ labels the excitations generated from the single KS excitation by the frequency-dependent kernel. This is an important condition for the approximation for $f_{\scriptscriptstyle\rm HXC}(\omega)$. We note that it is satisfied for any $\Omega(\omega)$ of the form

where $A,B,D$ are constants.

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 $\Psi_{I}$ contribute to $|G_{I}|^{2}$. Since $\sum_{I}G_{I}^{2}=1$, this means that $1-|G_{I}|^{2}$ contains the double-excitation (and any higher-excitation) contribution, in this limit.

While the “dressed TDDFT" kernel derived in Ref. Maitra et al. (2004) has the correct form of frequency-dependence to yield energies of double-excitations, it is based on the Tamm-Dancoff approximation, which reduces to the single-pole approximation (SPA) in the isolated KS excitation case, and cannot be used in a consistent way to obtain oscillator strengths. Instead, here we follow a similar approach to Ref. Maitra et al. (2004) but within the SMA. We begin by writing the many-body time-independent Schrödinger equation such that square of the frequency is the eigenvalue: $(H-E_{0})^{2}\Psi=\omega^{2}\Psi$, evaluating this eigenvalue problem in the KS basis, and working in the limit where a KS single-excitation $q$ and double-excitation $d$ are well-separated in energy from all the other excitations of the system. We define $\nu_{d}=\nu_{s1}+\nu_{s2}$ as the sum of the KS single-excitations $\nu_{s1}$ and $\nu_{s2}$ such that $\nu_{d}$ is close to $\nu_{q}$. Diagonalization then yields

where $H_{ij}$ are the matrix elements of the interacting Hamiltonian in the basis $\{q,d\}$ which spans the truncated Hilbert space containing the single and the double excitations, and $E_{0}$ is the ground state energy. However Eq. (9) is nothing but exact diagonalization in the KS basis; to take advantage of the correlation contained in ground-state TDDFT approximations and to extract a frequency-dependent kernel for an improved TDDFT approximation, we replace $H_{qq}-E_{0}$ with its adiabatic SMA value, and identify the remainder as defining a frequency-dependent kernel:

In Eq. (11) we also replaced $E_{0}$ in Eq. (9) by $H_{00}$ to better balance errors in the truncated matrix elements (as was done for the dressed SPA of Ref. Maitra et al. (2004)).

Eq. (11) is our central equation, and it becomes exact in the limit that the coupling between the single and double excitation induced by the electron-interaction is much stronger than the coupling with other single excitations. While the adiabatic approximation merely shifts the KS frequency, Eq. (11) creates a new one, yielding two positive frequencies $\omega_{1},\omega_{2}$ which approximate the mixed single- and double- excitations of the true system, similar to the dressed SPA kernel of Ref. Maitra et al. (2004). However, unlike that kernel, Eq.(11) preserves the oscillator strength:

It has the form of Eq. 8 and using Eq. (5) with Eq. (11), Eq. (3) explicitly finds that the eigenvectors $G_{1(2)}$ (one number for each eigenvalue in this $1\times 1$ case) are such that $G_{1}^{2}+G_{2}^{2}=1$, thus satisfying the sum-rule (Eq. 7),

The resulting $G_{1,2}$ from Eq. (3) gives the fraction of the KS oscillator strength that goes into the respective transition $I=1,2$ (see Eq. (6)), with $1-G_{I}^{2}$ giving the double (and higher) excitation contribution to that state in the limit that the ground-state is dominated by a single determinant (see earlier).

In the limit that the coupling between the single and double excitation $H_{qd}$ goes to zero, Eq. 11 correctly reduces to the adiabatic approximation chosen for $f_{\scriptscriptstyle\rm HXC}^{\rm A}$. When $H_{qd}$ is much smaller than all the other terms, then one frequency is a slightly corrected adiabatic value, while the other reduces to $\sqrt{(H_{dd}-H_{00})^{2}+H_{qd}^{2}}$. This motivates a variation of Eq. (11) where we replace the diagonal matrix element differences with KS values:

Another possibility is to replace them with their adiabatic TDDFT counterparts:

where $\Omega_{q,s1,s2}^{A}$ are the adiabatic TDDFT frequencies that correct the bare KS frequencies $\nu_{q,s1,s2}$. Eq. (13) and (14) have a numerical advantage in requiring less two-electron integrals than in Eq. (11) but all three flavours of DSMA capture excitation energies and oscillator strengths of double excitations as we will now demonstrate on explicit examples.

Figure 1 shows the performance on an exactly-solvable model system: two electrons in a one-dimensional harmonic plus linear potential, $v_{\rm ext}(x)=\frac{1}{2}x^{2}+\gamma|x|$ where $\gamma$ is a parameter in the range$[0,1]$; varying $\gamma$ 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: $\frac{1}{\sqrt{(x_{1}-x_{2})^{2}+1}}$. We observe in Fig. 1(a) that for small enough $\gamma$ 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 DSMA approximations 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_{\scriptscriptstyle\rm HXC}^{\rm 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 DSPA (see the Supplementary Material). However, the real and significant improvement of DSMA over DSPA lies in its capacity to accurately determine oscillator strengths. Fig. 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 $\gamma$ and the DSPA fractions shown in the Supplementary Material. All three flavors of DSMA respect the oscillator strength sum rule, Eq. (6) as shown in the figure, in contrast to the DSPA (see Supplementary Material for the analogous plot).

We next turn our attention to the Beryllium atom. Here the lowest two singlet ${}^{1}D$ states, have a mixed single and double excitation character. Ref. Loos et al. (2019) reports a roughly $30\%$ single-excitation character to the predominantly doubly-excited $1^{1}D$ state ($1^{1}S(2s^{2})\rightarrow 1^{1}D(2p^{2})$. As a result, adiabatic TDDFT fails to describe the states accurately. Figure 2 shows a plot of the two excitation frequencies of the $1^{1}D$ multiplet. Adiabatic TDDFT with the 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. Still, DSMA₀ and especially DSMA_S yield oscillator strengths which are close to that of the reference Ref. Loos et al. (2019) 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.

Finally we turn to the LiH molecule at stretched bond lengths, which we will use to illustrate the use of our DSMA to calculate coupling matrix elements between excited states. Adiabatic approximations to calculate these within quadratic response are plagued with unphysical divergences that appear when the difference between two excitation energies $\Omega_{a}$, and $\Omega_{c}$ coincides with another excitation frequency $\Omega_{b}=\Omega_{c}-\Omega_{a}$ Parker, Roy, and Furche (2016); Li and Liu (2014); Zhang and Herbert (2015); Ou et al. (2015) (see Supplementary Information for a figure of this). In Ref. Dar, Roy, and Maitra (2023) a frequency-dependent approximation for the second-order response kernel, $g_{\scriptscriptstyle\rm XC}^{\rm App}$, was derived in a truncated Hilbert space containing these states, that eliminated these divergences. This kernel was used to derive the transition density $n_{ca}({\bf r})=\langle\Psi_{c}|\hat{n}({\bf r})|\Psi_{a}\rangle$ for a case where $\Omega_{a}\approx\Omega_{b}\approx\Omega_{c}/2$:

Here $\nu_{1}$ and $\nu_{3}$ are the frequencies of the dominant KS states contributing to the interacting levels $\Omega_{a}$ and $\Omega_{c}$ respectively, $\langle f_{\scriptscriptstyle\rm HXC}(\Omega_{a})\rangle_{1}=\int\phi_{0}({\bf r})\phi_{1}({\bf r})f_{\scriptscriptstyle\rm HXC}({\bf r},{\bf r}^{\prime},\Omega_{a})\phi_{0}({\bf r}^{\prime})\phi_{1}({\bf r}^{\prime})d^{3}{\bf r}d^{3}{\bf r}^{\prime}$. In the case considered, state $a$ has predominantly single-excitation character, and a frequency-independent kernel suffices for this kernel. The third term involves $g^{\rm adia}_{{\scriptscriptstyle\rm XC}}({\bf r}_{1},{\bf r}_{2},{\bf r}_{3})=\left.\frac{\delta^{3}E_{\scriptscriptstyle\rm XC}[n]}{\delta n({\bf r}_{1})\delta n({\bf r}_{2})\delta n({\bf r}_{3})}\right|_{n=n_{0}}$, an adiabatic approximation to the second-order response kernel. The first term is the contribution from a double excitation, which should be weighted by the doubles contribution to state $c$, which, in the limit of weak enough interaction is $\sqrt{1-|G_{c}|^{2}}$, however until the present work, how to calculate this was unknown. We note that Eq. (15) differs a little from that given in Ref. Dar, Roy, and Maitra (2023): there we had instead proposed to approximately weigh the doubles contribution with $\sqrt{1-\alpha_{c}^{2}}$, where $\alpha_{c}$ is a spatially-independent approximation to $\sqrt{a_{{\scriptscriptstyle\rm S},3}^{-1}a_{c}}$ with $a_{c}({\bf r},{\bf r}^{\prime})=\langle\Psi_{0}|\hat{n}({\bf r})|\Psi_{c}\rangle\langle\Psi_{0}|\hat{n}({\bf r}^{\prime})|\Psi_{c}\rangle$ and likewise $a_{{\scriptscriptstyle\rm S},3}$ is the product of the KS transition densities where we use the subscript $3$ to denote the KS state that the TDDFT state $c$ is dominated by, but the weights in Eq. (15) are better justified in the weak-interaction limit. In the specific application to the LiH transition dipole in Ref. Dar, Roy, and Maitra (2023), reproduced in the Supplementary Material here, only the second term was computed (with $G_{c}=1$), since the underlying linear response TDDFT calculation was done in the adiabatic approximation over all frequencies, so would not detect any double-excitation contribution. In reality, a KS double-excitation lies near $\Omega_{c}$ (see Supplementary Material and also Fig. 3 shortly). The approximation of Ref. Dar, Roy, and Maitra (2023) has an overall trend that follows the exact results, but becomes an overestimate to the left of the resonance region while underestimating it to the right.

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 1st and 4th singlet excited states in LiH, $\mu_{14}=\langle\Psi_{1}|\hat{x}|\Psi_{4}\rangle=\int xn_{ca}^{\rm App}({\bf r})d{\bf r}$, where $x$ is along the internuclear axis. In Fig. 3 all calculations are done in the aug-cc-pVDZ basis set using PySCF Sun et al. (2018), while the pseudo-wavefunction calculation is performed with the Turbomole package tur  ¹¹1 The pseudo-wavefunction result was obtained using PBE0/aug-cc-pVDZ within the Turbomole package. To get smooth results in this basis, the computational parameters used were: SCF convergence threshold $<10^{-12}$; TDDFT residual threshold $<10^{-9}$. 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\nu_{1}^{\rm APBE0}$ with the 4th KS excitation frequency, justifying the need for the DSMA approach. The black lines demonstrate the “resonance condition" for FCI, $\Omega_{4}=2\Omega_{1}$, while near the intersection of the corresponding curves for APBE0 (blue lines), $\mu_{14}^{\rm APBE0}$ diverges (see Refs. Parker, Roy, and Furche (2016); Dar, Roy, and Maitra (2023) and Supplementary Material). The middle panel shows the 4th and 5th 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 4th state, we retrieve two states, $\Omega_{4}^{(1)}$ and $\Omega_{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 $\mu_{14}$ computed from Eq. (15) ($\Omega_{c}=\Omega_{4}^{(1)})$ using $G_{c}$ obtained from DSMA₀ and DSMA_S, against the reference FCI and the pseudo-wavefunction approximation Li and Liu (2014); Zhang and Herbert (2015); Ou, Alguire, and Subotnik (2015); Ou et al. (2015); Parker, Roy, and Furche (2016). 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\sim[2.6\mathrm{\mathring{A}}-3.4\mathrm{\mathring{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.

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. Also, 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 imply a frequency-dependence and an analogous formula for the generalized KS framework may be needed.