Exploring non-adiabatic approximations to the exchange-correlation functional of TDDFT

A decomposition of the exact exchange-correlation potential of time-dependent density functional theory into an interaction component and a kinetic component offers a new starting point for non- adiabatic approximations. The components are expressed in terms of the exchange-correlation hole and the difference between the one-body density matrix of the interacting and Kohn-Sham systems, which must be approximated in terms of quantities accessible from the Kohn-Sham evolution. We explore several preliminary approximations, evaluate their fulfillment of known exact conditions, and test their performance on simple model systems for which available exact solutions indicate the significance of going beyond the adiabatic approximation.


I. INTRODUCTION
Much effort has been put into the development of more accurate and reliable approximations to the ground state (gs) exchange-correlation (xc) functional of Density Functional Theory (DFT) [1][2][3].In order to make use of these approximations in the simulation of dynamical processes one invokes the adiabatic approximation: one utilizes a gs xc functional in place of the time-dependent one, and so implicitly assumes the many-body effects in the system are those of a groundstate even if the system is itself not near a ground state.Memory effects arising from the dependence of the xc functional on the history of the density, and on the initial states, are completely absent in these approximations.Adiabatic Time-dependent Density Functional Theory (TDDFT) is widely and successfully used to simulate excitation spectra and response properties in finite systems and solids that are initially in the gs [4][5][6][7] and probe the perturbative regime, where the system is only slightly disturbed from its ground state.The adiabatic approximation becomes less justified as we move away from the gs, as is the case when the perturbation is strong or resonant or when we start the simulation in an excited or arbitrary superposition state.In such situations, which have become of increasing interest in the past decade due to the development of experimental techniques and especially the growing interest in attosecond physics (e.g.Ref. [8]), the adiabatic approximation is at least questionable.How severely the adiabatic propagation deviates from the exact dynamics will depend on the system and the particular dynamics, and also on what the observables of interest are.For example, highly non-perturbative processes such as resonant dynamics and long-range charge transfer have been shown to be particularly challenging for adiabatic TDDFT [9][10][11][12][13][14][15][16][17][18][19][20][21][22] and even in the case of field-free evolution of a non-stationary state the ex-act xc functional develops non-trivial features that influence the dynamics but that adiabatic functionals cannot capture [23][24][25].On the other hand, for example, the adiabatic approximations give good predictions for photo-emission spectra and photo-angular distributions in clusters [26].Given the lack of alternative methods to simulate electron-dynamics in medium to large systems and the expanding range of time-and space-resolved experiments involving non-perturbative dynamics it is imperative to explore the capability of TDDFT beyond the adiabatic approximation.There have been only a few proposals for non-adiabatic xc functionals [27][28][29], but they are not commonly used for a variety of reasons including that they spuriously damp finite systems.Functionals that are explicit orbital-functionals also incorporate memory (e.g.time-dependent exact exchange [30,31]) but the orbital-dependence of these commonly only involves exchange which is inadequate to capture the most significant non-adiabatic effects.
In this work we present a new starting point for developing non-adiabatic approximations based on an exact decomposition of the xc functional.The exact decomposition, which is reviewed in section II, expresses the exact xc potential, in terms of an interaction component that depends on the xc hole of the interacting system, and a kinetic component that depends on the difference between the one-body density matrices of the interacting and Kohn-Sham systems [24,32].The idea is to approximate these terms in terms of quantities that are accessible from the evolving Kohn-Sham system and/or from the initial many-body wavefunction, and some preliminary approximations are proposed in sections II A. In section III we analyze these in terms of the fulfillment of exact conditions and in section IV we test their performance against the exact solution for some model systems.Finally, in Sec.V we conclude.
Eq. ( 1), with mixed separated boundary conditions (i.e.involving the potential and its gradient at each boundary point separately) has the form of a nonhomogeneous Sturm-Liouville problem for v XC (r, t).We can decompose this potential into an interaction (W) and a kinetic (T) component, where the interaction part (W) has an explicit dependence on the xc hole, and the kinetic (T) component, has an explicit dependence on spatial gradients of ∆ρ 1 (r , r, t) = ρ 1 (r , r, t) − ρ 1,S (r , r, t) evaluated at r = r .In Eq. ( 3) we have assumed that the right-hand-side approaches zero as r → ∞, so that no additional term is needed to fulfill the boundary condition.However, such a term, 1 n(r,t) ∇ × a(r, t), is in general needed in Eq. ( 4) to fulfill the Sturm-Liouville boundary condition; we will show in Appendix VI B that without this term, v XC would not be assured to satisfy the generalized translational invariance condition, which it must.There are cases in which this term is zero, including for finite systems in one-dimension; in many cases we will assume this is the case.
An analogous decomposition of the exact groundstate xc potential [33][34][35][36] has provided insight into features of the gs xc potential, while Ref. [24] explored how non-adiabatic features of the time-dependent xc potential manifest themselves in each of these terms.For example, it was found that the dynamical step features found in Refs.[15,23,37] predominantly appear in v T C , while v W XC appears generally much smoother.
In this work we take the first steps to move the decomposition beyond being simply a tool to analyze the exact potential to being a foundation on which to build non-adiabatic functional approximations.By virtue of the Runge-Gross theorem [4] the exact xc potential is a functional of the density n(r, t), the initial many-body state Ψ(r 1 σ 1 ...r N σ N ; 0), and the initial Kohn-Sham state 3) and ( 4) depend explicitly on instantaneous many-body quantities such as the xc hole density n XC (r, r , t) and the 1RDM ρ 1 (r , r, t), as well as the Kohn-Sham 1RDM ρ 1,S (r , r, t), and we do not know the functional form of these in terms of the density and the initial wavefunctions.Operating TDDFT via the time-dependent Kohn-Sham equations offers a partial simplification: one may search for explicit functionals of the time-evolving Kohn-Sham orbitals, which are themselves implicit functionals of n(r, t) and Φ(0), and, by putting more information in the variables that the functional depends on in this way, one hopes that the functional form itself can be simpler.For example, the Kohn-Sham 1RDM is directly accessible in terms of the instantaneous Kohn-Sham orbitals, ρ 1,S (r , r, t) = i,occ.φ * i (r , t)φ i (r, t), and has implicit dependence on n(r, t < t) and Φ(0), including memory.The terms depending on n XC and ρ 1 are however still unknown as functionals of φ i (r, t).
In the next section we propose some preliminary approximations to v W XC and v T C in terms of the instantaneous Kohn-Sham orbitals (orbital functionals) and many-body quantities evaluated at initial time (frozen approximations), that can be used in practical propagation schemes [38].We will assume for simplicity that the boundary term in Eq. ( 4) is zero, which is the case for the one-dimensional systems that we will test our approximations on.
A. Approximations arising from the exact decomposition Here we replace all the quantities from the interacting wavefunction with Kohn-Sham quantities on the RHS of Eq. ( 1), defining a "single-particle" approximation that we denote by v S XC .It is then immediately evident that the kinetic term goes to zero, and such an approximation is an approximation of the interaction component v W XC , making the following replacement in the RHS of Eq. ( 3): This approximation was first presented in Ref. [25], where the shape and size of v S XC was shown to be close to v W XC for a particular dynamics of a model system.This was field-free dynamics of an initial interacting state prepared in a 50:50 superposition of a ground and first excited-state.Ref. [25] found that for a wide range of choices of the initial Kohn-Sham state Φ(0), the approximate potential v S XC remained very close to v W XC throughout the dynamics.The approximation was also tested for time-resolved electron-atom scattering in Refs.[39,40], where it captured the approach of the electron to the target far better than the conventional approximations, although, like them, also ultimately severely underestimated the scattering probability.The integral in Eq. (3) tends to smooth out details of the xc hole somewhat, making v W XC less sensitive to the approximation for the xc hole than v T C is to the approximation made for ∆ρ 1 in Eq. ( 4) as we will see shortly.
From its definition, v S XC is an orbital functional, and thus generally is non-adiabatic in that it depends implicitly on the history of the density and on the Kohn-Sham initial state.When a Slater determinant is chosen for the Kohn-Sham initial state, it is highly plausible that v S XC reduces to time-dependent exact-exchange (TDEXX).For two electrons this follows directly from the formulae, but for more electrons, it is not as straightforward to show [41].For a more general Φ(0), v S XC goes beyond exact-exchange, and contains a portion of correlation [25].In either case, it has spatial-and time-nonlocal dependence on the density, and dependence on the KS initial state, through the dependence on the KS orbitals.
We now address approximations for the kinetic term to be used in conjunction with v S XC .Our first approximation that has a non-zero contribution from the kinetic part of the correlation potential is very simple: we freeze it to its initial value, that is, on the right-hand-side of Eq. ( 4) This is of a somewhat similar spirit to the frozen approximations that are popular in strong-field atomic and molecular physics, such as the single-active electron approximation [42] where it is assumed that only one electron is responsible for the dynamics while all the rest are frozen, such that if the approximation was self-interaction free, the potential that the dynamic electron sees would be static.It is also similar to the "instantaneous ground-state" approximation of Ref. [43] when any externally applied field is static.In addition to the non-adiabatic effects contained in v S XC this approximation also is dependent on the initial interacting and KS states through its dependence on ∆ρ 1 (0).
To "thaw" the frozen kinetic component somewhat, we consider now the next simplest approximation, which is to freeze only the 1RDM-difference in v T C , i.e. v T,∆ρ1(0) while self-consistently updating the density in the denominator on the RHS of Eq. ( 4).Like C inherits a degree of nonadiabatic density-dependence and Kohn-Sham initialstate dependence through the instantaneous orbitaldependence in v S XC , and a degree of both Kohn-Sham and interacting initial-state dependence through the kinetic component.
C Noting that ρ 1,S is accessible during the Kohn-Sham evolution, we define a further frozen approximation where only the interacting 1RDM is fixed to its initial value, while both the Kohn-Sham 1RDM and the density are self-consistently updated during the time-evolution: This approximation freezes the interacting 1RDM while the Kohn-Sham 1RDM may vary greatly.
In the next sections we will discuss the fulfillment of exact conditions and the performance in reproducing some model dynamics for each of the proposed approximations.

