Density Functional Theory for Electronic Excited States

This chapter provides a basic introduction to excited-state extensions of density functional theory (DFT), including time-dependent (TD-)DFT in both its linear-response and its explicitly time-dependent formulations. As applied to the Kohn-Sham DFT ground state, linear-response theory affords an eigenvalue-type problem for the excitation energies in a basis of singly-excited Slater determinants, and is widely known simply as"TDDFT"despite its frequency-domain formulation. This form of TDDFT is the mostly widely-used quantum-chemical method for excited states, due to a favorable combination of low cost and reasonable accuracy. The chapter surveys the accuracy of linear-response TDDFT, which is generally more sensitive to the details of the exchange-correlation functional as compared to ground-state DFT, and also describes some known systemic problems exhibited by this approach. Some of those problems can be corrected on a case-by-case basis using orbital-optimized, excited-state self-consistent field (SCF) calculations, in what is known as excited-state Kohn-Sham theory or a"Delta-SCF"procedure, a class of methods that includes restricted open-shell Kohn-Sham theory. Recent successes of these approaches are highlighted, including double excitations and core-level excitations. Finally, explicitly time-dependent (or"real-time") TDDFT involves propagation of the molecular orbitals in time following an external perturbation, according to the Kohn-Sham analogue of the time-dependent Schroedinger equation. The time-dependent approach has been used to model strong-field electron dynamics, and in the weak-field limit it provides a route to broadband spectra based on the time evolution of the dipole moment function. This is useful for describing high-energy excitations (as in x-ray spectroscopy) and in systems where the density of states is high, as demonstrated by a few examples.

Following its implementation in molecular quantum chemistry codes in the early 1990s, 1-6 density functional theory (DFT) quickly emerged as the most popular tool for ground-state electronic structure calculations due to its favorable balance of relatively low cost with reasonable accuracy for thermochemistry. The first excited-state implementations quickly followed, 7-12 based on a linear-response (LR) formalism [13][14][15] that mirrors much earlier work on time-dependent Hartree-Fock theory. 16 (The historical development of TDDFT has been summarized elsewhere. 17 ) The LR formulation is now known almost universally as "timedependent" (TD-)DFT, despite its frequency-domain formulation and implementation. In its most pedestrian applications, LR-TDDFT produces vertical excitation energies for closed-shell molecules at ground-state geometries to within a statistical accuracy of ∼ 0.3 eV, 18 at a cost that is often only a few times more than the cost of the ground-state self-consistent field (SCF) calculation and possessing the same formal scaling. 19 This is a useful accuracy for electronic absorption spectra. In view of its low cost, LR-TDDFT has therefore become the de facto method for computing electronic excitation spectra of finite molecular systems, although some fundamental problems remain in its application to periodic materials. [20][21][22] LR-TDDFT is also becoming increasingly popular for photochemical applications, [23][24][25] despite some problems with the description of conical intersections. [26][27][28] In part, this popularity is due to a growing recognition that complete active-space (CAS-)SCF methods cannot be considered quantitative approaches for excited-state dynamics, [29][30][31][32] due to a lack of dynamical electron correlation effects. This chapter provides an overview of TDDFT and other DFT-based methods for computing excitation spectra, excited-state properties, and for simulating photochemical reactions, emphasizing theory rather than applications but with some molecular examples to motivate the discussion. For those unfamiliar with the formal underpinnings of TDDFT, a natural question to ask is "what does time have to do with excitation energies?" In fact, one knows from basic quantum mechanics that the time evolution of a non-stationary wave function encodes the system's excitation energies via the Bohr frequencies, ω jk = (E j −E k )/ , therefore in principle the time evolution of a quantum system can be used to extract excitation energies. The existence of a time-dependent extension of DFT is formally justified by the Runge-Gross theorem, [33][34][35][36][37][38][39] which provides something akin to a time-dependent extension of the first Hohenberg-Kohn theorem for the ground state, 40 i.e., a density-to-potential mapping. In the time-dependent case there are important caveats about initialstate dependence and memory effects. [41][42][43] Those issues have yet to be fully resolved in a computationally feasible way, but this has not stymied the practical development and application of TDDFT.
Following a perturbation to the ground state, which creates a superposition of energy eigenstates, the Fourier components of the time-dependent dipole moment are precisely the Bohr frequencies. A Fourier transform of the dipole moment function is itself an absorption spectrum, 44 Excitation energies are also encoded in the isotropic frequency-dependent polarizability, α(ω), which has a sum-over-states expression α(ω) = e 2 m e n>0 f 0n ω 2 n0 − ω 2 (2) where m e is the electron mass, ω n0 = (E n − E 0 )/ , and is the dimensionless dipole oscillator strength for the |0 → |n transition. 44 The poles of response function α(ω) therefore encode excitation energies, with residues that encode oscillator strengths. 19 In the early days of quantum chemistry, Eq. (2) was actually used to compute excitation energies for atoms and atomic ions, 45,46 by computing α(ω) as the response to an applied field, and a version of this approach would eventually reappear in the form of "real-time" TDDFT. 47,48 The poles of the Kohn-Sham response function also serve this purpose, 9,49,50 and the LR formalism applied to the Kohn-Sham ground state turns this idea into a robust computational paradigm, in the form of an eigenvalue-type problem for the excitation energies. [13][14][15] Although the LR formulation exists strictly in the frequency or energy domain, the timedependent origins of the phenomenology have persisted in the name "TDDFT". Despite its overwhelming popularity, LR-TDDFT excitation energies do tend to be more sensitive to the details of the exchange-correlation (XC) functional as compared to ground-state properties computed with DFT. In some sense, the statistical accuracy of ∼ 0.3 eV that is quoted above should therefore be interpreted as representative of the best-case scenario with state-of-the-art functionals, and assuming that certain systemic pathologies can be avoided. LR-TDDFT may not be the theory that we want, but it remains the best theory that we have for excited states of large and even medium-size molecules. This theory is introduced formally in Section 2 and that discussion constitutes the most substantial part of this chapter, just as LR-TDDFT occupies the most significant place amongst excited-state DFT methods. It holds that position because it is easy to use, not significantly more expensive than ground-state DFT, and provides a slew of excited states in an automated way, starting from a ground-state SCF solution.
While the accuracy of LR-TDDFT is often quite reasonable, certain systemic problems have been identified and excited-state Kohn-Sham procedures have been developed to circumvent these. Rather than applying LR to the ground state, these methods look for an excited-state (non-aufbau) solution to the SCF equations, and for this reason the excited-state Kohn-Sham approach is often called a "∆SCF" method. Although not formally justified by the Runge-Gross theorem, the ∆SCF approach has an admirable record of rectifying the deficiencies of LR-TDDFT, again at a cost comparable to that of a ground-state DFT calculation. What is lost in the ∆SCF approach is the ability to compute a whole spectrum of states at once, making the state-specific ∆SCF procedure much more labor-intensive for the user. This approach is described in Section 3.
Finally, it is possible to take the time dependence in TDDFT at face value and to propagate Kohn-Sham molecular orbitals (MOs) in time, following a perturbation applied to the ground-state density. This is accomplished by solving the time-dependent Kohn-Sham (TDKS) equation, which is a one-electron analogue of the time-dependent Schrödinger equation. (Here, σ ∈ {α, β} is a spin index.) The one-electron effective Hamiltonian in Eq. (4) is the Fock operatorF σ that comes from the ground-state Kohn-Sham eigenvalue problem that determines the MOs: F σ ψ kσ (r) = ε k ψ kσ (r) .
The "real-time" approach to TDDFT, 51,52 which is described in Section 4, consists in solving Eq. (4) by propagating the MOs in time following a perturbation to the ground-state that creates a time-evolving density, expressed here for σ-spin electrons. (The total charge density is ρ = ρ α + ρ β .) This approach can be used to simulate strong-field electron dynamics, 53 which is a topic of contemporary interest in attosecond molecular science. [54][55][56][57] It also provides a route to broadband spectra via Fourier transform of the time-dependent dipole moment function, in a direct realization of Eq. (1). This chapter assumes a basis familiarity with ground-state DFT, as represented by the SCF eigenvalue problem in Eq. (5), which will serve as our starting point. It should therefore be familiar that the Fock operator takes the formF (7) in atomic units. The quantities v ext , v H , andv σ xc are known as the external, Hartree, and XC potentials, respectively. In the field-free case, the external potential is simply the interaction of the electrons with the nuclei, 58 v ext (r) = − More generally, v ext (r) might also contain a field-dependent contribution such as −E(r, t)·r in the presence of an electric field E(r, t). The Hartree (or Coulomb) potential v H (r) describes self-repulsion of the electrons, 58 equivalent to what is often called "J" in Hartree-Fock theory. 40,59 It is a functional of the density, given by v H [ρ](r) = ρ(r, r ) r − r dr .
The final component ofF σ isv σ xc = δE xc /δρ σ , the XC operator for σ-spin electrons. In "pure" Kohn-Sham theory, this quantity should be a local potential v σ xc (r) rather than an operator, but herein we allow the possibility for mixing some nonlocal Hartree-Fock exchange (HFX), as is done in the hybrid density functionals that are most useful in molecular DFT. For hybrid functionalsv σ xc is a nonlocal operator and this scenario is sometimes called generalized Kohn-Sham theory, 39 although the use of hybrid functionals can no longer be considered exotic in molecular DFT.
The textbook by Koch and Holthausen 40 is a good resource for ground-state DFT (though not for TDDFT), as are several literature overviews. 58,60 Updated ground-state benchmarks, relative to the rather dated ones in Ref. 40, can be found elsewhere. 61,62 For TDDFT, the textbook by Ullrich,63 or else overviews by Gross and co-workers, 36,64,65 provide the rigorous foundations of the theory, which are mostly omitted here. Several other reviews cover LR-TDDFT in a pedagogical way. [66][67][68] For overviews of molecular applications of LR-TDDFT, see reviews by Jacquemin and co-workers, [69][70][71][72][73] who have also reviewed the accuracy benchmarks 18 and functional selection. 74 2 Linear-Response ("Time-Dependent") DFT This section describes the formalism and application of LR-TDDFT, commonly known simply as "TDDFT". The starting point is the TDKS equation [Eq. (4)] that describes how the ground-state MOs ψ kσ evolve in time following a perturbation that is applied at t = 0. If that perturbation is taken to be a time-oscillating field at frequency ω, then in the weak-field limit (E → 0), the response of the ground state can be computed exactly using firstorder perturbation theory. 15 Formally, one ought to show that the poles of the frequency-dependent response function can be obtained from those of the independent-particle (Kohn-Sham) response function, 9 but for that exercise the reader is referred to reviews by Marques and Gross. 64,65 For a derivation of LR-TDDFT based on a variational principle, see Ref. 75.

Theoretical Formalism
The derivation from perturbation theory starts from the equivalent Liouville-von Neumann (LvN) form of the TDKS equation, which is is the time-evolving one-electron density operator for σ-spin electrons. Expanding Eq. (11) to first order in the perturbed Fock and density matrices, in the presence of the perturbation V (t), one obtains the unperturbed LvN equation at zeroth order. This is equivalent to the ground-state Kohn-Sham eigenvalue problem in Eq. (5). Working equations for LR-TDDFT are obtained by equating the first-order terms, 11,15,66 as elaborated below.

