Perspective : Fundamental Aspects of Time-Dependent Density Functional Theory

In the thirty-two years since the birth of the foundational theorems, time-dependent density functional theory has had a tremendous impact on calculations of electronic spectra and dynamics in chemistry, biology, solid-state physics, and materials science. Alongside the wide-ranging applications, there has been much progress in understanding fundamental aspects of the functionals and the theory itself. This Perspective looks back to some of these developments, reports on some recent progress and current challenges for functionals, and speculates on future directions to improve the accuracy of approximations used in this relatively young theory.

The wavefunction rose from the early days of quantum mechanics as the central player, the provider of all observable properties of atoms, molecules and solids.Within a few years it was however recognized that finding this complex-valued high-dimensional function was impossible for all but very few problems, even if we just focus on electrons alone.Almost 40 years later, the Hohenberg-Kohn theorem [1] proved that the groundstate density alone provides all observable properties of any static system.This is an astonishing result given the simplicity of the density, the probability of finding any one electron at a given point in space, compared with the wavefunction, a function of all electronic coordinates.Even more astounding, in 1984 Erich Runge and Hardy Gross proved that for time-dependent (TD) systems evolving from a given initial wavefunction, all TD properties can be extracted from the time-evolving density [2].Their theory, Time-Dependent Density Functional Theory (TDDFT) has become one of the most popular and successful methods for computing electronic excitation spectra today, with increasing applications for non-perturbative electron dynamics, and is the subject of this Perspective.I focus here mostly on fundamental aspects underlying the theory and its functionals, while applications enter throughout as a guide and inspiration.
The very first TD density functional calculation, to my knowledge, in fact predates not only the Runge-Gross theorem, but also the Hohenberg-Kohn theorem: Felix Bloch developed a time-dependent extension of Thomas-Fermi theory [3,4] to theoretically model the stopping power of a fast charged particle in matter.In the 1970's, the first time-dependent Kohn-Sham calculations were run [5][6][7], and in the early 1980's proofs of one-to-one density-potential mappings for certain cases began to appear, for time-periodic potentials [8] and adiabatic processes [9,10].The growing interest in running time-dependent Kohn-Sham calculations did not have to wait much longer before they could stand on a formal, and general, foundation.The Institute for Theoretical Physics at the Goethe University Frankfurt in the 1980's was especially known for nuclear physics and atomic scattering.Hardy Gross, a postdoc at the time, and stu-dent Erich Runge, were particularly motivated by the question of what is the time-dependent potential acting on an electron when it scatters from ions whose nuclei are treated classically.Given that Hohenberg, Kohn, and Sham had derived such a potential for an electron in a ground-state twenty years earlier, it was natural to wonder whether such a formulation could be extended to the time-dependent case, yielding the exact time-dependent potential acting on an electron.This led to the birth of the foundational theorem of TDDFT, the Runge-Gross theorem [2].
What Runge and Gross considered, were nonequilibrium systems of electrons treated nonrelativistically, described by the following universal Hamiltonian: We use atomic units (e 2 = = m e = 1) throughout.What distinguishes different N -electron systems is the "external" potential, v ext (r, t), that represents the potential the electrons experience due to the nuclear attraction and due to any field applied to the system.The time-dependence in this term can arise through nuclear motion (with nuclei treated as classical, but possibly dynamic, point particles) and/or through externally applied (scalar) fields.Of course, solving the TD Schr ödinger equation (TDSE), Ĥ(t)Ψ(t) = i∂ t Ψ(t) for the many-electron wavefunction Ψ(t) is not computationally feasible for more than a few electrons, as it runs into an exponential wall very quickly for nonperturbative fields.Moreover Ψ(t) provides usually far more information than we would be interested in.Typically, observables of interest involve integrated quantities, in particular one-and two-body probability densities.The idea of density-functional theories, broadly speaking, is to obtain such quantities directly without calculating Ψ(t).
Time-dependence opens up an ever-increasing wealth of interesting physical and chemical phenomena that tend to be classified into two groups: spectroscopy, which gives an electronic "fingerprint" of the atom, molecule, or solid (Sec.III), and real-time nonperturbative dynamics where the electronic system roams far from any ground-state (Secs.II and IV).In this Perspective, I cover only in passing a third group of phenomena arising from relaxation and dissipation in very large or extended systems, as induced by electronelectron viscosity.A fourth group of applications relevant to the ground-state correlation energy, which exploit the fluctuation-dissipation theorem to fold in longrange effects in a natural way [11], is not discussed.Further, I have brazenly ignored all-important computational concerns, numerical aspects, details, and algorithms.

I. DENSITY-POTENTIAL MAPPING
The Runge-Gross theorem proves the one-to-one mapping where Ψ(0) is the initial wavefunction and the density n(r, t) is the probability of finding any one electron of any spin σ at point r in space at time t: (3) As in the ground-state case, an elegant aspect of the original proof of this density-potential mapping, is that it can be completely followed from basic principles learned in the first year of a quantum mechanics course.Instead of the variational principle for energy minimization that leads to the Hohenberg-Kohn theorem, the TD case follows from Heisenberg equations of motion for the current-density, and the continuity equation.The basic idea of the proof is outlined here, and readers are referred to Refs.[2,12,13] for the details.
Consider two systems of N electrons, both initially in the same many-body state Ψ(0).Can two, physically distinct potentials, v ext (r, t) and v ext (r, t) be found such that the two systems evolve in a Hamiltonian of the form Eq. ( 1) with exactly identical densities n(r, t) = n (r, t) at all times?Physically distinct potentials are those that differ by more than just a purely TD function, i.e. v ext (r, t) and v ext (r, t) + α(t) are not physically distinct, because they lead to wavefunctions that differ only by a purely-TD phase factor, which cancels out when expectation values of Hermitian operators are taken in calculating observables.Now, if the answer to the above question can be shown to be "no", then it implies that the density n(r, t) can only be produced by one external potential, i.e. it proves the arrow pointing to the right in Eq. ( 2).(The arrow pointing to the left, that the potential points to a unique density follows from the usual assumption that the TDSE yields a unique Ψ(t), and hence via Eq.(3) a unique density).To prove the answer is in fact "no", Runge and Gross required that the considered potentials be time-analytic around the initial time, and used their Taylor expansion expressions in the equation of motion for time-derivatives of the current-density, evaluated about the initial time: Essentially, the difference ∂ l t j(r, t)| t=0 − ∂ l t j (r, t)| t=0 depends on the difference between the k<lth Taylor coefficients of the two potentials.If the two potentials are physically distinct, then there must exist a k for which their spatialdependences differ, from which the equation shows that the resulting currents will differ.In a second step, the equation of continuity, ∂ t n(r, t) = −∇ • j(r, t), is used to transform the one-to-one current-potential mapping to a one-to-one density-potential mapping.
The top part of Figure 1 illustrates this map.The initial states labelling the upper ellipses have the property that they have the same n(r, 0) and ∂ t n(r, 0) as each other, and the symbols inside the ellipses are examples of possible density-evolutions n(r, t) stemming from them.(The phase of the initial wavefunction determines ∂ t n(r, 0), via the equation of continuity).Different initial states can lead to the same density evolution in different potentials, so some density symbols appear in both.The Runge-Gross theorem states that no two lines from the same ellipse may lead to the same v ext in the central ellipse, which is the space of all one-body potentials, but lines from different densities in two different ellipses may point to either the same or different points in the central ellipse.The open symbols correspond to "non-v-representable" densities, meaning that these densities cannot be generated from evolving under any Hamiltonian of the form Eq. 1.We come back to the lower part of the figure shortly.