III. EXACT CONDITIONS
Exact conditions have been a key instrument guiding the development of physically-inspired and robust gs density functionals that work reliably over a wide range of systems [44,45].Except for the ones arising from the energy-minimization principle which cannot be extrapolated to the time-dependent case, many of the gs conditions should be also fulfilled by the exact td density functional [46] and pose a stringent test for approximations.In addition to these, are conditions that are particular to the time-dependent case, and involve aspects of memory-dependence [5].In this section we analyze the approximations proposed in the previous section in the light of the fulfillment of the exact conditions enumerated below.
(i) Zero Force Theorem (ZFT) The zero force theorem [47][48][49][50] ensures that the xc potential does not exert a net force, Since the net force exerted by the Hartree potential vanishes, Eq. ( 9) ensures that the inter-electron Coulomb interaction doesn't exert any force on the system [51].Linear response (around the gs) applied on Eq. ( 9) yields a link between the xc kernel f XC (r, r , ω) and the gs xc potential v 0 XC (r) [52] that shows that frequency-dependence in a memory-dependent kernel is incompatible with a concurrent local-in-space density-dependence.(The Gross-Kohn functional [53] thus violates ZFT).Violation of ZFT has been shown to lead to numerical instabilities [54] due to the system selfexciting.
(ii) Generalized Translational Invariance (GTI) Translational invariance requires the wavefunction in an accelerated frame (boost) to transform as, e −irj • ḃ(t) |Ψ(r 1 +b(t)...r N +b(t), t) (10) where b(t) is the position of the accelerated observer and b(0) = ḃ(0) = 0 such that accelerated and inertial systems coincide at initial time.The boosted density transports 'rigidly', Vignale proved in Ref. [50] that in order to fulfill GTI the xc potential must transform according to, A v XC that fulfills Eq. ( 12) automatically fulfills the ZFT Eq. ( 9) [50].A special case of GTI is the harmonic potential theorem (HPT), which states that for a system confined by a harmonic potential and subject to a uniform time-dependent electric field, the density transforms rigidly following Eq.( 11) where b(t) is the position of the center of mass [7,47,50,55].(It was shown that Gross-Kohn functional also violates GTI [50]).
(iii) Memory Condition The memory condition states that any instant in time can be regarded as the "initial" moment if the wavefunction is known at that time, and the xc potential should be invariant as to which previous time is used in its functional-dependence [56], i.e.
v XC [n t ; Ψ(t ), Φ(t )]( r, t) is independent of t for t > t , (13) where n t (r, t) = n(r, t) for t > t and is undefined for t < t .This is a very strict condition on the xc potential, tying together initial-state and history dependences, and has implications also when the initial-state is a ground-state.
(iv) Self-interaction Free (SI free) As in the ground-state case, the realization that an electron does not interact with itself yields the condition that for any one-electron system the correlation potential is zero and exchange cancels the Hartree potential, , evolving in a static Hamiltonian, is probed to determine the excitation frequencies, ω i , then the frequency predicted for a given transition should be independent of the state, i.e. of the C n .Within TDDFT this condition requires a subtle cancellation between time-dependence of the xc kernel and the KS response function in a generalized non-equilibrium response formalism [57,58]: the exact (non-equilibrium) Hartree-xc kernel must compensate time-dependence in the poles of the KS response function in order to maintain constant ω i .Given an approximation, this condition is difficult to test analytically; numerically it can be tested by turning off the external field at a given time T and kicking the system.The frequencies of the ensuing dipole moment should be independent of T if they satisfy this condition.Further, the fact that different choices of KS initial state give different peak positions is also an indication of CRC violation [58].