Linear Response Theory
To consider this in more detail, recognize that the perturbation V (t) in Eq. (10) is a one-electron operator whose spatial part can be expanded in the MO basis, leaving the time dependence to be carried by e ±iωt . Introducing a set of unknown coefficients z pqσ and z * qpσ , representing real and imaginary parts of the firstorder response, the first-order perturbed density matrix can be expressed as where P (0) pqσ is the unperturbed density matrix at t = 0. This change in the density matrix is accompanied by a corresponding change in the Fock matrix. Through first order, the Fock matrix is 11 where the unperturbed Fock operatorF (0) σ has the form given in Eq. (7). The first-order response of the density matrix is thus coupled to a term of the form 11 The first term, (p σ q σ |s τ r τ ) = (pq|sr), is a Coulomb integral expressed in Mulliken notation, 59 whilef στ xc = δ 2 E xc /δρ σ δρ τ . The latter quantity is discussed in more detail below.
So far, the MO indices p, q, r, s are arbitrary and could refer either to occupied or virtual orbitals. In fact, the idempotency conditionP 2 σ =P σ imposes restrictions. Through first order, the idempotency condition is σ . As a matrix, P σ contains only occupied-occupied and virtual-virtual blocks because the occupied-virtual block vanishes as a condition of SCF convergence. 76 Using i, j, . . . to index occupied MOs and a, b, . . . for virtual MOs, this means P (0) aiσ , so the constraint in Eq. (16) implies that the only non-vanishing coefficients in P (1) pqσ are z iaσ and z aiσ . 11,66 Conventional LR-TDDFT notation is obtained by relabeling these coefficients as Collecting these unknowns into vectors x and y, one may rewrite the first-order terms in the LvN equation in matrix form as [66][67][68] This represents a system of equations for the excitation energies ω n and the amplitudes x (n) iaσ and y (n) iaσ , and constitutes the basic working equations of LR-TDDFT. (The index n, which counts excited states, will usually be omitted for compactness.) The system in Eq. (18) is often called the Casida equations, 13,14 although these are formally identical to the equations of time-dependent Hartree-Fock theory. 16 The matrices A and B are known as orbital Hessians, 19 for reasons that are discussed below, and they derive from the derivative ofF σ with respect toP τ in Eq. (15). In the canonical MO basis, the matrix elements of A and B are 19,66 whereκ στ xc =f στ xc − α hfx (δ 2 E HFX /δρ σ δρ τ ). The quantity α hfx will be used throughout this chapter to mean the coefficient of HFX (often called "exact exchange") contained in the functional E xc [ρ], with 0 ≤ α hfx ≤ 1. For example, α hfx = 0.2 for B3LYP, meaning meaning 20% HFX and 80% semilocal exchange. We have chosen to separate the HFX terms in Eq. (19), which can be expressed in terms of electron repulsion integrals (ij|ab) and (ib|aj), leavingκ στ xc as the second functional derivative of the semilocal contribution, E xc − E HFX . The solution (x, y) of Eq. (18) parameterizes the transition density matrix for the excitation in question. In real space, this quantity is 10,15,19 The unknowns x and y satisfy a bi-orthogonal normalization condition, 15 which is also a feature of the much older time-dependent Hartree-Fock (TDHF) theory. 16 For historical reasons, TDHF is also known as the random phase approximation (RPA), 77,78 because it can be derived within an equation-of-motion formalism for the single-particle excitation operators, 68 similar to the historical RPA. 77 However, TDHF/RPA can also be considered to be a special case of LR-TDDFT corresponding to the Hartree-Fock functional, i.e., α hfx = 1 andκ στ xc ≡ 0. The number of unknown amplitudes in Eq. (18) is 2n occ n vir , hence to solve this equation for all of the excitation energies ω would incur sixth-order cost, O(n 3 occ n 3 vir ). Because matrix-vector products such as Ax or By can be computed with only fourth-order cost, in practical calculations Eq. (18) is solved iteratively for just the lowest few eigenvalues (n roots ). 7,10,[79][80][81][82] The cost of that calculation scales as n roots × O(n 2 occ n 2 vir ), 19 which is typically not significantly greater than the cost of the ground-state SCF calculation if n roots ∼ 10. Therefore if ground-state DFT is feasible then LR-TDDFT is probably tractable also, at least for the lowest few excited states. It is worth noting, however, that the memory footprint to solve Eq. (18) is n roots ×O(n occ n vir ), which is significantly more than the ground-state memory requirement. This can become a problem for large systems if a large number of excited states is desired, e.g., in models of semiconductors where a band structure is emerging. 83 For such applications, the real-time approach that is described in Section 4 offers a low-memory alternative to LR-TDDFT.
Some alternative forms of the basic LR-TDDFT equation are also worth considering. We first note that the matrices A and B in Eq. (19) can be rewritten as where K is a coupling matrix, 13,84 with a Hartree-XC kernel. 67 In practice, this looks like the energy-transfer coupling 85 between transition densities ρ iaσ (r) = ψ iσ (r) ψ aσ (r) and ρ bjτ (r ) = ψ bτ (r ) ψ jτ (r ). One can therefore consider that solution of the LR equations reveals how the zeroth-order, independent-particle excitations ψ iσ → ψ aσ are coupled to obtain excited states of the interacting system. Assuming that the orbitals are real, so that A * = A and B * = B, then Eq. (18) is equivalent to a pair of equations which makes it clear that Eq. (18) is not a conventional eigenvalue problem. However, upon introducing new variables which satisfy the more conventional normalization condition z † ± z ± = 1, 86,87 the LR-TDDFT equations can be transformed into either of two equivalent, Hermitian eigenvalue problems. 81,[86][87][88][89] These are where The Ω + version of Eq. (26) is especially convenient for semilocal functionals (α hfx = 0), because in that case A − B is diagonal and one obtains a Hermitian eigenvalue problem with half the dimension of the original pseudo-eigenvalue problem in Eq. (18). For a closed-shell (spin-restricted) ground state, another important transformation is which affords amplitudes for singlet (+) and triplet (−) spin functions. Making use of the unitary transfor- in addition to Eq. (28), one obtains singlet and triplet versions of A ± B that function as stability matrices. 7,8,86 In other words, these are Hessian matrices whose eigenvalues characterize whether the ground state is stable with respect to orbital rotations. For example, the singlet stability matrix is 86 A negative eigenvalue in A + + B + indicates an instability, which looks like a negative excitation energy from the standpoint of LR-TDDFT. This correspondence is a consequence of the Thouless theorem, 90 which states that orbital rotations (and therefore orbital relaxation) can always be parameterized as single excitations. Along similar lines, eigenvalues of the triplet instability matrix 86 indicate whether the ground-state solution is stable with respect to spin symmetry breaking (restricted → unrestricted). 91 In the presence of an unstable reference state, the transformation in Eq. (25) may become problematic, which can lead to failure of certain iterative LR-TDDFT algorithms.

Adiabatic Approximation
We have not yet discussed the key ingredient in the orbital Hessian matrices that makes LR-TDDFT different from TDHF/RPA, namely, the exchange-correlation kernel,f στ xc . A more careful application of LR theory would note that the quantities A(ω) and B(ω) are themselves functions of the excitation energy ω. [13][14][15] In wave function terms, this can be understood based on the fact that any exact theory of many-electron excitation energies that is formulated as an effective single-particle theory must ultimately invoke an effective Hamiltonian that is energy-dependent, in order to encapsulate the effects of higher-order excitations. 14,92,93 (In many-body theory this energy-dependent contribution is sometimes called the "self-energy". 94 ) Proofof-concept models for an energy-dependent kernel f στ xc (r, r , ω) have been put forward, 42,43,93,95 which have close connections to many-body perturbation theory and the Bethe-Salpeter equation. 92,93 However, there are no production models for molecular Hamiltonians at present.
To appreciate the nature of the approximation in neglecting the energy dependence off στ xc , consider that this quantity arises in Eq. (15) as the second functional derivative of the XC energy with respect to the density, or the first derivative of the XC potential. For a time-evolving density ρ σ (r, t), this means 41 (This expression leavesv σ xc in the form of an operator, which technically makes this an example of generalized Kohn-Sham theory. 39 ) The time dependence in f στ xc (r, r , t − t ) means that this quantity depends on the whole history of the time evolution of the density, 41-43 imparting a frequency dependence upon Fourier transformation: f στ xc (r, r , ω). For practical purposes, it is basically a requirement to invoke the adiabatic approximation, 13,68,96 which assumes locality in time and therefore differentiates with respect to the instantaneous density: 13 δv The "memory" of the kernel is thereby neglected, tantamount to assuming thatv xc [ρ](r, t) can be approximated using a conventional ground-state energy functional E xc [ρ], whose functional derivative is evaluated using the time-evolving density: 96v Time dependence is thus carried entirely by the density and not by the functional. The frequency dependence off στ xc disappears and conventional (ground-state) density functionals are all that is required for LR-TDDFT within the adiabatic approximation.
One immediate ramification of this approximation is that the LR-TDDFT equations have precisely 2n occ n vir solutions for the excitation energy ω, coinciding with the number of unknown amplitudes x iaσ and y iaσ . In wave function language, these are the "one-particle, one-hole" (1p1h) states, as in conventional configuration interaction with single excitations (CIS). States with significant double-excitation character (2p2h states) are either absent altogether, [95][96][97][98][99] or at best severely shifted. 98 The latter are therefore generally considered to be out of reach within the adiabatic approximation that is ubiquitous in practical LR-TDDFT calculations. 99

Tamm-Dancoff Approximation
Given a ground-state functional E xc [ρ], all that is required for LR-TDDFT within the adiabatic approximation are second functional derivativesf στ xc (r, r ) = from which the matrix elements of A and B can be evaluated. Upon solution of Eq. (18) or one of its equivalent forms, it is often found that the amplitudes y iaσ are 10 2 -10 3 times smaller than the largest x iaσ . Invoking the approximation y iaσ ≈ 0, one obtains a conventional Hermitian eigenvalue problem whose dimension is half that of the original LR-TDDFT pseudo-eigenvalue problem, and where the matrix B does not appear. The basis for this approximation can be understood from the fact that the matrix elements of A are typically much larger (at least along the diagonal) as compared to the matrix elements of B, because the leading contribution to A is a difference in one-particle energy levels (A iaσ,iaσ = ε aσ − ε iσ + · · · ). For historical reasons that are related to a similar approximation that is made in nuclear physics, 63 neglect of y is known as the Tamm-Dancoff approximation (TDA). For hybrid functionals, the reduction in dimension leads to a concomitant reduction in cost although for semilocal functionals the same reduction in dimensionality can be achieved using the Ω + version of Eq. (26). For the Hartree-Fock functional (α hfx = 1 and no correlation), Eq. (36) is equivalent to the CIS eigenvalue equation. 100 Excited-state wave functions in CIS are linear combinations of singly-excited Slater determinants |Ψ aσ iσ , and for this reason we identify the variables x iaσ as excitation amplitudes. The neglected amplitudes y iaσ represent de-excitation, insofar as TDHF/RPA was originally introduced in the nuclear and many-body physics literature as a means to add correlation to the ground state. 16,77 In fact, LR-TDDFT in the form of Eq. (18) was introduced in molecular quantum chemistry as the "DFT random phase approximation". 101,102 For that reason, solution of Eq. (18) or its equivalents, without invoking the TDA, is sometimes called "RPA". However, in view of a resurgence of interest in using the RPA formalism as a means for correlating the ground state, [103][104][105][106][107][108][109][110] it is better to refer to Eq. (18) as "full" LR-TDDFT, if there is a need to distinguish it from the TDA version in Eq. (36). Quantitatively, the impact of the TDA on excitation energies is often < 0.1 eV, 11 though a potentially detrimental impact is that oscillator strengths within the TDA no longer satisfy the Thomas-Reiche-Kuhn sum rule, 44 This constraint is satisfied by the full LR-TDDFT approach, at least in the complete-basis limit. 13,15,111 Some small-molecule tests suggest that errors incurred by the TDA are relatively mild (< 15%), 112 and perhaps not noticeable in band shapes once vibrational broadening is taken into account. 113 A more important consequence of the TDA is that it decouples the excitation energy problem from the ground-state stability problem, whereas for full LR-TDDFT a triplet instability in the ground state manifests as a negative excitation energy, or as an imaginary root of the Hermitian eigenvalue problem in Eq. (26). This may cause problems for eigensolvers that implicitly assume ω > 0. Note that triplet instabilities are associated with spin-symmetry breaking, i.e., with the emergence of an unrestricted solution that is lower in energy than the restricted solution. Where they appear, these instabilities cause significant artifacts in potential energy surfaces computed using LR-TDDFT, including the appearance of spurious cusps. 84,[114][115][116][117] In contrast, the variational nature of the CIS-type equation, as opposed to the pseudo-eigenvalue problem that characterizes full LR-TDDFT, prevents this from happening within the TDA. 68 Along similar lines, it has been appreciated for a long time that TDHF specifically is prone to triplet instabilities. 90,91,[118][119][120][121][122][123][124] In fact, the appearance of imaginary excitation energies at equilibrium geometries of small molecules led Furche and Ahlrichs to conclude that this method is "rather useless. . . for the investigation of excited potential energy surfaces". 125 In contrast, spin-symmetry breaking near the equilibrium geometry is often significantly mitigated when DFT is substituted for Hartree-Fock theory. 126 Because most molecular LR-TDDFT calculations use hybrid functionals that include some fraction of HFX, it can be expected that problems with triplet instabilities may increase as that fraction increases, which is precisely what is found in practice. [127][128][129][130][131] Similarly, in calculations using range-separated functionals, which incorporate HFX at long range in the electron-electron Coulomb potential, these instabilities are sensitive to the length scale on which that mixing is introduced. 132-137 Invoking the TDA thus improves the accuracy of triplet excitation energies. [134][135][136] For photochemical problems, where exploration of excited-state potential energy surfaces is paramount, Casida et al. suggest that the TDA is effectively mandatory, 138 in order to avoid excitation energies that drop to zero (and then become imaginary) as the system moves through a Coulson-Fischer point where spin-symmetry breaking occurs. Instabilities appear to proliferate as one moves away from the ground-state geometry on an excited potential surface. [138][139][140] For example, in the photochemical ringopening reaction of oxirane (C 2 H 4 O), 51% of configuration space is estimated to exhibit an instability with semilocal DFT, as compared to 93% of space with B3LYP. 138

