Progress towards understanding ultranonlocality through the wave-vector and frequency dependence of approximate exchange-correlation kernels

PHYSICAL REVIEW B 104, 125112 (2021) Progress towards understanding ultranonlocality through the wave-vector and frequency dependence of approximate exchange-correlation kernels Niraj K. Nepal,1Aaron D. Kaplan ,1J. M. Pitarke,2,3and Adrienn Ruzsinszky1,* 1Department of Physics, Temple University, Philadelphia, Pennsylvania 19122, USA 2CIC nanoGUNE BRTA and DIPC, E-20018 Donostia, Basque Country, Spain 3Materia Kondentsatuaren Fisika Saila and Centro Fisica Materiales CSIC-UPV/EHU, E-48080 Bilbao, Basque Country, Spain (Received 24 May 2021; revised 30 July 2021; accepted 23 August 2021; published 8 September 2021; corrected 16 December 2021) In the framework of time-dependent density functional theory (TDDFT), the exact exchange-correlation (xc) kernel fxc(n,q,ω) determines the ground-state energy, excited-state energies, lifetimes, and the time-dependent linear density response of any many-electron system. The recently developed MCP07 xc kernel fxc(n,q,ω)of Ruzsinszky et al. [Phys. Rev. B 101, 245135 (2020)] yields excellent uniform electron gas (UEG) ground-state energies and plausible plasmon lifetimes. As MCP07 is constructed to describe fxc of the UEG, it cannot capture optical properties of real materials. To verify this claim, we follow Nazarov et al. [Phys. Rev. Lett. 102, 113001 (2009)] to construct the long-range, dynamic xc kernel, limq→0fxc(n,q,ω)=−α(ω)e2/q2,ofa weakly inhomogeneous electron gas, using MCP07 and other common xc kernels. The strong wave-vector and frequency dependence of the “ultranonlocality” coefficient α(ω) is demonstrated for a variety of simple metals and semiconductors. We examine how imposing exact constraints on an approximate kernel shapes α(ω). Comparisons to kernels derived from correlated-wave-function calculations are drawn. DOI: 10.1103/PhysRevB.104.125112 I. INTRODUCTION By virtue of the Runge-Gross theorem [1], time-dependent density functional theory (TDDFT) is an extension of groundstate density functional theory (DFT) [2–5]. TDDFT is a more computationally feasible approach to compute excitation energies, compared to approaches based on many-body techniques [6,7]. The response of a many-electron system to a dynamic, external potential is characterized by a change in its charge density. If the amplitude of the external potential is small, then it exerts only a weak perturbation, and the density response can be taken to be linear. Through the linear density-density response function, linear response TDDFT (LR-TDDFT) [8] can predict the transition frequencies to electronic excited states, among other properties. The central equation of LR-TDDFT is a Dyson-like equation linking the true interacting density-response function of an arbitrary many-electron system χ(r,r;ω) to its noninteracting Kohn-Sham (KS) counterpart χ0(r,r;ω): χ(r,r;ω)=χ0(r,r;ω)+d3r d3rχ0(r,r;ω) ×[v(r,r)+fxc(r,r;ω)]χ(r,r;ω),(1) where v(r,r) represents the Coulomb interaction, and fxc(r,r;ω) is the exchange-correlation (xc) kernel. The TDDFT excitation energies are found as poles of χ(r,r;ω). *Corresponding author: [email protected] With fxc =0, χreduces to the random-phase approximation (RPA) [9]. Even if TDDFT is computationally more efficient than wave-function theories, its accuracy is restricted by the limitations of commonly used xc kernels [10–13]. The most widely used adiabatic local-density approximation (ALDA), for example, misses excitonic effects completely [14]. Recently developed model xc kernels describe the optical absorption spectra of small and medium band-gap semiconductors well [10,13]. Nonempirical xc kernels, which are not fitted to experimental or correlated-wave-function data, are typically insufficient to describe excited-state properties of semiconductors and insulators, especially large-gap materials. However, some modern kernels predict accurate properties of excitons in real materials [15–18]. In order to achieve this goal, the wave-vector and frequency dependence of the xc kernel fxc needs to be investigated further [19–21]. Recently, we proposed a parametrized fxc for a uniform electron gas (UEG) based upon satisfaction of exact constraints [22]. This model, called MCP07, depends upon both the wave-vector qand frequency ω. In the static (ω=0) limit, we modified the model of Constantin and Pitarke [12]as follows: fMCP07 xc (q,0) =4π q2B[e−kq2(1 +Eq4)−1] −4π k2 F C 1+(kq2)−2.(2) For definitions of the rs-dependent functions k,B,C, and E, refer to Ref. [22]. rs=(4πn/3)−1/3is the radius of a sphere containing, on average, one electron in a UEG of 2469-9950/2021/104(12)/125112(9) 125112-1 ©2021 American Physical Society NEPAL, KAPLAN, PITARKE, AND RUZSINSZKY PHYSICAL REVIEW B 104, 125112 (2021) density n.kF=(3π2n)1/3is the Fermi wave-vector. In the long-wavelength (q→0) limit, Eq. (2) yields the ALDA kernel, for which we use the Perdew-Zunger parametrization of the UEG correlation energy [23]. At finite frequencies, we use the long-wavelength model fxc(0,ω) of Gross, Kohn, and Iwamoto (GKI) [2,24] (analytically continued to complex frequencies when needed). This expression for fxc(0,ω) is then combined with the static limit of Eq. (2) to obtain [22] fMCP07 xc (q,ω)=1+e−kq2fxc(0,ω) fxc(0,0) −1fMCP07 xc (q,0). (3) In MCP07, one takes k=k, where kis the same function appearing in Eq. (2). Here we also discuss the effect of setting k=0. At long wavelengths (q→0), Eq. (3) yields the local GKI fxc(0,ω). At ω=0, Eq. (3) yields the nonlocal fMCP07 xc (q,0) of Eq. (2). The full MCP07 kernel of Eq. (3) provides a highly accurate description of ground-state correlation energies (up to rs=10, see Figs. 6and 7) and quasiparticle properties of the UEG. In particular, MCP07 predicts a finite plasmon lifetime that first decreases from infinity, and then increases as qgrows from 0 towards the electron-hole continuum. MCP07 also yields a static charge-density wave at rs69 that can be associated with a softening of the plasmon mode [25,26]. The exchange-only version of the static MCP07 xc kernel of Eq. (2) confirms Overhauser’s prediction [27] that correlation is essential for the creation of a charge-density wave. As MCP07 approximates fxc of the UEG, it cannot describe the 1/q2long-wavelength behavior of the exact fxc for nonuniform many-electron systems, called “ultranonlocality” [20]. Indeed, in reciprocal space, the kernel is a matrix characterized by reciprocal wave-vectors, whose spatial decay manifests in possible leading q-independent terms, called crystal local-field effects. The head and wings of adiabatic kernels, derived from semilocal density functional approximations, are independent of q, and thus are incorrectly nondivergent for q→0, as in the case of the UEG. The socalled “bootstrap” idea represents a very unique and effective route to account for ultranonlocality [28]. This work builds upon that of Nazarov et al. [29]by constructing, in the optical limit, the dynamic xc kernel of a weakly inhomogeneous electron gas using the MCP07 xc kernel (among other UEG-based kernels) as input. As MCP07 simultaneously describes the wave-vector and frequency dependence of the xc kernel, our approach serves as a basis for further investigations of the optical absorption of nonuniform systems and the real wave-vector and frequency dependence of xc kernels. II. WHAT IS KNOWN ABOUT THE OPTICAL LIMIT OF THE XC KERNEL? Giuliani and Vignale [30] provided a detailed discussion of exchange and correlation in uniform and nonuniform manyelectron systems. They discussed, in particular, the important difference between the short-range kernel of the uniform electron gas, where fxc(q,ω) tends to a finite constant as q→0, and the ultranonlocal kernel of nonuniform systems, where fxc is known to diverge in the long-wavelength limit as 1/q2. Concurrent works have generalized kernel development from density to current-density functionals [31]. Nazarov et al. [32] derived a general method for constructing a scalar TDDFT xc kernel using the tensorial kernel and KS currentdensity response function of time-dependent current-density functional theory (TDCDFT) [31,33]. This method proved to be particularly useful, as a local approximation to the xc kernel of TDCDFT results in a nonlocal approximation to the xc kernel of TDDFT. The resultant approximation is free of the contradictions that plague the standard local density approximation to TDDFT. This method also allowed the construction of the frequency-dependent xc kernel of a weakly inhomogeneous electron gas in the optical limit [29]: lim q→0fxc(q,q,ω)=− e2α(ω) q2,(4) where α(ω)=− 1 e2n2 0 G=0 (G·ˆ q)2fHL xc (G,ω)−fHL xc (G,0)|n(G)|2. (5) Here ˆ q=q/|q|,n(G) is the Fourier transform of the electron density evaluated at the reciprocal lattice vector G, and n0=n(G=0) is the average density. fHL xc represents the longitudinal component of the tensor xc kernel of a uniform (homogeneous) electron gas, which coincides with its scalar counterpart; this uniform-gas xc kernel is evaluated at the average electron density n0. By averaging over all ˆ qdirections, (G·ˆ q)2is replaced by G2/3. In the uniform limit, α(ω)→0. A negative sign is introduced in Eq. (4), as in Refs. [20,34], in such a way that a positive α(ω) could cancel the divergent Coulomb interaction 4π/q2as q→0. In Ref. [29], the xc kernel fHL xc (q,ω) entering Eq. (5) was approximated as fHL xc (q,ω)≈fHL xc (0,ω), taking the latter from Ref. [35]. Here we go a step further by using the nonlocal MCP07, which represents a reliable wave-vectorand frequency-dependent uniform-gas xc kernel. From this dependence on the wave-vector G, we expect that fxc(q,q,ω) can cover new spectral features that are simply absent when fHL xc (0,ω)isusedinEq.(5). At large wave-vectors q21/k, the frequency dependence of the MCP07 xc kernel is damped out significantly. In the limit q→∞, the MCP07 kernel approaches its static limit, therefore we expect a significant contribution to Eq. (5) for small values of Gonly. More details of the frequency dependence can be found in Ref. [22]. On the other hand, our calculations indicate that the wavevector dependence of the uniform-gas xc kernel—neglected in Ref. [29]—can largely affect ultranonlocality. Our calculation of the frequency-dependent coefficient α(ω)ofEq.(5) hinges upon the evaluation of: (i) the Fourier coefficients n(G) and (ii) the MCP07 xc kernel of Eq. (3), which is based on exact constraints. Consider a UEG of density n0=k3 F/(3π2) perturbed by a weak periodic external pseudopotential W(r)= G W(G)eiG·r(6) 125112-2 PROGRESS TOWARDS UNDERSTANDING … PHYSICAL REVIEW B 104, 125112 (2021) FIG. 1. α(ω) as a measure of ultranonlocality in fcc Al crystal. The upper (lower) figure shows the real (imaginary) part of the complex α(ω) for the GKI dynamic LDA (green, dot-dashed), MCP07 with k=0 (orange, dashed), the full MCP07 [22] (blue, solid), frequency-dependent Qian and Vignale [35] (QV, red solid), QV-MCP07-TD (purple, dashed), and 2p2h [21] (black, dot-dashed) xc kernels. We are using fHL xc (q,ω), as shown in the calculation of α. The labels of the curves indicate which fHL xc (q,ω) was used as input to Eq. (5). representing the actual crystal lattice. For W(G)= (n0/z)w(G) we apply the evanescent core pseudopotential w of an ion of valence zdescribed by Eq. (2.11) of Ref. [36]. The pseudopotential is designed to have a finite value at r=0, with vanishing first and third derivatives. This analytic behavior leads to a quick convergence of its Fourier transform in the limit G→∞. Linear response tells us how to find n(G): n(G)=χ(G)W(G)=χ0(G) (G)W(G),(7) with χ0(G)=−kF/π2F(y), where y=G/(2kF). Fis the static response function of the noninteracting UEG, known as the Lindhard function [37]: F(y)=1 2+1−y2 4yln  1+y 1−y .(8) The screening in the UEG is represented by the dielectric function: (G)=1−4π G2+fxc(G)χ0(G),|G|>0,(9) where fxc(G) is the static xc kernel of the UEG. The xc kernels fHL xc (G,0) and fHL xc (G,ω) entering Eq. (5) are taken from Eqs. (2) and (3), respectively. III. RESULTS AT THE OPTICAL LIMIT FROM LOCAL (LDA) AND NONLOCAL (MCP07) DYNAMIC KERNELS Figure 1displays the coefficient α(ω) of face-centered cubic (fcc) aluminum determined from either (i) Eq. (5) with FIG. 2. α(ω) as a measure of ultranonlocality in bcc Na crystal. The upper (lower) figure shows the real (imaginary) part of the complex α(ω) for the GKI dynamic LDA (green dot-dashed), MCP07 with k=0 (orange dashed), the full MCP07 (blue solid) [22], frequency-dependent Qian and Vignale [35] (QV, red solid), QV-MCP07-TD (purple, dashed), and 2p2h [21] (black, dot-dashed) xc kernels. the fully nonlocal MCP07 fxc of Eq. (3) as input, or (ii) the local density approximation (LDA) version of Eq. (5): αDLDA(ω)=− 1 3e2n2 0 G=0 G2fHL xc (0,ω)−fHL xc (0,0)|n(G)|2. (10) fHL xc (0,ω) is taken to be either the GKI dynamic LDA xc kernel [the long-wavelength limit of Eq. (3)] or the QianVignale (QV) dynamic LDA xc kernel of Ref. [35]. Note that the GKI and QV dynamic LDAs tend to distinct static limits; this is discussed further in Appendix B. Hereafter, we will use “the dynamic LDA” to refer to the GKI expressions. When using the fully nonlocal MCP07 xc kernel of Eq. (3), the parameter kis taken to be either equal to k, or equal to zero. We also introduce a hybrid kernel, which replaces the GKI frequency dependence in MCP07 with the QV model. As seen in Appendix B, ensuring that fxc(0,0) yields the ALDA yields more realistic correlation energies in the metallic range. This new kernel is called QV-MCP07-TD. All sums over G used a sufficiently large cutoff of G2 c/2<800 eV. We also present results for a correlated-wave-function-derived kernel, as discussed after Eq. (12). Our calculations clearly indicate that the coefficient α(ω) is particularly sensitive to the wave-vector dependence of the xc kernel, as seen by comparing MCP07 to its counterpart with the damping factor kof Eq. (3) set to zero. Physically, setting k=0 strengthens the frequency dependence of α(ω) significantly over MCP07 (where k=k). Indeed, |α(ω)|is reduced by more than 85% when kis increased from zero to its full MCP07 strength, k=k. These are consistent behaviors within the visible range of frequencies around 3 eV. The same physics is reported in Fig. 2for body-centered cubic (bcc) Na. The coefficient α(ω) is again particularly 125112-3 NEPAL, KAPLAN, PITARKE, AND RUZSINSZKY PHYSICAL REVIEW B 104, 125112 (2021) FIG. 3. α(ω) as a measure of ultranonlocality in cds Si crystal, calculated with the GKI dynamic LDA (green dot-dashed), MCP07 k=0 (orange dashed), MCP07 (blue solid), long-wavelength Qian and Vignale [35] (QV, red solid), QV-MCP07-TD (purple, dashed), and 2p2h [21] (black, dot-dashed) xc kernels. The horizontal gray dashed line is the static limit of the LRC ultranonlocality coefficient, α=0.13 [34]. The QV kernel comes closest, at low frequencies, to approximating this value. sensitive to whether the wave-vector dependence of the xc kernel is considered. Note that our values shown in Figs. 1 and 2for Re α(ω) at optical frequencies differ considerably from the metallic limit (α=−0.213) of Eq. (4) of Ref. [20], which was fitted to semiconductor data. Figures 3,4, and 5present our calculations of α(ω)for cubic diamond structure (cds) Si and C. In this case we obtain reliable ground-state valence electron densities using a planewave basis set. These calculations employed r2SCAN [38], a computationally efficient and highly accurate ground-state density functional, within the Vienna ab initio simulation package [39]. Refer to Appendix Afor definitions of the quantities n0and n(G) entering Eq. (5). To ensure that a large number of wave-vectors were used in the Hamiltonian, the calculation was performed on a -centered kmesh of spacing 0.08 Å−1, with an energy cutoff of 800 eV. The tetrahedron integration method was used to obtain reliable total energies converged within 10−7eV. The equilibrium volume was determined by fitting to the stabilized jellium equation of state [40]; r2SCAN predicts an equilibrium cubic lattice parameter of 5.440 Å for Si and 3.562 Å for C. The real-space density was then Fourier transformed to yield n(G) for reciprocal lattice vectors in the range |G|2/2<3200 eV. Our results for α(ω) in semiconductors are presented in Figs. 3and 4for cds Si, and Fig. 5for cds C. At optical frequencies in the range 3–5 eV, the coefficient α(ω) is positive, as expected, only when the frequency dependence of the xc kernel is taken to be that of Qian and Vignale. However, even the QV kernel considerably underestimates the expected values α=0.2[20]or0.28[34]. Note that the QV xc kernel uses the xc shear modulus μxc, which was tabulated only for rs=1,2,3,4,5inRef.[35]. To interpolate and extrapolate FIG. 4. α(ω) as a measure of ultranonlocality in cds Si crystal for the small omega regime only, calculated with the GKI dynamic LDA (green dot-dashed), MCP07 k=0 (orange dashed), MCP07 (blue solid), long-wavelength Qian and Vignale [35] (QV, red solid), QVMCP07-TD (purple, dashed), and 2p2h [21] (black, dot-dashed) xc kernels. The horizontal gray dashed line is the static limit of the LRC ultranonlocality coefficient, α=0.13 [34]. The inset shows the range 0<ω<10 eV. We emphasize that Ref. [29] used a different sign convention for α(ω), fxc =e2α/q2. Therefore, Fig. 2 of Ref. [29], which plots α(ω) for Si, appears to have the sign of α(ω)reversed. these values, we used the physically motivated [41]form μxc(rs) n=a rs +(b−a)rs r2 s+c,(11) with a=0.031152, b=0.011985, and c=2.267455 fitted to the values reported in Ref. [35] (in atomic units). The PerdewWang (PW92) [42] parametrization of the UEG correlation energy per electron εcwas used as input to the static compressibility of the QV and QV-MCP07-TD kernels. For the GKI dynamic LDA and MCP07 with k=0ork, we employed the Perdew-Zunger [23] parametrization of εc, consistent with Ref. [22]. For Si and C we also present (when possible) the static limit of the ultranonlocality coefficient determined by the empirical long-range contribution (LRC) xc kernel of Ref. [34]: fLRC xc (q,ω)=− α+βω2 q2.(12) Here αand βare material-dependent parameters that are fitted to spectroscopic data. For Si, α=0.13,β=0.00635 eV−2; for C, α=0.28,β=0.00135 eV−2. This gives a benchmark for the kernels presented here, and helps determine the validity of Eq. (5) for insulators. As an alternative to the LRC benchmark, we also consider the 2p2h kernel [21]. This kernel is determined directly from Fermi hypernetted chain calculations of the UEG, including two-particle, two-hole (2p2h) interactions. As these excitonic interactions are relevant for the optical regime, the 2p2h kernel may be the best point of reference for our work. 125112-4 PROGRESS TOWARDS UNDERSTANDING … PHYSICAL REVIEW B 104, 125112 (2021) FIG. 5. α(ω) as a measure of ultranonlocality in cds C crystal, calculated with the GKI dynamic LDA (green dot-dashed), MCP07 k=0 (orange dashed), MCP07 (blue solid), long-wavelength Qian and Vignale [35] (QV, red solid), QV-MCP07-TD (purple, dashed), and 2p2h [21] (black, dot-dashed) xc kernels. The static limit of the LRC ultranonlocality coefficient α=0.28 [34] is beyond the scale of the vertical axis. As cds C has a much larger band gap than cds Si, we begin to see the limited validity of applying Eq. (5) to insulators. See Appendix Cfor a discussion of the 2p2h kernel and its limitations. The real part of the coefficient α(ω), as defined in Eq. (4), is expected to vanish at ω=0 for metals and to be positive at ω=0 for semiconductors and insulators (e.g., +0.2 as in Ref. [20] for cds Si). Equation (5), which is formally exact in the limit of weak inhomogeneity and is, therefore, suitable for metals, always yields α(ω=0) =0, as expected for metals. For ωin the visible range, we expect Re α(ω) to be positive (and small for metals). The computed sign of Re α(ω) turns out, however, to be positive in the visible range only when the QV or 2p2h frequency dependence is used. The MCP07 xc kernel of Eq. (3) represents an extension, for finite wave-vectors, of the GKI dynamic xc kernel, which differs considerably at small ωfrom its QV counterpart: they satisfy distinct ω→0 limits, and thus have different Taylor expansions near ω=0. These important differences should be at the origin of the different behavior of α(ω), particularly at small ω, depending on whether the QV dynamic LDA xc kernel is used. While Re α(ω) is predicted to be negative (and not positive as expected) at small frequencies for the three GKI-based kernels under study, it does become (plausibly for metals) smaller in magnitude as we make the wave-vector dependence of this kernel increasingly sophisticated (from dynamic LDA to MCP07 k=0 to MCP07 k=k). Our work reveals the remarkable sensitivity of Eq. (5) to the wave-vector and frequency dependence of the UEG xc kernel. The long-range part of the xc kernel is expected to be highly nonmonotonic in its frequency dependence [43]. This was demonstrated in Ref. [43] for cds Si and cds C using the response function computed from both the Bethe-Salpeter equation and ground-state LDA eigenstates. While all kernels FIG. 6. The correlation energy per electron of the UEG as a function of the density parameter rs. Besides the PW92 ground-state LDA (black solid), the figure displays the ALDA (purple dotted), RPA (dark blue dashed), MCP07 (light blue solid), MCP07 k=0 (orange dashed), the GKI dynamic LDA (dark green dash-dotted), and the QV dynamic LDA (red solid) xc kernels. See Fig. 7and Table Iin Appendix Bfor a discussion of the relative errors made by these and other kernels. presented here show a nontrivial frequency dependence, none of them demonstrate the anticipated oscillatory behavior of α(ω). We do not advocate using MCP07 in its current form for general optical applications; however we recommend further testing of the QV-MCP07-TD kernel for the study of optical properties. This novel kernel demonstrates that introducing new constraints beyond the GKI-based frequency dependence could improve MCP07’s predictions of optical properties. Our aim here is to highlight the relevance of a proper description of both the frequency and wave-vector dependence of fxc in the optical limit. Future work could study the application of the QV-MCP07-TD kernel to describe the optical properties of a broader range of semiconductors. To complete the analysis of the coefficient α(ω), we have performed a calculation of the ground-state correlation energy of the UEG for electron-density parameters in the range from rs=1tors=10. The PW92 approximation [42] (black solid line in Fig. 6) is regarded as exact. In Fig. 6we plot the full MCP07 correlation energy (blue solid line in Fig. 6) together with the result of using various approximations to fxc: the fully nonlocal dynamic MCP07 with k=0, the GKI dynamic LDA, the QV dynamic LDA, the ALDA, and the RPA. The method for obtaining correlation energies through the fluctuation dissipation theorem [44], using the Cauchy integral formula for the frequency integral, is described in Ref. [22]. Here we see that while the ALDA is known to overestimate the correlation energy significantly for all values of rs, introducing frequency dependence (still within the LDA, q=0infxc) improves the correlation energy considerably. The fully nonlocal and dynamic MCP07 xc kernel yields excellent correlation energies, particularly with increasing rs. When the frequency dependence of the MCP07 xc kernel is undamped by taking k=0, the correlation energy is clearly worsened, but not to the same extent as in the case of α(ω). Note, however, the role of the wave-vector in determining accurate correlation energies. Both the QV kernel and 125112-5 NEPAL, KAPLAN, PITARKE, AND RUZSINSZKY PHYSICAL REVIEW B 104, 125112 (2021) dynamic LDA, which are wave-vector independent, tend to overcorrect the RPA. However, the MCP07 k=0 kernel, which adds a naive wave-vector dependence to the dynamic LDA, gives less reliable estimates of the UEG correlation energy. Using nonzero kand correctly interpolating between known limits, as in MCP07, improves upon the dynamic LDA. It is easy to show that the analytic continuation of fxc(ω)to imaginary frequency iu (with ωand upurely real), which is needed for the calculation of correlation energies, is Re fxc(iu)=f∞+1 2π∞ −∞ [ω2+u2]−1{[Re fxc(ω)−f∞]u +Im fxc(ω)ω}dω,(13) Im fxc(iu)=1 2π∞ −∞ [ω2+u2]−1{−[Re fxc(ω)−f∞]ω +Im fxc(ω)u}dω,(14) with f∞≡limω→∞ fxc(ω). The GKI and QV dynamic LDAs make Im fxc(ω) a purely odd function of real-valued frequency ω, and thus their corresponding Re fxc(ω)areeven functions of ω.FromEq.(14) we see that the integrand of Im fxc(iu)isoddinω, and thus fxc(iu) is purely real. IV. CONCLUSIONS We have used the dynamic MCP07 kernel to study the roles of the wave-vector and frequency dependence of fxc in the optical limit. Both ingredients have a significant impact in the visible region of light. The strong wave-vector dependence of α(ω) corroborates the presence of local-field effects that have been discussed earlier. Less is known about the role of the frequency dependence in the optical limit. To describe local-field effects, a kernel with the correct ultranonlocal limit is needed. Although there are a few efforts listed in the literature, no nonempirical xc kernels within TDDFT yield universal applicability for the optical absorption spectra of semiconductors and insulators. The TDCDFT framework allows the construction of, in the optical limit, a frequency-dependent xc kernel of a weakly inhomogeneous electron gas that exhibits ultranonlocality, as anticipated in Ref. [29]. In this work we rely on this formalism to evaluate the coefficient α(ω) as a measure of ultranonlocality in the optical limit. We have included three types of xc kernels in our analysis. As potentially better references, we have also included the dynamic LDA of Qian and Vignale [35] and the so-called 2p2h kernel [21]. Although there is no reference about what the correct ultranonlocality coefficient in the optical limit is, the completely local dynamic LDA sets a negative extreme with the GKI frequency model. Our work demonstrates the relevance of exact constraints. An approximate kernel can improve beyond the ALDA [which makes α(ω)=0] by enforcing known limiting behaviors of the exact fxc(q,ω), as the GKI, QV, and MCP07 kernels are constructed to do. The MCP07 kernel is an interpolation between the static and frequency-only limits. The impact of wave-vector and frequency have been independently investigated using the MCP07 model. Turning off the wave-vector and frequency dependence reduces all UEG-based xc kernels to the ALDA. Turning off the damping factor kin MCP07 results in a less-controlled frequency dependence, but the resultant kernel retains the correct wave-vector dependence of MCP07 in the static limit. The resulting α(ω) is significantly reduced compared to that of the ALDA, indicating a reduced degree of ultranonlocality. The full MCP07 kernel damps the frequency dependence of the GKI dynamic LDA, and reduces α(ω) further. We have tested all these kernels for fcc Al, bcc Na, cds Si, and cds C. These are metallic and semiconducting systems with small ultranonlocality for metals, and larger ultranonlocality for semiconductors. For the metallic systems, perturbation theory applies. We have obtained novel and informative results as an estimation of the ultranonlocality for metals. Our analysis confirms that even metals can have ultranonlocality, although to a much lesser extent than semiconductors. Neither the dynamic LDA nor MCP07 exhibit the estimated level and sign of ultranonlocality for Si and especially for diamond C. Our results are qualitative and indicate a strong sensitivity of α(ω) in the optical limit to the frequency and wave-vector dependence of the uniform gas kernel. This sensitivity was out of the reach in Ref. [29], as only the QV dynamic LDA kernel was used, in contrast with the various nonlocal kernels we have considered here. Furthermore, the nonlocal MCP07 and QV-MCP07-TD kernels considered here appear to be a reasonable basis for further improvements, especially when compared with the dynamic LDA. Our work clearly indicates the limitations of the GKI frequency model. Instead, the hybrid QV-MCP07-TD kernel introduced here unifies both wave-vector and frequency constraints, and shows a good promise within the limitations of the perturbation theory applied throughout this work. Our current work could (i) guide further modifications in the frequency-dependent MCP07 kernel for optical spectroscopy to build in more exact constraints on the frequency dependence and (ii) guide more efforts to extend the TDCDFT scheme to obtain α(ω) and the corresponding limq→0fxc(q,q;ω) for nonuniform systems. The Python code used to generate Figs. 1–5from VASP outputs is made freely available (without access restrictions) at [45] under the “code” directory. The raw data from VASP is also included there, under the “Si” and “C” directories. The processed data is located in the “code/data_files” subdirectory. ACKNOWLEDGMENTS The authors are grateful for fruitful discussions with Professor John P. Perdew. A.R. and N.K.N. acknowledge support from the U.S. National Science Foundation under Grant No. DMR-1553022. A.D.K. acknowledges support from the U.S. Department of Energy, Basic Energy Sciences, through the Energy Frontier Research Center for Complex Materials from First Principles Grant No. DE-SC0012575; and Temple University. N.K.N. and A.D.K. performed the calculations and revised the paper. J.M.P. and A.R. designed the research and revised the paper. A.R. wrote the paper. 125112-6 PROGRESS TOWARDS UNDERSTANDING … PHYSICAL REVIEW B 104, 125112 (2021) APPENDIX A: DEFINITIONS OF DENSITY VARIABLES Consider a density n(r) that varies weakly about an average uniform density n0. As is done in all standard electronic structure codes, we sample the (self-consistent) density in real space at Nrpoints R. To obtain the Fourier components n(G), we take the discrete Fourier transform n(G)=1 NR R n(R)e−iG·R,(A1) leveraging the fast Fourier transform for a suitable choice of R, R=n1 N1 a1+n2 N2 a2+n3 N3 a3,(A2) G=m1b1+m2b2+m3b3,(A3) G·R=2πn1m1 N1 +n2m2 N2 +n3m3 N3,(A4) where ni,mi=0,1,2,...,Ni−1, such that N1N2N3=NR. aiare the direct lattice vectors, and bjare the reciprocal lattice vectors such that ai·bj=2πδij. This convention is adopted by VASP; however, any standard plane-wave code will make similar choices for the discrete Fourier transform phase and normalization conventions, and will make similar choices for R. Real solids are typically not weakly varying; however, the valence densities of Si, and less so C, are approximately weakly varying about their average density, which we define as n0≡n(G=0).(A5) In pseudopotential codes like VASP, only valence electrons are considered, and the core electrons are replaced by a nonlocal effective potential inside a core radius. Therefore, the variables n(R) and n(G) represent the valence electron density and the Fourier transform of the valence electron density, respectively. Virtually all xc kernels based on the UEG paradigm require a real-density input; therefore, in Eq. (5) we evaluate fxc(q,ω)≡fxc(n0,q,ω).(A6) Let rs3=3/(4πn0). In our self-consistent calculations we found rs≈2.009 bohr for cds Si (which is within a reasonable metallic range, 2 rs5 bohrs), and rs≈1.315 bohrs for cds C (which is outside the normal metallic range). APPENDIX B: A MORE DETAILED DISCUSSION OF THE CORRELATION ENERGIES In this Appendix we compare a variety of xc kernels in predicting jellium correlation energies (per electron) for the physically relevant range of electron densities 1 ⩽rs⩽10. We take PW92 [42] to be essentially exact. The correlation energy per electron εccan be computed from the adiabatic-connection fluctuation-dissipation theorem [46]: εc=1 2d3q (2π)31 0 dλ λ∞ 0 dω4πλ q2[Sλ(q,ω)−S0(q,ω)], (B1) where fxc,λ(q,ω,rs)=λ−1fxc(λ−1q,λ −2ω, λrs)[47] and the spectral function Sλat coupling constant λis given by χλ(q,ω)=χ0(q,ω) 1−[4πλ/q2+fxc,λ(q,ω,rs)]χ0(q,ω),(B2) Sλ(q,ω)=− 1 πnIm χλ(q,ω).(B3) In the random phase approximation (RPA), fxc is taken to be zero. In this approximation, correlation energies are found to be too negative [47]. The use of appropriate approximations for fxc should build upon the RPA and improve its prediction of correlation energies. This section also presents a few new kernels. In addition to those presented previously, we also consider the static limit of the MCP07 kernel, and the dynamic kernel formed by replacing, within MCP07, the GKI dynamic LDA with the QV dynamic LDA. This combination can be done in two ways, as the TDDFT and TDCDFT static and q→0 limits are incompatible [41]: lim q→0lim ω→0fHL xc (q,ω)=fGKI xc (0) =fALDA xc ,TDDFT,(B4) lim ω→0lim q→0fHL xc (q,ω)=fQV xc (0) =fALDA xc +4 3 μxc(rs) n2,TDCDFT, (B5) where fHL xc (q,ω) is the longitudinal xc kernel of a homogeneous electron gas within TDCDFT and μxc is the xc shear modulus [see Eq. (11)]. For consistency, the two new QV-MCP07 kernels use the PW92 parametrization of the correlation energy, and have the form fxc(q,ω)=1+e−kq2fQV xc (ω) fxc(0,0) −1fMCP07 xc (q,0).(B6) To obtain the QV-MCP07-TD kernel, we take fxc(0,0) = fALDA xc in the MCP07 static kernel and the QV kernel, by setting μxc =0. To obtain the QV-MCP07-TDC kernel, we take fxc(0,0) =fALDA xc +4μxc(rs)/(3n2), which modifies k and fMCP07 xc (q,0), consistent with the construction principles of Ref. [22]. FIG. 7. Comparison of the relative errors εapprox. c−εPW92 cfor a variety of functionals. Error statistics are given in Table I. 125112-7 NEPAL, KAPLAN, PITARKE, AND RUZSINSZKY PHYSICAL REVIEW B 104, 125112 (2021) TABLE I. Comparison of approximate exchange-correlation kernels in predicting jellium correlation energies (per electron), for 91 values of rsin the range 1 ⩽rs⩽10. PW92 [42] is taken to be the reference energy. Mean error Mean absolute error Standard deviation Kernel ×10−2hartree ×10−2hartree ×10−3hartree RPA −1.4387 1.4387 1.8092 ALDA 1.2750 1.2750 1.3251 Dynamic LDA 0.1082 0.1289 1.1242 MCP07 static 0.0836 0.1277 1.1603 MCP07, k=0−0.4592 0.4592 1.2591 MCP07 −0.0496 0.1077 1.3382 QV 0.1455 0.1705 1.3580 QV-MCP07, TD −0.0047 0.1074 1.3104 QV-MCP07, TDC −0.1251 0.1491 1.7068 Figure 7plots the relative errors in the correlation energies per electron as a function of 1 ⩽rs⩽10. Table Ipresents the average errors and standard deviations in the UEG correlation energy for the same range of rs. It can be seen that the most accurate kernels in this range are the MCP07 and QV-MCP07-TD kernels, whose mean absolute errors are comparable, differing only by about 4×10−6hartree. APPENDIX C: TECHNICAL ASPECTS OF THE 2P2H KERNEL CALCULATION The 2p2h kernel, while capturing a broad range of manyelectron physics, is tabulated only for a limited range of rs,q, and ω. Furthermore, no analytic expression has been determined to interpolate it. Thus, we are forced to make two approximations to use the 2p2h kernel in a practical computation. To interpolate the kernel, we use a multivariate linear interpolation. For a one-dimensional function F(x), this amounts to F(x)≈F(xi)xi+1−x xi+1−xi +F(xi+1)x−xi xi+1−xi ,(C1) xi⩽x<xi+1,(C2) and the tabulated values of xiare sorted by increasing values. In three dimensions we perform a simple composition of one-dimensional linear interpolations. More sophisticated multidimensional interpolation schemes, like tricubic spline, assume the function to be interpolated is smooth in some sense: e.g., a cubic spline assumes continuity up to the second derivatives. We cannot make such an assumption about fxc(rs,q,ω). Second, the 2p2h kernel is tabulated only up to q=8kF and ω=3.98ωp(0), where ωp(0) =(3/r3 s)1/2is the semiclassical plasmon frequency. These values are too small for the 2p2h kernel to attain its zero separation and infinite frequency limits. While the limiting values fxc(q=0,ω →∞)[24] and fxc(q→∞,ω =0) [48] are known, we cannot, in general, extrapolate to q→∞with ω>0, or vice versa. Thus, we are forced to cut off the Gsum of Eq. (5)for|G|>8kF, and restrict α(ω)toω⩽3.98ωp(0). The numeric values of the cutoffs are given in Table II. Note that, for all solids considered here, the shortest reciprocal lattice vector is longer than the smallest wave-vector (0.1kF) for which the 2p2h kernel is tabulated. TABLE II. Summary of cutoff energies used in the computation of α(ω) for the 2p2h kernel only. For all other kernels, a cutoff of |G|2/2<800 eV was used. Note that qcut =8kFand ωcut = 3.98ωp(0). Solid (structure) rs(bohr) q2 cut/2 (eV) ωcut (eV) C (cds) 1.32 231.68 124.33 Al (fcc) 2.07 93.56 62.99 Si (cds) 2.01 99.35 65.89 Na (bcc) 3.93 25.96 24.08 [1] E. Runge and E. K. U. Gross, Phys.Rev.Lett.52, 997 (1984). [2] E. K. U. Gross and W. Kohn, Phys. Rev. Lett. 55, 2850 (1985); ibid.,57, 923(E) (1986). [3] C. Ullrich, Time-Dependent Density-Functional Theory: Concepts and Applications (Oxford University Press, Oxford, 2012). [4] C. Ullrich and Z. Yang, Braz. J. Phys. 44, 154 (2014). [5]N.T.Maitra,J. Chem. Phys. 144, 220901 (2016). [6] M. Rohlfing and S. G. Louie, Phys. Rev. B 62, 4927 (2000). [7] W. Hanke and L. J. Sham, Phys.Rev.Lett.43, 387 (1979). [8] M. Petersilka, U. J. Gossmann, and E. K. U. Gross, Phys. Rev. Lett. 76, 1212 (1996). [9] P. Nozières and D. Pines, Nuovo Cimento 9, 470 (1958). [10] C. E. Patrick and K. S. Thygesen, J. Chem. Phys. 143, 102802 (2015). [11] J. E. Bates, S. Laricchia, and A. Ruzsinszky, Phys. Rev. B 93, 045119 (2016). [12] L. A. Constantin and J. M. Pitarke, Phys. Rev. B 75, 245127 (2007). [13] A. V. Terentjev, L. A. Constantin, and J. M. Pitarke, Phys. Rev. B98, 085123 (2018). 125112-8 PROGRESS TOWARDS UNDERSTANDING … PHYSICAL REVIEW B 104, 125112 (2021) [14] E. K. U. Gross, J. F. Dobson, and M. Petersilka, in Density Functional Theory II: Relativistic and Time Dependent Extensions, edited by R. F. Nalewajski (Springer, Berlin, 1996), pp. 81–172. [15] P. E. Trevisanutto, A. Terentjevs, L. A. Constantin, V. Olevano, and F. D. Sala, Phys.Rev.B87, 205143 (2013). [16] S. Rigamonti, S. Botti, V. Veniard, C. Draxl, L. Reining, and F. Sottile, Phys.Rev.Lett.114, 146402 (2015). [17] J. A. Berger, Phys.Rev.Lett.115, 137402 (2015). [18] Z.-h. Yang, F. Sottile, and C. A. Ullrich, Phys. Rev. B 92, 035202 (2015). [19] F. Bruneval, F. Sottile, V. Olevano, and L. Reining, J. Chem. Phys. 124, 144113 (2006). [20] 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.B69, 155112 (2004). [21] M. Panholzer, M. Gatti, and L. Reining, Phys.Rev.Lett.120, 166402 (2018). [22] A. Ruzsinszky, N. K. Nepal, J. M. Pitarke, and J. P. Perdew, Phys.Rev.B101, 245135 (2020). [23] J. P. Perdew and A. Zunger, Phys. Rev. B 23, 5048 (1981). [24] N. Iwamoto and E. K. U. Gross, Phys.Rev.B35, 3003 (1987). [25] J. P. Perdew and T. Datta, Phys. Status Solidi (b) 102, 283 (1980). [26] J. P. Perdew, A. Ruzsinszky, J. Sun, N. K. Nepal, and A. D. Kaplan, Proc. Natl. Acad. Sci. USA 118, e2017850118 (2021). [27] A. W. Overhauser, Phys. Rev. 167, 691 (1968). [28] S. Sharma, J. K. Dewhurst, A. Sanna, and E. K. U. Gross, Phys.Rev.Lett. 107, 186401 (2011). [29] V. U. Nazarov, G. Vignale, and Y.-C. Chang, Phys. Rev. Lett. 102, 113001 (2009). [30] G. F. Giuliani and G. Vignale, Quantum Theory of the Electron Liquid (Cambridge University Press, Cambridge, 2005). [31] G. Vignale and W. Kohn, Phys.Rev.Lett.77, 2037 (1996). [32] V. U. Nazarov, J. M. Pitarke, Y. Takada, G. Vignale, and Y.-C. Chang, Phys.Rev.B76, 205103 (2007). [33] G. Vignale, C. A. Ullrich, and S. Conti, Phys.Rev.Lett.79, 4878 (1997). [34] S. Botti, A. Fourreau, F. Nguyen, Y.-O. Renault, F. Sottile, and L. Reining, Phys. Rev. B 72, 125203 (2005). [35] Z. Qian and G. Vignale, Phys. Rev. B 65, 235121 (2002);71, 169904(E) (2005). [36] C. Fiolhais, J. P. Perdew, S. Q. Armster, J. M. MacLaren, and M. Brajczewska, Phys.Rev.B51, 14001 (1995);53, 13193(E) (1996). [37] J. Lindhard, Mat. Fys. Medd. Dan. Vid. Selsk. 28, 1 (1954). [38] J. W. Furness, A. D. Kaplan, J. Ning, J. P. Perdew, and J. Sun, J. Phys. Chem. Lett. 11, 8208 (2020);11, 9248 (2020). [39] G. Kresse and J. Furthmüller, Phys. Rev. B 54, 11169 (1996); G. Kresse and J. Hafner, ibid. 49, 14251 (1994);47, 558 (1993); G. Kresse and J. Furthmüller, Comp. Mater. Sci. 6, 15 (1996). [40] V. N. Staroverov, G. E. Scuseria, J. Tao, and J. P. Perdew, Phys. Rev. B 69, 075102 (2004). [41] S. Conti and G. Vignale, Phys.Rev.B60, 7966 (1999). [42] J. P. Perdew and Y. Wang, Phys.Rev.B45, 13244 (1992). [43] R. Del Sole, G. Adragna, V. Olevano, and L. Reining, Phys. Rev. B 67, 045207 (2003). [44] F. Furche and T. Van Voorhis, J. Chem. Phys. 122, 164106 (2005). [45] https://github.com/esoteric-ephemera/ultranonlocal_ semiconductors. [46] D. C. Langreth and J. P. Perdew, Phys.Rev.B15, 2884 (1977). [47]M.Lein,E.K.U.Gross,andJ.P.Perdew,Phys. Rev. B 61, 13431 (2000). [48] M. Corradini, R. Del Sole, G. Onida, and M. Palummo, Phys. Rev. B 57, 14569 (1998). Correction: Minor errors in Eq. (B5) and in the text below Eq. (B6) have been fixed. 125112-9