A. Table of results
First, let us make some observations regarding the exact components of the xc potential; these are summarized in Table I.All the proofs may be found in Appendix VI.The exact interaction and kinetic components, Eqs. ( 3) and (4) independently fulfill ZFT, GTI, memory and are SI free.We can further break down the kinetic component v T C into a part v T int that depends on the interacting 1RDM and another v T S that depends on the KS 1RDM: v T int corresponds to replacing ∆ρ 1 → ρ 1 and v T S corresponds to replacing ∆ρ 1 → ρ 1,S in Eq. ( 4).
int and v T S fulfill ZFT independently.GTI is only fulfilled by the difference of the two terms.We therefore expect a strong violation of GTI when the kinetic component is approximated by a term that asymmetrically approximates the interacting and KS 1RDMs (as in the case of v S XC + v T,ρ1(0)

C
).None of the individual components of the exact v XC satisfy independently the CRC in general; the full sum of these terms is generally required.
In Table II we summarize the results for the approximations presented in section II A; the proofs can be found in Appendix VI for all conditions except the CRC.The latter is difficult to prove analytically, so we rely on numerical examples, and in all cases we found one example which explicitly showed peak shifting indicative of the violation of this condition.We separate the analysis of the v S XC component from the kinetic components, considering their satisfaction or violation of the exact conditions independently, although in practise we always evolve a kinetic component in conjunction with v S XC .The potential v S XC fulfills all exact conditions considered here, except for the CRC, and presents the more robust performance as we shall see in section IV, although it does not always give results that are the ones closest to the exact.On the other hand, v T C (0) violates all of them except SI-free.Freezing the difference in the 1RDMs, v T,∆ρ1(0) C fulfills ZFT and is SI free; it turns out to perform quite well for some of the dynamics when used in combination with v S XC as we will see in section IV.However GTI is violated since the system does not transform properly under a boost, but the violation is likely weaker than that for v T,ρ1(0) C : there are large cancellations between v T int and v T S that cannot occur when freezing ρ 1 while evolving ρ S 1 (see discussion later).v C also violates the memory condition and the CRC.On the other hand v C fulfills ZFT but, as mentioned before, violates GTI strongly due to the asymmetric treatment of the interacting and Kohn-Sham 1RDMs.Further, this approximation violates the memory condition and is not SI-free.In table II we also include the adiabatic local approximation (ALDA) and the adiabatically-exact (AE) approximation v AE XC , the later corresponds to propagating with the exact ground-state xc potential evaluated on the instantaneous density, v fulfills GTI and therefore also ZFT and so does ALDA [49,50,59], but ALDA is not SI-free.Because of the lack of dependence on the initial interacting and KS wavefunctions and on the history of the density, both adiabatic approximations trivially fulfill the memory condition [56].But both ALDA and AE violate the CRC as was shown in Refs.[57,58] respectively.

IV. ANALYSIS FOR NON-EQUILIBRIUM DYNAMICS
In this section we analyze the performance of the nonadiabatic approximations introduced in section II A to simulate electron dynamics.We focus on systems driven  well beyond linear response from their gs, i.e. nonequilibrium situations, where we expect memory (nonadiabaticity) to be quite relevant.In general, the performance of an xc functional approximation will likely depend on the system, the particular dynamics under study, and the observables of interest; as evident from the literature, even the simplest approximation, ALDA, has a performance range of abysmal to accurate enough to be very useful [5].This should be borne in mind when a limited range of studies is done, as here, however at the same time, these studies could be carefully used to indicate expected behavior of the functionals in more general situations.We choose a numerically exactly-solvable model system for which we can compute the exact evolution of the many-body wavefunction.We compare not only the performance of the approximation to simulate the density evolution but also the features of the potential as the dynamics evolves.The system is a one-dimensional model of the Helium atom: two soft-Coulomb interacting electrons in a soft-Coulomb well (atomic units are used throughout), where v app (x, t) is an applied potential, mimicking a laser field.All computations were performed using an in-house code.The spatial length of the simulation box is 40 a.u. with spacing 0.1 a.u. and has absorbing boundaries; the time-step used was ≤ 0.02 a.u.