Analytic Gradients
Photochemical simulations require analytic excited-state gradients. That formalism, which is closely connected to response theory for optical properties, 141 is not discussed here but can be found elsewhere. 19,125,[142][143][144][145] Nonadiabatic or "derivative coupling" vectors between excited states, 27,146 which are needed for nonadiabatic molecular dynamics simulations, [23][24][25][26][27] have also been formulated. [147][148][149][150][151][152][153][154] Evaluation of the nonadiabatic couplings has the same formal complexity as evaluation of the excited-state gradient. 25,147 Here, we comment briefly regarding how the gradient formalism bears on static properties of the excited state, such as the dipole moment or atomic population analysis. The density matrix for the excited state can be written which is sometimes called the "relaxed" density matrix, from which an "unrelaxed" contribution P unrlx = P 0 + ∆P is separated out. Here, P 0 represents the ground-state density matrix and the unrelaxed change upon excitation (∆P) can be obtained from the amplitudes x and y. The remaining contribution (Z) represents orbital relaxation. 19,144,145 The unrelaxed density change ∆P can be separated into particle (electron) and hole contributions, which are given by 125,145,155,156 These expressions can be rearranged to afford 155 Note that ∆P only contains occupied-occupied (∆P hole ) and virtual-virtual (∆P elec ) contributions to the excited-state density matrix, not occupied-virtual contributions. The latter are contained in Z, 19,125 the evaluation of which requires solution of so-called Z-vector equations. 157 This has the same formal complexity as an excited-state gradient calculation. Excited-state properties should be computed using the relaxed density matrix P rlx , because ∆P and Z make similar contributions. 100,158,159 Especially when the change in density is large, as for an excitation with significant charge-transfer (CT) character, the use of the unrelaxed density matrix may lead to unacceptable errors in excited-state properties. For example, excited states of p-nitroaniline involving intramolecular CT character exhibit relaxed and unrelaxed dipole moments (computed using P rlx versus P unrlx , respectively) that differ by more than 10 D in some cases! 158

Performance and Practice
As discussed above, the formal scaling of LR-TDDFT is n roots ×O(n 4 basis ) for hybrid functionals. 19 In practical terms, where only a few low-lying excited states are desired, this means that LR-TDDFT is generally feasible if the corresponding ground-state calculation is practical, perhaps up to about 400 atoms for single-point calculations or 150-250 atoms where excited-state geometry optimization is required, 70,71 with more severe limitations where excited-state frequency calculations are required. 71,160 These estimates are appropriate where basis sets of double-ζ quality are used, which is generally adequate. Triple-ζ basis sets may be considered to be converged. 67 Systems that possess a dense manifold of excited states (e.g., semiconductors) can create significant storage bottlenecks for the trial vectors that are required by the iterative eigensolver. In such cases, the "real-time" TDDFT approach of Section 4 may be advantageous, although Section 2.2.1 describes ways to reduce the cost while staying within the confines of LR-TDDFT. Following a discussion of XC functional approximations in Section 2.2.2, the accuracy of LR-TDDFT excitation energies is addressed in Section 2.2.3 by considering its performance for small-molecule benchmarks. Finally, techniques for visualizing and understanding the states that are computed are discussed in Section 2.2.4.

Restriction of the Excitation Manifold
Significant cost reduction in LR-TDDFT calculations for large molecules can be achieved by neglecting some of the amplitudes x iaσ , in addition to neglecting all of the amplitudes y iaσ . Figure 1 shows examples of excitation spectra computed for C 60 and for C 119 H 154 ClN 21 O 40 , performed by excluding over 70% of the virtual orbitals (based on orbital energies ε aσ ) without adverse effects on the spectral envelope. 161 Similar truncations of the excitation manifold can be used to access core-excited states. [162][163][164][165] There is significant interest in core excitations in contemporary quantum chemistry, 164-169 driven by the recent availability of tabletop laser sources with femtosecond time resolution. [170][171][172][173][174] However, core-to-valence excitations lie embedded in an ionization continuum and, at a practical level, lie above all of the valence-excited states, such that the use of iterative eigensolvers is prohibitively expensive if the spectrum must be computed starting from the lowest excitation energies. By retaining only those amplitudes x iaσ for which i is a core orbital on the atom of interest, core-excited states emerge as the smallest eigenvalues and can be computed directly. This "frozen valence occupied" approximation has historically been called "core-valence separation", [175][176][177] and it introduced negligible error for K-edge excitations where ψ iσ is a 1s orbital. 177 Another strategy to access core-level excitations is energy windowing, 162,178 in which the amplitudes x iaσ are excluded unless ε aσ − ε iσ lies within the window of interest.

Exchange-Correlation Functionals
Before considering the accuracy of LR-TDDFT it is useful to introduce a paradigm for classifying various density-functional approximations, for which we use the taxonomy of "Jacob's ladder". [179][180][181] At each rung on this metaphorical ladder, the functional dependence of E xc grows more intricate, incorporating more sophisticated functionality depending on the density, its gradients, the Laplacian, etc.: In a statistical sense (and only in a statistical sense), it is true that the best functionals on the higher rungs of the ladder outperform the best functionals on the lower rungs. 61,182 These rungs map onto various inputs ρ σ ,∇ρ σ ,∇ 2 ρ σ , . . . as follows.
• Rung 1: Local Density Approximation (LDA). The baseline LDA functional comes from the uniform electron gas model in which E xc is a functional of ρ(r) only, or of ρ α (r) and ρ β (r) if the system is spin-polarized. This approach does not afford useful accuracy for molecular quantum chemistry, with errors of 60-100 kcal/mol for atomization energies 61,183 and ∼ 20 kcal/mol for barrier heights. 61 • Rung 2: Generalized Gradient Approximations (GGAs). This class of functionals include a dependence on the density gradients∇ρ σ (r). These are often called "semilocal" approximations, to distinguish them from LDA while acknowledging that in their mathematical form, GGAs remain local in the sense that v σ xc (r) is a multiplicative potential. GGA functionals significantly improve thermochemistry relative to LDA; typical errors are 10-20 kcal/mol for atomization energies 61,183 and 5-15 kcal/mol for barrier heights. 61,182 • Rung 3: Meta-GGAs (mGGAs). These are also semilocal but incorporate additional derivatives including∇ 2 ρ σ and the kinetic energy density, The function τ σ (r) is related to the electron localization function, 184 and together with∇ 2 ρ σ it can be used to express the noninteracting kinetic energy. 185 The best mGGA functionals improve upon GGA thermochemistry, with errors of 5-10 kcal/mol for atomization energies 61 and 3-6 kcal/mol for barrier heights. 61,182 It is worth noting that some mGGAs introduce a considerable number of parameters, 61 and exhibit basis-set and grid sensitivities suggesting that they may be overfitted. [186][187][188][189] • Rung 4: Hyper-GGAs. As originally defined by Perdew et al., 179,180 this category consists of functionals that incorporate "exact exchange and compatible correlation". 179 A few genuine hyper-GGAs have been put forward, 190,191 but it has proven difficult to construct correlation functionals that work well with 100% HFX. As such, the fourth rung on Jacob's ladder has effectively been redefined to mean hybrid functionals, 181 which incorporate some fraction of HFX (0 < α hfx < 1), in conjunction with a fraction 1 − α hfx of semilocal exchange. These functionals are sometimes further categorized as either hybrid GGAs or hybrid mGGAs, depending on the nature of the semilocal contribution. The best hybrid functionals exhibit errors < 5 kcal/mol for atomization energies 61 and 2-4 kcal/mol in barrier heights. 61,182 This has made hybrids the de facto choice for molecular quantum chemistry.
• Rung 5: Double hybrid functionals. 192 These incorporate a fraction of the second-order Møller-Plesset (MP2) correlation energy in addition to fractional HFX, which introduces a dependence on the virtual MOs {ψ aσ } whereas functionals on the lower rungs depend only on the occupied MOs. Although double hybrids exhibit some of the best accuracy in contemporary DFT, with errors of 1-3 kcal/mol for barrier heights, 182,192 the introduction of MP2 correlation brings with it the much slower basis-set convergence of wave function methods, as well as a question of whether the orbital-dependent MP2 term should be self-consistently optimized. (Typically, it is not. 192 ) A comprehensive list of functionals, up-to-date as of 2017, can be found in Ref. 61. The aforementioned error statistics pertain to ground-state thermochemistry whereas accuracy for vertical excitation energies computed using LR-TDDFT is considered in the next section. For double hybrid functionals, the formulation of LR-TDDFT 193 begins to look more like CIS with perturbative double excitations, a method known as CIS(D), 194 which incurs a formal scaling of n roots × O(n 5 basis ) and often requires a form of quasi-degenerate perturbation theory. 195,196 For these reasons, the application of double hybrids to excited-state problems is still in its infancy and is not discussed in this chapter.
In the context of LR-TDDFT there is one further category of functionals that warrants mention, namely, range-separated hybrid (RSH) functionals. These partition the electron-electron Coulomb interaction (r −1 12 ) into a short-range (SR) component and a long-range (LR) background, typically using the error function (erf): 1 This partition introduces parameters α, β, and µ. The latter is the range separation parameter that determines the length scale (∼ µ −1 ) of the separation between the SR and LR components of r −1 12 ; see where E x,GGA is the semilocal GGA exchange functional and α ≡ α hfx is the coefficient of HFX. The RSH functional corresponding to Eq. (45) is Quantities labeled "SR" or "LR" in these equations are evaluated using the corresponding component of r −1 12 . The idea is to correct the asymptotic behavior of a semilocal potential v σ xc (r) at long range (using HFX), while attempting to have minimal impact on the SR behavior of the GGA or hybrid GGA in question, since that functional is responsible for the favorable thermochemical predictions in the ground state. RSH functionals have become popular enough that traditional hybrids such as PBE0 or B3LYP are often called "global hybrid" (GH) functionals in contemporary parlance, to emphasize that HFX is added at all length scales in these cases. The definition of E LR x,GGA in Eq. (47) is consistent with other literature, [197][198][199] although it is worth noting that this quantity is not truly long-ranged. In fact, the reason that RSH functionals were introduced in the first place is to address the fact that semilocal exchange falls off too rapidly with distance, [200][201][202] leading to an insufficiently attractive interaction potential between a well-separated electron and hole. [201][202][203] This is further discussed in Section 2.3.1.
Nomenclature and usage for RSH functionals has become somewhat muddled and the remainder of this section attempts to clarify it. An RSH functional is any that uses a partition of r −1 12 into SR and LR components, with Eq. (44) as the most common partition although other forms have been explored, [204][205][206][207][208][209] including variants with a three-way partition of r −1 12 into short-, middle-, and long-range contributions. 210-215 For many of these functionals, the range-separation parameter(s) are optimized or fitted alongside other parameters that define the functional and should not be modified. Examples include the "ωB97" class of functionals 216-220 and range-separated versions of the "Minnesota" functionals. 221-224 These functionals do not always afford correct the asymptotic behavior of the XC potential, however. For the ansatz in Eq. (44), the proper behavior requires α + β = 1 but that constraint is sometimes violated (in particular, by the popular CAM-B3LYP functional 225 ) in the interest of obtaining more accurate excitation energies for localized transitions.
In contrast to this empirical approach to range separation, long-range corrected (LRC) functionals represent a subset of RSH functionals that are constrained to include 100% HFX in the limit r 12 → ∞. 200,[226][227][228] For a given GGA (α hfx = 0) or hybrid GGA (α hfx > 0) functional, the corresponding LRC functional is The parameter µ in Eq. (44) still controls the length scale on which LR-HFX is activated, but α + β = 1 by construction and therefore v σ xc (r) ∼ −r −1 for any µ > 0. The LRC strategy is thus to graft correct asymptotic behavior onto an existing semilocal XC functional, while doing the least possible damage to that functional at short range. Non-empirical adjustment (or "tuning") of the parameter µ is often employed in this context, especially where CT states are involved, and this is discussed in Section 2.3.1.
LRC functionals require modification of the semilocal GGA exchange functional in order to use an attenuated Coulomb potential. (HFX integrals can be modified once and for all to separate them into LR and SR contributions. 229 ) There are several routes to modify E x,GGA . The first of these, originally introduced by Hirao and co-workers, [200][201][202] modifies the exchange inhomogeneity factor that multiplies the electron-gas exchange energy density. This author has suggested that these functionals should be denoted LRC-µGGA, 230,231 where "GGA" indicates the semilocal parent functional, e.g., GGA = BLYP or PBE. Note that "LC" is another common abbreviation for long-range correction so that functionals such as LC-BLYP 200 might more descriptively be called LRC-µBLYP, in order to emphasize which SR-GGA exchange function (µBLYP) is being used.
For semilocal exchange functional such as PBE that are based on a model for the exchange hole, 232,233 an alternative strategy is to combine that model with an attenuated Coulomb potential in order to obtain E SR x,GGA . 228,233 To distinguish this from the LRC-µPBE functional constructed using Hirao's approach, the author has suggested the nomenclature LRC-ωPBE for the model based on the PBE exchange hole, 230,231 which comports with the notation for the range-separation parameter (ω) that was introduced in Ref. 233. The term LC-ωPBE is synonymous with LRC-ωPBE, and LR-ωPBEh is sometimes used to indicate a shortrange hybrid (α hfx > 0). In contrast, ωPBE refers only to the modified exchange functional, E SR x,PBE , and should not be used to mean the LRC functional because E SR x,PBE is used in other capacities. For example the HSE functional 234 (sometimes called HSE06 235 ) uses ωPBE in conjunction with SR-HFX to construct a hybrid functional that is efficient for periodic calculations.

