Time-dependent density-functional theory for strongly interacting electrons

We consider an analytically solvable model of two interacting electrons that allows for the calculation of the exact exchange-correlation kernel of time-dependent density functional theory. This kernel, as well as the corresponding density response function, is studied in the limit of large repulsive interactions between the electrons and we give analytical results for these quantities as an asymptotic expansion in powers of the square root of the interaction strength. We find that in the strong interaction limit the three leading terms in the expansion of the kernel act instantaneously while memory terms only appear in the next orders. We further derive an alternative expansion for the kernel in the strong interaction limit on the basis of the theory developed in [Phys. Chem. Chem. Phys. 18, 21092 (2016)] using the formalism of strictly correlated electrons in the adiabatic approximation. We find that the first two leading terms in this series, corresponding to the strictly correlated limit and its zero-point vibration correction, coincide with the two leading terms of the exact expansion. We finally analyze the spatial non-locality of these terms and show when the adiabatic approximation breaks down. The ability to reproduce the exact kernel in the strong interaction limit indicates that the adiabatic strictly correlated electron formalism is useful for studying the density response and excitation properties of other systems with strong electronic interactions.


I. INTRODUCTION
Time-dependent density functional theory (TDDFT) [1][2][3][4][5][6][7] is a well-established approach to study the time-dependent and excitation properties of many-electron systems. One of the main reasons for its popularity is that within this formalism the time-dependent interacting many-body problem can be recast exactly into an equivalent one-particle framework, which is advantageous for numerical implementations. The corresponding one-particle equations, called the time-dependent Kohn-Sham equations, contain an effective potential, known as the Kohn-Sham potential, which is defined in such a way that the noninteracting system has the same time-dependent density as the original interacting many-body system. The Kohn-Sham potential is typically written as the sum of the external potential of the interacting system and of the Hartree and the exchange-correlation (xc) potential.
In practical applications of TDDFT, the xc potential v xc is approximated. The type of approximation employed crucially determines the quality of the results, and therefore a considerable amount of research has gone into the difficult task of finding reliable and accurate approximations for this quantity. This task is simpler in the linear response regime where we consider small variations in the density caused by applied perturbations. This regime is of interest as the knowledge of the density response function is sufficient to calculate the excitation energies and the absorption spectrum of the system [8]. For this purpose, it is enough to know v xc and its functional derivative δv xc /δn = f xc with respect to the density n, evaluated at the ground-state density. The quantity f xc is called the xc kernel and has been the subject of intense investigations.
The simplest possible approximation is the adiabatic local-density approximation (ALDA) in which the xc kernel is local in space and time. The ALDA has, however, a number of deficiencies [4] such as, for example, the inability to produce correct charge transfer excitations [9][10][11][12], Born-Oppenheimer surfaces of excited states in dissociating molecules [13][14][15] and semiconductor band gaps [16]. Some improvements have been made using hybrid functionals which contain mixtures of exact exchange and traditional local functionals. These methods are nonlocal in space but still adiabatic. However, they are not systematic and the optimal mixture of exact exchange is often system dependent [17][18][19][20]. Other, more systematic, approximations for f xc beyond the ALDA often rely on perturbative expansions [21][22][23][24][25][26][27][28][29][30] and many of them are restricted to the exchange-only approximation. Their perturbative nature makes these approaches questionable in the strong correlation regime which is relevant in various physical situations, notably the case of molecular dissociation, and hence it is highly desirable to develop new techniques to tackle this regime.
In a recent work [31], the so-called strictly correlated electrons (SCE) framework [32][33][34], a formalism well suited for the description of strong interactions, has been applied within the time-dependent domain. The authors derived an expression for the xc kernel in the so-called adiabatic approximation and established that the adiabatic SCE (ASCE) kernel satisfies the zero-force theorem [35], an exact property related to generalized translational invariance [36,37]. The kernel was furthermore studied for finite one-dimensional systems with different density profiles, some of which are prototypical of the dissociation of two-electron homonuclear molecules. It was found that the ASCE kernel is spatially nonlocal and exhibits a divergent behavior as the molecular bond is stretched. For adiabatic kernels, this diverging behavior is crucial [13] for describing bond-breaking excitations, which is a notoriously challenging problem in linear response TDDFT. Since the kernel was derived in the adiabatic approximation, not much could be concluded about the limitations of its adiabatic nature in the context of the physics of strong static correlation. The case of infinitely strong electron-electron correlation is quite peculiar and to date it is not known how accurate the adiabatic approximation can be in such a regime.
One way to shed light on this issue is to benchmark the ASCE kernel against an exact expression for f xc obtained from a model system where the density response function and thus the xc kernel can be calculated analytically. In this work, we consider such a model, namely two interacting electrons on a quantum ring [6,[38][39][40], for which we not only compute the exact density response function and xc kernel but also obtain these quantities for various two-body interaction strengths, including the infinitely strong one, and we compare these results with those given by the SCE theory in the adiabatic approximation. The leading order of the asymptotic expansion for the exact f xc and the expression for the ASCE kernel are found to be identical. We also derive the next-order correction term beyond the ASCE, called the adiabatic zero-point-energy approximation (AZPE), and show that also this adiabatic term is the same as the next order from the asymptotic expansion. The third order in the expansion for the kernel is still adiabatic, while a frequency dependence appears in the fourth order. In this order, an adiabatic approximation would break down. This is one of the central results of the paper and elucidates both the strengths and the weaknesses of the adiabatic approximation in the limit of strong electron-electron interaction.
The paper is organized as follows. In Sec. II we introduce the quantum ring model and, after computing the full spectrum of its Hamiltonian, we discuss asymptotic expansions for its eigenenergies and eigenstates in the case of strong interactions and analyze them. In Sec. III we study the density response function of strongly interacting systems, while in Sec. IV we focus on TDDFT in the same regime and give an asymptotic expansion for the xc kernel. Our conclusions are finally presented in Sec. V.

A. Two interacting electrons on a quantum ring
For our study of electron correlations, we consider an analytically solvable model, which we will refer to as the quantum ring model, of two electrons on a ring of length L which repel each other with a two-body interaction. The interaction strength can be adjusted using a parameter, which allows us to study the exact properties of the system ranging from weak to very strong interactions. The explicit form of the Hamiltonian of the quantum ring is given bŷ where λ 0 is a dimensionless parameter and V 0 has units of energy. The coordinates x 1 and x 2 are the coordinates of the electrons on the ring which run from 0 to L. The groundstate density n 0 = 2/L is spatially constant and independent of λ; for this reason the model can be used to illustrate several features of the coupling strength dependence in density functional theory. In order to calculate the properties of the system we have to determine the eigenfunctions which satisfy the stationary Schrödinger equationĤ = E , where E are the energy eigenvalues. For our two-particle system, these eigenfunctions can be written as a product of a spatial wave function and a spin function as follows: For the singlet case (which we will focus on), the normalized spin function is given by and is antisymmetric in the spin variables. For the triplet, there are three linearly independent symmetric spin functions which we, for simplicity, all denote by − (σ 1 ,σ 2 ). Since the two-electron wave function is antisymmetric under the simultaneous interchange of space and spin variables, it follows that the spatial wave functions ψ ± satisfy the symmetry relation Apart from these symmetry conditions, the Schrödinger equation needs to be solved with periodic boundary conditions on the variables x 1 and x 2 ; i.e., the wave function and its first spatial derivatives are invariant under the substitution x i → x i + L for i = 1,2. Using these conditions, we can solve the Schrödinger equation by a suitable coordinate transformation. Since these steps are carried out in detail in Ref. [41], here we just outline the main steps relevant for this work.
The Hamiltonian (1) becomes separable in the terms R = (x 1 + x 2 )/2, the center-of-mass coordinate, and z = π (x 1 − x 2 )/L, the dimensionless relative coordinate. This variable transformation giveŝ By inserting a product ansatz of the form ψ(R,z) = f (R)M(z) into the Schrödinger equation, we find that the spatial twoparticle eigenfunctions are of the form where k is an integer and the function M(z) satisfies the Mathieu equation, which we write in its standard form as [42,43] where the constants q and a are given by For a given value of q, the Mathieu equation (3)  the Mathieu cosines, this label starts at l = 0 and for the Mathieu sines at l = 1. The even Mathieu cosine function is denoted by C l (z; q) and its characteristic value by a + l (q), while the odd Mathieu sine function is denoted by S l (z; q) and its characteristic value by a − l (q). Since the center-of-mass wave functions are symmetric under the interchange of the spatial coordinates of the electrons, we see from Eq. (2) that the singlet wave functions must be described by even Mathieu functions whereas the triplet ones must be described by odd Mathieu functions. The final form of the normalized singlet and triplet wave functions therefore is in which the normalization of the Mathieu functions is chosen such that The Mathieu functions further have the periodicity property M l (z + π ) = (−1) l M l (z); i.e., they are periodic in π for even values of l and antiperiodic for odd values of l. Furthermore, the center-of-mass wave function in Eq. (2) changes with a prefactor (−1) k when x i → x i + L, for i = 1,2. Therefore, for the wave functions in Eqs. (6) and (7) to satisfy periodic boundary conditions, the labels k and l must be both even or both odd. Note that k runs over all integers while l only runs over the non-negative integers. From Eq. (5) we see that the energy eigenvalues are given by Since in the subsequent discussion of the density response function we focus on the singlet excitations in particular, we write the singlet wave functions in a slightly different form for the purpose of a better interpretation. By multiplying the spatial wave function of Eq. (6) with its singlet spin function, the full space-spin function can be written as where we defined the Slater determinant and we further defined the spatial normalized orbital by φ k (x) = e 2πikx/L / √ L which corresponds to a periodic singleparticle wave function of a free particle on the quantum ring. Let us consider the excitation from the ground state to another singlet state with l = 0, which requires that the excited state is characterized by an even k value. In that case, φ k/2 in the Slater determinant above is a proper periodic wave function as k/2 is an integer. According to Eq. (9) the excitation energy is which is independent of the interaction strength q as the excited state has the same relative wave function as the ground state. For the case that q = 0, we have C 0 (z; q = 0) = 1/ √ 2 and the ground and excited states both become pure Slater determinants. The excitation then represents a promotion of two electrons from a doubly occupied k = 0 state to a doubly occupied state with a one-particle quantum number k/2, which is commonly called a double excitation. When q is nonzero, this language is not accurate anymore as also the relative wave function becomes relevant. If the interaction strength becomes very large, the energy required to excite to a state with nonzero l becomes very large too and the excitations with energy E + k0 give the dominant contribution to the density response function, as we will see later.

B. The strong interaction expansion of the exact solution
As the interaction strength q increases, the electronic repulsion becomes more important and the electrons tend to stay in opposite positions on the ring. This physically intuitive picture can be analyzed in more detail using the Mathieu equation. According to Eqs. (6) and (7), the square of the spatial wave function is given by where M l (z; q) is either a Mathieu cosine C l or a Mathieu sine S l depending on whether the wave function is a singlet or a triplet one. We therefore see that the probability to find a given electron at x 2 given an electron at x 1 only depends on the relative coordinate x 1 − x 2 , as one would expect on the basis of the symmetry of the system. This probability distribution is given by the square of the Mathieu function M l . Let us analyze the properties of this function in the large interaction limit which according to Eq. (3) satisfies a singleparticle Schrödinger type of equation in a potential of the form V (z) = 2q cos(2z). For large values of q, we can see that the relative wave function described by M l becomes localized in the minimum of the potential at z = π/2, which corresponds to a relative distance of the particles of L/2. We can expand the potential around this minimum to obtain This potential describes (apart from a shift of the minimum) a harmonic oscillator with frequency = 2 √ q. The eigenfunctions of the harmonic oscillator are well known to consist of Gaussians of width proportional to 1/ √ = 1/( √ 2 q 1/4 ). In the limit of large q, the harmonic frequency increases and the wave functions become localized around z = π/2. This behavior is illustrated in Fig. 1. The eigenenergies l of the harmonic oscillator are well known and given by l = −q + (l + 1/2) = a/2. This also immediately provides an asymptotic formula for the characteristic value of the Mathieu equation for large values of q: and consequently also an asymptotic expansion for the eigenenergies of the quantum ring from Eq. (9).
FIG. 1. The squared ground-state wave function |ψ 00 | 2 for two values of the interaction strength q, and the interaction cos 2 z using a suitable scaling for showing it in the same plot. For large q, the wave function localizes around z = π/2 where cos 2 z is almost parabolic and |ψ 00 | 2 then tends to a sharp Gaussian.
A more rigorous connection to the harmonic oscillator wave functions can be made on the basis of the substitution u(z) = √ 2 q 1/4 cos z which transforms the Mathieu equation (3) to the new form where we defined = (a + 2q)/(4 √ q) and M(u(z)) = M(z; q). In the large q limit, this equation attains the form of the Schrödinger equation for the harmonic oscillator. Its eigenfunctions are well known and, apart from a normalization, are given by the parabolic cylinder functions D m (u) defined by where H m (u) are the Hermite polynomials. On the basis of this analysis, we may suspect that it is possible to find an asymptotic large-q expansion of the Mathieu functions in terms of harmonic oscillator functions D m of argument u. Sips [44][45][46][47] already derived such an expansion on the basis of the transformed Mathieu equation. For reference, in the next section we briefly outline its main features for the case of the Mathieu cosine, which is relevant for the discussion of singlet states. The general form of the Sips expansion is given by in which we defined D m<0 = 0. The specific form of the coefficients c 2n,l (q) is given in the work of Sips [44][45][46], who outlined a systematic procedure to obtain them. In general they can be obtained from a recursion relation [47] and we refer to Appendix B for a more detailed discussion. In Eq. (13)  even values of m. This follows directly from the derivation by Sips [44] but we see with hindsight that this condition is necessary to make the Mathieu cosine satisfy C l (z + π ; q) = (−1) l C l (z; q). Namely, if we replace z by z + π then the variable u changes to −u, yielding this desired property for a series of the form (13) Sips expansion, Eq. (13), will be used in the next section to determine the large interaction expansion of the density response function.
We conclude the section with a remark on the eigenenergies of the quantum ring. We can obtain an asymptotic expansion for them as the work of Sips also derives the large q behavior of the Mathieu characteristic values in terms of an asymptotic series expansion in power of q 1/2 (see Appendix B). Taking the first few leading orders, we obtain the following expression for the eigenenergies of the quantum ring: The asymptotic expansion is the same for the singlet and triplet energies as their difference becomes exponentially small in the large q limit (see Appendix B). To illustrate the q dependence of the eigenenergies, we present in Fig. 2 some of the lowest eigenvalues and their asymptotic expansion from Eq. (14) as a function of q. We see that the asymptotic expansion converges more slowly for higher values of l, and for these l we need high values of q in order to have a reliable estimate.

III. DENSITY RESPONSE OF STRONGLY INTERACTING ELECTRONS
After having discussed the two-particle wave function and energy spectrum of the system, let us now move to its response properties. Particularly relevant to TDDFT is the induced density change δn(r,t) when a small time-dependent external potential δv(r,t) is applied. They are related by the retarded 042505-4 density response function χ (rt,r t ) as follows: where χ is defined via (16) wheren H is the density operator in the Heisenberg picture and 0 is the ground state of the system. Since the unperturbed system is time independent, the density response function is a function of the relative time τ = t − t only and we can Fourier transform it with respect to τ : Before addressing in greater detail the properties of χ of the quantum ring in the large interaction limit, let us first make some considerations about the static density response function in a more general context.

A. Static density response in the strong interaction limit
Let us consider an interacting many-electron system in its ground state. The Hamiltonian consists of a kinetic energy operator, an external potential v(r), and a two-body interaction. If we consider a small variation δv(r) in the static external potential, the ground-state density will vary by an amount δn(r), which can be expressed as where χ (r,r ) = χ (r,r ,ω = 0) is the static density response function. Let us now consider a shifted potential v (r) = v(r + R). The ground-state density for this new potential is given by n (r) = n(r + R). For small translations, we can write that δn(r) = n (r) − n(r) = R · ∇n(r) and similarly δv(r) = v (r) − v(r) = R · ∇v(r). Since this is valid for all small vectors R, we find from Eq. (17) that This equation relates the gradient of the external potential to the gradient of the ground-state density, and amounts to the static limit of an equation derived for the dynamic density response function by Vignale [48].
Let us now consider a system in which we scale the twobody interaction with a parameter λ and let us choose the external potential v λ (r) in such a way that the density n(r) is the same for all values of λ. According to the Hohenberg-Kohn theorem [49], such a potential is unique when it exists. For such a system, the density response function will depend on λ as well and Eq. (18) becomes Let us now consider the limit of very large values of λ. One can show, for a general inhomogeneous system, that asymptotically v λ (r) = λu(r) + · · · , where u(r) is the socalled strictly correlated electron potential [32,33,50]. The result is intuitively clear as the linearly growing repulsive twobody interaction must be compensated by a linearly growing attractive one-body potential in order to keep the density profile constant. This has consequences for the behavior of χ λ . We consider two cases. Let us first assume that for λ → ∞ the response function χ λ attains a finite value α(r,r ). Because the left-hand side of Eq. (19) is independent of λ, this implies that which means that the three vector components of ∇u must be eigenfunctions of α with zero eigenvalue. Since the density response function reaches a finite limit, the system does not become rigid even when the interaction becomes infinitely large. This is a possible situation in systems in which the Hamiltonian by a coordinate transformation can be separated in two parts, in which one of the parts is weakly dependent on the interaction strength. That such inhomogeneous systems exist is demonstrated for the harmonic model system described in Appendix A and for which we demonstrate that Eq. (20) is indeed valid.
A probably more common situation is that such a separation is either not possible or that both parts of such a Hamiltonian are still strongly λ dependent. In this case, one would expect that the energy required to excite the system grows with λ, and as such the density response function would vanish for large λ. If this is the case, it is to be expected that χ λ in Eq. (19) asymptotically behaves as where β is a λ-independent function and the terms that follow decay faster than 1/λ. This means that for a given perturbation δv(r), the density response δn(r) decays as 1/λ and therefore the strong interaction makes the system more rigid and suppresses density variations. In such a case, Eq. (19) reduces to which is an exact equation for the leading order in λ.
When we finally consider systems in which the groundstate density and external potential are spatially constant, the reasoning that we carried out does not apply anymore since the gradients in Eq. (19) are identically zero. However, such systems are homogeneous, which implies that the center of mass can be separated off and we can therefore expect the response function to attain a finite value in the large interaction limit. This is exactly the case of our quantum ring model. Indeed, we saw in Eq. (11) that the quantum ring admits excitation energies that are independent of the interaction strength and these correspond to excitations that only change the center-of-mass wave function and do not affect the relative probability distribution of the particles. As we will see in more detail below, such excitations give a contribution to the density response function that survives in the large interaction limit, while the remaining excitations give a contribution which behaves as in Eq. (21).
It is interesting to connect this analysis to the f -sum rule for the dynamic density response function. In a system where the density is kept independent of λ with v λ (r), the f -sum rule attains the form [51] 1 π dω ω χ λ (r,r ,ω) = ∇[n(r)∇δ(r − r )].
We therefore see that the frequency integration removes the λ dependence. This is not in contradiction with Eq. (21). Although the density response function itself can become very small for large λ, the integrand in Eq. (23) can remain finite as it is weighted by the frequency ω. As a consequence the integral gets contributions proportional to the excitation energies, which grow with increasing interaction strength.

B. Exact density response of the quantum ring
After having discussed the general static case, we now turn our attention to the exact dynamical density response function of the quantum ring. Inserting a complete set of eigenstates of the HamiltonianĤ from Eq. (1) into the one-dimensional analog of Eq. (16), we find the Lehmann representation [51] of the retarded response function where we defined the excitation energies as E p kl (q) = E p kl (q) − E + 00 (q). The expression contains an infinitesimal parameter η > 0 that arises from the Fourier transform of the Heaviside function and the limit η → 0 is implied after the evaluation of all terms. Furthermore,n(x) is the density operator in the Schrödinger picture and p = ± labels the singlet or triplet eigenstates. The label k runs over all positive and negative integers while l runs over non-negative integers, with the condition that both are even or both are odd. The expression in Eq. (24) is simplified by the fact that the triplet terms vanish because the triplet spin function is orthogonal to the singlet spin function of the ground state, which yields + 00 |n(x)| − kl = 0. The remaining nonzero terms can be evaluated as where ψ + kl denotes the spatial part of the singlet wave function of Eq. (6) expressed in the original coordinates and the excitation amplitudes D kl (q) read The amplitude D kl (q) has a number of properties directly related to properties of the Mathieu functions. Since C l (z; q) is real, D * kl (q) = D (−k)l (q), and as a consequence of the orthogonality of the Mathieu functions, D 0l (q) = δ l0 . Moreover, the fact that C l (z + π ; q) = (−1) l C l (z; q) and that these functions are even in z leads to D kl (q) = (−1) k+l D * kl (q). Making use of the symmetry properties of D kl (q) described, combined with E + kl (q) = E + (−k)l (q), yields the following expansion of the response function: where in which the sum runs over even values of l for k even and over odd values of l for k odd. We see from Eq. (27) that χ (k,ω) can be regarded as the discrete Fourier transform of χ (x,x ,ω) with respect to the relative spatial coordinate x − x as was to be expected on the basis of the symmetry of the system. It will be now convenient to define the spatially discrete and temporally continuous Fourier transform of a function f (x,t) and its inverse as By using this Fourier transformation in Eq. (15) we rewrite the density response as in which k is an integer and ω is a continuous variable. We will make use of this relation below. We have now obtained an explicit form of the density response function that allows for an analytical analysis in the strong interaction limit. However, before moving to that, we briefly give the form of the response function for the noninteracting system, i.e., q = 0, which in the density functional context will be the same as the Kohn-Sham response function, since the system has the same density for all values of q. For the noninteracting case, the Mathieu characteristic value is a + l (0) = l 2 and the Mathieu cosine functions are given by C 0 (z; 0) = 1/ √ 2 and C l (z; 0) = cos(lz) for l 1. The excitation energies are given by Eq. (9), while the eigenstates are given by Eq. (6) as in which k ± l is always even. Note that l = |k| yields a doubly excited state. The corresponding excitation amplitude can be calculated from Eq. (26). Apart from the amplitude D 00 (0) = 1, which does not contribute to the Lehmann sum 042505-6 since E + 00 = 0, for (kl) = (00) we have that Since only terms where k + l is even contribute, we see that the only nonzero excitation amplitudes are the ones with l = |k|. This implies the absence of double excitations in the densityresponse function, a well-known property of noninteracting systems [4,52]. Inserting Eq. (33) into Eq. (28) we find that the noninteracting response function χ s (k,ω) is given by with E + kk (0) = 2(πk/L) 2 . We note that in k space, χ s has only a single pole for ω > 0 and no zeros.
Having determined the noninteracting response function, it now remains to study the density response function in the complementary limit of very strong interactions. For this purpose, we need to study the excitation energies E + kl (q) and excitation amplitudes D kl (q) in the limit of large q. This is the topic of the next section.

C. Strong interaction expansion of the dynamic density response function
Let us focus on the excitation energies E + kl (q) and the excitation amplitudes D kl (q) for large q. Because for E kl (q), explicit asymptotic expansions are known (see Appendix B), and this leaves us with the determination of D kl (q) defined by Eq. (26). We start by inserting the Sips expansion of Eq. (13) into Eq. (26), which gives where we defined where u(z) = √ 2 q 1/4 cos z. Since the coefficents c 2n,l (q) are known (see Appendix B for explicit expressions), it remains to evaluate J n 1 n 2 kl (q). Changing the integration variable to u and defining b = √ 2q 1/4 gives the expression where we defined the function and its Taylor coefficients a r (k). Inserting this Taylor series into Eq. (37) then gives the expansion where the interchange of integral and sum is allowed as we have an absolutely convergent series. Due to the Gaussian decay of the functions D n (u), in the limit q → ∞, we make an error which, as a function of q, decays faster than any polynomial function if we replace b in the limits of the integral by infinity. We therefore obtain the asymptotic expansion where we introduced coefficients of the form This integral can be computed analytically, and the explicit expression is given in Appendix C. Also note that due to the parity properties of the integrand I l n 1 n 2 ,r vanishes unless r and l are both even or both odd. Therefore, depending on whether l is even or odd, the summation index r in Eq. (39) can be taken to run only over even or only over odd values.
Expression (39) together with Eq. (35) gives an explicit procedure to calculate the large q expansion of the excitation amplitudes. The asymptotic expansions of D kl (q) and of |D kl (q)| 2 are given in Appendix C in Eqs. (C5) and (C6) respectively. Together with the asymptotic expansion for the excitation energies, Eq. (B4), inserted into Eq. (28), we find that the asymptotic expansion of the density response function is given by From this expression, we can draw a number of interesting conclusions. We find that the response function behaves quite differently for even and odd Fourier coefficients in the strong interaction regime. If we apply a potential with general coefficients δv(k,ω) to the system, the density change δn(k,ω) is strongly suppressed for odd k as χ (k,ω) becomes very small. In this limit it will therefore mainly have even Fourier coefficients, which implies that δn(x,t) = δn(x + L 2 ,t), i.e., the density change at antipodal points of the ring is the same. If we, however, apply a potential with only odd Fourier coefficients, the density change has the symmetry δn(x,t) = −δn(x + L 2 ,t) and is therefore opposite in antipodal points of the ring. From Eq. (41) we find that the leading term in real space in this case is given by We can understand the dependence on q as follows. The generation of an antisymmetric antipodal density requires excitation to states with an odd number of nodes in the relative wave function, which requires a large energy in the strong interaction limit and therefore the density response is suppressed for large interaction strength. From Eq. (42) we also see that the density increases instantaneously around the points where the potential has positive curvature. The instantaneous nature of the response has a simple explanation. If we perturb the system with a potential δv(k,ω) which is only nonzero for frequencies ω well below the first excitation energy, the temporal variation of the perturbation is much slower than a typical time scale of the free evolution of the system and the density response can be regarded as instantaneous. Since the excitation energies E + kl for odd k (which must have odd l as well) increase proportionally to √ q, the density response function in this case is well approximated by a frequency-independent function for ω √ q, which explains the instantaneous dependence of the density variation on the perturbation in Eq. (42).
For even values of k, the density response function has a more interesting frequency dependence. In the strong interaction limit The poles of this response function correspond to the center-of-mass excitations of Eq. (11). Being independent of q, they are not shifted towards infinity when we increase the interaction strength. This is a peculiarity of the quantum ring system as the Hamiltonian is separable in a λ-dependent and a λ-independent part. This happens also for some other homogeneous systems such as the three-dimensional electron gas with periodic boundary conditions or for electrons restricted to the surface of a sphere [53]. The analysis based on Eq. (19) shows that such a separation is usually not possible in inhomogeneous systems. When it is possible, both parts will generally still depend on λ (see Appendix A for an example).
To illustrate the accuracy of the expansion in Eq. (41) we display the exact response function and the expanded one in the top panels of Fig. 3 for k = 2 and for some values of the interaction strength. For small interactions (q = 1/3), the exact response function has two poles; one is approximately at the same location as the Kohn-Sham response function (ω = 2(πk/L) 2 ), while a new pole with a small weight appears at the center-of-mass excitation energy at ω = (πk/L) 2 . The expansion captures this pole, albeit with a very different weight. When we increase the interaction strength, the pole originally at the Kohn-Sham energy will shift to higher energies to a position proportional to √ q, while the pole corresponding to the center-of-mass excitation stays fixed and increases in weight. Already at q = 5, the asymptotic expansion yields good results for this k value. We thus see that the expansion is accurate for frequencies that are small compared to √ q. This result was to be expected since in the expansion of Eq. (41) we treated the frequency ω as a constant that is small compared √ q. Having obtained the exact response function, we have obtained all information needed to study the xc kernel of TDDFT, which will be the topic of the next section.

A. The exchange-correlation kernel
In TDDFT an effective noninteracting system, known as the Kohn-Sham system, is constructed in such a way as to have exactly the same density as the interacting many-particle system of interest. The external potential in this system, v s ([n]; rt), is a functional of the density [6,41] and is often written as follows v s (rt) = v(rt) + dr w(r,r )n(r ,t) + v xc (rt). (44) Here v(rt) is the external potential of the interacting system of interest and w(r,r ) is the two-particle interaction of that system. The second term in Eq. (44) is the Hartree potential and the last one is the exchange-correlation (xc) potential.
Taking the functional derivative of Eq. (44) with respect to the density, one obtains 042505-8 Here χ −1 s is the inverse of the Kohn-Sham density response function whereas χ −1 is the inverse of the density response function of the interacting system and f Hxc , the Hartree-xc kernel, is defined as where v Hxc is the sum of the Hartree and the xc potential. Equation (45) is commonly used to calculate χ from the knowledge of χ s at the price of approximating f Hxc . For the discussion in the next section, it is important to note that the functional derivative is not a uniquely defined function [41,54] due to the fact that for a system with a fixed number of particles the density change must integrate to zero at any time, i.e., 0 = dr δn(rt).
Let us define a new functioñ f Hxc (rt,r t ) = f Hxc (rt,r t ) + g(r,t,t ) + h(r ,t,t ) (48) with g and h arbitrary functions. The change in the Hartree-xc potential produced by the kernel of Eq. (48) due to a density change δn is given by δṽ Hxc (rt) = dr dt f Hxc (rt,r t )δn(r t ) = dr dt f Hxc (rt,r t )δn(r t )

+ dr dt [g(r,t,t ) + h(r ,t,t )]δn(r t )
= δv Hxc (rt) + C(t), (49) where the integral over g integrates to zero as a consequence of Eq. (47) and the integral over h yields a function C(t) of time t only, which is merely a gauge of the potential.
We therefore see thatf Hxc and f Hxc are physically equivalent integral kernels. The quantity that is defined unambiguously 1 is the mixed spatial derivative a property that will be used below.
Let us turn to the specific case of the quantum ring. The density response function is diagonal in the momentum-energy representation and for the Fourier components the following relation holds: By Fourier transforming the kernel f Hxc we impose a dependence on the relative coordinate in real space, which reduces the ambiguity of Eq. (48) to that of adding an arbitrary spatially constant function. In Eq. (51) this freedom is reflected in the fact that the kernel is well defined for all k values except for k = 0, since in this case both the response functions vanish. For the homogeneous quantum ring, the Kohn-Sham response function coincides with the response function of truly noninteracting electrons of Eq. (34). Using this equation, together with the expansion of Eq. (41), we obtain an explicit expression for f Hxc in the strong interaction limit: where we have reintroduced the variable λ rather than q as the λ notation is commonly used in the density functional context, which will be central for the discussion in the next section.
Since these quantities only differ by a numerical prefactor [see Eq. (4)] we will refer to the large interaction regime as the regime in which either of these two variables tends to infinity.
To illustrate the accuracy of the expansion in Eq. (52) we display the exact f Hxc kernel and the expanded one in the bottom panels of Fig. 3 for k = 2 and some values of the interaction strength. We see that the corresponding asymptotic expansion for f Hxc is accurate up to the lowest excitation energy corresponding to a change in the relative wave function. In this energy region f Hxc has a pole at an energy corresponding to a zero in χ , as a consequence of Eq. (51). It is worth noticing that since the noninteracting response function has no zeros, there is no pole in f Hxc originating from the first term in Eq. (51). This is peculiar to our quantum ring system for which the noninteracting response function has only a single pole for ω > 0. Instead, for a general system, the Kohn-Sham response function will have multiple poles and zeros which implies that, to order λ 0 , f Hxc is frequency dependent, causing the adiabatic approximation to fail in this order. In our system, f Hxc for even k tends to a frequency-independent function for all ω when the interaction strength approaches infinity. Its static value, given by f is the value needed to shift the Kohn-Sham pole to the q → ∞ pole.
For the odd k values (not shown here), all poles in the response function shift to infinity as q → ∞. The asymptotic expansion for f Hxc captures this, and the kernel becomes frequency independent in this limit. In fact, for odd values of k in Eq. (52), all the leading terms up to order λ 0 are frequency independent. However, frequency-dependent terms will appear to order λ −1/2 (not presented here) as is also the case for even k.
For the discussion in the next section it is useful to recast the kernel in real space. The expressions in Eq. (52) are sufficient to calculate this quantity to order λ 0 in real space using the and r = x − x is the relative distance between the points x and x , which are the one-dimensional counterparts of the spatial points in Eq. (46). Since x and x are both in the interval from 0 to L we have that r ∈ [−L,L]. We find that The leading term is given explicitly by in which we choose the arbitrary constant function [see the discussion below Eq. (51)] such that the Fourier coefficient of f 1 becomes zero for k = 0. For f 2 we find (up to a constant) that Since in the next section we will focus mostly on f 1 and f 2 , we do not report here the real space representation of f 3 . In the next section, we will show how f 1 and f 2 can be calculated in an alternative manner using the SCE theory in the adiabatic approximation.

B. Expanding the xc kernel in the theory of strictly correlated electrons
To date, no good and reliable approximations for electrons in the strong correlation regime have been developed within TDDFT. A ground-state theory of so-called strictly correlated electrons (SCE) [32,50] has been constructed and applied within the adiabatic approximation to calculate the xc kernel to the leading order in the interaction strength. In this section, we will first benchmark this approximation against the exact solution for the quantum ring model and then derive and compare the next order.
Let us begin with a brief overview of the ingredients of SCE theory that we will use. The Hartree-xc energy for a system with interaction strength λ can be written as [55] where we defined In this expression,Ŵ is the two-particle interaction and λ [n] is the ground-state wave function of a system with a local external potential, interaction λŴ , and ground-state density n.
In the strong interaction limit, W λ [n] can be expanded as [33] where the first term is the SCE energy, which has the explicit form where w is the two-particle interaction which we assume to depend only on the distance between the particles. The where D is the spatial dimensionality of the system and ω n are the harmonic frequencies. Inserting Eq. (59) into Eq. (57), we find the large λ expansion of the Hartree-xc energy to be The ground-state theory can be used in the adiabatic approximation [31] to find an approximate exchange-correlation kernel from where we added the superscript A to indicate that we make the adiabatic approximation. This approximation yields a frequency-independent Hartree-xc kernel when transformed to frequency space, which we denote as f A Hxc (r,r ). The second-order variation of Eq. (62) with respect to the density gives an expansion in orders of √ λ for the Hartree-xc kernel: where we defined the adiabatic SCE and ZPE kernels as Let us now turn again to the case of the quantum ring. In this case, there are two electrons and just one simple comotion function f : [0,L] → [0,L] given by If one electron is at x, this function simply puts the other electron at the antipodal point of the quantum ring. From this comotion function, it is straightforward to calculate the SCE energy. Since |x − f(x)| = L/2, we have from our interaction w(x) = V 0 cos 2 (πx/L) in Hamiltonian (1) that w(|x − f(x)|) = 0 and therefore V SCE = 0 for our density. Physically this means that the electrons are simply localized at the bottom of a potential well with zero energy.
To next order, oscillations around these equilibrium positions start to appear. These zero-point oscillations give an energy contribution which can be calculated using Eq. (61). For a one-dimensional two-electron system there is only one nonzero harmonic frequency given by [56] where w (x) = ∂ 2 x w(x). If we calculate this frequency for our quantum ring, we find ω 1 We can verify that this is in accordance with the exact strong interaction expansion of the Hartree-xc energy. Since the Kohn-Sham kinetic energy as well as the external potential is zero for the quantum ring system, E λ Hxc simply coincides with the total energy, which is known in the strong interaction limit from the large-q expansion of the lowest Mathieu characteristic value (see Appendix B). This gives where the leading coefficient indeed exactly gives 2V ZPE . The expansion of W λ can be calculated from Eq. (57) to give This result agrees with a direct calculation of W λ using the Sips expansion of the Mathieu functions and indeed has the structure of the expansion in Eq. (59) in which we also included the term to order λ −3/2 .
Let us turn to the calculation of the kernels of Eqs. (65) and (66).
The ASCE kernel is obtained from the expression derived in Ref. [31] and reads where f 1 is the function given in Eq. (55). Because of the freedom in Eq. (48) f ASCE is physically equivalent to the kernel f 1 (x − x ) and agrees with the leading term in the expansion of Eq. (54) in the strong interaction limit. Both kernels are shown in Fig. 4. Since f 1 (x − x ) has a simpler shape, we will restrict ourselves to this function. The kernel describes how a By taking the second derivative of this equation, we obtain where we stress that δn(x,t) is periodic with L. Since f 1 (x − x ) is linear everywhere except at the kinks at x − x = 0, ± L/2, the second derivatives yield δ functions at these points. Also note that the sign of the density change yields the curvature of the induced potential. Equation (72) has some interesting consequences. If we make a localized density variation δn(x,t) in a very small interval of the ring, there will not only be a change δv Hxc in the same interval, but at the same time a similar change in the potential with opposite sign in an antipodal interval. This shows very clearly that the Hxc-potential depends nonlocally, but instantaneously, on the density. Let us analyze the next orders in the strong interaction expansion. The calculation of f AZPE allows for a comparison with the next leading term in Eq. (54), which is proportional to √ λ. The kernel f AZPE was obtained by taking the second functional derivative of the one-dimensional counterpart of Eq. (61), using Eq. (68) and the functional derivative of the comotion functions (derived already in Ref. [31]). We also find agreement between this expression and that of the next leading term in Eq. (56), i.e., f AZPE (x,x ) = f 2 (x − x ) modulo the addition of arbitrary functions of x and x separately [see again Eq. (48)]. We observe that the first two leading terms of the expansion of the Hartree-xc kernel from the adiabatic SCE theory agree with the exact results for the quantum ring. An interpretation of this fact will be presented below.
We have thus seen that, in this model, the ASCE and AZPE terms agree with the terms f 1 and f 2 respectively of the exact asymptotic expansion. To better elucidate their role in the strong interaction limit, we show in Fig. 5 the first three terms contributing to f Hxc , all scaled by λ, and compare them with the expression for the exact kernel, in Fourier space for k = 3 and k = 5. As was pointed out earlier, the accuracy of the expansion depends on the value of k: high k values require higher λ values to achieve better accuracy. The first term, that is the ASCE, is constant, while the second one, that is the AZPE, only improves on it for large λ values and worsens it for smaller ones, as one would expect for an asymptotic expansion. The third term, beyond the AZPE, also exhibits a non-negligible contribution in the small-λ regime.
We will now offer a physical interpretation of the above terms and make some considerations about their properties in the case of more general systems than the quantum ring model. In the (infinitely) strong interaction limit, a given system behaves very rigidly, since the position of the reference electron determines the positions of all the remaining electrons. Upon application of a perturbation, the response of the system is instantaneous, or adiabatic, while maintaining its rigidity, unless special symmetries are present. This behavior is likely to apply to a wider class of systems, both with a uniform (such as the quantum ring) and a nonuniform density. On the other hand, it is unclear whether the frequency independence of f 2 is equally general as we already move away from the strictly correlated electron limit by introducing zero-point vibrations: Thus the adiabaticity of f 2 for general systems is still an open issue. Finally, as already noted before, the third term f 3 of the expansion will be nonadiabatic for general systems.

V. CONCLUSIONS
In this work, we have considered an exactly solvable model consisting of two interacting electrons on a quantum ring. We focused on the response properties and calculated the energy spectrum, the excitation amplitudes, the density response function, and the exchange-correlation kernel of time-dependent density functional theory. In the limit of strong interaction, we derived the asymptotic expansion in powers of the square root of the interaction strength for the response function and kernel. For the kernel, we found that its leading terms are local in time but nonlocal in space. This already shows that the commonly used adiabatic local-density, or semilocal, approximations will fail for such strongly correlated systems, since they are local both in time and space. We compared the expansion for the kernel to a similar one obtained from the adiabatic-SCE formalism [31] which has the spatial nonlocality built in. The leading term of the exact expansion was found to coincide with the adiabatic-SCE kernel derived in Ref. [31]. After working out the next order term, the so-called zero-point energy contribution, we found that it also coincided with the exact next-to-leading term. For our model, the subsequent term in the expansion is still adiabatic, but we showed that in general systems this term will be nonadiabatic.
The agreement with our exact results puts the adiabatic-SCE and the adiabatic-ZPE approximations on firmer ground and gives confidence in employing the formalism of strictly correlated electrons in the adiabatic approximation for calculating response properties of strongly correlated systems. To illustrate some properties of the response function in the large interaction limit for an inhomogeneous system, let us analyze a model system of two harmonically confined electrons in three dimensions with a harmonic repulsion [57]. The Hamiltonian of the system is given bŷ in which the harmonic frequency ω λ is chosen in such a way that the density is independent of λ. Using the coordinate transformation s = (r 1 + r 2 )/ √ 2 and r = (r 1 − r 2 )/ √ 2, the Hamiltonian can be written as that of two independent harmonic oscillatorŝ where |n| = n 1 + n 2 + n 3 and H n (x) is the Hermite polynomial of order n. The energy eigenvalues are given by E nm = ω λ |n| + 3 2 + ν λ |m| + 3 2 .
From this function, the density is readily obtained as where we defined β = 2ω λ ν λ ω λ + ν λ . (A8) If we insert ν λ = ω 2 λ − 2λ into this relation, we can determine the λ dependence of ω λ as β is independent of λ. We find where y solves the quartic equation In the limit of large interaction, we find that y = 1 + β/(2 √ 2λ) + O(λ −1 ) such that We see that the harmonic frequency of the center-of-mass mode approaches infinity, whereas the one of the relative mode approaches a finite value. This has interesting consequences for the excitation spectrum. For the excitation energies of the relative mode, it implies that while all other excitation energies diverge to infinite values at large interaction strength. The latter correspond to excitations of the center-of-mass mode. In contrast to the quantum ring, only the excitation energies of the relative mode remain finite in the large interaction limit, which is due to the very different nature of the two-body interaction. Let us turn our attention to the density response function. In the response function, only the singlet excitations contribute. This means that the spatial wave functions that we need to consider are symmetric in the interchange of the particle positions. For this to be true, the relative wave functions need to be even and we have to require that |m| only attains even values. The response function therefore has the form χ (r 1 ,r 2 ,ω) = n,m D nm (r 1 )D * nm (r 2 ) ω − E nm + iη − D nm (r 2 )D * nm (r 1 ) where E nm = E nm − E 00 are the excitation energies and we further put the restriction that we sum over all m such that |m| is even. The excitation amplitudes corresponding to these excitations are given by D nm (r 1 ) = 0 |n(r 1 )| nm = 2 dr 2 * 00 (r 1 ,r 2 ) nm (r 1 ,r 2 ), where we rewrote the eigenfunctions in terms of the original coordinates. Let us now consider the large interaction limit λ → ∞ of the response function. Since only the excitation energies of the form E 0m remain finite in this limit, we only need to consider the excitation amplitudes of the form D 0m (r 1 ) = 2 ω λ ν λ π 2 3/2 dr 2 H m ( ν λ /2(r 1 − r 2 )) × e − ω λ 2 (r 1 +r 2 ) 2 − ν λ 2 (r 1 −r 2 ) 2 . (A17) If we now use that lim λ→∞ ω λ 2π 3/2 e − ω λ 2 (r 1 +r 2 ) 2 = δ(r 1 + r 2 ) (A18) is a limit representation for the delta distribution and the fact that ν λ → β/2 in the large interaction limit, we find that lim λ→∞ D 0m (r 1 ) = 2 β π For the response function in the large interaction limit, we therefore obtain α(r 1 ,r 2 ,ω) = lim λ→∞ χ (r 1 ,r 2 ,ω) = n(r 1 )n(r 2 ) m H m ( β r 1 )H m ( β r 2 ) where the sum runs over m values such that |m| is even. We see that in the large interaction limit it is still possible to excite the relative modes of the system.
Let us now have a look at the external potential in the strong interaction limit v λ (r) = 1 2 ω 2 λ |r| 2 = λ|r| 2 + O( which implies ∇v λ (r) = 2λr + O( √ λ). Since the response function does not vanish in the large interaction limit, Eq. (20) tells us that 0 = dr 2 α(r 1 ,r 2 ; ω = 0) r 2 must hold, where we took the static limit of the response function of Eq. (A20). Since the polynomial functions in Eq. (A20) are all even (since |m| is even), and r 2 is an odd function, we find that this relation is indeed satisfied. more localized and thereby makes the expansion of Eq. (38) used in the integrand more accurate.
We finally note that the frequency-sum rule, Eq. (23), can be used to check the validity of some of the terms of Eq. (C6). If we insert the explicit form of Eq. (28) into the one-dimensional equivalent of Eq. (23) for the frequency sum rule, we can derive that l (k 2 + [a + l (q) − a + 0 (q)])|D kl (q)| 2 = k 2 , where the sum runs over even l if k is even and over odd l when k is odd. The right-hand side is independent of the interaction strength q; thus in the sum on the left-hand side the q dependence in the excitation amplitudes has to be compensated by the q dependence of the Mathieu characteristic values to give a result independent of the interaction strength.