Time-dependent density-functional theory for strongly interacting electrons

TL;DR

This paper analytically investigates the exact exchange-correlation kernel in time-dependent density functional theory for strongly interacting electrons, revealing the dominant instantaneous terms and the conditions under which the adiabatic approximation fails.

Contribution

It provides an analytical expansion of the exchange-correlation kernel in the strong interaction limit and compares it with the strictly correlated electron approach, highlighting the approximation's validity.

Findings

Leading terms in the kernel act instantaneously in the strong interaction limit.

Memory effects appear only in higher-order terms of the expansion.

The adiabatic strictly correlated electron approximation matches the exact leading terms.

Abstract

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.

Figures

Equations

$$\hat{H}=-\frac{1}{2}({\partial }_{x_1}^2+{\partial }_{x_2}^2)+\lambda V_0{\cos }^2\left[\frac{\pi }{L}(x_1-x_2)\right]$$

$$\Psi (x_1{\sigma }_1,x_2{\sigma }_2)={\psi }^{\pm }(x_1,x_2){\Xi }^{\pm }({\sigma }_1,{\sigma }_2).$$

$${\Xi }^{+}({\sigma }_1,{\sigma }_2)=\frac{1}{\sqrt{2}}\left({\delta }_{{\sigma }_1\uparrow }{\delta }_{{\sigma }_2\downarrow }-{\delta }_{{\sigma }_1\downarrow }{\delta }_{{\sigma }_2\uparrow }\right)$$

$${\psi }^{\pm }(x_1,x_2)=\pm {\psi }^{\pm }(x_2,x_1).$$

$$\hat{H}=-\frac{1}{4}{\partial }_R^2-\frac{{\pi }^2}{L^2}{\partial }_z^2+\lambda V_0{\cos }^2(z).$$

$$\psi (R,z)=\exp \left(\frac{2\pi ikR}{L}\right)M(z)$$

$$\left[-{\partial }_z^2+2q\cos (2z)\right]M(z)=aM(z)$$

$$q$$

$$a$$

$${\psi }_{kl}^{+}(R,z;q)$$

$${\psi }_{kl}^{-}(R,z;q)$$

$$\underset{0}{\overset{\pi }{\int }}dz{\left|C_l(z;q)\right|}^2=\underset{0}{\overset{\pi }{\int }}dz{\left|S_l(z;q)\right|}^2=\frac{\pi }{2}.$$

$$E_{kl}^{\pm }(q)={\left(\frac{\pi }{L}\right)}^2\left[k^2+a_l^{\pm }(q)+2q\right].$$

$${\Psi }_{kl}(x_1{\sigma }_1,x_2{\sigma }_2)={\Phi }_k(x_1{\sigma }_1,x_2{\sigma }_2)\sqrt{2}{C}_l\left[\frac{\pi }{L}\left(x_1-x_2\right)\right]$$

$$\Phi_{k}(x_{1}\sigma_{1},x_{2}\sigma_{2})=\frac{1}{\sqrt{2}}\left|\begin{array}[]{cc}\phi_{k/2}(x_{1})\delta_{\sigma_{1}\uparrow}&\phi_{k/2}(x_{1})\delta_{\sigma_{1}\downarrow}\\ \phi_{k/2}(x_{2})\delta_{\sigma_{2}\uparrow}&\phi_{k/2}(x_{2})\delta_{\sigma_{2}\downarrow}\end{array}\right|$$

$$\Delta E_{k0}^{+}=E_{k0}^{+}-E_{00}^{+}={\left(\frac{\pi k}{L}\right)}^2,$$

$$|{\psi }_{kl}^{\pm }(R,z;q){|}^2=\frac{2}{L^2}{M}_l^2(z;q)$$

$$2q\cos (2z)=-2q+4q(z-\frac{\pi }{2}{)}^2+\dots$$

$$a_l^{\pm }(q)=-2q+2\sqrt{q}(2l+1)+\dots$$

$$\left[-\frac{1}{2}(1-\frac{u^2}{2\sqrt{q}})\frac{d^2}{d{u}^2}+\frac{u}{4\sqrt{q}}\frac{d}{du}+\frac{u^2}{2}\right]\mathcal{M}(u)=ϵ\mathcal{M}(u)$$

$${\mathcal{D}}_m(u)$$

$$=\frac{1}{2^{m/2}}{e}^{-u^2/2}{H}_m(u)$$

$$C_l(z;q)=\sum _{n=-\infty }^{\infty }{c}_{2n,l}(q){\mathcal{D}}_{l+2n}\left[u(z)\right]$$

$$E_{kl}^{\pm }$$

$$+\frac{(2l+1)}{128\sqrt{q}}\left((2l+1{)}^2+3\right)]+\mathcal{O}\left(q^{-1}\right).$$

$$\delta n(\mathbf{r},t)=\int d{\mathbf{r}}^{\prime }\int d{t}^{\prime }\chi (\mathbf{r}t,{\mathbf{r}}^{\prime }{t}^{\prime })\delta v({\mathbf{r}}^{\prime },t^{\prime }),$$

$$\chi (\mathbf{r}t,{\mathbf{r}}^{\prime }{t}^{\prime })$$

$$-i\theta (t-t^{\prime })⟨{\Psi }_0|[{\hat{n}}_H(\mathbf{r},t),{\hat{n}}_H({\mathbf{r}}^{\prime },t^{\prime })]|{\Psi }_0⟩$$

$$\chi (\mathbf{r},{\mathbf{r}}^{\prime },\omega )=\int d\tau \chi (\mathbf{r},{\mathbf{r}}^{\prime },\tau )e^{i\omega \tau }.$$