Accuracy for Vertical Excitation Energies
There have been numerous systematic surveys of the accuracy of various XC functionals for use in LR-TDDFT, [236][237][238][239][240][241][242][243][244][245] enough to have spawned a meta-review of the benchmark studies themselves. 18 Two of these studies are highlighted here, to provide some sense for how various categories of functionals can be expected to perform. As usual in DFT (and even more so in LR-TDDFT), for any given molecule it is likely that one could find some XC functional that outperforms the statistically-best approach. As such, it is only by understanding trends amongst functionals (and likely trying the same calculation with more than one functional) that results can be taken seriously. Both of the studies highlighted here compare vertical excitation energies to experimental data, and while that has the advantage of being a direct comparison against numbers that one might hope to simulate, it has the disadvantage that vertical excitation energies are not strictly measurable quantities and various effects including solvatochromatic shifts and vibrational averaging are folded into the comparison. Other studies have compared LR-TDDFT excitation energies against correlated wave function benchmarks, [236][237][238] which makes for a much more straightforward test of the theory although unfortunately such comparisons often tend to get little traction outside of quantum chemistry circles, where comparison of theory against theory is often viewed with derision. Fortunately, the trends that are highlighted here are reasonably similar to those obtained by comparing against ab initio benchmarks. In assessing the performance of various functionals, we will use the taxonomy of Jacob's ladder (Section 2.2.2) as an organizing principle.
A first set of benchmarks is depicted in Fig. 3, for a set of 101 transitions in 14 gas-phase molecules. 243 Error statistics are grouped and color-coded by category, including (GH-)GGAs and (GH-)mGGAs but not RSH functionals. Errors are further separated into singlet excitations, triplet excitations, nπ * , ππ * , and Rydberg excitations. Examining these data, it quickly becomes apparent that the GH functionals significantly outperform the semilocal ones, across all types of data, although it is less clear whether GH-mGGA functionals are categorically superior to GH-GGAs. Perhaps surprisingly, the PBE0 and B3LYP functionals outperform most other functionals, including much newer mGGAs and some GH-mGGAs of the Minnesota type, 246 although M06-2X does exhibit slightly smaller errors. The B3LYP and PBE0 functionals, which for many years have served as the closest there is to a "default" setting in molecular DFT, continue to outshine many other functionals for vertical excitation energies.
The best-performing functionals (PBE0, B3LYP, and M06-2X) exhibit mean absolute errors (MAEs) of ∼ 0.3 eV for the entire data set. Unlike other functionals examined in Fig. 3, they do not seem to be systematically worse for nπ * states as compared to ππ * states. In contrast, none of the GGA functionals has a MAE below 0.5 eV and and the semilocal mGGAs also have MAEs 0.4 eV, with M06-L as the bestperformer amongst the latter. All of the semilocal functionals perform significantly worse for nπ * excitations than they do for ππ * excitations.
The comparison between Rydberg and valence excitations in Fig. 3(d) warrants special attention. With few exceptions, errors are significantly larger for the Rydberg excitations. Significant errors in Rydberg excitation energies were noted in early molecular applications of LR-TDDFT, 247 leading to the understanding that these excitation energies are quite sensitive to the long-range behavior of the XC potential. That behavior is incorrect for almost all of the functionals evaluated in Fig. 3. Later this analysis was extended to CT excitation energies in general, 203 of which Rydberg excitations can be considered a special case insofar as these states involve excitation into a diffuse orbital, relatively far from the molecular core. As will be discussed in Section 2.3.1, this observation led to the understanding that HFX is the only component of modern functional construction that exhibits the proper asymptotic behavior for a charge-separated state, whereas semilocal XC potentials fall off much too rapidly with distance and thus significantly underestimate Rydberg and CT excitation energies. It is therefore no accident that the only functionals in Fig. 3(d) for which the valence excitation error is larger than the Rydberg excitation error are precisely the ones with the largest fractions of HFX: M06-2X (α hfx = 0.54), 248 M06-HF (α hfx = 1.0), 248 PBE0 (α hfx = 0.25), 249 and BH&HLYP (α hfx = 0.5). 250 A second statistical survey is presented in Fig. 4, taken from one of the largest statistical assessments of LR-TDDFT to date: 239 614 singlet excitation energies in 483 solution-phase organic molecules. Vertical excitation energies have been corrected for solvent effects and compared to experimental band maxima. (For a discussion of dielectric continuum solvation models and their application to LR-TDDFT, see Ref. 251.) Functionals are once again grouped by category and this larger data set makes it clear that the GH functionals generally outperform the semilocal mGGA functionals, which themselves outperform the semilocal GGAs. For most of the GH functionals, MAEs are 0.2-0.3 eV as compared to 0.4-0.5 eV for the semilocal functionals, but the mean signed errors [ Fig. 4(a)] are much smaller for the GH functionals. Signed errors are nearly zero for PBE0 and B3LYP, indicating no systematic error in these cases. In contrast, errors are much larger for GH functionals with a large fraction of HFX, including BMK (α hfx = 0.42), 252 M05-2X (α hfx = 0.56), 253 and BH&HLYP (α hfx = 0.5). 250 These large-α hfx functionals exhibit bias towards overestimation of excitation energies, whereas semilocal functionals consistently underestimate them. Also included in Fig. 4 are error statistics for a set of RSH functionals. MAEs for these functionals span a wide range from 0.2-0.5 eV and in that sense are not better than the GH functionals. Furthermore, whereas semilocal functionals systematically underestimate excitation energies, GH functionals are largely free of that bias except when α hfx > 0.4. Finally, RSH functionals systematically overestimate excitation energies, which has also been observed in more recent benchmarks for biochromophores. 236 Putting these observations together, it seems that some HFX is optimal, perhaps α hfx ≈ 0.20-0.25, with excitation energies   that are too low for smaller values and too high for larger ones. Included the the latter category are many LRC functionals that use α hfx = 1 in the asymptotic limit. With that in mind, it is interesting to compare error statistics for LC-ωPBE and LC-ωPBE (20) in Fig. 4. The former uses a range-separation parameter ω = 0.4 bohr −1 that was optimized for ground-state properties, 254 whereas in LC-ωPBE (20) this parameter is set to ω = 0.2 bohr −1 , leading to significant reduction in the errors. Attempts to fit ω using both groundstate properties as well as excitation energies typically lead to values in the range ω = 0.2-0.3 bohr −1 , depending on whether short-range HFX is present or not. 226,228,255 This is consistent with the revised choice in LC-ωPBE(20).

Visualization
Having computed an excitation energy, there are a variety of tools available to visualize the excited state in question. One could simply examine each pair of occupied and virtual MOs for which the coefficient x iaσ is large, but this is often tedious due to significant configuration mixing, especially in the virtual space. At the CIS level, it is easy to understand why the canonical MOs are not a good basis for visualization purposes, since Koopmans' theorem implies that the virtual MOs are reasonable orbitals for electron attachment, not excitation. 59 The Hartree-Fock virtual MOs feel the full repulsive potential of the N -electron charge density, whereas the occupied MOs feel only N − 1 electrons, which has the effect of rendering the virtual MOs significantly more diffuse than the occupied MOs. Often the Hartree-Fock virtual MOs are simply unbound and therefore represent discretized continuum states, 256 whose shapes are sensitive to small changes in basis set. 257 Significant configuration mixing is therefore necessary to obtain a localized valence excitation. In principle, exact Kohn-Sham MOs are a much better basis for excitations, 257-259 since both occupied and virtual MOs are subject to the same N -electron potential, and in practice it is often the case that the first few Kohn-Sham virtual orbitals are bound (ε aσ < 0). Hybrid functionals, however, push the virtual orbitals and their eigenvalues back towards the Hartree-Fock picture and even 20-25% HFX can be enough to engender significant configuration mixing due simply to the diffuseness of the virtual MOs.
This type of configuration mixing is artificial, in the sense that it can be removed via orbital rotation and therefore does not represent true multiconfigurational character in the excited state. The relevant transformations of the canonical occupied MOs is a unitary matrix U o that diagonalizes ∆P elec in Eq. (42a): The n occ × n occ diagonal matrix Λ 2 contains the eigenvalues, which are strictly non-negative and are normalized such that i λ 2 i = 1. The corresponding transformation of the canonical virtual MOs is These two transformations define the natural transition orbitals (NTOs), [260][261][262][263] which are the natural orbitals (eigenfunctions) of the excited-state density matrix. 263 They can equivalently be defined based on a singular value decomposition of the n occ × n vir matrix of amplitudes, x + y: 261,262 This form demonstrates that no more than n occ of the singular values {λ i } are non-zero. These eigenvalues appear in pairs (±λ 2 i ) when ∆P elec and ∆P hole are diagonalized, because the natural occupation numbers of the excited-state density matrix are 263 The matrices U o and U v transform the canonical occupied and virtual MOs, respectively, into a set of "hole" orbitals {ψ hole i } along with corresponding "particle" (or "electron") orbitals, {ψ elec i }. These are the NTOs, and their diminishing importance for describing the excitation in question is quantified by the values λ 2 1 > λ 2 2 > λ 2 3 > · · · . NTOs often provide a much more compact description of the wave function as compared to the canonical MOs, yet one that preserves the phase and nodal structure that can be helpful in qualitatively characterizing the nature of the excitation. This is illustrated in Fig. 5 for the S 0 → S 2 excitation of a five-unit polyfluorene molecule in which a carbonyl defect in one of the terminal fluorene monomers serves to localize the excitation. That localization, however, is not obvious from the canonical MOs, which are delocalized across the length of the molecule, but arises from a coherent superposition of four different occupied MOs. Upon transformation to the NTO basis there is only one significant singular value, with λ 2 1 = 0.96. The principle NTO pair, ψ hole 1 → ψ elec 1 in Fig. 5(b), thus paints a picture that is 96% complete. There is an unfortunate tendency in the literature to refer to ψ hole 1 as the "highest-occupied NTO" (HONTO), with ψ elec 1 then deemed the "lowest-unoccupied NTO" (LUNTO). This terminology is incorrect insofar as "highest" and "lowest" (as in HOMO and LUMO) refer to orbital energies, which are not welldefined for the NTOs because the Fock matrix is not diagonal in that representation. The HONTO/LUNTO terminology should therefore be avoided so as not to conflate visual depictions of NTOs with qualitative arguments that might be based on one-electron energy levels, which are only well-defined in the canonical MO basis. The term principle transition orbitals (or perhaps principle NTOs) is suggested instead, to describe the pair of orbitals corresponding to the largest λ i . One might therefore describe these in sequence of principle NTOs (pNTOs) as pNTO, pNTO − 1, pNTO − 2, . . . for λ 2 1 > λ 2 2 > λ 2 3 > · · · . Another common tool to visualize an excitation is the density difference as compared to the ground state. The unrelaxed density difference ∆ρ(r) = ∆ρ elec (r) + ∆ρ hole (r) has particle and hole components that are the real-space analogues of the density matrices ∆P elec and ∆P hole in Eq. (41). Using the NTOs, the particle and hole densities may be expressed as Note that ∆ρ elec (r) is positive definite and ∆ρ hole (r) is negative definite, consistent with Eq. (41). Because the NTOs are defined by a singular value decomposition, which distills the n occ × n vir matrix x + y into the fewest number of non-zero parameters, the densities in Eq. (54) are often dominated by the principle NTO pair. Although it is not widely recognized, the quantities ∆ρ elec (r) and ∆ρ hole (r) are precisely the attachment density and the detachment density, respectively, which have long been used to visualize excited states. [264][265][266] (These were originally introduced in a different way, 264 based on eigenvectors of ∆P that afford positive or negative eigenvalues, respectively.) In the author's view, NTOs are still the preferable description since phase information is lost upon squaring the orbitals in Eq. (54). Figure 6 illustrates these densities for the same S 0 → S 2 excitation of polyfluorene that was examined in Fig. 5. Because λ 2 1 ≈ 1, the particle and hole (or attachment and detachment) densities have the same information content as the principle NTO pair in Fig. 5(b). Also shown in Fig. 6 is the transition density T (r) ≡ T (r, r) where T (r, r ) is defined in Eq. (20). For an excitation |Ψ 0 → |Ψ , the general definition of this quantity is 267 T (r, r ) = N Ψ * 0 (r , r 2 , . . . , r N ) Ψ(r, r 2 , . . . , r N ) dr 2 · · · dr N , hole (detachment) density particle (attachment) density transition density difference density Figure 6: Visualization of the S0 → S2 excitation of the fluorenone-terminated polyfluorene whose orbital depiction is given in Fig. 5, represented here in terms of difference densities. These include the particle density (or attachment density) ∆ρ elec (r), the hole density (or detachment density) ∆ρ hole (r), the unrelaxed difference density ∆ρ(r) = ∆ρ elec (r)+∆ρ hole (r), and the transition density T (r). Each isosurface encompasses 95-97% of the density in question.
and for LR-TDDFT in the NTO representation it is Thus the NTOs distill the content of the transition density into the smallest possible number of particle-hole pairs. In that well-defined sense, the NTOs are the best orbitals for visualization purposes, and detection of more than one significant singular value λ i indicates unresolvable multideterminant character in the excited state. For the excitation depicted in Fig. 6 there is little of that character, and T (r) ≈ ψ elec 1 (r) ψ hole 1 (r) is well described by the principle NTO pair. The nature of this product accounts for the somewhat more complicated nodal structure as compared to ∆ρ elec (r) ≈ |ψ elec 1 (r)| 2 or ∆ρ hole (r) ≈ |ψ hole 1 (r)| 2 .