A. Field-free evolution of a superposition state
Our first case study is a superposition state evolving freely, i.e. v app (x, t) = 0. We prepare the He atom in a 50:50 superposition of the ground and first-excited singlet state and we let it evolve freely, where ) denote the many-body eigenstates of Hamiltonian Eq. ( 15) with v app (x, t) = 0.The frequency of the dynamics corresponds to the energy difference between ground and first excited state, Our proposed approximate functionals depend on the choice for the KS initial state Φ(0).Note that this is a separate statement than simply the fact that different initial states yield different dynamics under the same potential: here the potential itself depends on the initial state.The theorems of TDDFT state that any KS initial state may be chosen provided it has the same initial density and initial first time-derivative of the density as the true interacting state.The exact xc potential is different for each choice, but yields the same density-dynamics for each choice.The potentials in an adiabatic approximation are invariant to the choice and yield different density-dynamics for each (and some investigations have been made whether one can judiciously pick an initial Kohn-Sham state for a given physical state in which an adiabatic approximations performs best [25,39,60]).Like the exact xc potential, our approximate potentials are sensitive to the choice of initial KS state, but unlike the exact xc potential, the dynamics differs greatly depending on this choice.
We explore the following choices: 1. single Slater Determinant (SSD): where ϕ(x) = n(x, 0)/2, with n(x, 0) the initial density of the interacting system, 2. non-interacting 50:50 superposition: where Φ 0 is the ground state and Φ 1 is the first noninteracting singlet single excitation with ϕ 0 and ϕ 1 the ground and first excited orbitals of a non-interacting potential such that the one-body density of Φ 50:50 (0), n Φ 50:50 (0 To obtain ϕ 0 (x) and ϕ 1 (x) we use an iterative procedure that targets the initial density.The procedure is close to the ones found in Ref. [61,62].At every iteration the eigenstates ϕ i (x) are found for a guess of the potential v S (x) and a density n(x) is obtained from Φ 0 or Φ 50:50 (0).Then a correction in the potential is added, proportional to n(x) − n target (x).Here we actually use λ(n(x) − n target (x))/(n(x) + ε) where λ and ε are small numerical parameters.The convergence criterium is the difference between the two densities.
We compare the densities and potentials resulting from propagating these initial KS states under our approximations, with the exact ones.
To find the exact time-dependent potential, a variant of the global fixed point iteration method of [63,64] was used.While the method of Ref. [63] has been designed to also have a high accuracy even in the low density regions, which requires a careful treatment of the boundary conditions and low density regions, we are not so much focused on low densities and tail regions in this work.We therefore adopted a more pragmatic approach, introducing a cut-off that flattens the potential in the small-density regions.The cut-off takes the form g where c is the cutoff and α is a large constant.For a potential of the form v(y) = y f (x)dx applying the cut-off corresponds to v(y) = y g(x)f (x)dx.Since the method of [63] in each time-step tries to get the right density, and thus within each time-step to compensate the error in the density caused by for example the cut-off, it causes an overly oscillatory potential in the low density regions.We therefore further limited the number of iterations to three to avoid this issue.As initial guess we used the analytic formulas, Eqs. ( 3) and ( 4), since we also propagate the correlated wave function.(Though analytically exact, these formulas alone do usually not lead to a stable numerical propagation: the numerical error created by the multiple derivatives and integrations lead to a density that drifts away from the reference density).
1. SSD choice for Φ(0) In Fig. 1 we show the dipole d(t) = xn(x, t)dx and norm N (t) = n(x, t)dx for the field-free evolution of superposition state Eq. ( 16) with first KS initial state choice of the previous subsection, the SSD Eq. ( 17).The exact dynamics is shown in black.Propagating using v S XC (in green) which, for two electrons and SSD choice, corresponds to adiabatic exact-exchange (AEXX) yields a dipole that appears to oscillate with more than one frequency, giving an envelope to the oscillations, and neither the dominant frequency nor the amplitude of oscillations are captured well.The ALDA propagation (blue) is also shown, and it approximates the frequency of the oscillation a little better than AEXX, but the amplitude is poor and again there is a beating over long time.The propagations that include the various approximate v T FIG.1: Field-free dynamics of the 50:50 superposition state Eq. ( 16) where the approximate TDDFT calculations propagate the SSD choice for initial KS state Φ(0), Eq. (17).Norm (upper panel) and dipole moment (lower panel) for exact (black C (light blue) and v ALDA XC (blue).Time is given in a.u.here and in all following graphs.improve over the traditional approximations for only a very short time but soon deviate dramatically with qualitatively wrong behavior even before half a cycle of the evolution.We will now discuss what goes wrong with each of these, with the help of the density and potential time-snapshots within the first optical cycle shown in Fig. 2.
We first note that when a SSD is chosen for the KS system, with one doubly-occupied orbital, the features of the exact v XC (in black) have the periodicity of the density dynamics: initially the density has a minimum on the left side coincident with a large dynamical step and peak in v T C , which decreases and disappears as the density becomes symmetric around t = 3a.u., before building up again on the other side, such that at half-cycle (t ≈ 6 a.u.) v XC and the density are a mirror images of their initial values.
Consider now the potentials of the traditional approximations, v S XC =AEXX and ALDA (in green and blue in Figs.1-2).Although both the AEXX and ALDA dipoles are close to the exact for t < 3a.u. the shape of the density at this time has already started to deviate.The v S XC potential is generally a very good approximation to the exact interaction part v W XC of the xc potential [25] but it completely misses the step and peak features present in the exact v T C , simply smoothly cradling the density as it evolves, and ALDA also does not capture these structures.These missing structures are important to get the details of the dynamics right, including the period, as evident in the plots.
C potential is scaled by a factor of 10 −5 at t = 6a.u. and t = 9a.u..

Consider now propagation with v S
XC + v T C (0) (magenta curves in Figs.1-2).As v T C (0) is static, the large initial step in the kinetic correlation part of the potential remains for all time, and as a consequence the density leaks to the left and is unable to slosh back to the right (see magenta density in left panel Fig. 2).The density eventually gets absorbed by the boundary around t ≈ 30a.u. as can be seen in the evolution of the norm in the upper panel of Fig. 1.

We turn now to propagation under
C , where again the dipole swings far too much to the left (Fig. 1), resulting in a large unphysical absorption evident in the norm in the top panel.The reason for this becomes evident from the right panel of Fig. 2: Freez- S (in black) for field-free dynamics evolution of superposition state Eq. ( 16) at initial time and at about 1/4 cycle (t = 3 a.u.).Top panel: SSD choice Eq. ( 17) for KS initial state.Lower panel: Superposition choice Eq. ( 18) for the KS initial state.
ing ∆ρ 1 in v T C results in a step that, although varying in size, is always on the left, pulling the density consistently over to the left.

On the other hand, propagation under v
C yields a dipole that displays a most peculiar flattening after half a cycle (light blue in Figs. 1 and 2).The freezing of the interacting 1RDM while keeping the KS 1RDM dynamical has drastic consequences, as can be seen from first recalling the exact v T C : Both the exact v T int and v T S display large dynamical steps that tend to counteract each other, yielding a smaller step in the exact v T C , at times even canceling each other.(see Fig. 3 and movie 1 in supplementary material).Freezing only the interacting 1RDM as in v T,ρ1(0) C leads to the development of a disproportionately large step at times, resulting in instabilities and the poor dynamics as shown here.
Satisfaction of the ZFT is demonstrated in Fig. 4 which plots the left-hand-side of Eq. ( 9) for the different approximations.From Table II, all approximations except for v S XC + v T C (0) are expected to satisfy the ZFT, and give zero.The deviation of v S XC + v C is a numerical error related to the very large steps that appear in this approximation as just discussed.The spike in v S XC + v C at around 25 a.u. is a numerical artifact caused by the step in the potential at that time being very large and sharp.Notice the relatively large violation of the ZFT for v S XC + v T C (0) in this dynamics.We conclude that except for quite short times, the new approximations do not improve the dynamics of a fieldfree superposition state when the initial KS state is cho- FIG.4: Violation of ZFT: The left-hand-side of Eq. ( 9) in field-free evolution of superposition state Eq. ( 16), in propagation under v S XC (green C (light blue) and v ALDA XC (blue).Upper panel: SSD choice Eq. ( 17) for KS initial state.Lower panel: Superposition choice Eq. (18) for the KS initial state.
sen to be a SSD; in fact the traditional approximations ALDA and AEXX (coinciding with v S XC in this case) are more reliable, although inaccurate.
One might be tempted to argue that these results have limited relevance, since starting the KS system in an SSD when the interacting state is a 50:50 superposition state would be a bad choice, making a challenging job for approximate functionals to do well.However, this could well be the situation reached when the system began in the ground-state, with some field applied that brought the interacting state to a 50:50 superposition state, at which point the field is turned off.Then the natural choice for the initial KS state is the KS ground-state, and we are stuck with a SSD for the whole evolution.

Superposition choice for Φ(0)
Fig. 5 shows the field-free superposition state dynamics for the different approximations when the KS initial state is chosen to be Eq.(18).Considering first the traditional approximations ALDA and AEXX = −v H /2, we observe a great improvement in their performance with this choice of KS initial state, compared to the Slater determinant choice.In fact ALDA is remarkably good: choosing a KS initial state structure that is close to that of the true system seems to override the errors from the ground-state xc assumption inherent in the adiabatic nature of ALDA, its self-interaction error, and its locality in space.It is the configuration of the initial KS state that is key: if instead we used LDA orbitals but still within the  16).Initial KS state Φ(0) is chosen as in Eq. ( 18).Shown are results of propagation under the exact (black), v S XC (green C (light blue), AEXX= − vH 2 (yellow) and ALDA (blue).
first-excited Slater determinant configuration of Eq. ( 18), the dynamics is similar although not as accurate (see shortly).In fact ALDA outperforms the spatially nonlocal AEXX, suggesting the importance of an adequate accounting of correlation for this example, with even the local-density ground-state correlation approximation doing a good job here.We will see that freezing the kinetic component of the correlation in various ways in our new approximations, together with the correlation contained in v S XC , do not provide a sufficiently good correlation component to compete with the simple ALDA for more than short to intermediate times.
Turning then to the new approximations, we first note that, unlike the SSD choice, none of the approximations lead to the unphysical situation of the density reaching the boundaries and being absorbed, at least until much later (for v S XC + v T C (0)), as shown by the norm.Two of the approximations, namely v S XC (which gives dynamics close to AEXX) and v S XC + v C do a reasonable job however none of them beats the remarkable performance of ALDA for the dipole dynamics as noted above.Inclusion of v T,∆ρ1(0) C tends to yield a worse amplitude but a better dominant frequency than v S XC .Performing a Fourier transform of the dipole yields the frequencies, ω = 0.515 ± 0.006 a.u.for v S XC + v C as compared with ω = 0.584 ± 0.006 a.u.for v S XC while the exact frequency is ω 0 = 0.534 a.u, best approximated by ALDA's frequency of ω = 0.527 ± 0.006a.u.
These two new approximations, v S XC and v S XC + v T,∆ρ1(0) C yield a good density evolution, as can be seen in Fig. 6 in green and red respectively, where they can be compared with the exact density and potential (black).The exact v T C for this choice of KS initial state is initially much smaller than for the SSD choice (c.f.Fig. 2), and v S XC , which, as we noted before is a very good approximation to v W XC , then plays a relatively more important role than in the SSD case.We note that v S XC in this case is not equivalent to TDEXX: it contains some correlation due to the KS state configuration, and has memorydependence in that it is an implicit functional of the density at earlier times and of the KS initial state.Including v T,∆ρ1(0) C now partially captures the step and peak features of the exact potential and the step "switches sides" unlike for the SSD case.Notice that from approximately t = 6a.u.onwards the exact v XC is no longer periodic while the density is (unlike the case in the previous section where a SSD is used).The dynamical features of the exact v T C become more and more complex as the system evolves and at later times they become eventually as large as the ones observed for the SSD choice (see also Ref. [25]); the movie 2 in supplementary materials displays the densities and xc potentials from exact propagation, ALDA, v S XC , and C .Turning now to the other two approximations, C , the performances are poor, although the release of density then absorbed by the boundaries is less compared to the SSD case.This is because the step in the initial kinetic components is less, so even if frozen as in v S XC + v T C (0), it results in less density moving out to the left.Propagation with v S XC + v C suffers from the same problem as in the SSD choice: the uneven treatment of the interacting and KS density matrices results in a huge step in the potential that prevents the density from sloshing back and forth as it should (right panel of Fig. 6).Finally, since typically orbitals that reproduce the exact density are not easily obtained, we show in Fig. 7 the results of propagating a state of the form Eq. ( 18) where the LDA ground and first-excited orbitals are inserted into the KS states.We show the best of our new approximations, along with ALDA.We notice that the densitydynamics is generally significantly worsened for both propagation with ALDA and v S XC + v C compared to when initial orbitals that reproduce the exact initial density are used.

B. Arbitrary field starting in the ground state
Next we study the dynamics of system Eq.( 15) initially in the gs and driven by a non-resonant, fairly strong, laser field, v app (x, t) = E(t)x, where E(t) = 0.1 sin(0.4t).Since we start in the gs the natural choice is a SSD for the KS initial state, and this is the only choice we consider here.In Fig. 8, the dipoles and norm for the exact and approximate KS evolutions, beginning in the initially exact KS orbital ϕ(0) = n0 2 are shown.x (a.u.) FIG. 6: Time snap-shots of the density (upper panels) and v XC potential (lower panels) for the dynamics shown in Fig. 5 for: exact (black), v S XC (green), C (red), and C (light blue).The latter has been scaled by 10 −4 for t = 9au and t = 15au.
The propagation with v S XC =AEXX alone (in green) is excellent till almost t = 12 a.u., after which it begins to significantly deviate from the exact.The other traditional approximation ALDA (in blue) is not as good as AEXX in this case for the first few cycles but stays closer to the exact than AEXX at later times.Fig. 9 shows that in fact both these approximations and the new ones, except for v S XC + v C which is not shown (see below), capture the structure of the density and the potential quite well in the central region; the noticeable deviations of the dipole moments of the approximations in Fig. 8 arise from errors in the density further out where the density itself is not very large, but these errors get enhanced by the multiplication by x in the calculation of the dipole moment (which, additionally, highlights asymmetric structures in the density).
Adding the correlation potential frozen to its initial value, v S XC + v T C (0) (in magenta) does not improve the  16) and using KS initial state of the form Eq. ( 18) but with the gs and first excited LDA orbitals.Propagation with C (red) and v ALDA XC (blue).In black the exact dynamics.results over v S XC (except at short times) as can be seen by comparing the green and magenta curves in Fig. 8 and in the density time snapshots of Fig. 9.In this case the system starts off weakly correlated, and the exact potential is a smooth well initially, but as the laser field kicks in and begins to drive the interacting system far away from approximately a single Slater determinant, correlation increases, as is reflected in the increased size of the structures in the exact v T C evident in Fig. 9; these structures appear in v T C rather than in v W XC which remains quite smooth.As v S XC + v T C (0) freezes the v T C to its initial value, it cannot capture these.
The approximation v C does include partial step structures, although not always quite in phase with the exact (see also the movie 3 in the supplementary material).However, these do not apparently have a significant effect on the ensuing dynamics, neither in the region where the density is appreciable, nor in the dipole moment that weighs more the low-density regions further out, and the dipole moment is in fact worse than the traditional approximations except at short times.One aspect of the largest of these step structures for this dynamics, is that they can change very rapidly, for example the large step switches sides from t = 43 a.u. to 44 a.u.Also, at t ≈ 60 a.u. the exact dynamics begins to show some absorption (see upper panel Fig. 8) which is underestimated in C .This leads to some discrepancy in the dipole moments, but the density in the central region is very well reproduced still (see lower panel in red in Fig. 9).
The propagation under v S XC + v C develops very large features near the boundaries (no compensation of the large step in v S C , as in the discussion in the previous section) that progressively approach the central region, growing in size, and by t ≈ 5 a.u. the step has become so large and sharp, that the propagation results in a freezing of the density on the left (in light blue in Fig. 8) .Finally, in practice we rarely have access to the exact initial KS orbitals, and approximate DFT orbitals are used for the initial state.In Fig. 10 LDA orbitals are propagated using v S XC + v C and ALDA, showing a noticeable deviation from the results when the exact initial orbitals are used, especially, in this case for the propagation under v S XC + v C .

V. CONCLUSIONS
Building approximations based on the exact expression for the time-dependent xc potential, Eq. ( 1), allows us to immediately break free of the adiabatic approximation in TDDFT.The expression requires approximations for the time-evolving xc hole and for the 1RDM of the interacting system, quantities that are not accessible in a KS calculation.In this work, we have proposed a number of practical approximations to these terms, evaluated their satisfaction of some known exact constraints in TDDFT, and explored their performance on a few model systems.The approximations have memorydependence, in particular initial-state dependence, and how well they perform in turn depends on the choice of KS state in which to begin the propagation.The functionals can be viewed as orbital-functionals, however they do not require TDOEP in their operation, and have an explicit dependence on the initial interacting manybody state and the initial choice of KS state.Two of these newly proposed approximations stand out.One is v S XC in which all quantities in the exact expression for v XC are replaced by their KS counterparts.This approximation has been explored in recent literature where it was found to be an excellent approximation to the interaction component v W XC of the exact potential [25,39] and to yield a significant improvement in dynamics of an electron approaching a target in scattering compared to ALDA and AEXX [39,40] albeit ultimately failing to actually scatter even qualitatively.In the examples studied here, it provides a reasonable approximation, although ALDA was in fact generally better.A dynamical accounting of correlation beyond what is present in v S XC is required.The approximation satisfies the most exact conditions out of the new approxi- mations considered, including ZFT, GTI, SI-free, and the memory condition.
The second approximation that worked reasonably well was v S XC + v T,∆ρ1(0)

C
, where the kinetic component is approximated by freezing the difference in the 1RDMs to the initial value.In the examples considered here, it did not perform particularly better than the traditional approximations; a better accounting of correlation is required.The approximation satisfies the ZFT and is SIfree, but violates the other exact conditions we considered.A key factor in determining the accuracy of the performance of any approximation, is the choice of the KS initial state; the field-free propagation of the superposition state demonstrated clearly the vast improvement of all approximations once the KS initial state was chosen with a configuration close to that of the interacting state.
It is clear that better approximations are needed for the kinetic term, that involve a dynamical approximation to the 1RDM of the interacting system, and such developments are underway [65].This term is more challenging as it has much more structure and a stronger non-adiabatic dependence than the interaction term [24].How well these approximations work for realistic many-electron three-dimensional systems remain to be seen.For robust performance, we expect that certainly the approximation should satisfy the ZFT; satisfaction of GTI, SI-free, memory-condition, and CRC will likely improve its accuracy.As in the ground-state case, it is not always clear which exact conditions are most important for predicting a functional approximation's suc-cess.The approach pointed to in this work, based on approximations to the 1RDM-difference and xc hole appearing in Eq. ( 1), provides a starting point to develop such practical memory-dependent approximations.r.To obtain the third line, the definition of the 1RDM in terms of the wavefunction is inserted, and we have noted that ∇ = ∇ r acts only on Ψ * (r , r 2 ..r N ) while ∇ = ∇ r acts only on Ψ(r, r 2 ..r N ).Noting that the integrand in the final line is in fact a total gradient, we see that upon performing the integral over r, the integral vanishes for finite or periodic systems where the integrand vanishes at infinity or is identical at either end of the periodic system.Hence, the ZFT is satisfied by v T int .An identical argument can be made using KS wavefunctions, to show that v T S also satisfies the ZFT.We conclude that both n(r, t)∇v T int (r, t)d 3 r and n(r, t)∇v T S (r, t)d 3 r vanish independently, and therefore v T C fulfills ZFT independently of v W XC .

Approximations
C That v S XC satisfies the ZFT follows from the same argument as applied to v W XC in the previous section, noting that the KS pair density P S (r, r , t) obeys the same symmetry relations as for the true pair density.
Turning to the frozen kinetic components, consider first v T C (0), for which At the initial time, the right-hand-side is zero, as follows from the previous section, but at later times, due to the ratio of the densities factor, the arguments of the previous section cannot be made, and there is no guarantee that the integral will be zero.In fact it is generally non-zero, and therefore v T C (0) violates the ZFT.This is also evident in the numerics as can be seen in Fig. 4 and the violation leads to numerical instabilities from "selfexcitation"of the system [54].
Fortunately, the other two frozen approximations do satisfy the ZFT.By keeping the density in the denominator of Eq. ( 4) dynamical and equal to the evolving density, the density factor in the zero force expression n(r, t)∇v C d 3 r is exactly cancelled and the argument follows closely to that of Eq. ( 21) and subsequent discussion in the previous subsection.Likewise, the proof there that the two components of v T C independently satisfy the ZFT, ensures that v T,∆ρ1(0) To study the translational invariance of the kinetic correlation potential we again separate it into interacting and KS contributions, ∇v T C = ∇v T int − ∇v T S , and study the behavior of each term separately under the boost.That is, where Working through the derivatives eventually leads to where r b (t) ≡ r + b(t).
Eq. (26) shows that ∇v T int does not satisfy GTI by itself, since the last three terms on the right are not generally zero.However, once we subtract the corresponding equation for ∇v T S , noting that the KS system has the same density as the interacting system, we find ) where j S (r, t) is the current-density of the KS system.In 1D, the second term vanishes as the KS and interacting currents are the same, but in 3D they are generally not the same.Yet, we know that the exact full v XC satisfies GTI, and also that the exact v W XC does, and therefore v T C = v XC − v W XC must also.The resolution lies in the additional term ∇×a(r,t) n(r,t) in Eq. ( 4) that we assumed zero throughout most of this paper: in 3D this term is not generally zero, and it plays a crucial role here in restoring GTI.
In summary, although ∇v T int and ∇v T S separately are in general not translationally invariant, their difference, does satisfy GTI.

Approximations
However none of the frozen approximations to v T C satisfy the GTI.In v T C (0) the entire kinetic potential remains constant in time, so does not transport rigidly under a boost.Likewise it can be seen that freezing ∆ρ 1 in v T,∆ρ1(0) C or ρ 1 alone in v T,ρ1(0) C will not transform correctly under a boost.

C. Memory Condition Satisfaction
Adiabatic xc functionals trivially fulfill the memory condition, Eq. ( 13) because of the lack of dependence on the initial interacting and KS states Ψ(0) and Φ(0).When we introduce partial memory as in the approximations proposed here, this condition becomes susceptible to being broken.which is independent of t < t, because it depends explicitly on the exact xc hole, n XC (r, r , t), at time t, which, although an implicit functional of the history of the density and the initial interacting state, it is an explicit and instantaneous functional of Ψ(t), from which it is obtained directly.Likewise, the memory condition is fulfilled by the two components of the exact correlation potential independently: v T int explicitly depends on the instantaneous interacting 1RDM, an implicit functional of the history of the density and the initial interacting state, but is obtained explicitly from the instantaneous interacting state, and so again is independent of t < t.An analogous argument can be made for v T S , and so the kinetic component of the correlation potential satisfies the memory condition, v T C [∆ρ 1 [n t , Ψ(t ), Φ(t )](t)](r, t) because it is independent of t < t.Note that the boundary term (∇×a(r, t))/n(r, t) does not affect this analysis, since it depends on instantaneous quantities.In the case of v S XC the proof follows a similar path as for the exact interaction component Eq. (30).The potential v S XC is an explicit functional of the KS xc hole n S XC (t) which is obtained explicitly from the instantaneous KS wavefunction Φ(t), hence v S XC [n t , Ψ(t ), Φ(t )](rt) = v S XC [n S XC [Φ(t)](t)](r, t) (33) is independent of t < t.
In the case of the frozen approximations to the kinetic component, we first note that these approximations are defined with freezing parts of the kinetic component to their values at the time that is considered the initial time, i.e. the time t , in v XC [n t , Ψ(t ), Φ(t )](rt).Therefore, freezing the entire kinetic potential as in v T C (0) means freezing it to v T C (t ).This clearly is t -dependent, and so violates the memory condition.The approximation v T,∆ρ1(0) C (t) also violates the condition, since since for N = 1, Ψ(t) = Φ(t) and thus ρ 1 (t) = ρ 1,S (t).
XC results if we make the exact same approximation, Φ = Ψ = Slater determinant in this equation.The latter corresponds to v Loc X in Ref. [9], which the authors compare against exchange-only-KLI for some particular dynamics and find the two yield similar results.It is extremely difficult to show directly from the formulae that the two expressions are equal [67], but the argument above suggests it is likely that they are.

FIG. 2 :
FIG. 2: Density (upper panel) and v XC potential (lower panel) for the same dynamics as Fig.1 shown at initial time, and at three other times indicated within the first cycle.Left panel: exact (black), v S XC =AEXX (green), v S XC + v T C (0) (magenta) and v ALDA XC (blue).Right panel: exact (black), v S XC + v T,∆ρ1(0) C

FIG. 3 :
FIG. 3: Exactcomponents v T int (in red) and −v T S (in blue) of the exact kinetic component v T C = v T int − v T S (in black)for field-free dynamics evolution of superposition state Eq.(16) at initial time and at about 1/4 cycle (t = 3 a.u.).Top panel: SSD choice Eq. (17) for KS initial state.Lower panel: Superposition choice Eq. (18) for the KS initial state.

FIG. 7 :
FIG.7: Time snap-shots of the density (upper panel) and v XC potential (lower panel) for field-free dynamics starting in 50:50 superposition state Eq.(16) and using KS initial state of the form Eq. (18) but with the gs and first excited LDA orbitals.Propagation with v S XC + v

FIG. 8 :
FIG. 8: Dynamics of the laser-driven gs.Norm (upper panel) and dipole (lower panel) for exact (black), v S XC

FIG. 10 :
FIG. 10: Time snapshots of the density (upper panels) and v XC (lower panels) for laser-driven dynamics starting in the LDA ground-state: propagation under exact (black), ALDA (blue), v S XC + v T,∆ρ1(0) C

1 .
Exact components: v W XC , v T C , v T int , and v T SFrom the expression, Eq. (12), for the boosted wavefunction discussed in section III it follows that the xc hole transforms asn b XC (r 1 , r 2 , t) = n XC (r 1 + b(t),r 2 + b(t), t) under a boost.Given the invariance of w(|r 1 − r 2 |) under such a transformation we conclude that the ∇vT,∆ρ1(0) XC [n t , Ψ(t ), Φ(t )](t) = D∆ρ 1 (r , r, t )| r =r 4n(r, t) (34) where D = (∇ − ∇) ∇ 2 − ∇ 2 will not in general be independent of t .Likewise v T,ρ1(0) C (t) will violate it in general, due to the t -dependence of the interacting 1RDM in this approximation.D. Self-interaction1.Exact components: v W XC , v T C , v T int , and v T SThe exact interaction component cancels the Hartree component for one electron as it should,∇v W XC (r, t) = − dr n(r , t)∇w(|r − r |) = −v H (r, t)(35) since n XC (r , r, t) = −n(r , t) for N = 1.The two exact components of the kinetic potential are not independently SI-free, v T int = 0 and v T S = 0 for N = 1 (36) but it is the subtraction of the two what gives the exact kinetic correlation component, which properly vanishes for one electron, v T C (r, t) = 0 (N = 1)

TABLE II :
Fulfillment of exact conditions by the proposed approximations to the interaction and kinetic components.Also included are the adiabatically-exact (AE) and ALDA approximations.
XC fulfills GTI can be carried over for the single-particle approximation v S XC since the KS xc hole transforms properly under a boost, namely n S,b XC (r 1 , r 2 , t) = n S XC (r 1 + b(t), r 2 + b(t), t), v S,b XC [n S XC ](r, t) = v S XC [n S XC ](r + b(t), t).
T,ρ 1 (0)CThe same argument used in previous section to prove the exact v W