$$\delta n(\mathbf{r})=\int d{\mathbf{r}}^{\prime }\chi (\mathbf{r},{\mathbf{r}}^{\prime })\delta v({\mathbf{r}}^{\prime })$$

$$\nabla n(\mathbf{r})=\int d{\mathbf{r}}^{\prime }\chi (\mathbf{r},{\mathbf{r}}^{\prime }){\nabla }^{\prime }v({\mathbf{r}}^{\prime }).$$

$$\nabla n(\mathbf{r})=\int d{\mathbf{r}}^{\prime }{\chi }_{\lambda }(\mathbf{r},{\mathbf{r}}^{\prime }){\nabla }^{\prime }{v}_{\lambda }({\mathbf{r}}^{\prime }).$$

$$0=\int d{\mathbf{r}}^{\prime }\alpha (\mathbf{r},{\mathbf{r}}^{\prime }){\nabla }^{\prime }u({\mathbf{r}}^{\prime })$$

Full text

I Introduction

Time-dependent density functional theory (TDDFT) Gross et al. (1996); Burke et al. (2005); van Leeuwen (2001); Ullrich (2012); Marques et al. (2012); Ruggenthaler et al. (2015); Maitra (2016) 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, 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 non-interacting 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_{\textrm{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 Petersilka et al. (1996). For this purpose it is enough to know $v_{\textrm{xc}}$ and its functional derivative $\delta v_{\textrm{xc}}/\delta n=f_{\textrm{xc}}$ with respect to the density $n$, evaluated at the ground-state density. The quantity $f_{\textrm{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 Ullrich (2012) such as, for example, the inability to produce correct charge transfer excitations Stein et al. (2009); Gritsenko and Baerends (2004); Maitra (2005); Autschbach (2009) and Born-Oppenheimer surfaces of excited states in dissociating molecules Gritsenko et al. (2000); Hieringer and Görling (2006); Giesbertz and Baerends (2008) and semiconductor band gaps Botti et al. (2004). Some improvements have been made using hybrid functionals which contain mixtures of exact exchange and traditional local functionals. These methods are non-local in space but still adiabatic. However, they are not systematic and the optimal mixture of exact-exchange is often system-dependent Dreuw et al. (2003); Hofmann et al. (2012); Dierksen and Grimme (2004); Peach et al. (2008). Other, more systematic, approximations for $f_{\textrm{xc}}$ beyond the ALDA often rely on perturbative expansions Görling (1998); Tokatly and Pankratov (2001); von Barth et al. (2005); Fuchs et al. (2005); Hellgren and von Barth (2009, 2010); Heßelmann and Görling (2011a, b); Bates and Furche (2013); Erhard et al. (2016) 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 Lani et al. (2016) the so-called strictly correlated electrons (SCE) framework Seidl (1999); Gori-Giorgi et al. (2009); Malet and Gori-Giorgi (2012), 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 Mundt et al. (2007), an exact property related to generalized translational invariance Vignale (1995a); Dobson (1994). 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 non-local and exhibits a divergent behavior as the molecular bond is stretched. For adiabatic kernels this diverging behavior is crucial Gritsenko et al. (2000) 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_{\textrm{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 Viefers et al. (2004); Ruggenthaler et al. (2012); Loos and Gill (2012); Ruggenthaler et al. (2015), for which not only we compute the exact density response function and xc-kernel, but we 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_{\textrm{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.

II An exactly solvable system

II.1 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 by:

[TABLE]

where $\lambda\geq 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 [math] to $L$. The ground-state density $n_{0}=2/L$ is spatially constant and independent of $\lambda$; 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 $\Psi$ which satisfy the stationary Schrödinger equation $\hat{H}\Psi=E\Psi$ 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

[TABLE]

For the singlet case (which we will focus on) the normalized spin function is given by

[TABLE]

and is anti-symmetric in the spin variables. For the triplet there are three linearly independent symmetric spin functions which we, for simplicity, all denote by $\Xi^{-}(\sigma_{1},\sigma_{2})$. Since the two-electron wave function $\Psi$ is anti-symmetric under the simultaneous interchange of space and spin variables it follows that the spatial wave functions $\psi^{\pm}$ satisfy the symmetry relation

[TABLE]

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}\rightarrow 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. Ruggenthaler et al. (2013) 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=\pi(x_{1}-x_{2})/L$, the dimensionless relative coordinate. This variable transformation gives

[TABLE]

By inserting a product Ansatz of the form $\psi(R,z)=f(R)M(z)$ into the Schrödinger equation we find that the spatial two-particle eigenfunctions are of the form

[TABLE]

where $k$ is an integer and the function $M(z)$ satisfies the Mathieu equation which we write in its standard form as Arscott et al. (1964); DLMF

[TABLE]

where the constants $q$ and $a$ are given by:

[TABLE]

For a given value of $q$ the Mathieu equation (3) has only periodic solutions for particular values $a(q)$ which are called the Mathieu characteristic values. Moreover this equation has either even or odd periodic solutions which are called the Mathieu cosine and Mathieu sine functions respectively. Both sets of functions form a countable set, therefore its members can be labeled by a non-negative integer $l$. For 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

[TABLE]

in which the normalization of the Mathieu functions is chosen such that

[TABLE]

The Mathieu functions further have the periodicity property $M_{l}(z+\pi)=(-1)^{l}M_{l}(z)$, i.e. they are periodic in $\pi$ for even values of $l$ and anti-periodic for odd values of $l$. Furthermore the center-of-mass wave function in Eq.(2) changes with a prefactor $(-1)^{k}$ when $x_{i}\rightarrow 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

[TABLE]

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. Multiplying the spatial wave function of Eq.(6) with its singlet spin function, the full space-spin function can be written as

[TABLE]

where we defined the Slater determinant

[TABLE]

and we further defined the spatial normalized orbital by $\phi_{k}(x)=e^{2\pi ikx/L}/\sqrt{L}$ which corresponds to a periodic single-particle 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 $\phi_{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

[TABLE]

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/\sqrt{2}$ and the ground and excited state 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 non-zero 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 non-zero $l$ becomes very large too and the excitations with energy $\Delta E_{k0}^{+}$ give the dominant contribution to the density response function, as we will see later.

II.2 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

[TABLE]

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 single-particle 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=\pi/2$, which corresponds to a relative distance of the particles of $L/2$. We can expand the potential around this minimum to obtain

[TABLE]

This potential describes (apart from a shift of the minimum) a harmonic oscillator with frequency $\Omega=2\sqrt{q}$. The eigenfunctions of the harmonic oscillator are well known to consist of Gaussians of width proportional to $1/\sqrt{\Omega}=1/(\sqrt{2}\,q^{1/4})$. In the limit of large $q$ the harmonic frequency increases and the wave functions become localized around $z=\pi/2$. This behaviour is illustrated in Fig.(1). The eigenenergies $\epsilon_{l}$ of the harmonic oscillator are well-known and given by $\epsilon_{l}=-q+\Omega(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$:

[TABLE]

and consequently also an asymptotic expansion for the eigenenergies of the quantum ring from Eq.(9).

A more rigorous connection to the harmonic oscillator wave functions can be made on the basis of the substitution $u(z)=\sqrt{2}\,q^{1/4}\cos z$ which transforms the Mathieu equation (3) to the new form

[TABLE]

where we defined $\epsilon=(a+2q)/(4\sqrt{q})$ and $\mathcal{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 $\mathcal{D}_{m}(u)$ defined by

[TABLE]

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 $\mathcal{D}_{m}$ of argument $u$. Sips Sips (1949, 1959); Blanch (1960); Frenkel and Portugal (2001) 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

[TABLE]

in which we defined $\mathcal{D}_{m<0}=0$. The specific form of the coefficients $c_{2n,l}(q)$ is given in the work of Sips Sips (1949, 1959); Blanch (1960) who outlined a systematic procedure to obtain them. In general they can be obtained from a recursion relation Frenkel and Portugal (2001) and we refer to Appendix B for a more detailed discussion.

In Eq.(13) we see that for odd values of $l$ the Mathieu cosine is expanded in functions $\mathcal{D}_{m}$ with only odd values of $m$ while for even $l$ it is expanded in functions $\mathcal{D}_{m}$ with only even values of $m$. This follows directly from the derivation by Sips Sips (1949) but we see with hindsight that this condition is necessary to make the Mathieu cosine satisfy $C_{l}(z+\pi;q)=(-1)^{l}C_{l}(z;q)$. Namely, if we replace $z$ by $z+\pi$ then the variable $u$ changes to $-u$ yielding this desired property for a series of the form (13) since $\mathcal{D}_{m}(-u)=(-1)^{m}\mathcal{D}_{m}(u)$. The 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

[TABLE]

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 $\delta n(\mathbf{r},t)$ when a small time-dependent external potential $\delta v(\mathbf{r},t)$ is applied. They are related by the retarded density response function $\chi(\mathbf{r}t,\mathbf{r}^{\prime}t^{\prime})$ as follows

[TABLE]

where $\chi$ is defined via

[TABLE]

where $\hat{n}_{H}$ is the density operator in the Heisenberg picture and $\Psi_{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 $\tau=t-t^{\prime}$ only and we can Fourier transform it with respect to $\tau$:

[TABLE]

Before addressing in greater detail the properties of $\chi$ 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.

III.1 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(\mathbf{r})$ and a two-body interaction. If we consider a small variation $\delta v(\mathbf{r})$ in the static external potential, the ground-state density will vary by an amount $\delta n(\mathbf{r})$ which can be expressed as

[TABLE]

where $\chi(\mathbf{r},\mathbf{r}^{\prime})=\chi(\mathbf{r},\mathbf{r}^{\prime},\omega=0)$ is the static density response function. Let us now consider a shifted potential $v^{\prime}(\mathbf{r})=v(\mathbf{r}+\mathbf{R})$. The ground-state density for this new potential is given by $n^{\prime}(\mathbf{r})=n(\mathbf{r}+\mathbf{R})$. For small translations we can write that $\delta n(\mathbf{r})=n^{\prime}(\mathbf{r})-n(\mathbf{r})=\mathbf{R}\cdot\nabla n(\mathbf{r})$ and similarly $\delta v(\mathbf{r})=v^{\prime}(\mathbf{r})-v(\mathbf{r})=\mathbf{R}\cdot\nabla v(\mathbf{r})$. Since this is valid for all small vectors $\mathbf{R}$ we find from Eq.(17) that

[TABLE]

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 Vignale (1995b).

Let us now consider a system in which we scale the two-body interaction with a parameter $\lambda$ and let us choose the external potential $v_{\lambda}(\mathbf{r})$ in such a way that the density $n(\mathbf{r})$ is the same for all values of $\lambda$. According to the Hohenberg-Kohn theorem Hohenberg and Kohn (1964) such a potential is unique when it exists. For such a system the density response function will depend on $\lambda$ as well and Eq.(18) becomes

[TABLE]

Let us now consider the limit of very large values of $\lambda$. One can show, for a general inhomogeneous system, that asymptotically $v_{\lambda}(\mathbf{r})=\lambda u(\mathbf{r})+\ldots$ where $u(\mathbf{r})$ is the so-called strictly-correlated electron potential Seidl (1999); Seidl et al. (2007); Gori-Giorgi et al. (2009). The result is intuitively clear as the linearly growing repulsive two-body interaction must be compensated by a linearly growing attactive one-body potential in order to keep the density profile constant. This has consequences for the behavior of $\chi_{\lambda}$. We consider two cases. Let us first assume that for $\lambda\rightarrow\infty$ the response function $\chi_{\lambda}$ attains a finite value $\alpha(\mathbf{r},\mathbf{r}^{\prime})$. Because the left-hand side of Eq.(19) is independent of $\lambda$ this implies that

[TABLE]

which means that the three vector components of $\nabla u$ must be eigenfunctions of $\alpha$ 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 $\lambda$-dependent. In this case one would expect that the energy required to excite the system grows with $\lambda$, and as such the density response function would vanish for large $\lambda$. If this is the case it is to be expected that $\chi_{\lambda}$ in Eq.(19) asymptotically behaves as

[TABLE]

where $\beta$ is a $\lambda$-independent function and the terms that follow decay faster than $1/\lambda$. This means that for a given perturbation $\delta v(\mathbf{r})$ the density response $\delta n(\mathbf{r})$ decays as $1/\lambda$ and therefore the strong interaction makes the system more rigid and suppresses density variations. In such a case Eq.(19) reduces to

[TABLE]

which is an exact equation for the leading order in $\lambda$.

When we finally consider systems in which the ground-state 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 $\lambda$ with $v_{\lambda}(\mathbf{r})$ the $f$-sum rule attains the form Stefanucci and van Leeuwen (2013)

[TABLE]

We therefore see that the frequency integration removes the $\lambda$-dependence. This is not in contradiction with Eq.(21). Although the density response function itself can become very small for large $\lambda$ the integrand in Eq.(23) can remain finite as it is weighted by the frequency $\omega$. As a consequence the integral gets contributions proportional to the excitation energies which grow with increasing interaction strength.

III.2 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 $\hat{H}$ from Eq.(1) into the one-dimensional analogue of Eq.(16) we find the Lehmann representation Stefanucci and van Leeuwen (2013) of the retarded response function

[TABLE]

where we defined the excitation energies as $\Delta E_{kl}^{p}(q)=E_{kl}^{p}(q)-E_{00}^{+}(q)$. The expression contains an infinitesimal parameter $\eta>0$ that arises from the Fourier transform of the Heaviside function and the limit $\eta\rightarrow 0$ is implied after the evaluation of all terms. Furthermore $\hat{n}(x)$ is the density operator in the Schrödinger picture and $p=\pm$ 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 $\langle\Psi_{00}^{+}|\hat{n}(x)|\Psi_{kl}^{-}\rangle=0$. The remaining non-zero terms can be evaluated as

[TABLE]

where $\psi_{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

[TABLE]

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)=\delta_{l0}$. Moreover, the fact that $C_{l}(z+\pi;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 $\Delta E_{kl}^{+}(q)=\Delta E_{(-k)l}^{+}(q)$, yields the following expansion of the response function

[TABLE]

where

[TABLE]

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 $\chi(k,\omega)$ can be regarded as the discrete Fourier transform of $\chi(x,x^{\prime},\omega)$ with respect to the relative spatial coordinate $x-x^{\prime}$ 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:

[TABLE]

By using this Fourier transformation in Eq.(15) we rewrite the density response as

[TABLE]

in which $k$ is an integer and $\omega$ 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 non-interacting 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/\sqrt{2}$ and $C_{l}(z;0)=\cos(lz)$ for $l\geq 1$. The excitation energies are given by Eq.(9),

[TABLE]

while the eigenstates are given by Eq.(6) as

[TABLE]

in which $k\pm l$ is always even. Note that $l\neq|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 since $\Delta E_{00}^{+}=0$, for $(kl)\neq(00)$ we have that

[TABLE]

Since only terms where $k+l$ is even contribute, we see that the only non-zero excitation amplitudes are the ones with $l=|k|$. This implies the absence of double excitations in the density-response function, a well-known property of non-interacting systems Maitra et al. (2004); Ullrich (2012). Inserting Eq.(33) into Eq.(28) we find that the non-interacting response function $\chi_{s}(k,\omega)$ is given by

[TABLE]

with $\Delta E_{kk}^{+}(0)=2(\pi k/L)^{2}$. We note that in $k$-space, $\chi_{s}$ has only a single pole for $\omega>0$, and no zeroes.

Having determined the non-interacting 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 $\Delta E_{kl}^{+}(q)$ and excitation amplitudes $D_{kl}(q)$ in the limit of large $q$. This is the topic of the next section.

III.3 Strong interaction expansion of the dynamic density response function

Let us focus on the excitation energies $\Delta E_{kl}^{+}(q)$ and the excitation amplitudes $D_{kl}(q)$ for large $q$. Because for $\Delta E_{kl}^{+}(q)$ explicit asymptotic expansions are known (see Appendix B), 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

[TABLE]

where we defined

[TABLE]

where $u(z)=\sqrt{2}\,q^{1/4}\cos z$. Since the coefficients $c_{2n,l}(q)$ are known (see Appendix B for explicit expressions) it remains to evaluate $\mathcal{J}_{kl}^{n_{1}n_{2}}(q)$. Changing the integration variable to $u$ and defining $b=\sqrt{2}q^{1/4}$ gives the expression

[TABLE]

where we defined the function

[TABLE]

and its Taylor coefficients $a_{r}(k)$. Inserting this Taylor series into Eq.(37) then gives the expansion

[TABLE]

where the interchange of integral and sum is allowed as we have an absolutely convergent series. Due to the Gaussian decay of the functions $\mathcal{D}_{n}(u)$, in the limit $q\rightarrow\infty$, 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

[TABLE]

where we introduced coefficients of the form

[TABLE]

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_{n_{1}n_{2},r}^{l}$ 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 $\left|D_{kl}(q)\right|^{2}$ are given in Appendix C in Eq.(105) and Eq.(106) respectively. Together with the asymptotic expansion for the excitation energies, Eq.(97), inserted into Eq.(28) we find that the asymptotic expansion of the density response function is given by

[TABLE]

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 $\delta v(k,\omega)$ to the system, the density change $\delta n(k,\omega)$ is strongly suppressed for odd $k$ as $\chi(k,\omega)$ becomes very small. In this limit it will therefore mainly have even Fourier coefficients which implies that $\delta n(x,t)=\delta n(x+\frac{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 $\delta n(x,t)=-\delta n(x+\frac{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

[TABLE]

We can understand the dependence on $q$ as follows. The generation of an anti-symmetric 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 were the potential has positive curvature. The instantaneous nature of the response has a simple explanation. If we perturb the system with a potential $\delta v(k,\omega)$ which is only non-zero for frequencies $\omega$ well below the first excitation energy, the temporal variation of the perturbation is much slower than a typical timescale of the free evolution of the system and the density response can be regarded as instantaneous. Since the excitation energies $\Delta E_{kl}^{+}$ for odd $k$ (which must have odd $l$ as well) increase proportionally to $\sqrt{q}$ the density response function in this case is well approximated by a frequency independent function for $\omega\ll\sqrt{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

[TABLE]

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 $\lambda$-dependent and a $\lambda$-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 Seidl (2007). 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 $\lambda$ (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 $(\omega=2(\pi k/L)^{2})$, while a new pole with a small weight appears at the center-of-mass excitation energy at $\omega=(\pi 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 $\sqrt{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 $\sqrt{q}$. This result was to be expected since in the expansion of Eq.(41) we treated the frequency $\omega$ as a constant that is small compared $\sqrt{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.

IV TDDFT in the strong interaction limit

IV.1 The exchange-correlation kernel

In TDDFT an effective non-interacting 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];\mathbf{r}t)$, is a functional of the density Ruggenthaler et al. (2013, 2015) and is often written as follows

[TABLE]

Here $v(\mathbf{r}t)$ is the external potential of the interacting system of interest and $w(\mathbf{r},\mathbf{r}^{\prime})$ 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

[TABLE]

Here $\chi_{s}^{-1}$ is the inverse of the Kohn-Sham density response function whereas $\chi^{-1}$ is the inverse of the density response function of the interacting system and $f_{\textrm{Hxc}}$, the Hartree-xc kernel, is defined as

[TABLE]

where $v_{\textrm{Hxc}}$ is the sum of the Hartree and the xc-potential. Equation (45) is commonly used to calculate $\chi$ from the knowledge of $\chi_{s}$ at the price of approximating $f_{\textrm{Hxc}}$.

For the discussion in the next section it is important to note that the functional derivative is not a uniquely defined function Hellgren and Gross (2012); Ruggenthaler et al. (2013) 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.

[TABLE]

Let us define a new function

[TABLE]

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 $\delta n$ is given by

[TABLE]

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 that $\tilde{f}_{\textrm{Hxc}}$ and $f_{\textrm{Hxc}}$ are physically equivalent integral kernels. The quantity that is defined unambiguously 111Note that Ref. Lani et al. (2016) did not discuss the possible addition of arbitrary functions of one spatial variable to the kernel. is the mixed spatial derivative

[TABLE]

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:

[TABLE]

By Fourier transforming the kernel $f_{\textrm{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 non-interacting electrons of Eq.(34). Using this equation, together with the expansion of Eq.(41), we obtain an explicit expression for $f_{\textrm{Hxc}}$ in the strong interaction limit:

[TABLE]

where we have reintroduced the variable $\lambda$ rather than $q$ as the $\lambda$ 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_{\textrm{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_{\textrm{Hxc}}$ is accurate up to the lowest excitation energy corresponding to a change in the relative wave function. In this energy region $f_{\textrm{Hxc}}$ has a pole at an energy corresponding to a zero in $\chi$, as a consequence of Eq.(51). It is worth noticing that, since the non-interacting response function has no zeroes, there is no pole in $f_{\textrm{Hxc}}$ originating from the first term in Eq.(51). This is peculiar to our quantum ring system for which the non-interacting response function has only a single pole for $\omega>0$. Instead, for a general system, the Kohn-Sham response function will have multiple poles and zeroes which implies that, to order $\lambda^{0}$, $f_{\textrm{Hxc}}$ is frequency dependent causing the adiabatic approximation to fail in this order. In our system $f_{\textrm{Hxc}}$ for even $k$ tends to a frequency-independent function for all $\omega$ when the interaction strength approaches infinity. Its static value, given by $f_{\textrm{Hxc}}=-\frac{3L}{8}\left(\frac{\pi k}{L}\right)^{2}=\chi_{s}^{-1}\left[k,\omega=\left(\frac{\pi k}{L}\right)^{2}\right]$, is the value needed to shift the Kohn-Sham pole to the $q\to\infty$ pole.

For the odd $k-$values (not shown here), all poles in the response function shift to infinity as $q\to\infty$. The asymptotic expansion for $f_{\textrm{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 $\lambda^{0}$ are frequency independent. However, frequency-dependent terms will appear to order $\lambda^{-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 $\lambda^{0}$ in real space using the Fourier transform of Eq.(53)

[TABLE]

and $r=x-x^{\prime}$ is the relative distance between the points $x$ and $x^{\prime}$ which are the one-dimensional counterparts of the spatial points in Eq.(46). Since $x$ and $x^{\prime}$ are both in the interval from [math] to $L$ we have that $r\in[-L,L]$. We find that

[TABLE]

The leading term is given explicitly by

[TABLE]

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

[TABLE]

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.

IV.2 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) Seidl (1999); Seidl et al. (2007) 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 $\lambda$ can be written as Langreth and Perdew (1975):

[TABLE]

where we defined:

[TABLE]

In this expression $\hat{W}$ is the two-particle interaction and $\Psi_{\lambda}[n]$ is the ground-state wave function of a system with a local external potential, interaction $\lambda\hat{W}$ and ground-state density $n$. In the strong interaction limit $W_{\lambda}[n]$ can be expanded asGori-Giorgi et al. (2009),222Strictly speaking the absence of a term of order $\lambda^{-1}$ in this expansion has only been shown for a Coulomb system but we find it to be true in our model as well:

[TABLE]

where the first term is the SCE energy which has the explicit form:

[TABLE]

where $w$ is the two-particle interaction which we assume to depend only on the distance between the particles. The functions $\mathbf{f}_{i}([n];\mathbf{r})$ are the so-called co-motion functions, which specify the position of $N-1$ electrons given the position of one electron at $\mathbf{r}$. The second term in Eq.(59) contains the so-called zero-point energy (ZPE) which describes the vibrations of the electrons around their equilibrium positions and is given explicitly as:

[TABLE]

where $D$ is the spatial dimensionality of the system and $\omega_{n}$ are the harmonic frequencies. Inserting Eq.(59) into Eq.(57) we find the large $\lambda$ expansion of the Hartree-xc energy to be:

[TABLE]

The ground-state theory can be used in the adiabatic approximation Lani et al. (2016) to find an approximate exchange-correlation kernel from

[TABLE]

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^{\textrm{A}}_{\textrm{Hxc}}(\mathbf{r},\mathbf{r}^{\prime})$. The second order variation of Eq.(62) with respect to the density gives an expansion in orders of $\sqrt{\lambda}$ for the Hartree-xc kernel:

[TABLE]

where we defined the adiabatic SCE and ZPE kernels as

[TABLE]

Let us now turn again to the case of the quantum ring. In this case there are two electrons and just one simple co-motion function $\textrm{f}:[0,L]\rightarrow[0,L]$ given by

[TABLE]

If one electron is at $x$ this function simply puts the other electron at the antipodal point of the quantum ring. From this co-motion function it is straightforward to calculate the SCE energy. Since $|x-\textrm{f}(x)|=L/2$ we have from our interaction $w(x)=V_{0}\cos^{2}(\pi x/L)$ in Hamiltonian (1) that $w(|x-\textrm{f}(x)|)=0$ and therefore $V_{\textrm{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 non-zero harmonic frequency given by Malet et al. (2014)

[TABLE]

where $w^{\prime\prime}(x)=\partial^{2}_{x}w(x)$. If we calculate this frequency for our quantum ring we find $\omega_{1}(x)=2\pi\sqrt{V_{0}}/L$ and $V_{\textrm{ZPE}}=\pi\sqrt{V_{0}}/(2L)$. 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_{\textrm{Hxc}}^{\lambda}$ 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

[TABLE]

where the leading coefficient indeed exactly gives $2V_{\textrm{ZPE}}$. The expansion of $W_{\lambda}$ can be calculated from Eq.(57) to give

[TABLE]

This result agrees with a direct calculation of $W_{\lambda}$ 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 $\lambda^{-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. Lani et al. (2016) and reads

[TABLE]

where $f_{1}$ is the function given in Eq.(55). Because of the freedom in Eq.(48) $f^{\textrm{ASCE}}$ is physically equivalent to the kernel $f_{1}(x-x^{\prime})$ 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^{\prime})$ has a simpler shape, we will restrict ourselves to this function. The kernel describes how a density variation induces a change in $v_{\textrm{Hxc}}$. To leading order in $\lambda$ we have

[TABLE]

By taking the second derivative of this equation, we obtain

[TABLE]

where we stress that $\delta n(x,t)$ is periodic with $L$. Since $f_{1}(x-x^{\prime})$ is linear everywhere except at the kinks at $x-x^{\prime}=0,\pm L/2$, the second derivatives yield delta functions at these points. Also note that the sign of the density change yields the curvature of the induced potential. Eq.(72 ) has some interesting consequences. If we make a localized density variation $\delta n(x,t)$ in a very small interval of the ring, there will not only be a change $\delta v_{\textrm{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 non-locally, but instantaneously, on the density.

Let us analyze the next orders in the strong interaction expansion. The calculation of $f^{\textrm{AZPE}}$ allows for a comparison with the next leading term in Eq.(54) which is proportional to $\sqrt{\lambda}$. The kernel $f^{\textrm{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 co-motion functions (derived already in Ref. Lani et al. (2016)). We also find agreement between this expression and that of the next leading term in Eq.(56), i.e. $f^{\textrm{AZPE}}(x,x^{\prime})=f_{2}(x-x^{\prime})$ modulo the addition of arbitrary functions of $x$ and $x^{\prime}$ 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_{\textrm{Hxc}}$, all scaled by $\lambda$, 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 $\lambda$ 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 $\lambda$-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-$\lambda$ 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 non-uniform 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 non-adiabatic 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 developed an 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 non-local in space. This already shows that the commonly used adiabatic local-density, or semi-local, 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 Lani et al. (2016) which has the spatial non-locality built in. The leading term of the exact expansion was found to coincide with the adiabatic-SCE kernel derived in Ref. Lani et al. (2016). 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 non-adiabatic.

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.

Acknowledgements.

LC, DK and RvL acknowledge the Academy of Finland for support under project no. 267839. GL acknowledges the Université Pierre et Marie Curie through an ATER position, which granted her considerable flexibility and freedom to contribute to this work.

Appendices:

Appendix A Response function of a harmonic model system

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 repulsionMoshinsky et al. (1970). The Hamiltonian of the system is given by

[TABLE]

in which the harmonic frequency $\omega_{\lambda}$ is chosen in such a way that the density is independent of $\lambda$. Using the coordinate transformation $\mathbf{s}=(\mathbf{r}_{1}+\mathbf{r}_{2})/\sqrt{2}$ and $\mathbf{r}=(\mathbf{r}_{1}-\mathbf{r}_{2})/\sqrt{2}$ the Hamiltonian can be written as that of two independent harmonic oscillators

[TABLE]

where $\nu_{\lambda}^{2}=\omega_{\lambda}^{2}-2\lambda$. The eigenfunctions and eigenvalues of this Hamiltonian are well-known. The normalized eigenfunctions are given by

[TABLE]

where we defined the triplet of non-negative integers $\mathbf{n}=(n_{1},n_{2},n_{3})$ and $\mathbf{m}=(m_{1},m_{2},m_{3})$ and the functions

[TABLE]

in which we denoted

[TABLE]

where $|\mathbf{n}|=n_{1}+n_{2}+n_{3}$ and $H_{n}(x)$ is the Hermite polynomial of order $n$. The energy eigenvalues are given by

[TABLE]

The ground-state wave function $\Psi_{0}=\Psi_{\mathbf{0}\mathbf{0}}$ has the explicit form

[TABLE]

From this function the density is readily obtained as

[TABLE]

where we defined

[TABLE]

If we insert $\nu_{\lambda}=\sqrt{\omega_{\lambda}^{2}-2\lambda}$ into this relation, we can determine the $\lambda$-dependence of $\omega_{\lambda}$ as $\beta$ is independent of $\lambda$. We find

[TABLE]

where $y$ solves the quartic equation

[TABLE]

In the limit of large interaction we find that $y=1+\beta/(2\sqrt{2\lambda})+\mathcal{O}\left(\lambda^{-1}\right)$ such that

[TABLE]

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

[TABLE]

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 $|\mathbf{m}|$ only attains even values. The response function therefore has the form

[TABLE]

where $\Delta E_{\mathbf{n}\mathbf{m}}=E_{\mathbf{n}\mathbf{m}}-E_{\mathbf{00}}$ are the excitation energies and we further put the restriction that we sum over all $\mathbf{m}$ such that $|\mathbf{m}|$ is even. The excitation amplitudes corresponding to these excitations are given by

[TABLE]

where we rewrote the eigenfunctions in terms of the original coordinates. Let us now consider the large interaction limit $\lambda\rightarrow\infty$ of the response function. Since only the excitation energies of the form $\Delta E_{\mathbf{0}\mathbf{m}}$ remain finite in this limit we only need to consider the excitation amplitudes of the form

[TABLE]

If we now use that

[TABLE]

is a limit representation for the delta distribution and the fact that $\nu_{\lambda}\rightarrow\beta/2$ in the large interaction limit we find that

[TABLE]

For the response function in the large interaction limit we therefore obtain

[TABLE]

where the sum runs over $\mathbf{m}$ values such that $|\mathbf{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

[TABLE]

which implies $\nabla v_{\lambda}(\mathbf{r})=2\lambda\mathbf{r}+O(\sqrt{\lambda})$. Since the response function does not vanish in the large interaction limit, Eq.(20) tells us that

[TABLE]

must hold, where we took the static limit of the response function of Eq.(92). Since the polynomial functions in Eq.(92) are all even (since $|\mathbf{m}|$ is even), and $\mathbf{r}_{2}$ is an odd function, we find that this relation is indeed satisfied.

Appendix B Sips’ expansion of the Mathieu function

In this Appendix we describe the details of the expansion of the Mathieu functions for large $q$. A recent general discussion is given by Frenkel and Portugal Frenkel and Portugal (2001) who give an overview of various expansions for the Mathieu functions for large and small values of $q$ in different regions of their domain and recursion formulas to determine the expansion coefficients. The expansion that we are interested in for this work is the large $q$ expansion for the Mathieu cosine function $C_{l}(z;q)$ in the region enclosing the value $z=\pi/2$. Such an expansion was derived originally by Sips Sips (1949, 1959); Blanch (1960); DLMF who developed a systematic theory. For the Mathieu cosine this expansion is of the form of Eq.(13). We here give the explicit expressions for the coefficients $c_{2n,l}(q)$ in Eq.(13), and for this it will be convenient to define new coefficients $g_{2n,l}(q)$

[TABLE]

which only differ from the coefficients $c_{2n,l}(q)$ by a prefactor $\mathcal{C}_{l}(q)$. This is done to ensure that $g_{0,l}(q)=1$, which is convenient for a recursive calculation of the remaining coefficients $g_{2n,l}(q)$ as is done in Blanch (1960); Frenkel and Portugal (2001). The prefactor $\mathcal{C}_{l}(q)$ is chosen such that the Mathieu cosine satisfies the normalization of Eq.(8). It has the explicit asymptotic expansion Sips (1949); Blanch (1960); Frenkel and Portugal (2001)

[TABLE]

In terms of the new coefficients $g_{2n,l}(q)$ the Sips expansion of Eq.(13) becomes Blanch (1960)

[TABLE]

where the functions $\mathcal{D}_{l}(u)$ are defined in Eq.(12) and $u=\sqrt{2}\,q^{1/4}\cos z$. Note that the expansion here differs in the choice of argument $u$ compared to that of Ref. DLMF by a factor of $\sqrt{2}$ as we preferred to use the physicist’s convention for the Hermite polynomials appearing in the harmonic oscillator functions $\mathcal{D}_{l}$. We prefer here to avoid a general discussion on the determination of $g_{2n,l}(q)$ and refer the interested reader to Refs. Blanch (1960); Frenkel and Portugal (2001) for details. Instead we give the explicit forms of the coefficients $g_{2n,l}(q)$ up to order $1/q$ which is sufficient for this work:

[TABLE]

where we define $\binom{n}{m}=0$ if $n<m$. From the formulas given in Ref. Frenkel and Portugal (2001) it is readily seen that to know the Sips expansion to order $q^{-k/2}$ we need to know the coefficients $g_{2n,l}(q)$ for $n\in\{-2k,\ldots,2k\}$. For example, the knowledge of the next order $q^{-3/2}$ would need the knowledge of the coefficients from $g_{-12,l}(q)$ up to $g_{12,l}(q)$, see for example Blanch (1960).

To show the accuracy of the Sips expansion we display in Fig.6 the Mathieu cosine $C_{l}(q;z)$ for $l=0$ and $l=3$ for $q=1$ and $q=10$ and compare them to the Sips expansion using the coefficients $g_{2n,l}(q)$ given explicitly above. We see that already at $q=10$ the Sips expansion performs very well. In general we find that for higher values of $l$ more terms in the expansion need to be taken into account for high accuracy.

It remains to give the asymptotic values of the Mathieu characteristic value. The characteristic values $a_{l}^{+}(q)$ and $a_{l}^{-}(q)$ have the same asymptotic expansion of the form Frenkel and Portugal (2001)

[TABLE]

while the difference $a_{l}^{+}(q)-a_{l}^{-}(q)$ is exponentially small in the large $q$ limit DLMF .

Appendix C Expansion of the excitation amplitudes

Here we further outline some general features of the expansion of the excitation amplitudes $D_{kl}(q)$. If we rewrite the expansion of Eq.(35) using the coefficients of Eq.(94) we have the expression

[TABLE]

For the products of the pre-factors we can write

[TABLE]

where $F_{l}(q)$ has an expansion in powers of $q^{-1/2}$, i.e.

[TABLE]

and higher powers can be calculated from expression (95). With these definitions and Eq.(39) we can rewrite the expansion of Eq.(98) as

[TABLE]

where the factor $q^{1/4}$ from Eq.(99) has been absorbed in the last sum. The function $F_{l}(q)$ and the coefficients $g_{2n,l}(q)$ have an expansion in powers of $q^{-1/2}$. Now since the coefficients $I_{n_{1}n_{2},r}^{l}$ vanish unless $k$ and $l$ are both even or both odd, we conclude that $D_{kl}(q)$ has an expansion only in odd powers of $q^{-1/4}$ if $l$ is odd and only in even powers of $q^{-1/4}$ otherwise. The explicit expression for $I_{n_{1}n_{2},r}^{l}$ is given by

[TABLE]

where we defined $s=r/2-n_{1}-n_{2}-l/2$ Weisstein . Finally the coefficients $a_{r}(k)$ can be obtained from a Taylor expansion of the function $f_{k}(x)$ defined in Eq. (38). Taken all these terms together we find the following expansion of the excitation amplitude:

[TABLE]

The excitation amplitude indeed has an expansion in odd powers of $q^{-1/4}$ for odd $l$ and in even powers of $q^{-1/4}$ for even $l$, as we demonstrated above. This implies that $|D_{kl}(q)|^{2}$ has an expansion in powers of $q^{-1/2}$. The expansion derived here has been checked numerically to ensure the correctness of the derivations. As a further check on the result we see that the excitation amplitude satisfies $D_{kl}(q)=(-1)^{k+l}D_{kl}^{*}(q)$ as well as $D_{0l}(q)=\delta_{l0}$.

The order for $D_{kl}(q)$ given in Eq.(105) is enough to yield the following asymptotic expansion for the absolute value squared of the excitation amplitudes to order $q^{-3/2}$

[TABLE]

We compare the expansion of $|D_{kl}(q)|^{2}$ to the exact values obtained by numerically integrating Eq.(26) in Fig.7. We see from this figure that higher values of $k$ require higher values of $q$ to make the expansion accurate. This was to be expected since larger $k$-values imply a more oscillatory integrand in Eq.(26) while a larger value of $q$ makes the Mathieu functions more localized and thereby making 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.(106). 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

[TABLE]

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.