A. The TD Kohn-Sham system
The statement in Eq. ( 2) implies that all observables are functionals of the initial state Ψ(0) and the TD density n(r, t): O[n; Ψ(0)].This is because it states that, for a given initial state, knowledge of the density identifies the external potential, and therefore the entire Hamiltonian Eq. ( 1), and hence the wavefunction Ψ(t) (up to a physically irrelevant purely TD phase).But, as in the ground-state case, it is difficult to find these functionals, and so one resorts to mapping the system to a non-interacting system of orbitals, the TD Kohn-Sham (KS) system, which reproduces the exact density.Having orbitals eases the job of the functionals in many cases; in particular, the kinetic energy, approximated as the orbital functional 2 |φ i does a far better job than available explicit kinetic-energy densityfunctionals.The Runge-Gross argument holds for noninteracting electrons (in fact, for any particle-particle interaction), as illustrated in the lower part of Figure 1.Again, no two lines from the same lower ellipse can lead to the same one-body potential in the central ellipse, representing now v S (r, t), and again, lines from the same symbol appearing in both lower ellipses leading to two different potentials, reflects the fact that two different initial states can yield the same density evolution in different potentials.A density evolution in both an upper ellipse and a lower ellipse means that the density is "non-interacting v-representable" (more in Sec I B).
To run a TDDFT calculation for a given problem then, one first picks an initial state to be used for the noninteracting KS propagation, Φ(0), that is compatible with the true initial state of the physical problem Ψ(0) in that it has the same density and first time-derivative of the density.Typically a Slater determinant of singleparticle orbitals φ i (r, 0) is chosen.Even when the initial density is stationary, there are infinitely many Slater determinants to choose from that reproduce a given density [14,15].The KS orbitals evolve via a single-particle TDSE whose potential, v S [n; Φ(0)] is defined such that it guarantees the KS one-body density reproduces the exact density: The KS potential v S (r, where v H [n](r, t) = d 3 r n(r ,t) |r−r | is the Hartree potential, and v XC [n; Ψ 0 , Φ 0 ](r, t) is the exchange-correlation (xc) potential.In practical calculations, approximations enter first when selecting an initial KS state (the exact initial density is generally only known approximately, e.g.via a ground-state DFT calculation, or a correlated many-body calculation), and second in approximating the xc functional.
The initial-state dependence of the 1-1 mapping for the interacting and non-interacting systems endows the xc potential with a functional-dependence on the initial states Ψ(0) and Φ(0), as manifest in Eq. ( 5), a feature that has no precedent in ground-state DFT.In a large class of applications (e.g.calculations of response and spectra), the initial true and KS states are ground states, which, by virtue of the Hohenberg-Kohn theorem, are themselves functionals of their density; initial-state dependence is then redundant.But it is generally unknown what the functional dependence of the xc potential is on the initial states.The exact xc potential is known to depend on the past also through the history of the density, n(t < t), and this dependence together with the dependence on Ψ(0) and Φ(0) is called memory-dependence [16,17].In the vast majority of applications today, an adiabatic approximation is used, where the TD density is inserted into a ground-state approximation, audaciously neglecting all memory-dependence.
Note that the TD KS equations do not arise from a variational principle: all that was needed above is the reassurance that a given density evolution beginning from a given initial state can be produced by at most one potential, and the assumption of non-interacting v−representability (see more shortly).Still, thinking of the ground-state case, where many approximations are derived directly for the xc energy, whose functional derivative then gives the ground-state potential, it might be similarly useful if there is an analog for the TD case.TD variational principles are based on the action, but, if one attempted to write v XC [n, Ψ(0), Φ(0)](r, t) as δA XC [n, Ψ(0), Φ(0)]/δn(r, t) one encounters a predicament: δv XC [n, Ψ(0), Φ(0)](r, t)/δn(r , t ) would equal δ 2 A XC [n, Ψ(0), Φ(0)]/δn(r, t)δn(r , t ) which is symmetric in t and t , while in reality density-changes at times t > t cannot affect the xc potential at earlier times.This causality violation was first pointed out in Ref. [18] and several different possible solutions have been proposed: defining the action on the Keldysh contour [19], using Liouville space pathways [20], or including boundary terms with the usual real-time contour [21].

B. Solving an Existential Crisis Leads to Alternative Mapping Proofs
How do we even know that the time-evolving density of an interacting electronic system can be reproduced by a non-interacting system?The Runge-Gross proof ensures the uniqueness of the potential that yields a given density evolution from a given initial state but it does not prove its existence.The question of existence of a KS system for a given density is called "non-interacting v-representability".
In Ref. [22], Robert van Leeuwen showed that a unique TD KS potential exists for a prescribed TD density generated in an interacting system, under some conditions, and provided a method for its construction.His method also provides an alternate proof of the densitypotential mapping theorem of the statement (2).Like the Runge-Gross theorem, the van Leeuwen proof requires that the potential be time-analytic, but additionally the latter assumes the density is time-analytic.Now the time-analyticity of the potential covers most physical cases of interest (especially with the piece-wise analytic extension [16]), but time-analyticity of the density is a more restrictive condition.The coupling of space and time-derivatives in the TDSE means that spatial singularities in the potential can lead to non-analyticities in time in the wavefunction and density, even when the potential is time-analytic.Technically the proof in Ref. [22] unfortunately does not apply to Coulomb potentials.Limiting consideration to the linear response regime, alternate proofs of the density-potential mapping however have no time-analyticity requirements [23,24].
The search for a general proof free of the timeanalyticity requirements has seen a resurgence of interest in recent years.Ref. [25] reformulated the question of the density-potential mapping and v−representability in terms of uniqueness and existence, respectively, of a particular TD nonlinear Schr ödinger equation (NLSE).However the nature of the particular NLSE involves fourth-order derivatives of the wavefunction and because of this together with the singularity in the Coulomb potential, existence and uniqueness results are not known.Refs.[26][27][28] proposed a proof based on global fixed-point iteration, which further yields an algorithm to numerically construct the exact KS potential for a prescribed density-evolution from a given initial state.The convergence is technically not assured for Coulombic systems, but demonstrations on systems with smooth potentials have already led to interesting insights on features that exact xc potentials have [29][30][31].
Ref. [32] provides a detailed review of the mathematical aspects of the density-potential mapping, delving also into uniqueness and existence of solutions of the TDSE when Coulomb interactions are present.
Discretizing space, conditions for existence and uniqueness of the density-potential mapping and KS system in "lattice-TDDFT" have been proven in Refs.[33,34].For practical computations on a real-space grid, one might be tempted argue that the lattice-TDDFT is what is relevant, however without well-posedness of the continuum limit, there is no guaranteed convergence with respect to grid size.

C. TD Current Density Functional Theory
The first part of the RG theorem actually proves that no two scalar potentials yield the same one-body current, for a fixed initial state, implying a one-to-one mapping between the external potential and the current.However, currents of interacting systems are generally not non-interacting v-representable, as shown in Ref. [35]: even if no external magnetic field is applied to the true system, one needs a magnetic field in a noninteracting system for it to reproduce the current.That is, the KS current in TDDFT is, generally, not equal to the true physical curent.The longitudinal component of the TDDFT KS current is equal to that of the physical current, thanks to the continuity equation, but they may differ in their transverse component.
But even stepping back to the very first stage, the question of whether one can find an initial Slater determinant that reproduces a given initial density and current-density of an interacting system is highly nontrivial, unless the velocity j(r, 0)/n(r, 0) is curl-free; recent results can be found in Ref. [36].
A more useful one-to-one map exists between currentdensities and vector potentials [18,[37][38][39][40], yielding what is known today by TD current-density functional theory (TDCDFT).Here non-interacting Arepresentability was proven under conditions of timeanalyticity [41], similar to the van Leeuwen proof of v−representability in the TDDFT case.Closely related is TD deformation functional theory [42], which can be viewed as a Lagrangian version of TDCDFT, with the deformation metric tensor as basic variable.In fact, Ilya Tokatly has illuminated the relationships between the NLSEs of TDDFT, TDCDFT, and that underlying TD deformation functional theory [43].
Even when the external potential is purely scalar, there can be advantages to use TDCDFT instead of TDDFT when approximations are involved.Going beyond the adiabatic approximation was a motivating factor for Giovanni Vignale and Walter Kohn, in deriving the "Vignale-Kohn" (VK) approximation in TDCDFT.Consider trying to incorporate memory-dependence in TDDFT.The xc potential at time t depends on the density at earlier times t < t, when a small volume element which is now at r was, at that earlier time, at r [44], suggesting that memory resides with the fluid element rather than with the local density.This entanglement of spatial-and time-nonlocality was first brought to light by John Dobson in Ref. [45], where he showed that simply modifying ALDA by using the finite-frequency response of the uniform gas while retaining only local density-dependence as in the Gross-Kohn approximation [46], violated exact conditions, such as the harmonic potential theorem.Vignale showed this approximation also violates the zero-force theorem (Sec.II) [47].Any nonadiabatic xc functional must depend in a strongly spatially-non-local way on the density, even when in the limit of spatially-slowly-varying densities.A theory local in space and non-local in time is instead possible in terms of the current-density [40,[48][49][50], and TDCDFT with the VK approximation has proven successful for a number of applications where the standard approximations in TDDFT perform poorly.Due to its inclusion of retardation effects, it has had some success in capturing relaxation due to electron-electron interaction in large or periodic systems [51,52], viscous forces responsible for the stopping power in elec-tron liquids [53], spin-Coulomb drag [54], the frequency and linewidth of intersubband transitions in quantum wells [55], and many-body dynamical corrections to conductance in molecular junctions [56].Spatially local functionals of the current incorporate non-local densitydependence, which allows VK to correct the overestimation in ALDA/GGA of polarizabilities for many polymer chains [57] and to improve the dielectric function in semiconductors to some extent [58].However its general application has been limited, partly because it spuriously damps finite systems [52], and, even if only applied to extended systems, the electron gas response functions that enter the VK functional are not known with certainty at general frequencies and wavevectors.

II. EXCHANGE-CORRELATION POTENTIALS: WHAT DO WE KNOW EXACTLY?
The TDDFT theorems say that the dynamics of the electrons in any atom, molecule, or solid can be exactly extracted from completely removing the electron repulsion while distorting the attractive potential of the nuclei.What kind of crazy distortion is required for this to happen?Here, we discuss some general aspects about what is known about the exact potential.
The exact TD KS potential By inverting the TD KS equation, we find an expression for the exact KS potential in terms of the density n i (r, t) and phase α i (r, t) of any one of the occupied KS orbitals φ i (r, t) = n i (r, t)e iαi(r,t) : expresses the continuity equation for orbital φ i .In onedimension (1D), this simplifies in terms of n i , the orbitalvelocity u i = j i /n i and the orbital-acceleration ∂ t u i , of any one of the occupied orbitals: As an illustration, in Fig. 2, we plot the exact KS potential for a model problem: field-free evolution of the initial state Ψ(0) = (Ψ 0 + Ψ 1 )/ √ 2 where Ψ 0,1 are the ground and first excited singlet states in a 1D He atom.The initial KS state Φ(0) has the same configuration but in terms of non-interacting Slater determinants [31], with the KS orbitals φ 0 and φ 1 that form these determinants being eigenstates of some non-interacting Hamiltonian.The top panels show the KS potential and its decomposition into the three terms on the right-handside of Eq. ( 7), labelled as adia, veloc, and acc, respectively, evaluated using φ 0 (φ 1 ) in the upper(lower) panels, at the initial time, and about half a period later.We notice that the components from each orbital can individually be very different, but they sum to the same v S .It is clear that the KS potential has a complicated structure coming from all contributions, going beyond just a static screening of the "nuclear charge" in v ext .FIG. 2. Exact KS potential for field-free evolution of the 50:50 superposition state example (see text).The upper plots show the three contributions to Eq. ( 7) from one orbital at two different times, while the lower plots show those of the other orbital.
Although Eq. ( 6) and Eq. ( 7) give an explicit orbitalfunctional for the KS potential, the functional itself is utterly useless!Without input (knowledge) of the externally applied potential, it clearly cannot be used to propagate the orbital.(Nevertheless Eq. ( 6) and Eq. ( 7) can be instructive in analysing/understanding features of the exact xc potential, e.g.step features arising from the integral of the local acceleration [59,60], see Sec IV).
In fact, if one did try to propagate using Eq. ( 6) or ( 7), one immediately runs into a problem: At the initial time, the choice of the initial KS state specifies the orbitals, so that, on the right-hand-side of Eq. ( 7) for example, n i and u i are known, however the orbital-acceleration ∂ t u i is not determined by the initial state.In fact, the external potential determines the acceleration, and is hiding in the functional-dependence in this term.It is essential to realize that the external potential is given by the problem at hand, and is not to be treated as a densityfunctional [61]; otherwise one runs into misunderstandings concerning the predictivity of the TD KS evolution.With the external potential v ext (r, t), initial many-body state Ψ(0) and initial choice of KS state Φ(0), as inputs, the exact TD KS propagation is predictive: extracting an exact expression for the xc potential alone shows that it depends only on times past [12,61] as we next show.
The exact xc potential If we differentiate the continuity equation with respect to time for each of the interacting systems and KS systems, requiring the density to be identical, we obtain [22,[62][63][64] Here T is the kinetic energy operator (first term on the right in Eq. ( 1)), and Ŵ is the electron-electron interaction (third term in Eq. ( 1)).At the initial time, the right-hand-side involves expectation values in the given initial states, showing that v XC (r, t = 0) can be constructed as a functional of the initial states alone (unlike v S and v ext ), and therefore together with the v ext prescribed by the problem at hand, and v H computed from the density, v S can be constructed and the propagation can proceed.Now, Eq. ( 8) is an exact expression for the xc potential, which is instructive to rewrite as: where is the spin-summed one-body density-matrix of the true system of electrons with two-body interaction potential w(|r − r |), ρ 1,S (r , r, t) is the one-body densitymatrix for the KS system, and n XC (r , r, t) is the xc hole, defined via the pair density, Eq. ( 9) is a Sturm-Liouville equation for v XC that has a unique solution for a given boundary condition [22].The first term represents a kinetic contribution to the potential, and the second term is an interaction component.The problem of approximating the xc potential as a functional of the density and initial-states is therefore here re-cast as a problem of modeling the interacting density-matrix and xc hole as functionals of the density and initial-states.
The adiabatically-exact approximation In most practical applications, an adiabatic approximation is used for the xc potential, which inserts the instantaneous density into a ground-state approximation, Even if the instantaneous density happens to be a density of some ground state, the underlying interacting and KS states typically are not ground states, and so use of v A XC can lead to significant error.There are two sources of error in v A XC ; one is in neglecting all memory-dependence, and the other is in the choice of the ground-state functional approximation.To disentangle the two sources, the "adiabatically-exact" (AE) approximation is a useful analytic tool, defined by using the exact ground-state xc potential [65,66]: Obtaining such an approximation is clearly only possible for simple model systems where the exact ground-state xc potential corresponding to the instantaneous density may be found.In this way, the AE isolates errors due to the adiabatic approximation alone, and is the "best" that an adiabatic approximation could do.If a ground-state orbital functional, such as exactexchange, or self-interaction corrected LDA (SIC-LDA), is used as v g.s.
XC then it contains some memory: a functional that depends on instantaneous orbitals depends on the history of the density and the initial KS state since the orbitals are themselves implicit functionals of these.The xc potential is technically given by the timedependent optimized effective potential (OEP) equations [19,67].Although not always [68,69], these are most often solved within the KLI approximation, at the sacrifice of the zero-force theorem [70].
Known exact conditions The satisfaction of exact conditions has played, and continues to play, an important role in the development of ground-state functionals (see Ref. [71] for a recent advance).When such functionals are adopted in an adiabatic functional for TDDFT, these exact conditions can endow v A XC with desirable properties also for dynamics and excitations.An important example of this is the asymptotic behavior of the potential far from an atom or molecule: v XC (r) → −1/r [72].This is perhaps more important in the TD case than the ground-state case, because the density is more likely to "sample" regions far from the nuclei over time, and certainly the accuracy of energies of higher-lying excitations depend on this decay behavior (see also Sec.III).Relatedly, via Koopman's theorem, H = −I, the highest occupied molecular orbital (HOMO) energy determines the onset of continuous absorption in the spectrum.
Several exact conditions that are specific to the TD case are known, and some have shaped the development of non-adiabatic functional approximations.For example, the harmonic potential theorem [45] and the zero-force theorem [47] were instrumental in the discovery that a spatially-local but time-non-local approximation is not possible in terms of the density, but is possible in terms of the current-density (Sec.I C).The harmonic potential theorem, or generalized Kohn's theorem, states that the electron density in a harmonic potential subject to a uniform TD selectric field rigidly follows the classical oscillation of the center of mass.The zeroforce theorem is an expression of Newton's third law, requiring the total force exerted on the system by the xc and Hartree potentials to vanish.For TDCDFT, the analogous theorem for the xc torque is also satisfied [49].Imposing these conditions were key in leading to the nonadiabatic VK approximation of TDCDFT [40,48] (Sec II.C).
Exact conditions unique to the TD problem, such as the zero-force and harmonic potential theorems, naturally involve the memory-dependence (time-nonlocality) of the xc functional.Another such condition is the "memory-condition" which entangles the initialstate dependence and history-dependence of the xc potential [16].Essentially it states that one can consider any fixed time t along an evolution as the initial time, and so v XC [n t ; Ψ(t ), Φ(t )](r, t > t ), where n t is the density function from time t onwards, must be the same for all choices of such t .This is a very hard condition to satisfy for any memory-dependent functional, and is actually violated by VK.
Due to its complete lack of memory, the adiabatic approximation satisfies the memory-condition, the zeroforce and harmonic potential theorems.But it violates the "resonance condition" [73], with severe consequences for resonantly-driven processes, and possibly for general non-equilibrium dynamics.This is a condition on the density-response function of a system in a general non-stationary state, and expresses the fact that the excitation frequency of a given transition is independent of the state.This involves a delicate cancellation of time-dependences between the KS response function around the KS non-stationary state and shifts induced by the generalized xc kernel such that the poles of the true density-response function remain at the same frequency, and is generally hard to satisfy.

III. EXCITATION SPECTRA: LINEAR RESPONSE
The vast majority of applications of TDDFT today lie in spectroscopy, to compute the linear response and excitation energies of atoms, molecules, and solids.TDDFT's rise to fame is largely due to its unprecedented balance between accuracy and efficiency for these applications.Figure 3 gives a pictorial representation of a few of the most recent successes; these are chosen to illustrate the wide range of applications in which TDDFT today is giving useful results.I refer the reader to Refs.[74][75][76][77][78][79][80] for recent reviews that give more detail on the state-of-the-art performance for energies and excited state properties and instead here I highlight challenges for the field after first giving a brief outline of the linear response TDDFT formalism.Consider perturbing a system initially in its groundstate Ψ 0 with density n 0 , with a weak field: ext represents the field-free unperturbed external potential, and δv ext (r, t) the perturbation, turned on at t = 0.The Fourier transform of the first-order response n (1) (r, t) to the δv ext (r, t) yields excitation energies and oscillator strengths of the unperturbed system.Polarizabilities can be obtained from the linear response of the dipole to electric fields while higher-order response yields hyperpolarizabilities.The density-density response function is defined as χ[n 0 ](r, r , t − t ) = δn(r, t)/δv ext (r , t ), and, for finite systems, has the Lehmann (or spectral) representation in frequency-space: where n(r) = i δ(r − ri ) is the density operator, Ω I = E I − E 0 is the excitation frequency of the Ith excited state, and the sum goes over all interacting states Ψ I .The notation c.c.(−ω) means the complex conjugate of the first term evaluated at −ω.Here the experiment is from Ref. [83], while the TDDFT results are from: the nanoquanta kernel [84], the dynamical long-range corrected kernel (dynLRC) [85], the bootstrap kernel on top of an exact-exchange OEP groundstate (BO oep) [86], the jellium-with-gap kernel on top of the scissors-corrected LDA ground-state (JGM lda sx) [87], the guided iteration bootstrap kernel on top of scissors-corrected LDA (GIBO lda sx) [88], screened exchange on top of scissorscorrected LDA (SXX lda sx) [89], and the optimally-tuned range-separated hybrid (OTSRSH) [90].(d) Table of errors from different functionals for 103 vertical excitation energies of singlet valence excitations in 28 small and medium-size organic molecules from data of Refs.[75,91], as compiled in Ref. [92].
Now, because the TD KS system gives the same TD density when perturbed via the corresponding KS perturbation δv S (r, t), all response quantities can be directly obtained from the KS system.The KS density-density response function, χ S [n 0 ](r, r , t−t ) = δn(r, t)/δv S (r , t ) naturally has a Lehmann representation in frequencyspace for finite systems, but the poles lie at KS singleexcitations, i.e. orbital-energy differences between orbitals that are occupied and unoccupied in the groundstate.Further, the residue of the pole ω ia = a − i reduces to simple orbital product: φ i (r)φ a (r)φ i (r )φ a (r ).
A Dyson-like equation relates χ S to the interacting χ [93,94]: using to denote spatial convolution integrals such as χ S (ω) f HXC (ω) = d 3 r 1 χ S (r, r 1 , ω)f HXC (r 1 , r , ω), we have This equation features the indispensable Hartree-xc kernel, The kernel plays a crucial role in dressing the KS response function, shifting the KS poles to the true excitation energies, and mixing the KS oscillator strengths to form the true oscillator strengths.With the ubiquitous adiabatic approximation, since v A XC [n](r, t) depends only on n(t) at the present time, the right-hand-side of Eq. 13 becomes proportional to δ(t − t ) and the Fourier transform frequency-independent, f A XC (r, r , ω) = f XC (r, r ).Equation ( 12) is the basis of TDDFT linear response.There are evidently two ingredients: first the elements of the KS density-response function, involving all occupied and unoccupied KS orbitals living in the groundstate KS potential v (0) S .Second, the xc kernel, which introduces so-called dynamical effects.Technically, the approximation used for the ground-state xc potential should be linked to that used for the xc kernel since they are generated from one and the same TD xc potential, v XC [n](r, t) = v (0) XC (r) + δv XC .But in practise, the two approximations are often treated independently.
Typically in quantum chemistry codes, the equation is re-cast in the form of matrix equations derived by Mark Casida, based on considering poles and residues of the frequency-dependent polarizability [95,96].The eigenvalues of the matrix R qq (ω) = ω 2 q δ qq + 4 √ ω q ω q f qq HXC (ω) , where f qq HXC (ω) = d 3 rd 3 r φ i (r)φ a (r)f HXC (r, r , ω)φ i (r )φ a (r ), lie at the squares of the excitation energies, and oscillator strengths can be obtained from the eigenvectors.Here q = (i, a) represents a double-index, with i labelling an occupied orbital and a an unoccupied one.Other equivalent derivations and formulations appear in Refs.[97,98]; note that the derivation in Ref. [97] is valid only when an adiabatic approximation is assumed from the outset.
The matrix formulation is only valid for discrete spectra, so used mostly for atoms and molecules, while Eq. ( 12) is usually solved when dealing with continuous spectra as in solids.
To obtain the continuous part of spectra of atoms and molecules, in particular resonance widths and positions, the Sternheimer approach is often used [99][100][101].But even for bound excitations, the Sternheimer approach, also known as density functional perturbation theory, or coupled perturbed KS, can be more efficient, if just a few low-lying excitations are of interest; it avoids solving for many unoccupied states by considering instead perturbations of the occupied KS orbitals in frequency-space.Response can also be obtained from real-time propagation of occupied orbitals [102] under a weak perturbation; often a δ-kick is applied so to uniformly stimulate all excitations, and the system is then evolved freely for a time long enough such that Fourier transform of the dipole yields well-enough resolved peaks at the resonant frequencies.This is, for example, how the absorption spectrum of such a large system as that in Fig 3b could be computed [82].

A. Challenging Excitations
The successes of linear response TDDFT have been celebrated in a number of reviews cited earlier.Its efficiency and overall reasonable reliability have led to it becoming an easy-to-use "go to" method for simulating spectra that, for other methods of comparable accuracy, are not easily computationally accessible.However, one must be cautious in treating it completely as a black-box: there are some notorious problems where either a careful choice of functional is needed or, where the results with currently available functional approximations are simply unacceptable.We discuss some of these here.
Rydberg excitations Here the problem lies primarily with the ground-state approximation, rather than the kernel: Because LDA and GGA ground-state KS potentials depend locally or semilocally on the density, they asymptotically approach zero exponentially with r, much faster than the −1/r fall-off of the exact KS potential.Without the −1/r tail, they cannot generate a Rydberg series.The frequencies of the highestlying excitations that are bound within these approximations tend to be underestimated, while even higher excitations are absent, with their oscillator strengths being squeezed into the continuum albeit with accurate optical intensity [103][104][105][106][107]. Different solutions include asymptotically-corrected functionals [104,106,108], exact-exchange methods [109,110], and hybrids [111][112][113], including range-separated hybrids [114][115][116].Conventional hybrid functionals, formally justified by generalized KS theory [117][118][119], mix in a fraction C X of Hartree-Fock exchange into otherwise semi-local functionals, which means their potentials decay asymptotically as −C X /r.On the other hand, range-separated hy-brids produce the requisite −1/r: the essential idea is to separate the Coulomb repulsion using a parameter γ as, e.g., where the first term dominates at short-range and is treated via a semi-local functional, and the second term dominates at longer ranges and is treated via exact (Hartree-Fock) exchange.The result yields a selfinteraction-free description at large electron-electron separations, while retaining the balanced description of exchange and correlation at short range.Double excitations It was recognized soon after the derivation of the linear response equations, that states of double-excitation character were missing from the approximate TDDFT spectra; this time the fault lies in the adiabatic approximation to the kernel [105,120].It is clear they are absent in the non-interacting response function χ S : considering Slater determinants in the numerator of Eq. ( 11), doubly-excited KS states give a zero contribution to χ S because the matrix elements of the density-operator vanish between Slater determinants that differ by more than one orbital due to its one-body property.This is just the mathematical statement of the fact that exciting two electrons in a non-interacting system requires two photons, so is beyond linear response.Excited states, as well as the ground-state, of the interacting system, on the other hand, are superpositions of single, double, and higher-electron-number excitations, and so the numerator Ψ 0 |n(r)|Ψ I remains finite from contributions between components of the two manybody wavefunctions that differ by one orbital.A simple number-counting argument [121,122] shows that no adiabatic kernel can introduce doubly-excited configurations.The strong frequency-dependence of the exact kernel near a state of double-excitation character was deduced [121,123] by in effect inverting the Dyson equation Eq (12) for a model where one KS single and one double excitation interact, from which a practical approximation was proposed.More rigorous derivations using many-body theory subsequently essentially verified this form [124][125][126][127] (see Ref. [78] for a recent review), and the approximate kernel, dubbed "dressed LR-TDDFT", has been implemented in some codes (Niedoida and deMon2k), enabling extensive tests [125,128,129].In particular, Ref. [125] found that the non-adiabatic dressing has best performance when used in conjunction with a global hybrid functional.
Long-range charge-transfer excitations As illustrated in Eqs. ( 12)-( 14), TDDFT excitations are obtained from shifting and mixing single-excitations of the KS system via matrix elements of f XC , often referred to as dynamical corrections.When there is little spatial overlap between the occupied and unoccupied orbitals relevant in the transition, the dynamical corrections from conventional xc-kernels vanish, and the TDDFT excitation energy reduces to KS orbital energy difference.Excita-tions of charge-transfer character over long range are a notorious example of this, and was first brought to light in Ref. [112].In the large separation limit, the exact lowest charge-transfer excitation energy tends to ω exact CT → I D − A A − 1/R, where I D is the ionization energy of the donor, A A is the electron affinity of the acceptor, and −1/R is the first electrostatic correction between the species after the charge-transfer.As analyzed in Refs.[130][131][132][133] however, TDDFT with LDA or GGA yields the orbital energy difference between the LUMO of the acceptor and the HOMO of the donor, L (A) − H (D). With approximate ground-state functionals H underestimates the true ionization energy, but even with the exact ground-state functional L lacks relaxation contributions to the electron affinity, and further the −1/R tail is missing.A useful diagnostic of the performance of a given functional based on considering the overlaps was given in Ref. [115].Significant effort has been devoted to solving this problem, e.g.Refs.[113,116,130,[134][135][136][137][138][139], especially because charge-transfer excitations play an important role in many topical applications today that involve systems large enough that few practical alternatives to TDDFT exist, e.g.biomolecules, photovoltaics, conductance.Most modify the ground-state functional to correct the underestimation of the ionization potential, and mix in some degree of Hartree-Fock exchange, especially via range-separated hybrids [117][118][119], since they display the −1/R asymptotic behavior.The amount of long-range exchange needed to give reasonable answers varies widely depending on the system, which motivated the development of the non-empirical "optimallytuned" range-separated hybrids of Refs.[137,138].The range-separation parameter is tuned so that H is as close as possible to the negative ionization potential computed by the same functional from total energy differences, typically for both the neutral and anionic system.Size-consistency violation and other problems [140] arising from system-dependence of these functionals means caution must be applied for general usage; the range-separation parameter can have very different values for the individual subsystems and the combined system, and varies along the internuclear separation, leading to poor binding energies and sometimes erratic potential energy surfaces [140].
The xc-kernel in pure (i.e.not generalized) KS DFT must diverge exponentially with internuclear separation [133], and it has been argued that discontinuities of the xc-kernel with respect to particle number play a crucial role [141] whenever the KS orbital overlap is too small.The exact-exchange approximation (within OEP) captures these features when the chargetransfer occurs between closed-shell fragments, when used in its fully non-adiabatic form [139,141].The selfinteraction-corrected LDA, applied within a generalized OEP scheme, has also been shown to work [69].
However, one case where none of the above suggestions works is that of charge-transfer between open-shell fragments [142,143].At the root of the difficulty in this case is the KS ground-state description: the HOMO (and LUMO) straddle both fragments, delocalized over the whole molecule, quite distinct from the Heitler-London-like nature of the physical wavefunction.Although in this case the HOMO and LUMO have finite overlap, their orbital energy difference tends to zero as the molecule stretches (for either the exact functional or approximations), a sign of static correlation in the KS system, and f HXC must be responsible for the entire charge-transfer energy.It can be shown that the exact kernel not only diverges with interatomic separation but also is strongly frequency-dependent.
Conical intersections Photoinduced dynamics is heavily influenced by conical intersections, seams of exact degeneracy between electronic potential energy surfaces that are ubiquitous in the landscape of nuclear configurations.They funnel nuclear wavepackets between surfaces, providing a rapid pathway for non-radiative decay, and even affect nuclear trajectories that encircle the intersection at large distances, by generating a geometric phase.Unfortunately, current approximations in TDDFT typically provide a poor description of the surfaces near conical intersections with the ground-state.In some cases the shape is dramatically exaggerated (e.g.H 3 ), while in other cases, the dimensionality of the seam is wrong (e.g.H 2 0) [144].There is trouble here for both the ground-state approximation and the xc kernel approximation, and, conceptually, the issues are not unlike the case of charge-transfer excitations between open-shell fragments.Static correlation resulting from nearly degenerate determinants near the conical intersection creates unusual features in the exact groundstate xc functional in order to compensate for the KS Slater determinant description being so far from the true correlated ground-state.The xc kernel then must be strongly frequency-dependent because of the doubleexcitation to the near-lying determinant.In practise, the very small HOMO-LUMO orbital energy difference leads to an instability in the solution of the LR equations, that can be tamed by making the Tamm-Dancoff approximation [145].In fact, in Ref. [145], where ringopening in oxirane was studied, it was argued that although TDDFT gets the dimension of seam wrong, giving slightly interpenetrating cones rather than touching ones, that the resulting dynamics (using surfacehopping) still yielded useful results regarding reaction pathways, although rate constants were likely not accurate.Also, spin-flip TDDFT gets around the problem of the ground-state degeneracy by "viewing" the conical intersection from a high-spin reference state [146].
Optical response of solids Simple ALDA accurately captures plasmon dispersion curves in simple metals, with the xc kernel providing an improved description of the plasmon dispersion and linewidth compared to the random phase approximation that uses the Hartree kernel only [147].For non-metallic systems, ALDA provides a good description of electron energy loss spec-tra (EELS) [80,148], and for the optical response at finite photon wavevector q [149], but does a poor job for the spectrum at vanishing q.There are two main problems: one is that the optical gap, defined by the onset of continuous absorption, can be underestimated by as much as 30-50%, and the second is that the excitonic structure is missing from the spectrum.For both of these problems, a long-ranged xc kernel is necessary [50,80,[150][151][152]: the xc kernel must have a long-wavelength behavior of 1/q 2 as q → 0. In fact, this behavior can be deduced from exact conditions satisfied by the xc kernel, such as the zero-force theorem, that inextricably link nonlocality in time and space [50,153].ALDA, on the other hand, being local in space, is constant in q.Further, to open the gap, the xc kernel must have a non-vanishing imaginary part [12,50] which, due to Kramers-Kronig relations, implies that it must have a frequency-dependence, i.e. be non-adiabatic.Instead, today a common workaround is to use LDA bandstructure with a scissors operator [154,155].
A non-adiabatic kernel is not necessary to capture the excitonic peak, but the 1/q 2 behavior is essential [150].(The excitonic Rydberg series however does appear to require a frequency-dependent kernel [156]).Interestingly, although 1/q 2 is the Fourier-transform of −1/r, this required 1/q 2 -dependence of the kernel is not a consequence of the long-rangedness of the Coulomb interaction, and it has been dubbed an "ultranonlocal" effect since it is even present in systems arbitrarily close to the homogeneous liquid [50,153].
Much progress has been made to correct this deficiency of ALDA, along several different avenues, as evidenced by the success in Fig. 3c of capturing the strongly bound exciton in LiF from a variety of different modern functional approximations.Exact exchange initially seemed to give promising results [157] until it was realized that the calculations only worked because of an artificial screening induced by a long-wavelength cut-off in the numerics [158].Kernels derived from the Bethe-Salpeter equation, such as the nanoquanta kernel, have seen much success [84,150,159].The long-range corrected (LRC) kernel [85,150,160], α/q 2 , with α determined empirically from the static dielectric constant of the crystal, can be viewed as a simple scalar approximation to the nanoquanta kernel; a dynamical version uses also the plasma frequency.The "bootstrap" kernel proposed in Ref. [86] uses instead self-consistent iterations of Eq. ( 12) with an approximate expression for the (static) xc kernel in terms of the dielectric function; going beyond the scalar approximation, it includes localfield-effects.A guided-iteration version has been presented in Ref. [88] where the sensitivity of the bootstrap approach to details such as choice of groundstate orbitals and energies from which to build χ S (usually scissors-shifted LDA, sometimes instead with GWcorrected energies, sometimes OEP), and treatment of the local field effects, has been stressed.The VK approximation has recently been shown to provide a parameter-free broadening of the spectra [161].In yet another direction, a kernel built on a jellium-with-gap system [87] has been developed.Meta-GGA functionals also display the requisite α/q 2 behavior [162], as do hybrid functionals [163], and optimally-tuned screened rangeseparated hybrids have recently been applied [90].Relatedly, screened exact-exchange [89], which formally falls into generalized KS theory, has been proposed as a bridge between TDDFT and BSE; it is computationally much cheaper than BSE yet of comparable accuracy.Finally, Ref. [164,165] showed how to directly obtain the excitonic binding energies instead of the full spectra by adapting the Casida-equation formalism.
Molecular Conductance An area of great interest is in simulating charge transport at the nanoscale.Most ab initio descriptions are in the steady-state using DFT within a Landauer formalism, but give poor agreement with experiment due to problems with both the groundstate xc potential and the dynamical f XC corrections.Regarding the former, the lack of derivative discontinuity [166] in LDA and GGA leads to a large overestimate of conductance [167,168] and artificial level broadening, and is intimately related to Coulomb blockade, a fully dynamical picture of which was provided in Ref. [169].The discontinuity also needs to be accounted for in the xc kernel, where, in addition, viscous effects can be important [56,170].An exactly solvable impurity model has demonstrated the importance of long-range nonadiabatic terms in the xc potential [171].The exact theoretical framework for transport in TDDFT for strong and weak bias, requiring a real-time description beyond Landauer, was formulated in Refs.[172,173].To include dissipation and decoherence induced by the environment, several approaches to extend TDDFT to open quantum-systems were developed [174][175][176][177][178], relevant both in and out of the transport context.
The lack of derivative discontinuity in (semi-)local functionals is sometimes called the "gap problem", however it is worth realizing that several distinct gaps can be defined, and functionals have no problem with some of these.The KS HOMO-LUMO energy difference approximates neutral excitations such that, once the TDDFT f XC corrections are applied, one obtains the lowest optical excitation, or optical gap; (semi-)local functionals give pretty accurate results for molecules.The fundamental gap, on the other hand, I−A, involves electron addition and removal, and in particular L should not be interpreted as an affinity level.When a fraction of electron is added to a neutral system, the exact KS potential jumps up by a positive constant equal to the derivative discontinuity, which then shifts the orbital energy up to the affinity level [166]; such a shift is not captured by semi-local functional approximations that are smooth functions of particle number [179,180].(That it is in fact possible to extract derivative discontinuities from LDA and GGA has been pointed out in Refs.[181,182], but this is not how these functionals are usually used.)Because Hartree-Fock virtual orbitals rep-resent addition-energy levels rather than neutral excitations, hybrid functionals evaluated within the generalized Kohn-Sham scheme effectively open the gap, especially the range-separated ones [183].Lack of the discontinuity is partly responsible for the so-called level alignment problems of standard functionals at moleculemetal, molecule-semiconductor, and molecule-molecule interfaces in applications relevant to transport as well as to photovoltaic modeling and to photocatalysis [184][185][186][187].For metal-molecule interfaces non-local polarization effects also contribute significantly to the renormalization of the frontier energy levels, for which "imagecharge" models have been helpful to capture [184,185].

IV. FULLY NON-LINEAR DYNAMICS
Notwithstanding the successes of TDDFT for linear response, TDDFT has particular relevance in the regime of non-perturbative time-resolved dynamics.There is really no alternative practical method for accurately describing correlated electron dynamics when more than a few electrons are involved, yet many phenomena which are fascinating in themselves and also have technological applications, lie in this regime.Although no strong fields are present during simulations of photoinduced coupled electron-ion dynamics (e.g. in light-harvesting and artificial photosynthesis) and other photo-chemical or -physical processes, the electrons do not remain close to any ground state, so also lie beyond the perturbative regime.
TDDFT has provided useful predictions for a variety of strong-field processes and below I select some to briefly describe.The examples are chosen to illustrate the wide range of applications of TDDFT in this regime, and emphasise recent achievements over older ones; inevitably such a short description is a little random.
In solids, coherent phonon generation produced by intense laser pulses has been simulated [188], capturing the transition from the impulsively stimulated Raman scattering mechanism to the displacive excitation mechanism and the associated phase change of the generated phonon.Dielectric breakdown in crystalline silicon has been studied, where the importance of including both electric and magnetic fields self-consistently has been stressed; the authors developed a multiscale model solving Maxwell's equations alongside TDKS [189].The same method showed that TDDFT calculations of laserinduced ablation by very short, intense laser pulses in α-quartz gives reasonable agreement with experiment [190].TDDFT has elucidated the possible mechanisms of the experimentally-observed laser-induced ultrafast demagnetization and details of spin-dynamics in bulk ferromagnets in Ref. [191].
Multiphoton processes have been studied in atoms and molecules, with the first calculations of highharmonic generation spectra and ionization rates in atoms emerging in the mid-1990's [192] capturing cor-relation effects missed in the single-active-electron approximation.Because the electrons make large excursions away from the nuclei, exact-exchange or selfinteraction-corrected-LDA are often used, implemented within a KLI approximation to the OEP [193].TDDFT revealed interesting interference effects between contributions from different molecular orbitals in highharmonic generation and multiphoton ionization in molecules, e.g.Refs.[194][195][196][197], where, additionally, orientation effects have been studied.Time-resolved Coulomb explosion imaging of molecular torsion dynamics in complex molecules has been shown to give results consistent with experiment [198].TDDFT has begun to elucidate the mechanism of charge-resonance enhanced-ionization in triatomic molecules [199], identifying which molecular orbitals are responsible; together with Ehrenfest dynamics for the nuclei, enhanced ionization and orientation dependence has been analyzed in acetylene and ethylene [200].The timedependent electron localization function introduced in Ref. [201] has been used to understand time-resolved dynamics of chemical bonds in scattering and excitation processes.Recently attosecond "self-streaking" has been proposed and demonstrated on the endohedral fullerene Ne@C 60 , from which one can directly observe electronic density oscillations in real-time [202].TDDFT has facilitated the analysis of photoelectron spectra and angular distributions in molecules and clusters after irradiation by strong pulses; a review can be found in Ref. [203].
Optimal control theory has been formulated with TDDFT [204], and used already to obtain desired targets in an array of different scenarios, from selectively cleaving a molecular bond [205] to enhancing a particular harmonic in high-harmonic generation [206].Very recently, the possibility for control of spin-orbit mediated ultrafast demagnetization via lasers in solids has been demonstrated [207].
Alongside this rosy picture, however, studies on small molecules where numerically exact or high-level wavefunction methods can be applied, have shown that approximate TDDFT functionals can yield significant errors in their predictions of the dynamics, sometimes even failing qualitatively.We next discuss these problems.

A. Challenges
The problems fall into two main categories: the inaccuracy of the approximate xc potential functional, and the need for additional functionals to extract observables of interest from the KS system.For the first, we focus here mostly on problems associated with making the adiabatic approximation, rather than the specific choice of ground-state functional that goes into such an approximation.Exactly-solvable model systems are useful to explore the impact of these problems, and to compare approximate functionals, especially the AE, with the exact one.
Memory appears to generally play a more important role in non-perturbative time-resolved dynamics compared to the linear response regime [52,59,63,65,68,73,[208][209][210][211][212].The problems appear to require a memorydependence beyond what is contained in the presentlyavailable orbital-functionals.To investigate the effects of the adiabatic approximation, comparison of exact xc functionals with the AE of Eq. 10 has proven quite instructive [59, 64-66, 211, 213-218].It was shown that for strong field ionization, the xc potential appears not to have significant non-adiabatic dependence in some cases, but this depends on the details of the field.Sometimes the observables of interest are not so much affected by non-adiabatic features present in the xc potential, e.g.integrating over regions of space to find total ionization probabilites can wash over differences in the density variations.In fact, the famous knee in the double-ionization curve that signals a switch from a correlated to sequential ionization mechanism, can be captured by a non-adiabatic approximation [219] provided it has a derivative discontinuity.
But in the general case, it appears that neglect of memory-dependence in the xc-potential can lead to large errors, and part of the problem appears to come from the lack of dynamical step features that the exact xc functional possesses, that have a non-adiabatic dependence on the density [59,60,64,211,220].In the dynamics of charge-transfer beginning in the ground-state, the lack of these features have a drastic effect [211,217,[221][222][223]; even the AE fails to transfer charge even though its predictions for the charge-transfer excitations themselves can be accurate [213,214].The non-adiabatic step structures have been shown to appear quite generically in all sorts of dynamics, and are not associated to ionization or fractional charge as other observed timedependent steps are [141,169,219].They appear in the presence or absence of a field, for few electrons [59] or many electrons [60], whether one uses a Slater determinant or not [31] for the KS system (see also Fig. 2).Unlike step features that appear in ground-state potentials, they are not a signature of static correlation, and their relationship with time-dependent natural orbital occupation numbers is non-trivial [64].The AE potential completely lacks these steps, as illustrated in Figure 4 (Multimedia view) which, together with a movie, plots the xc potential in a model He atom driven by a weak resonant field that causes Rabi oscillations between the ground and first-excited state.
Violation of the resonance condition mentioned at the end of Sec.II can wreak havoc and confusion in analysing time-resolved spectra as is now experimentally accessible via ultrafast pump-probe spectroscopy [224,225].Several researchers had noted the peaks in their spectra would unphysically shift when evolving far from the ground-state [221][222][223][226][227][228]; having peak positions that are somewhat off are one thing, but having them continuously drift in time is more serious.
Turning now to the second category of problem: the need for "observable-functionals" when the observable of interest is not directly related to the density.Although in principle all properties of interest can be extracted from the KS system, it is not generally known how to extract them.Only in some cases can one use the usual quantum mechanical operators acting directly on the KS wavefunction, e.g.high-harmonic generation spectra are measured from the dipole moment of the system.But for general observables, e.g.crosssections in atomic collisions [229], ionization probabilities [230][231][232][233], or momentum-distributions [231], substituting the KS wavefunction into the usual quantum mechanical operator typically gives poor results.For example, in ion-recoil upon ionization of a model He atom, the KS momentum distributions were found to be drastically wrong, displaying a single maximum instead of the characteristic two-hump structure, and with a significantly overestimated magnitude [231].More sophisticated functionals have been developed for a few particular cases; momentum-densities for the above problem in Ref. [231], and double-ionization probabilities, that "lower" the knee, in Ref. [234].
The one-body nature of the KS Hamiltonian is somewhat of a double-edged sword in the non-perturbative regime.On the one hand, it makes impressive calculations on large systems doable, but on the other hand, an initial Slater determinant remains a Slater determinant throughout the evolution, even when the interacting wavefunction significantly changes its degree of correlation [235].Even if initially the physical system is weakly correlated, and adiabatic approximations with the usual ground-state approximations initially work well, over time excited states get populated with the degree of correlation possibly significantly changing, as reflected in the natural orbital occupation numbers.If one perfects a time-dependent non-adiabatic functional approximation (for either the xc functional or an observable) to work in the weakly-correlated regime, it will not work well for general dynamics.

B. Coupled electron-ion dynamics
Dynamics of electrons in molecules usually cannot be considered in isolation from the dynamics of the ions to which they are coupled and the phenomena arising from the correlated dynamics of electrons and ions even without worrying about electron-electron interaction themselves form a very rich field.In simulating photo-induced processes, it is often the case that the nuclear motion drives the electronic dynamics after the initial excitation.
Typically, although not always, a classical treatment of nuclear motion is coupled to the electronic dynamics.The most common methods are Ehrenfest dynamics and trajectory surface-hopping, and here we restrict our discussion to an analysis of their use specifically in conjunction with TDDFT, rather than giving a more extensive discussion of their relative merits and pitfalls.
It is well known that the mean-field nature of Ehrenfest dynamics is inadequate when branching of nuclear trajectories occurs, e.g. in electron-transfer reactions, reactions near metal surfaces, and photochemistry generally speaking.However, as a formalism, Ehrenfest is well-suited to be used in conjunction with TDDFT: the classical force on the Ith nucleus is obtained via where V nn is the classical nuclear-nuclear potential, and H el is the Hamiltonian from Eq. ( 1) with R appearing in v ext evaluated at the instantaneous nuclear positions.The only term that the gradient in Eq. ( 16) acts on is the electron-nuclear attraction V en = − Nn I=1 N i=1 Z I /|r i − R I | (where N n is the number of nuclei and Z I their charge), and so the first term, the quantum "backreaction" on the nuclear trajectories, becomes − Ψ(t)|∇ I H el |Ψ(t) = d 3 rn(r, t)∇ I V en ({r i }, {R I }).Ehrenfest therefore needs nothing but the electronic density, which can be directly obtained from the TD KS evolution.Notably, Ehrenfest was recently used to demonstrate the relevance of coherent coupled electronic and vibronic motion in driving the chargetransfer in a prototypical light-harvesting system [236,237].
Surface-hopping on the other hand, is able to capture wavepacket branching but its marriage with TDDFT is a little more stormy due to its need for non-adiabatic couplings.Trajectories evolve on Born-Oppenheimer surfaces obtained from TDDFT linear response [238][239][240][241][242] (although approximate methods that use the bare KS surfaces have also been used [243], with notable applications in photovoltaics [244]), and make jumps between them based on Tully's fewest-switches surface-hopping algorithm [245].Such calculations have shed light into the photochemical reactions going on in vitamin D synthesis [246], for example, and combined with local control theory to design pulses to enhance proton transfer reaction in 4-hydroxyacridine [247].The algorithm requires non-adiabatic couplings (NACs) between surfaces, and it has been shown that these are rigorously available from linear response TDDFT for NACs involving the ground-state [248][249][250][251][252].But NACs between excited states, which are clearly most relevant for photoinduced dynamics, where typically the simulation starts in an electronic excited state, are strictly not available from linear response.They can be obtained instead from quadratic response theory, however it was recently shown that the adiabatic approximation yields unphysical divergences in the couplings whenever the energydifference between two excited states equals another excitation energy [253][254][255][256].
Besides this, the accuracy of these methods are clearly limited by the accuracy of the forces obtained from the TDDFT surfaces, so the problems of Sec III enter.In particular, the dynamics probes nuclear configurations away from equilibrium, so even in situations where these problems are not so important at equilibrium configurations, they can plague the dynamics, e.g. the case has been made that without adequate treatment of doubly-excited states, global potential energy surfaces are unlikely to be accurate [144].

V. OUTLOOK
From some of the discussion above, the features of the exact xc functional and the subtleties of extracting observables, one might be tempted to conclude that some magic is required to make the ingredients and goals of TDDFT get along, and that TDDFT actually stands for Thaumaturge-Dependent DysFunctional Theory.But one should bear in mind that TDDFT is a relatively young field; indeed the success it has had already is remarkable.The variety and size of the systems studied these days in TDDFT is impressive, as both Fig. 3 and the non-linear calculations described at the beginning of Sec.IV suggest.Its impact in materials science and quantum chemistry is incontrovertible, and evidenced by its increasing use in an ever-widening range of applications, in both the response regime and the fully nonperturbative regime.
Nevertheless, the results on the model systems where exact potentials are available more often than not suggest that, in the non-perturbative regime, the currently used functional approximations may have significant errors.The model systems do probe perhaps the worstcase scenarios for approximate TDDFT, and the system-size scaling of the xc features missing in approximations should be explored.Characterizing when and why the current approximations, particularly the adiabatic ones, work in a given situation remains an open problem today, but an important one in establishing the reliability of the results in different situations.In the linear response regime, such characterization has greatly improved the predictability of TDDFT for excitation spectra, and our understanding of when the results can be trusted.The adiabatic approximation trivially satisfies many of the exact conditions, some of which were mentioned in Section II, which might be behind the unexpected fact that an approximation based on ground-states can give useful results for problems far beyond the ground-state, as the applications particularly in Sec.IV have shown.
New developments for functionals might take advantage of the initial-state dependence of the Runge-Gross theorem.For a given initial state of the physical interacting system, there is an infinite choice of allowed KS initial state, and whether one can make a judicious choice for the best one to use when stuck with an adiabatic xc functional has only begun to be explored [29,31,215].Even in the linear response regime, using excited states as alternative references around which perturbations yield excitation frequencies could be very useful; for example, a double-excitation is a single-excitation of a singly-excited state, so the need for non-adiabaticity could be avoided.Conical intersections involving the ground-state would appear simply as a degenerate deexcitation from an excited state.Already, calculations of excited state absorption in molecules have appeared by kicking excited-state densities obtained from linear response [257].
It may be time to consider cutting the apron strings from ground-state DFT, and try starting points other than the adiabatic approximation for building approx-imations.The expression for the exact xc potential in Eq. ( 8), approximated to yield an orbital-dependent functional, may be a useful starting point.
Alongside new functional approaches, I have no doubt that the applications of TDDFT will continue to grow, taking on ever-larger systems possibly with the help of embedding-and fragment-based methods [258][259][260][261][262], and adventuring into new regimes such as the recently-introduced quantum electrodynamical density functional theory that includes coupling to quantized radiation fields [263][264][265][266], and to quantum ions in a multicomponent formalism [267].It will be exciting to see where TDDFT will be in another thirty years.

FIG. 1 .
FIG. 1. Illustration of the one-to-one density-potential map in TDDFT.The upper (lower) ellipses contain density evolutions of interacting (non-interacting) systems arising from the initial state labelling the ellipse.The four initial states shown all have the same density as each other, and first time-derivative of the density as each other.The potentials in the central ellipse are one-body potentials, representing vext (vS) when pointed to from the upper (lower) ellipses.

FIG. 3 .
FIG.3."Splash"-plot illustrating the range of TDDFT applications today for excitations and response: (a) In Ref.[81], the rate of yellowing of Leonardo da Vinci's self-portrait began to be ascertained by comparing non-destructive optical reflectance measurements with TDDFT calculations of the relevant chromophoric spectra.The image here is from the Wikimedia Commons page for Leonardo da Vinci [Public domain].(b) In Ref.[82], Published by the PCCP Owner Societies, the absorption spectrum of the complete chlorophyll network of the light-harvesting complex from green plants, was computed using TDDFT, yielding implications for intrinsic energy transfer pathways within the LHC-II complex.Reproduced under a Creative Commons Attribution 3.0 Unported Licence.(c) The bound exciton in the absorption spectrum of LiF was successfully captured by several approximate functionals developed from quite different approaches (see also text).Here the experiment is from Ref.[83], while the TDDFT results are from: the nanoquanta kernel[84], the dynamical long-range corrected kernel (dynLRC)[85], the bootstrap kernel on top of an exact-exchange OEP groundstate (BO oep)[86], the jellium-with-gap kernel on top of the scissors-corrected LDA ground-state (JGM lda sx)[87], the guided iteration bootstrap kernel on top of scissors-corrected LDA (GIBO lda sx)[88], screened exchange on top of scissorscorrected LDA (SXX lda sx)[89], and the optimally-tuned range-separated hybrid (OTSRSH)[90].(d) Table of errors from different functionals for 103 vertical excitation energies of singlet valence excitations in 28 small and medium-size organic molecules from data of Refs.[75,91], as compiled in Ref.[92].

FIG. 4 . 2 + 1 +
FIG. 4. The xc potential (black), plotted at a quarter of a Rabi period, in a model He atom driven by a weak resonant field that causes Rabi oscillations between the ground and firstexcited state [59].The two soft-Coulomb interacting electrons live in the potential vext(x, t) = −2/ √ x 2 + 1 + Axsin(ωt) where A = 0.0067a.u., and ω = 0.533a.u., resonant with the lowest singlet excitation.The density is shown in red, while the AE potential is in blue.Linked here (Multimedia view) is a movie, directed by Peter Elliott, showing the xc potential evolve in time, over the duration of about one Rabi period, TR = 854.86au.