Exact exchange-correlation kernels for optical spectra of model systems
EExact exchange-correlation kernels for optical spectra of model systems
M. T. Entwistle and R. W. Godby
Department of Physics, University of York, and European TheoreticalSpectroscopy Facility, Heslington, York YO10 5DD, United Kingdom (Dated: March 22, 2019)For two prototype systems, we calculate the exact exchange-correlation kernels f xc ( x, x (cid:48) , ω ) oftime-dependent density functional theory. f xc , the key quantity for optical absorption spectra ofelectronic systems, is normally subject to uncontrolled approximation. We find that, up to the firstexcitation energy, the exact f xc has weak frequency-dependence and a simple, though non-local,spatial form. For higher excitations, the spatial behavior and frequency-dependence become morecomplex. The accuracy of the underlying exchange-correlation potential is of crucial importance. Time-dependent Kohn-Sham density functional theory[1, 2] (TDDFT) is in principle an exact and efficienttheory of the excited-state properties of many-electronsystems, including a wide variety of important spectro-scopies such as optical absorption spectra of moleculesand solids. However its application is restricted by thelimitations of the available approximate functionals forelectron exchange and correlation – in particular theexchange-correlation kernel, f xc , the functional deriva-tive of the exchange-correlation potential with respect tothe electron density. To assist the construction of morepowerful approximations for f xc , we calculate the exact f xc for small prototype systems, and analyze its char-acter, including key aspects in which it differs from thecommon approximations.In the Runge-Gross formulation [1] of TDDFT the realsystem of interacting electrons is mapped onto an auxil-iary system of noninteracting electrons moving in an ef-fective local Kohn-Sham (KS) potential v KS = v ext + v H + v xc , with both systems having the same electron density n at all points in space and time. Many TDDFT calcu-lations are done within the framework of linear responsetheory, which describes how a system responds upon ap-plication of a weak, time-dependent external perturba-tion. The induced density is described by the interact-ing density-response function, the functional derivative χ = δn/δv ext . χ is related to the non-interacting density-response function of the KS system, χ = δn/δv KS , bythe Dyson equation [3] [4] χ = χ + χ ( u + f xc ) χ , where u is the bare Coulomb interaction. χ is to be obtainedfrom a ground-state DFT calculation. χ can then be usedto compute, for example, the optical absorption spectrumof the system, σ ( ω ) = − πωc (cid:90) (cid:90) Im (cid:0) χ ( x, x (cid:48) , ω ) (cid:1) x x (cid:48) dx dx (cid:48) . (1)In practice, both v xc and its functional derivative f xc must be approximated. While there have been somesuccesses, the commonly used adiabatic TDDFT func-tionals, such as the adiabatic local density approxima-tion [5, 6] (ALDA), often fail in extended systems. Forexample, the optical absorption spectra of many semi-conductors and insulators are not even qualitatively de- scribed, with excitonic effects and many-electron excita-tions omitted [7, 8], and the optical gap underestimated.Here, approximations for f xc achieve little improvementover the random phase approximation (RPA), in which f xc is neglected entirely [9]. Attempts to improve approx-imations for f xc include exact-exchange methods [10–15], diagrammatic expansions using perturbative meth-ods [16, 17], and adding long-range contributions [18–21].Another approach involves calculations of the homoge-neous electron gas (HEG) [22–26]. Kernels derived fromthe Bethe-Salpeter equation [27–31] have had some suc-cess, but require a relatively expensive many-body per-turbation theory calculation as their input, and are out-side the KS TDDFT framework.There have been a limited number of studies conductedon analyzing the character of the exact f xc , all of whichfocus on its frequency-dependence. One approach hasbeen to calculate the exact adiabatic f xc for model sys-tems [32], in order to investigate its performance uponapplication and deduce when memory effects become im-portant. This approach has been used in simple Hub-bard systems [33, 34] and extended by analyzing addi-tional properties, such as the frequency-dependence ofthe full f xc around double excitations. Other researchhas explored how this frequency-dependence of f xc turnsthe single-particle quantities of exact KS TDDFT intomany-body excitations [35] and its behaviour for long-range excitations has been analyzed in order to developapproximate kernels [36].In this paper, we explore the properties of exact xckernels, including full spatial and frequency-dependence,in order to inform the development of improved approx-imate functionals. We employ our iDEA code [37] whichsolves the many-electron Schr¨odinger equation exactlyfor small, one-dimensional prototype systems [38] [39].From the many-electron eigenstates of the system we cal-culate the exact χ using the Lehmann representation, χ ( x, x (cid:48) , ω ) = (cid:88) m (cid:20) (cid:104) | ˆ n ( x ) | m (cid:105)(cid:104) m | ˆ n ( x (cid:48) ) | (cid:105) ω − ( E m − E ) + iη + c.c. ( − ω ) (cid:21) , (2)where | (cid:105) , E , | m (cid:105) and E m are the ground state andits energy, and the m − th excited state and its energy, a r X i v : . [ c ond - m a t . o t h e r] M a r respectively, ˆ n is the density operator in the Heisenbergpicture and η is a positive infinitesimal. χ has poles atthe excitation energies of the system, E m − E . It isconvenient to calculate Im( χ ), with Re( χ ) following fromthe Kramers-Kronig relations [40, 41]. As is customary,we replace η with a small positive number, to broadenthe absorption peaks for ease of viewing.We then determine the exact KS potential through ourreverse-engineering algorithm [42]. From the exact KSorbitals, we calculate the exact non-interacting density-response function, χ ( x, x (cid:48) , ω ) = (cid:88) i,j ( f i − f j ) φ ∗ i ( x ) φ j ( x ) φ ∗ j ( x (cid:48) ) φ i ( x (cid:48) ) ω − ( ε j − ε i ) + iη , (3)where the φ i , ε i are the exact solutions to the Kohn-Sham equations of ground-state DFT, and f i is the Fermioccupation (0 or 1) of φ i .The Dyson equation may be manipulated into an ex-pression for f xc , f xc = χ − − χ − − u, (4)but the inverses of χ and χ are not well defined. Forinstance, a spatially uniform perturbation of any angu-lar frequency induces no change in density, so both χ and χ have a zero eigenvalue and therefore a zero de-terminant. To overcome this, we find a pseudo-inverseof χ using truncated singular-value decomposition, dis-carding those eigenvectors with eigenvalues smallest inmagnitude, which we term the eigenvalue cutoff. Thisprocedure is repeated for χ , discarding the same numberof eigenvectors. From the modified response functions akernel f xc is now well defined. We confirm the validityof this procedure by verifying that the calculated f xc ,together with the unmodified χ , closely reproduces theunmodified χ via the Dyson equation. Additionally, weensure that the zero-force sum rule is obeyed – a wellknown property of the exact f xc [43].We begin by considering a system of two interactingelectrons confined to a harmonic well potential ( ω =0 .
25 a.u.), where ω is the angular frequency of the well(inset of Fig. 1). We compute the exact optical absorp-tion spectrum of the system [44], with the first excita-tion at ω = ω (Fig. 1). Additionally, we compute theabsorption spectrum of the exact Kohn-Sham system, inwhich the absorption frequency is slightly too low ( ≈ . exact χ . This last point providesa strong reminder of the challenge of f xc : starting fromthe exact Kohn-Sham orbitals, a much better absorptionpeak is obtained by ignoring the induced changes in theHartree and xc potentials ( χ ) than by accounting for thefirst exactly and either neglecting (RPA) or approximat-ing (ALDA) the second. This highlights the importance of obtaining a good approximation to the ground-statexc potential v xc , which leads to χ . .
00 0 .
05 0 .
10 0 .
15 0 .
20 0 .
25 0 .
30 0 . ! (a.u.) ( a . u . ) ExactKS RPAALDA x (a.u.) . . . . . n , V e x t , V K S ( a . u . ) FIG. 1. Two interacting electrons in a harmonic potential.The inset shows the electron density (solid blue), along withthe external (dotted-dashed green) and exact Kohn-Sham(dashed purple) potentials. In the main panel, the absorp-tion spectra (detailing the first excitation) of the exact andKohn-Sham systems, along with the RPA and ALDA approxi-mations. We check that the calculated f xc is correct by solvingthe Dyson equation and comparing the resultant absorptionspectrum (short-dashed black) with the exact. We now turn to the spatial characteristics of f xc (Fig. 2). Typically, several different choices of the eigen-value cutoff yield kernels f xc with varying degrees of spa-tial structure, all of which essentially yield the correct χ from the exact χ as set out above. Of these, we se-lect the eigenvalue cutoff with the largest magnitude, re-sulting in the smoothest possible spatial structure with-out detriment to the exact absorption spectrum (Fig. 1).We observe that while f xc has real and imaginary parts(see later), the real part alone is sufficient to reproducethe position and weight of the first excitation ( ω = ω ).Fig. 2(a) and Fig. 2(b) show Re( f xc ) at ω = 0 and ω ,respectively. The behavior of f xc away from the diago-nal, x (cid:54) = x (cid:48) , represents the kernel’s non-locality, and itis evident that this is fairly simple in nature; analysis(Fig 2(c)) shows it to be similar to the negative of theCoulomb interaction, with which it therefore tends tocancel in the expression for χ . The ω -dependence of f xc up to the first excitation is seen to be extremely weak,as observed in other model systems [33–35]. We analyzethis more closely in the inset to Fig. 2(c).To gain insight into these observations, we analyze theexact χ and χ . Fig. 3 shows Re( χ ) and Re( χ ); upto the first excitation, these exhibit strong – but closelysimilar – ω -dependence. The similarity arises in partfrom the exact many-electron wavefunction being well ap-proximated by the exact Kohn-Sham wavefunction [46],which reflects the dominance of exchange (including self- x (a.u.) x ( a . u . ) ω = 0 x (a.u.) ω = 0.25 . . . . . x x (a.u.) . . . . . . f x c ( a . u . ) . . . . ω (a.u.) − . − . − . (a) (b) (c) x = 0 x = 2 x = 4 − u Re( f xc ) FIG. 2. The real part of the exact f xc of the harmonic wellsystem: (a) In the adiabatic limit ( ω = 0), and (b) at the firstexcitation ( ω = 0 . f xc has a rather simple non-localdependence, which is similar to the negative of the Coulombinteraction u . Here we focus on f xc at ω = 0 .
25. Inset: Weobserve f xc to have strikingly weak ω -dependence up to thefirst excitation (vertical line). We illustrate this by plotting itsvalue (solid red) at a point along the main diagonal ( x = x (cid:48) )which corresponds to the peak in electron density in the insetof Fig. 1 ( x = 1 . interaction correction) in the harmonic potential system[47]. Therefore χ − and χ − largely cancel, so that f xc is similar to − u , with weak ω -dependence.This can be demonstrated succinctly through a simplemodel, in which we take Eq. (2) and Eq. (3), and re-place the spatially-dependent numerators (the oscillatorstrengths) with scalars. Specifically, we consider a sys-tem with a single excitation at ω = 1, set the numeratorequal to 1, and let η = 0 .
05 (top inset of Fig. 4). Wedo the same for the Kohn-Sham system, but choose theexcitation to occur at ω = 0 .
9. By taking their inverses,we calculate Re( χ − − χ − ), which is the ω -dependentpart of Re( f xc ) in Eq. (4), and find this to be small inamplitude and have a fairly weak ω -dependence up to thefirst excitations (bottom inset of Fig. 4). The inclusionof higher excitations, and taking the limit η →
0, changelittle at these low frequencies.Including higher excitations in the model χ causesRe( χ ) to pass through zero between excitations. At thesepoints Re( χ − ) also passes through zero, and Im( χ − ) . . . . . x ( a . u . ) ω = 0 . . . . . x ( a . u . ) ω = 0.245 . . . . . x (a.u.) . . . . . x ( a . u . ) ω = 0.255 ω = 0 ω = 0.236 . . . . . x (a.u.) ω = 0.246 Re( χ ) Re( χ ) . . . . . . FIG. 3. The exact χ and χ in the harmonic well system: Left:Re( χ ) at ω = 0 and on either side of the transition at 0.250.Right: Re( χ ) at ω = 0 and on either side of the transition at0.241. Re( χ ) and Re( χ ) exhibit remarkably similar spatialstructure. peaks sharply (Fig. 4). As Im( f xc ) = Im( χ − − χ − ),we find that the f xc in our simple model only has animaginary component when χ or χ passes through zerobetween excitations, and hence is completely real up tothe first excitations (as η → f xc ) wasvery small up to the first excitations, and Re( f xc ) wassufficient to reproduce the peak in the absorption spec-trum.We now consider a system whose absorption spectrumincludes higher excitations – two interacting electrons ina softened atomic-like potential (top inset of Fig. 5). Asin the harmonic well system, the absorption spectrumof the exact Kohn-Sham system is slightly too low forthe first excitation (Fig. 5). Again, we find f xc to bedominated by its real part and nearly ω -independent,while exhibiting relatively simple spatial structure, up toand including the first excitation (Fig. 6(a)). The secondexcitation does not appear in the absorption spectrum,and so we move to the third excitation, which is muchsmaller in amplitude than the first, and once more ob-serve the peak in the Kohn-Sham system to be slightly ! (a.u.) ( a . u . ) ReIm ! (a.u.) . . . . ! (a.u.) . . . . R e ( , ) FIG. 4. The top inset shows χ for a simple model withoutspatial dependence and a single excitation at ω = 1; the bot-tom inset shows the near cancellation (solid green) betweenRe( χ − ) (dashed red) and Re( χ − ) (dotted dark red), where χ has an excitation at 0.9, causing f xc to exhibit weak ω -dependence. In the main panel, two further excitations at ω = 2 and 3 have been included, to show that Re( χ − )passes through zero between excitations, which leads to anon-zero Im( f xc ), as does the corresponding feature in χ − (not shown). below but still very close to the exact (bottom inset ofFig. 5). Again, the closeness between the two peaks arisesfrom the strong similarity between χ and χ . In order toreproduce this excitation, a smaller eigenvalue cutoff isrequired, leading to higher spatial frequencies in f xc [48](Fig. 6(b)).For the atom-like system, we have investigated the ex-tent to which local kernel approximations for f xc maybe meaningful. As we have observed f xc to largely can-cel with u at low ω , we choose to focus on the Hartreeexchange-correlation kernel f Hxc = f xc + u which is morelocal. We incorporate the nonlocal parts of f Hxc by pro-jecting them onto a local kernel [49]. We find this largelycorrects the difference χ − χ , and hence the position ofthe peak in the absorption spectrum, for the first exci-tation, but fails to correct the height of the peak. Sucha local kernel is completely inadequate to describe thethird excitation.In summary, we have calculated the exact f xc ( x, x (cid:48) , ω )for two prototype systems. At low ω , we find the imag-inary component of f xc to be small, with the real partalone sufficient to reproduce the first excitation. Up toand including the first excitation, Re( f xc ) exhibits strik-ingly weak ω -dependence, stemming from strong, butclosely similar ω -dependence between the interacting andnon-interacting density-response functions – boding wellfor the applicability of adiabatic kernels. Additionally,Re( f xc ) here has a rather simple spatial form, which issimilar to the negative of the Coulomb interaction u , in- .
00 0 .
02 0 .
04 0 .
06 0 .
08 0 .
10 0 . ! (a.u.) ( a . u . ) ExactKS .
10 0 .
11 0 .
12 0 . ! (a.u.) . . . ( a . u . ) x (a.u.) . . . . n ( a . u . ) . . . V e x t , V K S ( a . u . ) FIG. 5. Two interacting electrons in an atomic-like poten-tial. The top inset shows the electron density (solid blue),along with the external (dotted-dashed green) and Kohn-Sham (dashed purple) potentials. In the main panel, theabsorption spectra of the exact and Kohn-Sham systems; thebottom inset shows the third excitation (fourth in the KS sys-tem) in more detail, which is the next to appear after the firstexcitation.
10 0 10 20 x (a.u.) x ( a . u . ) ω = 0.045 . . . . .
10 0 10 20 x (a.u.) ω = 0.119 (a) (b) Re( f xc ) FIG. 6. The real part of the exact f xc of the atom-like sys-tem: (a) At the first excitation ( ω = 0 . ω = 0 . f xc to be nearly ω -independent and exhibit a rela-tively simple spatial form up to the first excitation. However,more complex spatial structure is needed to capture higherexcitations. dicating that approximations to f Hxc may be more appro-priate than those for f xc alone. For higher excitations, f xc exhibits both additional spatial structure and stronger ω -dependence, indicating that more sophisticated approx-imations are needed. Throughout, the absorption spec-trum of the exact Kohn-Sham system provides a verygood starting point, signifying the crucial importance ofan accurate approximation for the ground-state v xc .We thank Phil Hasnip for helpful discussions. Datacreated during this research is available from the YorkResearch Database [50]. [1] E. Runge and E. K. U. Gross, Phys. Rev. Lett. , 997(1984).[2] E. K. U. Gross, J. F. Dobson, and M. Petersilka, “Den-sity functional theory of time-dependent phenomena,” in Density Functional Theory II: Relativistic and Time De-pendent Extensions , edited by R. F. Nalewajski (SpringerBerlin Heidelberg, Berlin, Heidelberg, 1996) pp. 81–172.[3] Matrix multiplication for the spatially non-local quanti-ties χ , χ and f xc , and ω -dependence, are implied.[4] M. Petersilka, U. J. Gossmann, and E. K. U. Gross,Phys. Rev. Lett. , 1212 (1996).[5] A. Zangwill and P. Soven, Phys. Rev. A , 1561 (1980).[6] E. K. U. Gross and W. Kohn, Phys. Rev. Lett. , 2850(1985).[7] C. Jamorski, M. E. Casida, and D. R. Salahub,The Journal of Chemical Physics , 5134 (1996),https://doi.org/10.1063/1.471140.[8] N. T. Maitra, F. Zhang, R. J. Cave, and K. Burke,The Journal of Chemical Physics , 5932 (2004),https://doi.org/10.1063/1.1651060.[9] V. I. Gavrilenko and F. Bechstedt, Phys. Rev. B , 4343(1997).[10] Y.-H. Kim and A. G¨orling, Phys. Rev. Lett. , 096402(2002).[11] Y.-H. Kim and A. G¨orling, Phys. Rev. B , 035114(2002).[12] A. G¨orling, Phys. Rev. A , 3433 (1998).[13] M. Hellgren and U. von Barth, The Jour-nal of Chemical Physics , 044110 (2009),https://doi.org/10.1063/1.3179756.[14] A. Ipatov, A. Heßelmann, and A. G¨orling, Interna-tional Journal of Quantum Chemistry , 2202 (2010),https://onlinelibrary.wiley.com/doi/pdf/10.1002/qua.22561.[15] A. G¨orling, Phys. Rev. Lett. , 5459 (1999).[16] I. V. Tokatly and O. Pankratov, Phys. Rev. Lett. ,2078 (2001).[17] I. V. Tokatly, R. Stubner, and O. Pankratov, Phys. Rev.B , 113107 (2002).[18] S. Botti, F. Sottile, N. Vast, V. Olevano, L. Reining, H.-C. Weissker, A. Rubio, G. Onida, R. Del Sole, and R. W.Godby, Phys. Rev. B , 155112 (2004).[19] S. Sharma, J. K. Dewhurst, A. Sanna, and E. K. U.Gross, Phys. Rev. Lett. , 186401 (2011).[20] P. E. Trevisanutto, A. Terentjevs, L. A. Constantin,V. Olevano, and F. D. Sala, Phys. Rev. B , 205143(2013).[21] S. Rigamonti, S. Botti, V. Veniard, C. Draxl, L. Reining,and F. Sottile, Phys. Rev. Lett. , 146402 (2015).[22] V. U. Nazarov, I. V. Tokatly, S. Pittalis, and G. Vignale,Phys. Rev. B , 245101 (2010).[23] Z. Qian and G. Vignale, Phys. Rev. B , 235121 (2002).[24] C. F. Richardson and N. W. Ashcroft, Phys. Rev. B ,8170 (1994).[25] S. Conti, R. Nifos`ı, and M. P. Tosi, Journal of Physics: Condensed Matter , L475 (1997).[26] M. Panholzer, M. Gatti, and L. Reining, Phys. Rev.Lett. , 166402 (2018).[27] G. Onida, L. Reining, and A. Rubio, Rev. Mod. Phys. , 601 (2002).[28] L. Reining, V. Olevano, A. Rubio, and G. Onida, Phys.Rev. Lett. , 066404 (2002).[29] F. Sottile, V. Olevano, and L. Reining, Phys. Rev. Lett. , 056402 (2003).[30] G. Adragna, R. Del Sole, and A. Marini, Phys. Rev. B , 165108 (2003).[31] A. Marini, R. Del Sole, and A. Rubio, Phys. Rev. Lett. , 256402 (2003).[32] M. Thiele and S. K¨ummel, Phys. Rev. A , 012514(2009).[33] F. Aryasetiawan and O. Gunnarsson, Phys. Rev. B ,165119 (2002).[34] D. J. Carrascal, J. Ferrer, N. Maitra, and K. Burke, TheEuropean Physical Journal B , 142 (2018).[35] M. Thiele and S. K¨ummel, Phys. Rev. Lett. , 083001(2014).[36] N. T. Maitra and D. G. Tempel, The Jour-nal of Chemical Physics , 184111 (2006),https://doi.org/10.1063/1.2387951.[37] M. J. P. Hodgson, J. D. Ramsden, J. B. J. Chapman,P. Lillystone, and R. W. Godby, Phys. Rev. B , 241102(2013).[38] We perform calculations for systems of two spinless elec-trons interacting via the appropriately softened Coulombrepulsion [51] u ( x, x (cid:48) ) = ( | x − x (cid:48) | + 1) − , and work inHartree atomic units: m e = (cid:126) = e = 4 πε = 1.[39] See Supplemental Material at ‘Link’ for the parametersof the model systems, and details on our calculations toobtain converged results.[40] D. Pines, Elementary excitations in solids: lectures onphonons, electrons, and plasmons , Lecture notes and sup-plements in physics (W. A. Benjamin, 1964).[41] M. Marder,
Condensed Matter Physics (Wiley, 2010).[42] J. D. Ramsden and R. W. Godby, Phys. Rev. Lett. ,036402 (2012).[43] See Supplemental Material for more details.[44] For this harmonic well system, at the level of linear re-sponse theory, only one excitation appears in the absorp-tion spectrum.[45] M. T. Entwistle, M. Casula, and R. W. Godby, Phys.Rev. B , 235143 (2018).[46] We define this as a Slater determinant of the occupiedKS orbitals.[47] E. Richardson, private communication.[48] As expected for higher energy excited states.[49] We fold the nonlocal f Hxc with an envelope function thatsuppresses the more distant nonlocal parts and projectsthe remainder onto the diagonal x = x (cid:48) .[50] M. T. Entwistle and R. W. Godby,http://dx.doi.org/10.15124/56f576ee-b2de-4ca5-9251-831bfc3cae6f (2019).[51] A. Gordon, R. Santra, and F. X. K¨artner, Phys. Rev. A72