Systemic Problems
The utility of LR-TDDFT lies in its combination of low cost, which facilitates calculations on systems such as C 119 H 154 ClN 21 O 40 [Fig 1(b)] or conjugated polymers (Figs. 5 and 6), along with an accuracy of ∼ 0.3 eV for localized valence excitations. That level of accuracy requires a treatment of dynamical correlation effects, as seen from the CIS errors in Fig. 4 that exceed 0.8 eV, comparable to the ∼ 1 eV of correlation energy for a pair of electrons. That is the good news. In this section we discuss some of the bad news, namely, systematic errors that make certain types of problems extremely challenging for LR-TDDFT. Of these, the most widely discussed is severe underestimation of excitation energies for states with substantial CT character, ultimately manifesting in spurious CT states in a sufficiently large system (Section 2.3.1). The second is a problem with the topology of conical intersections that involve the ground state, which presents a problem for direct ab initio simulation of internal conversion following photoexcitation (Section 2.3.2).

Description of Charge Transfer
Problems with the description of long-range CT excitations manifests in small, gas-phase molecules as Rydberg excitation energies that are systematically too low, 247 even when reasonable accuracy is obtained for valence excitations. This was noticed in early studies of LR-TDDFT and was quickly diagnosed as a symptom of incorrect asymptotic decay of the XC potential in GGA functionals that existed up to that point. [268][269][270][271] The same problem was quickly recognized to affect CT excitation energies. 84,115 Both CT and Rydberg excitations are sensitive to the long-range behavior of the potential, which ought to be v σ xc (r) ∼ −r −1 for a charge-neutral molecule. [272][273][274][275][276] This asymptotic behavior ought to be borne by the exchange potential because correlation dies off more quickly, 273,277 but in practice so does semilocal exchange.
Consider the form of the LR-TDDFT pseudo-eigenvalue problem for an excitation between MOs ψ iσ and ψ aσ that are well-separated in space, such that ψ iσ (r) ψ aσ (r) ≈ 0 everywhere. A semilocal expression for v σ xc (r) affords a semilocal XC kernel, such that the matrix elements (ia|κ στ xc |jb) in A and B vanish in such a situation, for all j and b. Ignoring spin by setting σ = τ in Eq. (19), this leaves Only the integral (ij|ab), which comes from the HFX term, survives to provide distance dependence for the i → a excitation. A pictorial illustration is provided in Fig. 7, which plots the distance dependence of the lowest CT excitation energy (ω CT ) in the (C 2 H 4 ) · · · (C 2 F 4 ) dimer as a function of intermolecular separation. 203 Only Hartree-Fock theory affords the correct distance dependence for ω CT (R), which varies according to the Mulliken formula, 66,278,279 ω CT (R) = IE donor + EA acceptor − 1 R in atomic units. For hybrid functionals the last term becomes −α hfx /R rather than −1/R, leading to a too-small penalty for long-range CT. For semilocal functionals where α hfx = 0, the CT excitation energy has no distance dependence whatsoever, once the donor and acceptor moieties are sufficiently far apart such that their orbitals do not overlap. This is reflected in the flat ω CT (R) profile for the LDA functional in Fig. 7. As a result, long-range CT excitation energies are almost invariably too small in LR-TDDFT unless the functional contains 100% HFX, which it usually does not because fully nonlocal exchange is somewhat unbalanced given the local nature of existing correlation functionals. The M06-HF functional is an example that does use 100% HFX, leading to reasonable performance for Rydberg states but larger errors for valence excitations [ Fig. 3(d)].
Where small-molecule benchmarks are available, errors in CT excitation energies can exceed 3 eV, 228 but the problem gets worse in larger molecules so that error is likely limited only by the size of the benchmark systems for which reliable ab initio results are available. A consequence of this severe underestimation of CT excitation energies is the appearance of completely spurious CT excited states in large systems, especially  solvated chromophores 280-287 but also large molecules. 159,288,289 When the system size is sufficiently large, there are inevitably well-separated occupied and virtual MOs such that the orbital energy gap ε a − ε i is small. For α hfx ≈ 0, the electron-hole interaction vanishes and the diagonally-dominant A matrix consists of weakly-coupled blocks corresponding to these spurious CT transitions. The kernel f σσ xc (r, r ) lacks the long-range exchange (or a derivative discontinuity, 63,68,279 or frequency dependence 290 ) that is needed to provide an energetic penalty for CT and an upshift to ω CT ≈ ε a − ε i in Eq. (57).
A physical example is shown in Fig. 8 for a model of aqueously-solvated uracil. 286 Whereas this system ought to have only a 1 nπ * and a 1 ππ * state below 6 eV, 227 a hybrid LR-TDDFT calculation using the PBE0 functional results in numerous low-energy solvent-to-chromophore CT states, e.g., 27 states below 6 eV for the (uracil)(H 2 O) 25 cluster that is shown in Fig. 8 and additional states as the size of the water cluster grows. 227,286 Many of these states are accidentally near-degenerate with the optically-bright 1 ππ * state and as a result these nominally dark CT states can acquire intensity from the bright state, which diminishes the intensity of the latter because total oscillator strength is conserved by the Thomas-Reiche-Kuhn sum rule, Eq. (38). The state ω 9 in Fig. 8 exhibits the largest coefficient with ππ * character, 286 yet due to spurious intensity borrowing it does not exhibit the largest oscillator strength and itself contains some contribution from solvent-to-chromophore CT. Fortunately, the same sum rule can be used to argue that the overall spectral envelope may still be valid upon ensemble averaging and broadening, even if some fraction of the oscillator strength has been ported onto spurious CT excitations.
In large chromophores, such as conjugated polymers, spurious low-energy CT excitations can manifest as artificial delocalization of the excitation across the length of the chromophore, 159,291,292 whereas the CIS method predicts that exciton size eventually saturates even as conjugation length increases. 292 As such, there is a need to develop a metric for whether a particular excited state has too much CT character for its excitation energy to be trusted. The first such CT metric to see widespread use was the quantity Λ defined by 293 where measures the overlap of |ψ iσ (r)| and |ψ aσ (r)|. (Absolute values are required since the occupied and virtual MOs are orthogonal.) This overlap is then weighted by the LR-TDDFT amplitudes and normalized such that 0 ≤ Λ ≤ 1. For calculations that do not invoke the TDA, however, the denominator in Eq. (59) is an odd choice, given the normalization condition in Eq. (21), and this inconsistency has propagated into other CT metrics used in LR-TDDFT. 294,295 Regarding the metric in Eq. (59), an early benchmark study concluded that 0.45 ≤ Λ ≤ 0.89 for localized valence excitations, making values in this range "safe" for LR-TDDFT, whereas 0.08 ≤ Λ ≤ 0.27 for Rydberg excitations, which are unsafe. 293 It was suggested that excitation energies for which Λ 0.3-0.4 (depending on the functional) should not be trusted. Various LR-TDDFT errors have been rationalized by appeal to Λ or similar metrics. 135,243,296,297 The point at which CT character becomes a problem is dependent on the manner in which it is quantified, 298 and several alternative CT metrics have been suggested. [293][294][295][299][300][301][302][303][304][305] Ciofini and co-workers introduced a widely-used "D CT metric", 299 originally defined in a rather complicated way but which ultimately measures the distance between the centroids of ∆ρ elec (r) and ∆ρ hole (r). The centroid of ∆ρ elec (r) is with an analogous definition for ∆ρ hole (r). If one defines then the distance between centroids of the electron and the hole is d − elec/hole , whereas d + elec/hole is the average position of the center of mass of the exciton. The quantity d − elec/hole is equivalent to the D CT metric but is more directly connected to the physics of the excitation. Other similar descriptors can be envisaged. 306,307 For example, by defining the root-mean-square size of the electron and the hole, one may define a charge-displacement distance, The quantity d CD connects directly to the properties of the exciton and is a more physically-motivated version of the "electron displacement" metric introduced by Adamo and co-workers, 295 and one that avoids the incorrect normalization in Eq. (59) and is thus rigorously invariant to orbital rotations even when the TDA is not invoked. To the best of our knowledge, d CD is introduced here for the first time but we suggest that d − elec/hole and d CD should replace alternative CT metrics that serve essentially the same purpose. To combat the long-range CT problem without going beyond the adiabatic approximation, LRC functionals are used to provide an XC potential with correct asymptotic behavior for an electron-hole pair. The LRC modification was introduced in Eq. (48) and contains an additional parameter that controls the separation between semilocal GGA exchange at short range and nonlocal HFX at long range. Figure 9 shows an example of how these functionals can be used to mitigate the growth in spurious CT states around a chromophore in aqueous solution. 287 Whereas the number of CT states increases extremely rapidly as water molecules are added around the system, and hybrid functionals such as B3LYP only partially mitigate this increase, the functionals LRC-ωPBE 228 and LRC-ωPBEh 255 control this growth completely. The LRC-ωPBEh functional is a short-range hybrid with α hfx = 0.2 whereas LRC-ωPBE is semilocal at short range (α hfx = 0), but both functionals employ 100% HFX in the long-range limit. This should be contrasted with functionals such as CAM-B3LYP, 225 which use range separation but sacrifice proper asymptotic behavior in an effort to obtain more accurate excitation energies for localized valence transitions. Although RSH functionals such as CAM-B3LYP and ωB97X-D are good choices in many respects for valence excitations, neither improves the accuracy of LR-TDDFT for CT excitations. 308 Standard double-hybrid functionals contain only a fraction of HFX and thus do not improve the situation for CT states, 309 unless the LRC scheme employed. 310 In the early development of LRC functionals, the range-separation parameter was often fit to minimize error in some benchmark thermochemical or excitation energy data. 200,202,226,228,254 However, excitation energies were found to be quite sensitive to this parameter, 226,227,255 especially for states where CT character is involved. 227,255 More recently, the community has increasingly turned to a more theoretically wellgrounded "optimal tuning" strategy, 231,311-314 which is grounded in the ionization energy (IE) theorem of exact DFT. 315,316 That theorem simply states that IE = −ε HOMO for the exact Kohn-Sham functional, consistent with the fact that the IE is set by the asymptotic decay of the wave function. 256 This condition is violated badly by common GGA and even hybrid functionals, often by several electron volts. 317,318 The optimal tuning (or "IE tuning") procedure consists in enforcing this condition for an approximate XC functional, by adjusting the range-separation parameter µ such that Here, IE(N, µ) represents the ∆SCF value of the IE for the N -electron molecule, computed using a LRC functional with range-separation parameter µ. Alternatively, one might try to find the value of µ that comes closest to satisfying Eq. (65) for both the N -electron molecule and its (N + 1)-electron anion, representing donor and acceptor for electron transfer. That procedure has been shown to reproduce not only CT excitation energies but also to afford Kohn-Sham gaps (ε LUMO − ε HOMO ) in good agreement with fundamental gaps (IE−EA). 197 The optimally-tuned value of µ does exhibit a strong dependence on system size, however. [319][320][321][322] Strategies to mitigate this dependence have been suggested.   55 60 65 θ (deg) 55 60 65 θ (deg) 55 60 65 θ (deg) 55 60 65 (c) SF-CIS (d) SF-TDDFT (BH&HLYP)

Conical Intersections
A different systemic problem with LR-TDDFT, which is relevant in the context of computational photochemistry, is that it predicts the wrong topology around any conical intersection that involves the ground state. 146,325 The TDHF method suffers from the same deficiency, which is not a DFT artifact per se but rather a LR artifact, arising from an unbalanced description of the ground (reference) state and the excited (response) states. 27 The result is that the branching space around a conical seam that involves the two lowest electronic states is necessarily one-dimensional rather than two-dimensional.
(For examples, see Refs. 146 or 325.) Even the CIS method can exhibit erratic behavior when the ground state becomes quasi-degenerate with the first excited state, 26,147 because in the absence of double excitations the ground-and excited-state eigenvalue problems are decoupled (according to Brillouin's theorem), 59 leading to an unbalanced description. 27,325 This is not a problem for conical intersections between two excited states because those states are coupled by the matrix A, in both CIS and LR-TDDFT. An example of a conical intersection involving the ground state is Jahn-Teller symmetry lowering from D 3h to C 2v , which is illustrated for the H 3 radical in Fig. 10. 147 In the vicinity of the D 3h conical intersection, the upper-state potential surface exhibits erratic behavior at both CIS and LR-TDDFT levels of theory. This warping of the potential surface around a conical intersection has consequences in nonadiabatic molecular dynamics simulations, including SCF convergence difficulties 326 and incorrect internal conversion timescales. 26 As a result, nonadiabatic trajectory surface-hopping calculations based on LR-TDDFT should probably be terminated prior to internal conversion to the ground state. 26 The "spin-flip" (SF) variant of LR-TDDFT 327 has been suggested as a way to overcome this problem, as discussed in detail in Ref. 27. Briefly, SF-TDDFT uses a sacrificial reference state that is not the ground state of interest, but rather a state with higher spin multiplicity S + 1, for target states with total spin S. By combining single excitations with a single α → β spin flip, SF-TDDFT generates both ground and excited states of multiplicity 2S + 1 as spin-flipping excitations, meaning that both are obtained as solutions to a common eigenvalue problem. This eliminates the imbalance and restores correct topology to conical intersections involving the ground state, as seen for H 3 in Figs. 10(c) and 10(d). Functionals with ≈ 50% HFX perform well in the context of SF-TDDFT, 27,328 and the Becke "half-and-half" functional BH&HLYP (with α hfx = 0.5) has become the de facto standard for SF-TDDFT. 27 An unfortunate side effect of SF-TDDFT is that it tends to exacerbate spin contamination, 27,329 especially as one moves away from the Franck-Condon point on the potential surface and starts to enter regions of photochemical interest. This necessitates the use of state-tracking algorithms to maintain a consistent spin multiplicity. 329,330 There have been various attempts to find a more theoretically appealing solution to this conundrum by adding additional determinants to the excitation space in order to restoreŜ 2 symmetry. 27 Methods developed along these lines include a fully spin-complete version of SF-TDDFT, 329 which adds the minimal number of additional determinants needed to obtainŜ 2 eigenstates (based on the equation-ofmotion formalism), 77 and also a "mixed-reference" spin-flip (MRSF) approach, which uses a combination of high-spin and low-spin S + 1 reference states to generate target states with spin S. [331][332][333][334][335][336] Although the MRSF-TDDFT excitation manifold is not formally spin-complete, in practice the spin contamination is very small. 331 The analytic gradient 334 and nonadiabatic derivative couplings 335 for MRSF-TDDFT have recently been formulated, facilitating nonadiabatic molecular dynamics simulations.

Excited-State Kohn-Sham Theory: The ∆SCF Approach
For periodic DFT calculations, LR-TDDFT is theoretically ill-posed if semilocal functionals are used within the adiabatic approximation. [20][21][22] Specifically, the too-rapid asymptotic decay of v σ xc (r) causes the lowest LR-TDDFT excitation energy to collapse to the Kohn-Sham gap, ω = ε LUMO − ε HOMO . 20,21 Semilocal LR-TDDFT also does not produce bound excitons in periodic systems, 22 and in large (but finite) conjugated polymers, the exciton delocalization length typically extends to the length of the entire molecule. 292,337 This observation can be conceptualized as incomplete cancellation of self-interaction that grows worse with system size, and infinitely worse under periodic boundary conditions. 20 Equivalently, it is a manifestation of the systematic underestimation of CT energies that was discussed in Section 2.3.1.
In recognition of these and other systemic problems exhibited by LR-TDDFT, there has been growing interest in "∆SCF" approaches that attempt to determine excited-state solutions to the Kohn-Sham SCF equation. 338,339 Having found such a solution, the excitation energy is computed simply as the difference relative to the ground-state energy, hence "∆SCF". In contrast to the well-automated machinery of LR-TDDFT, these methods are less "black-box", involving more effort and finesse on the part of the user, because each excited state requires a separate calculation. On the other hand, the ∆SCF approach can exploit ground-state gradient technology for geometry optimizations and vibrational frequency calculations. 340 For this reason, the ∆SCF procedure has sometimes been called excited-state Kohn-Sham theory. 338,340 In cases where LR-TDDFT exhibits known deficiencies, the ∆SCF approach may be more accurate and more reliable even if the formal justification (based on the Hohenberg-Kohn theorems 40,58 ) is absent because the system is not in its ground state. The method therefore rests upon the assumption that the description of short-range dynamical correlation depends upon the local environment of an electron and can be ported to a "non-aufbau" solution of the SCF equations, in which an electron has been promoted into a virtual MO. Such a state does not formally satisfy the non-interacting v-representability requirement of ground-state DFT. 17,58,68,341 Excited-state SCF solutions do contain full orbital relaxation, yet these solutions are inherently unstable because they are saddle points rather than local minima in the space of orbital rotations. Attempts to locate these non-aufbau solutions, each characterized by a virtual (empty) level that is lower in energy than the HOMO level, may suffer "variational collapse" to the ground state or to a lower-lying SCF solution. It is up to the user to determine that the SCF solution corresponds to the state of interest; if not, then the search must be begun anew, using a different SCF convergence algorithm or a different initial guess. Several modified SCF algorithms have been developed to try to locate non-aufbau solutions, based on overlaps with a set of user-specified MOs [342][343][344][345] or else based on direct search. [346][347][348] These algorithms are described in

General Considerations
A flowchart of the SCF procedure is illustrated in Fig. 11. At each iteration, the occupied MOs {ψ iσ } are used to construct the Fock matrix F σ , which is then diagonalized to obtain new MOs. Notably, diagonalization results in n basis = n occ + n vir MOs and one must decide how to choose the occupied set. Ordinarily, the lowest eigenvalues ε rσ are selected (aufbau principle), resulting in the ground-state determinant upon SCF convergence. To locate an excited-state SCF solution instead, one seeds the procedure with initial-guess MOs from a ground-state calculation but with non-aufbau occupancies, promoting an electron from HOMO to LUMO, for example. This makes the LUMO into an occupied level and the HOMO into a virtual level, resulting in a "hole", i.e., a virtual level whose energy lies below that of the highest occupied level. When the Fock matrix is constructed from this new set of occupied MOs and then diagonalized, the question becomes which of the new MOs should be the occupied ones, since energy levels may have shifted. The SCF procedure therefore deviates from the usual one only when it comes to selecting the occupied subset from among the n basis MOs. Several different options have been explored, as discussed in the next section. Before reviewing algorithms for locating non-aufbau SCF solutions, however, it is important to note some properties of those solutions that are different from ground-state properties. First, because the effective HamiltonianF σ [{ψ iσ }] depends on the MOs themselves, the ground-and excited-state Slater determinants are eigenfunctions of different Hamiltonians and are therefore not orthogonal. One consequence is that the formula for oscillator strengths in terms of transition dipole matrix elements [Eq. (3)] is not strictly valid, 349 as that formula is derived using the assumption that the eigenfunctions of the Hamiltonian form a complete orthonormal set. 44 In small-molecule tests, however, overlap integrals between ground-and excited-state determinants are found to be 0.1. 342 Another general concern is that excited states are always open-shell species, even if the ground state is closed-shell, so any single-determinant approximation is certain to be spin-contaminated, perhaps badly so. Indeed, single-determinant approximations for open-shell singlet states are often characterized by Ŝ 2 ≈ 1 (in atomic units of 2 ), which is equal to the average of pure-state singlet and triplet values. A similar phenomenon occurs, for similar reasons, in the case of the spin-unrestricted Hartree-Fock wave function in the separated-atom limit, 59 because a spin-pure state with two half-filled orbitals can be described using a minimum of two Slater determinants. The same two determinants (with different relative signs) are needed to describe both singlet spin-coupling (total S = 0) as well as the M S = 0 component of triplet spin-coupling (S = 1).
In practice, ∆SCF excitation energies for open-shell singlet excited states are often surprisingly accurate despite significant spin contamination, 342,344,350 although there are exceptions. One such exception is the 1 B 1u state of ethylene, whose underestimation by almost 2 eV is attributed to severe spin contamination. 342 Yamaguchi and co-workers have developed spin-projection techniques that can be used to recover spin-pure states in such cases, [351][352][353][354] and approximate spin purification is often used as a practical workaround in "broken-symmetry" DFT calculations of transition metal complexes. [353][354][355] For an open-shell singlet, the most common approach is to approximate the singlet energy as Here, E mix is the energy of the contaminated (mixed-spin or broken-symmetry) state that is obtained in searching for a singlet solution, whereas E triplet is the triplet energy for the same system, for which spin contamination is typically less severe. This procedure has a long history, 355,356 and Eq. (66) can be viewed as an approximate form of spin projection, generalizable to cases where the target state has spin S > 0. 352 The formula in Eq. (66)

Orbital-Optimized Non-Aufbau SCF Solutions
The simplest means to construct a non-aufbau occupied set is known as the maximum overlap method (MOM). [342][343][344][345] Starting from an initial guess corresponding to non-aufbau occupation of the ground-state MOs, this approach uses an overlap criterion to identify the new MOs at each subsequent SCF iteration. To do this, one must compute the projections p rσ of the MOs ψ (n) rσ at the nth iteration onto a reference set of MOs. The reference set might be the MOs at the previous iteration, in which case or else it could be the initial set of ground-state MOs, {ψ iσ }: The first choice [Eq. (67)] represents the original version of the algorithm, 342 whereas Eq. (68) has been called the "initial MOM" (IMOM) algorithm and tends to have better success at converging orbital-relaxed non-aufbau states. 344 The signature of success is a "hole below the Fermi level", i.e., a virtual MO whose energy is lower than the HOMO energy. The MOM and IMOM algorithms consist simply in replacing the aufbau selection of occupied MOs with a selection based on the n occ largest values of the overlaps p rσ . All other aspects of the SCF algorithm remain the same. This approach exhibits the same cost per SCF iteration as the ground-state algorithm and when it succeeds, the rate of convergence (measured in number of SCF cycles) is typically on par with a conventional ground-state calculation. There are certainly cases where MOM and IMOM fail, 347,348 however, typically resulting in variational collapse to the ground-state SCF solution. In such cases, more robust SCF convergence algorithms are required.
One such approach is the "σ-SCF" method, 346 which is based on minimizing the functional for a specified energy ω. This idea stems from recognizing that eigenstatesF |Ψ = ω|Ψ satisfy the zerovariance condition F 2 = F 2 . The σ-SCF approach avoids variational collapse by solving a proper minimization problem, but the appearance ofF 2 means that four-particle operators are required and the requisite transformations endow this method with O(n 5 basis ) scaling. 346 This makes the σ-SCF approach much more expensive than conventional SCF theory.
An alternative approach with the same formal scaling as the ground-state SCF problem is squared-gradient minimization (SGM). 347 Here, the idea is to convert an inherently unstable saddle-point optimization into a search for a local minimum by optimizing an objective function equal to the squared gradient of the energy with respect to orbital rotations. A local minimum can always be converged (if slowly), whereas a saddle point can be missed, and this makes SGM more robust as compared to MOM or IMOM. While the cost remains O(n 3 basis ), it is 2-3 times more expensive per SCF iteration as compared to a conventional SCF calculation, due to the cost of constructing the objective function. 347 It is also known that the squared gradient ∇ V (x) 2 of a function V (x) may contain minima that do not correspond to stationary points of the original function. [367][368][369][370] From the standpoint of trying to locate an orbital-relaxed excited-state Slater determinant, these are spurious solutions.
A middle way between MOM and SGM is state-targeted energy projection (STEP), 348 which does not increase the cost per SCF iteration yet shows much more robust convergence behavior as compared to MOM or IMOM. The STEP approach constructs a projection operator onto the virtual space, where the summation runs over some or all of the virtual MOs. The matrix representation ofQ σ is Q σ = C σ C † σ , where C σ consists of column vectors corresponding to whichever MOs are included in Eq. (70). The Fock matrix is then modified according to where S is the atomic orbital (AO) overlap matrix. The effect of the additional term is to shift all of the orbitals that are included in Eq. (70) by an energy η. By pre-selecting a virtual MO from the ground-state calculation that will be occupied in the first iteration of STEP, one can modify the Fock matrix to shift other virtual orbitals (including a lower-energy one that was occupied in the ground state but whose electron was promoted) to energies above the non-aufbau orbital. For example, upon HOMO → LUMO promotion, the original HOMO is unoccupied and should be included in Eq. (70), whereas the LUMO becomes occupied and should be excluded fromQ σ . The STEP algorithm is a form of level-shifting that tends to ensure that the SCF algorithm converges to the "closest" stationary point in the space of MO coefficients, which therefore resembles the initial guess. 348 Like MOM and IMOM, STEP can be used in conjunction with ground-state gradient technology to perform geometry optimizations and vibrational frequency calculations.

Transition Potential Methods
The methods described in Section 3.1.2 involve state-specific orbital optimization, meaning that the SCF procedure must be iterated to convergence separately for each excited state of interest. This has the advantage of including full orbital relaxation effects (beyond LR theory), but the disadvantage that there is no guarantee that an excited state resembling the one of interest can actually be found. A simpler (if cruder) approach was devised long ago by Slater,371,372 and forms the basis of several popular techniques for estimating x-ray excitation energies from Kohn-Sham eigenvalues. [373][374][375][376][377] To understand Slater's method, imagine that E({n i }) is the energy of a single-determinant wave function with orbital occupation numbers {n i }, some of which might be fractional. Expanding the energy as a Taylor series around a reference energy E 0 = E({n 0 i }), keeping the orbitals fixed, one obtains According to the Slater-Janak theorem, 378 the first derivative is an orbital eigenvalue: ε i = ∂E/∂n i . Now consider promotion of one electron from an occupied MO to a virtual MO. It suffices to deal with just a pair of occupancies (n i , n a ), in terms of which the transition in question can be abbreviated as (1, 0) → (0, 1). If a fractional-occupancy state with n i = 1/2 = n a is used for the reference {n 0 i }, then using Eq. (72) to compute the excitation energy ∆E = E(0, 1) − E(1, 0) leads to an estimate where the approximation neglects terms of order (n j − n 0 j ) 3 . 374 This forms the basis of the Slater transition method (STM), wherein an SCF calculation is carried out for the fractional-occupancy state (n i = 1/2, n a = 1/2) and then the energy difference ∆E = ε a − ε i affords an estimate of the excitation energy.
Variants of STM have historically been popular for x-ray spectroscopy, [373][374][375][376][377] particularly in the context of periodic DFT calculations where LR-TDDFT with semilocal functionals is problematic. 21 In principle this method requires a separate SCF calculation for each excitation of interest, and while it is generally easy to converge the x-ray "edge" in this way (i.e., a core → LUMO transition), higher-lying states will require a convergence algorithm that can avoid variational collapse. Moreover, this state-by-state approach leads to nonorthogonal MOs and therefore exhibits the same ambiguities regarding oscillator strengths as the ∆SCF method. 349 For these reasons, it is common to omit the 1/2 electron in the virtual space (with only pragmatic justification), leaving n i = 1/2 in the core-excited MO. This variant of the procedure has been called the transition potential (TP) method. 377,[379][380][381][382] By neglecting to occupy any core-excited virtual states at all, this approach sidesteps the issue of nonorthogonality, at least for a given choice of n i . Oscillator strengths can be computed in a straightforward way from matrix elements | Ψ 0 |μ|Ψ a i | 2 constructed from the orbitals obtained from the fractional-occupancy SCF calculation.
Modifications to the formula in Eq. (73) have also been proposed, 377,383-385 sometimes involving more than one SCF calculation with differing fractional occupancies, or by combining eigenvalues of both the neutral molecule and its cation or anion. 386,387 These modifications represent attempts to eliminate higherorder errors in Eq. (72). An example is the "generalized STM" (gSTM) method, 377,383 which replaces Eq. (73) with This is based on an alternative approximation for the integral The original STM in Eq. (73) corresponds to a midpoint approximation for this integral. 377,383 The gSTM approach requires two separate SCF calculations, one with (n i = 1, n a = 0) and the other with (n i = 1/3, n a = 2/3). Variants that set n i = 0 (removing the entirety of the core electron) have also been suggested and are sometimes called "full core hole" (FCH) methods. [373][374][375][376][377] The TP approach is then a "half core hole" (HCH) method. Although the FCH approach deviates significantly from Slater's original idea, it can be conceptualized as an attempt to restore charge balance, once the 1/2 electron in the virtual space has been abandoned for reasons of convenience. The excited core hole (XCH) approach 388 is yet another variant that creates a charge-neutral state (which is important for periodic DFT calculations) by placing the excited electron in the LUMO and using the full virtual spectrum from that calculation: 377,388 ∆E XCH = ε a (n i = 0, n LUMO = 1) − ε i (n i = 0, n LUMO = 1) . Together, these STM-and TP-type procedures are known as occupancy-constrained ∆SCF methods. In that context there has been some discussion of "many-electron" effects on oscillator strengths for x-ray transitions. 389,390 What "many-electron" means in this context is multi-determinant character in the final state, which is of course included automatically in a LR-TDDFT calculation.

Examples
The primary purpose of this chapter is to survey methods rather than applications but we will highlight a few recent applications of the ∆SCF approach in order demonstrate that it can be an elegant and lowcost alternative in cases where LR-TDDFT performs poorly, such as CT states. 338,391 Whereas LR-TDDFT systematically (and sometimes dramatically) underestimates CT excitation energies, the same excitation energies are systematically overestimated by the uncorrelated CIS method. 392 At the CIS level, a longrange excitation uses up the one occupied → virtual excitation that is included in the ansatz and leaves no excitations to facilitate orbital relaxation around either the electron or the hole, hence the overestimation. LR-TDDFT and CIS may therefore bracket the correct answer for a CT state but these upper and lower bounds can be several electron volts apart! 255 The ∆SCF approach includes full orbital relaxation and is also less sensitive to the asymptotic behavior of the XC potential. There has also been some preliminary work on the description of conical intersections and nonadiabatic dynamics using ∆SCF methods. 339,393 States with double-excitation character represent another categorical failure of LR-TDDFT within the adiabatic approximation, 99 with the most famous example being the optically-dark S 1 (2 1 A − g ) state in carotenoids, [395][396][397] or the analogous 2 1 A − g state in butadiene and other conjugated polyenes. [398][399][400][401][402] Doublyexcited states can be captured accurately using ∆SCF methods, 338,347,348 as shown for a few examples in Fig. 12. For these challenging cases, taken from a benchmark data set of double excitations, 394 several mGGA and hybrid functionals prove to be significantly more accurate than the CC3 method, which includes triple excitations and is generally close to CCSD(T) in quality, 403 with similar scaling. 404 For the full data set from Ref. 394, the hybrid GGA functional ωB97X-V achieves a mean absolute error (MAE) of 0.6 eV and a maximum error of 1.1 eV, whereas for CC3 the MAE is 1.0 eV and the maximum error is 1.8 eV. 348 The mGGA functional B97M-V does even better, with a MAE of 0.15 eV and a maximum error of 0.46 eV. 348 The ∆SCF methodology can also be used to compute an electronic absorption spectrum, although this must be done one state at a time by converging a sequence of non-aufbau determinants representing each excited state, and there is no guarantee that some states are not accidentally omitted. A successful example is shown in Fig. 13, reproducing the absorption spectrum of the chlorophyll a molecule that was only recently measured in the gas phase. 405,406 Using a STEP-based ∆SCF procedure, the major peaks in that spectrum  can be identified with transitions amongst the frontier MOs, 348 confirming the basic picture of Gouterman's four-orbital model. 407 LR-TDDFT calculations of the same molecule require twice as many states in order to resolve the spectrum up to 300 nm. Many of these states have near-zero oscillator strengths, 405 suggesting possible contamination by spurious CT states. Core-valence excitation energies are fertile ground for ∆SCF techniques. These states appear at photon energies ω > 200 eV and therefore it is not feasible to reach them by iterative solution of an eigenvalue problem starting from the lowest excitation energies. The frozen-valence approximation is one way to reach these states in LR-TDDFT, which is very accurate for K-edge transitions 177 but may be questionable for Lor M-edge excitations. Fortunately, core-to-LUMO excitations are relatively easy to locate using MOM. 343 Table 1 shows some error statistics for a benchmark set of K-and L-edge transitions. 408 Except for the LDA functional, all of the errors are < 2 eV and several functionals achieve errors < 0.5 eV, in excitation energies that are hundreds of electron volts. Notably, some of these same functionals also afford accurate L-edge transition energies, if spin-orbit interactions are included in order to describe the splitting of the 2p subshell into 2p 1/2 and 2p 3/2 states. 164 This splitting can be quite large, e.g., ≈ 13 eV for Fe(II). 409,410 Errors of < 0.5 eV are also possible for heavier elements using ROKS with relativistic corrections. 411 This excellent performance is perhaps somewhat surprising due to the substantially different self-interaction errors in core versus valence orbitals. 412,413 The "optimal tuning" of LRC functionals that was described in Section 2.3.1, in which the range-separation parameter is adjusted to set ε HOMO = −IE, can be understood as an attempt to cancel the self-interaction error associated with the HOMO, but that is likely to leave residual self-interaction in the much more compact core orbitals. These errors are exposed in ∆SCF calculations of core-level electron binding energies (for x-ray photoelectron spectroscopy), where many functionals afford errors 10 eV for transition metals. 414 Even the SCAN functional, which performs well for core-excited states of second-row atoms (Table 1), affords errors of ∼ 1 eV for core-level binding energies. 415 In several cases, Hartree-Fock theory proves to be more accurate than standard functionals that include correlation, even upon accounting for relativistic corrections. 414,416 This is consistent with other results indicating that the restricted open-shell (RO-)CIS method is a reasonable level of theory for M-and L-edge spectra of solidstate transition metal oxides, despite its lack of correlation effects, provided that spin-orbit corrections are included. 417 Resolution of this apparent paradox remains an open question.
For second-row atoms, errors in both core-level IEs 343,416,418 and also core-level excitation energies 343,348,408 are comparatively small when using the ∆SCF approach, although even for these elements Hartree-Fock theory is competitive with DFT, 418 suggesting that orbital relaxation is much more important than correlation. These rather small errors should be contrasted with much larger ones encountered when LR-TDDFT is applied to the same states using the frozen-valence approximation. Table 2 shows results for two different molecules (HF and CH 4 ) using two different functionals (SCAN and ωB97X-V) that both perform well in ∆SCF benchmarks. In contrast to the sub-eV errors obtained using the ∆SCF approach, LR-TDDFT calculations exhibit errors in excess of 10 eV for the carbon K-edge transition and ∼ 20 eV for the fluorine K-edge transition. For K-edge transitions of Mn(II) at ω ∼ 6540 eV, LR-TDDFT errors of ≈ 32 eV are obtained using B3LYP, and errors using the GGA functional BP86 are ≈ 62 eV. 419 Errors are even larger for heavier elements. 420 Notably, the sign of the errors in Table 2 points to underestimation of the excitation energy, consistent with too-soft asymptotic decay of the potential for a transition with CT character from a very compact 1s orbital to a radially-diffuse LUMO. LRC functionals perform much better in this capacity. 421,422 That said, the precision of core excitation energies computed using LR-TDDFT is rather good even if the accuracy is not, 422 meaning that chemical shifts can be obtained even if absolute excitation energies must be shifted to match experiment. It is the magnitude of the required shifts that is somewhat unnerving. This inspired work on RSH functionals that partition r −1 12 into short-, middle-, and long-range components, 211-214 with the intention to use a larger fraction of exact exchange (α hfx ≈ 0.87) at a length scale µ −1 ≈ 0.24Å. 423 These "short-range corrected" (SRC) functionals work rather well for LR-TDDFT x-ray calculations, 213 although they are empirically parameterized specifically for that purpose and may not be good functionals for other applications such as ground-state thermochemistry or valence excitation energies. Furthermore, the SRC functionals are not LRC functionals in the sense of Eq. (48), because they do not go to a limit of 100% HFX as r 12 → ∞. Proper LRC functionals have been developed to describe core-level excitations. 211 As illustrated by the porphyrin spectrum in Fig. 13 and other examples, 338 accurate electronic absorption spectra can be computed via ∆SCF calculations that are carried out in a state-by-state manner. However, this does present something of a nuisance as compared to automatic generation of numerous excitation energies in a LR-TDDFT calculation. Some of the TP or occupancy-constrained ∆SCF approaches that were described in Section 3.1.3 bypass this annoyance by using orbital eigenvalues from one or two fractional-occupancy SCF calculations to obtain the entire spectrum of excitation energies. These methods are assessed (versus experiment) in Fig. 14, for a data set of K-edge transitions. 377 Following the notation of Ref. 377, these methods are characterized as either "explicit" or "implicit", with explicit methods involving occupancy constraints that are applied state-by-state, as suggested by the original Slater method, whereas the implicit methods place no electron at all into the virtual space, which makes for a very simple computational scheme. The "neutral implicit" methods follow the paradigm of the XCH approach, 388 placing an electron (or a fraction of an electron) into the LUMO only; it is often possible to optimize a core-to-LUMO SCF solution without specialized algorithms. Finally, the "ground state" results in Fig. 14 represent a control experiment in which the ground-state eigenvalue difference ε LUMO −ε 1s is used to approximate the K-edge transition energy. That method fares surprisingly well, or perhaps the other approaches should be said to fare surprisingly poorly. If the user desires to avoid state-by-state optimization of an occupancy-constrained determinant, then the XCH approach would seem to be the best option. The cost of the ∆SCF approaches described in Section 3 is typically no more than a few times the cost of ground-state DFT, depending somewhat on the algorithm that is used to converge the non-aufbau SCF solution. With the exception of the σ-SCF method, the algorithms described in that section each possess the same O(n 3 basis ) scaling as the ground-state calculation. However, these methods must be applied in a state-by-state manner, constructing a different initial guess for each state of interest, which limits their applicability to problems where only a small number of states is desired or required. In contrast, using LR-TDDFT it is possible to obtain a large number of states in an automated fashion, at least for mediumsized molecules. For large molecules, the cost of computing a very large number of excited states can become prohibitively expensive, especially in systems where the density of states is large. The appearance of spurious CT states in the spectrum exacerbates the cost, even if their affect on the overall spectral envelope is nil. 286,424 In situations such as these, where proliferation of states (whether real or spurious) makes the iterative diagonalization cost-prohibitive, the TDKS approach may be advantageous. Using this method, a broadband spectrum can be computed from the oscillating dipole moment function obtained from time-dependent electron dynamics.

Theory
The TDKS approach is also known as "real-time" (RT-)TDDFT, 51,52 to distinguish it from LR-TDDFT. Starting from the ground-state Kohn-Sham determinant, an external electric field E(r, t) is turned on at t = 0, either as an impulse or as a continuous wave, and the resulting perturbation creates a time-evolving superposition state whose Fourier components encode the excitation energies. This is analogous to propagation of a non-stationary wave packet according to the time-dependent Schrödinger equation. Similar to that situation, the time-dependent MOs ψ kσ (r, t) are complex-valued for t > 0. Unlike the many-electron time-dependent Schrödinger equation, the effective HamiltonianF in Kohn-Sham theory depends on its own time-evolving eigenfunctions, andF (t) does not commute withF (t ) for t = t . Therefore the time evolution operatorÛ for time-independent HamiltonianĤ must be generalized tô for TDKS calculations, whereT is a time-ordering operator. This leads to so-called Magnus expansion that generalizes the time-independent Baker-Campbell-Hausdorff expansion. [425][426][427] A variety of time-propagation algorithms have been developed for TDKS simulations, 52,427-431 the simplest of which is a "modified-midpoint" algorithm, 429 a type of explicit Euler integration scheme that requires only one Fock matrix construction per time step. Denoting the Fock and density matrices in the orthonormal MO basis as F and P, respectively, the Fock matrix F is first propagated forward in time by a half step, (∆t)/2. At each subsequent instant in time, t n = n∆t, the matrix representation of the propagator is constructed according to again in the MO basis. 427 The MO density matrix is then propagated from t n−1/2 to t n+1/2 according to (Spin indices are omitted as these equations are valid for either spin.) Construction of U n requires diagonalization of the Fock matrix, which is not a problem in Gaussian basis sets but is not feasible in plane-wave basis sets or on a grid. See Ref. 428 for a discussion of alternatives when F is too large to diagonalize. For Gaussian basis calculations, self-consistent propagators based on predictor-corrector algorithms have also been developed. 427,432 These may require more than one Fock build per time step but allow for the use of somewhat larger time steps as well as for automatic detection of time steps that are too long, which is not always obvious from the usual criterion of checking to make sure that fluctuations in the total energy are bounded. 427 For the modified-midpoint algorithm that is encapsulated by Eqs. (79) and (80), the cost of a single time step is comparable to the cost of a single SCF cycle of the ground-state calculation. The storage requirement is also modest, amounting to a few complex-valued matrices of dimension n basis × n basis . This should be contrasted with the storage requirement for iterative solution of the LR-TDDFT pseudo-eigenvalue problem, which is O(n roots n occ n vir ) with a prefactor that reflects the number of iterations and therefore the size of the iterative subspace. The modified-midpoint approach works well provided that the time step ∆t is sufficiently small; values ranging from 0.01-0.50 a.u. are typical, where 1 a.u. ≈ 2.42 × 10 −17 s = 24.2 attoseconds. (The time for one orbit in the Bohr model of the hydrogen atom is 2π times this value, which establishes a timescale for electronic motion.) The maximum acceptable value of ∆t is limited not only by stability of the time integration but also by the excitation energies that one desires to access, as discussed below.
Unlike LR-TDDFT, which operates by definition in the limit of a vanishingly weak external field, the TDKS approach is non-perturbative and in principle can be used to simulate electron dynamics in strong laser fields, e.g., to simulate nonlinear optical properties of materials, 433 or to make contact with emerging attosecond spectroscopies 57 that create electronic wave packets that are out of equilibrium with the nuclei and thus outside of the Born-Oppenheimer approximation. [53][54][55][56] In practice, there are various issues related to the use of the adiabatic approximation (Section 2.1.2), 290,434,435 meaning the use of ground-state functionals with no memory, such that the time dependence is carried solely by the time-evolving density, E xc [ρ(r, t)]. On the other hand, within the adiabatic approximation the initial-state dependence vanishes since the XC kernel is fully specified in terms of the instantaneous time-evolving density. 436 The topic of strong-field electron dynamics and how it can be described using TDKS calculations is not considered here, except to note that there have been successful TDKS simulations of strong-field photoionization, [437][438][439][440][441][442][443][444][445][446][447][448] and also of high harmonic generation, [449][450][451][452][453][454][455][456][457] both in Gaussian-orbital representations of the density. Unlike the gridbased treatments that are common in atomic physics, Gaussian-based methods are scalable to molecules. However, self-interaction error is known to significantly suppress strong-field ionization rates, 458 therefore much of the aforementioned work has been performed at the TD-CIS level, where self-interaction is not a concern and the exchange potential has the correct asymptotic form.
This section focuses on the use of TDKS simulations to obtain broadband spectra. Within the electric dipole approximation (which is also invoked in LR-TDDFT insofar as oscillator strengths are proportional to transition dipole matrix elements), the absorption spectrum corresponds to the dipole strength function 163,459 S(ω) = 4πω 3c Im α xx (ω) + α yy (ω) + α zz (ω) (81) where for example is an element of the dynamic polarizability tensor, α(ω). The quantities µ x (ω) and E y (ω) are the Fourier transforms of the time-dependent dipole moment and the external electric field, respectively, although for an impulsive δ-function pulse the denominator in Eq. (82) can be replaced by the field amplitude while the numerator is replaced by Here, w(t) is a windowing or padding function, as in standard signal processing. 52 To obtain a linear absorption spectrum the external field E(r, t) should be weak, impulsive, and off-resonant, and should contain components in all three Cartesian directions in order to excite states of all symmetries, effectively averaging over molecular orientations. With an appropriately chosen integration time step, this procedure reproduces the same spectrum that is obtained using LR-TDDFT, if all of the LR-TDDFT excitations within the energy window of interest are included. 427,459 The TDKS approach affords the entire broadband spectrum from a single Fourier transform [Eq. (83)] following sufficient time propagation, but it is not straightforward to assign the features in the TDKS spectrum to transitions between specific MOs. Techniques to do so have been developed, based on identifying individual Fourier modes in the dipole moment matrix (expressed in the MO basis) at a specific transition frequency. [460][461][462][463][464] This does require some post-processing and some insight regarding the important MOs. Time-dependent generalization of the NTO basis (Section 2.2.4) have also been proposed. 465 The cost of the time steps needed to propagate the Kohn-Sham MOs in time, and thus to obtain the time-evolving density and dipole moment function, is comparable to the cost of a single ground-state SCF cycle but many time steps are required. A typical simulation time might be 30 fs to obtain a fully-converged spectrum, 427 but with ∆t = 0.1 a.u. = 2.4 as (1 as = 10 −18 s), this represents > 10 5 time steps. Recently, Padé approximant techniques have been introduced in order to obtain µ(ω) based on a short time series of input data µ(t). 52,462 Using this approach, spectra that are well-converged (with respect to LR-TDDFT results) can be obtained in < 10 fs of time propagation and rough spectra can be obtained with as little as 3-5 fs. 163 The time step ∆t dictates the spectral window that can be accessed, via the usual time-energy Figure 15: Experimental x-ray absorption spectra of α-quartz at the L2,3 edge, 466 along with LR-and RT-TDDFT calculations using a cluster model Si5O16H12 of the bulk material (as shown), with modified hydrogen charges to enforce charge neutrality. 467 The spectrum labeled "RT-TDDFT + iΓ" includes phenomenological lifetime parameters for the virtual orbitals. All calculations were performed with an optimally-tuned version of LRC-ωPBEh. The dashed line labeled E0 is the experimental ionization energy of Si(2p). Reproduced from Ref. 467; copyright 2015 American Chemical Society. uncertainty relationship (∆E)(∆t) that comes from the Fourier transform. In practice, the spectrum is only reliably converged up the Nyquist frequency f Ny = π/(∆t), and perhaps only a fraction of that value. 163,427 This implies that especially small time steps are required for x-ray applications. A time step ∆t = 0.1 a.u., for example, corresponds to a Nyquist frequency f Ny = 854 eV, which is well above the K edge for second-row elements (C, N, O, etc.) but not for third-row elements. The K edge for elements Al-Cl lies above 1,500 eV.

Examples
As an example of a broadband spectrum of interest in materials science, Fig. 15 presents x-ray spectra at the L 2,3 edge of α-quartz, computed using several different TDDFT methods. 467 Starting from the lowest valence excitations, a LR-TDDFT calculation with n roots = 300 reproduces the first two features in the experimental spectrum (labeled "A" and "B") but is unable to resolve the higher-energy features. (The calculations do not include spin-orbit coupling and thus do not reproduce the doublet for peak A, which arises due to the 0.6 eV splitting of the 2p 3/2 and 2p 1/2 levels. 466 ) Although the most intense feature (peak B) is evident in the RT-TDDFT spectrum, at a peak position that precisely matches the corresponding LR-TDDFT spectrum, the RT-TDDFT spectrum is quite noisy and other features in that spectrum are obscured by this noise.
That noise is actually a basis-set artifact arising from the absence of proper continuum states (or the inability to describe ionization within a finite-basis approximation), which has the effect of artificially trapping metastable excitations that lie above the ionization threshold. 467,468 A solution to this problem is to incorporate phenomenological lifetimes for the unbound MOs, meaning those with ε pσ > 0. This procedure is described in Ref. 468 and corresponds to a modification ε pσ → ε pσ + iΓ in the MO basis, where Γ −1 is a phenomenological lifetime that is modeled as a function of energy, decreasing exponentially above the vacuum level. As seen in Fig. 15, this "iΓ" modification removes the noise from the RT-TDDFT spectrum, such that most of the experimental features become evident even at energies far above what can feasibly be reached with LR-TDDFT. This does require some phenomenological modeling, however.
When computing K-edge x-ray spectra using high-quality basis sets, similar artifacts can manifest as spurious pre-edge features that are not seen in compact basis sets where the density of levels ε pσ is more sparse. In a molecule that contains both nitrogen and oxygen, for example, excitations from N(1s) core orbitals to the highest-energy virtual MOs can manifest as spurious pre-edge features at the oxygen Kedge, 163,469 as shown in Fig. 16. These artifacts appear despite the fact that the nitrogen K-edge lies > 100 eV below the oxygen K-edge! These intruder peaks could potentially be mitigated via heuristic lifetime models for the unbound states, as described above, although a simpler fix is to modify the timedependent dipole moment matrix in the MO basis. That matrix is D jkσ,x (t) = −e ψ jσ (t)|x|ψ kσ (t) (84) for the x component. The time-dependent dipole moment function that is needed in Eq. (83) is then µ x (t) = tr(P α D α,x + P β D β,x ) .
By eliminating the rows and columns of D σ,x that correspond to occupied MOs other than the ones of interest, prior to computing the Fourier transform in Eq. (83), the undesired resonances can be removed from the spectrum. 469 For the oxygen K edge, this means retaining only those rows and columns where j and k refer to the O(1s) orbitals. This is precisely analogous to the frozen-valence truncation of the LR-TDDFT excitation manifold that is used to obtain core-level spectra (Section 2.2.1), and may have similar limitations for L-and M-edge spectra. For the oxygen K edge, Fig. 16 shows that this procedure affords a spectrum in good agreement with LR-TDDFT, free of contamination by N(1s) excitations. This filtering procedure does require the user to decide in advance which edges are of interest, so that multiple edges can no longer be computed in a single calculation unless the entire trajectory of dipole moment matrices is stored.
The TDKS approach has been extended to compute excited-state absorption spectra, 436,470 using an excited-state density prepared via LR-TDDFT as the initial density at t = 0. Because the initial state is non-stationary, this requires that the "field on" simulation be referenced to a time-evolving "field-off" simulation. 436 This approach has recently been applied to simulate emerging transient x-ray experiments, [471][472][473][474] carried out at free-electron laser facilities using x-ray pulses with femtosecond time resolution. As an example, it is possible to follow metal-to-metal CT dynamics in the mixed-valence [(CN) 5 Fe II CNRu III (NH 3 ) 5 ] − compound, which occur on a ∼ 60 fs timescale following excitation at 800 nm, using time-resolved x-ray emission spectroscopy at the iron K edge (2p 3/2 → 1s transition at 7,114 eV). 474