{'fraction_non_alphanumeric': 0.032818532818532815, 'fraction_numerical': 0.011728709841917389, 'mean_word_length': 4.732929630403007, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Understanding plasmon dispersion in nearly-free-electron metals: the relevance of exact constraints for novel exchange-correlation kernels within time-dependent density functional theory

Small-wavevector excitations in Coulomb-interacting systems can be decomposed into the high-energy collective longitudinal plasmon and the low-energy single-electron excitations. At the critical wavevector and corresponding frequency where the plasmon branch merges with the single-electron excitation region, the collective energy of the plasmon dissipates into single electron-hole excitations. The jellium model provides a reasonable description of the electron-energy-loss spectrum (EELS) of metals close to the free-electron limit. The random phase approximation (RPA) is exact in the high-density limit but can capture the plasmonic dispersion reasonably even for densities with rs>1. RPA and all beyond-RPA methods investigated here, result in a wrong infinite plasmon lifetime for a wavevector smaller than the critical one where the plasmon dispersion curve runs into particle-hole excitations. Exchange-correlation kernel corrections to RPA modify the plasmon dispersion curve. There is however a large difference in the construction and form of the kernels investigated earlier. Our current work introduces recent model exchange-only and exchange-correlation kernels and discusses the relevance of some exact constraints in the construction of the kernel. We show that, because the plasmon dispersion samples a range of wavevectors smaller than the range sampled by the correlation energy, different kernels can make a strong difference for the correlation energy and a weak difference for the plasmon dispersion. This work completes our understanding about the plasmon dispersion in realistic metals, such as Cs, where a negative plasmon dispersion has been observed. We find only positive plasmon dispersion in jellium at the density for Cs.


I. INTRODUCTION
Due to its computational feasibility and relatively high accuracy, approximate Kohn-Sham density-functional theory (KS-DFT) simulations are the basis of present-day first-principles computational materials science. Particle-particle interactions can require treatment beyond semilocal KS-DFT [1]. Experimental applications heavily rely on understanding and guiding electron-electron interaction. A relevant example of electron-electron interaction is scattering resulting in electron energy loss. Experimentally the electron loss spectrum can be realized by electron energy loss spectroscopy or inelastic scattering [2,3].
Excited states can be accurately characterized by expensive Greens function techniques.
By virtue of the Runge-Gross theorem [4], DFT can be extended to time-dependent processes. Time-dependent density-functional theory (TDDFT) is becoming an attractive alternative to many-body perturbation theory and can offer, in principle, an unbiased and independent framework complementary to experimental observations, enabling the interpretation of specific experimental observations and predictions of new materials with targeted properties [5].
Plasmon excitations are collective oscillations of electrons in the absence of an external electric field, that incorporate Coulomb interaction between electrons [6,7]. Due to the electron-electron interactions, plasmon excitations establish a high barrier when testing ab initio theories. When the external perturbation is weak, linear response TDDFT [8] is a useful tool to describe optical excitation energies. In TDDFT the electron energy loss is quantified by the imaginary part of χ , the spatially nonlocal and dynamic density-density response function. The poles of the interacting density-density response function contain information about the optical excitation energies. The same density-density response function can deliver further information about the plasmons for a range of wavevectors.
Plasmon dispersion in nearly-free-electron alkali metals has attracted a great interest among experimentalists and theorists. The negative dispersion in the volume plasmon of the low-density alkali metals such as Rb and Cs has triggered a debate about the origin of the anomaly observed by most theoretical approximations within TDDFT and Fermi liquid theories [9]. A strong failure of approximations based on TDDFT with a static exchange-correlation kernel is the lack of a damping mechanism which results in an infinite lifetime of plasmons for a region of wavevectors smaller than the critical wavevector that separates plasmonic and particle-hole excitations [9,10]. The negative plasmon dispersion in low-density alkali metals can be attributed to correlation effects or band structure [11][12][13][14].
There are pros and contras for both explanations in the literature [13,14]. In heavier alkali metals such as Cs, electron transition to the near-Fermi-level d bands can occur [11,12].
The transition energy of these electrons is comparable to the plasmon energy, potentially causing a negative dispersion. Additional corrections to the interacting response function or dielectric function can originate in a weak lattice potential and core polarization effects [9].
The random phase approximation (RPA) is a Greens function-based method that is often used to obtain the ground-state correlation energy of bulk and two-dimensional materials [15][16][17][18][19]. Although RPA relies on the linear response TDDFT framework, its excitation energies are inaccurate because of the overestimated short-range exchange-correlation effects [20]. The exact exchange-correlation kernel f xc that would provide these effects is a functional derivative of the exchange-correlation potential. RPA is often interpreted within DFT to have roots in the adiabatic connection fluctuation dissipation theorem [21].
The bare RPA f xc = 0 without a band structure does not yield a negative dispersion in heavy alkali metals [11]. Theoretical predictions from the late 80's in polycrystalline metals showed negative dispersion for Cs beyond RPA when the band structure was included [9]. Some investigations of the role of correlation in the plasmon dispersion discuss kernels, but most tests consider the exchange-correlation effects at the level of adiabatic local density functional approximation (ALDA) only. The work of Tatarczyk et al. [22] steps beyond this limitation to some extent by considering some more model kernels based either on the uniform electron gas paradigm or on other constraints.
In our current work we aim to fill the gap in analyzing the plasmon dispersion with recent exchange-correlation kernels, beyond the early ones developed in the 90's. Our work aims to go beyond a simple analysis of exchange-correlation effects with kernels. The major goal here is to give an a priori numerical analysis why kernels by themselves (without the band structure effects) can not predict negative dispersion in low density alkali metals. In this work, we rely on the jellium model. The dimensionless density parameters (r s ) for the jellium model corresponding to different metals are taken from Ref. [23]. We demonstrate that the exact constraints can lead to kernels that correctly predict a positive plasmon dispersion in jellium.

II. METHODOLOGY: EXCHANGE-CORRELATION KERNELS WITHIN LIN-EAR RESPONSE TDDFT
Nonempirical construction of density functionals has allowed widespread and successful applications of these approximations for the ground state [24]. Various exact constraints such as the uniform electron gas limit, Lieb-Oxford bound or the one-electron limit are known for the ground state [25]. According to linear response theory the interacting and noninteracting density-density response functions are coupled by the Dyson equation: χ λ (q, ω) and χ 0 (q, ω) are the interacting and noninteracting response functions, respectively, v c (q) = 4πλ q 2 and f λ xc (q, ω) are the Coulomb and exchange-correlation kernels. λ is the coupling constant that provides the adiabatic connection between a noninteracting Kohn-Sham (λ = 0) and the interacting real system (λ = 1) response. When Eq. 1 is applied to the uniform electron gas, χ 0 (q, ω) becomes the Lindhard function with complex frequencies [26], a basic input to our current research. In the adiabatic approximation, the exact kernel is a second functional derivative of the ground state exchange-correlation energy. In practice the exact kernel is unknown but can be modeled by satisfying exact physical constraints.
Many real systems have densities close to the paradigm uniform electron gas, as in alkali metals. The uniform electron gas is therefore a simple model system with physical relevance. All exchange-correlation kernels in this work model the uniform electron gas with known limiting behavior at q → 0 and q → ∞. The simplest approximation is known as the adiabatic local density approximation (ALDA) kernel for λ = 1 [27]: dn 2 , where k F = (3π 2 n) 1/3 is the Fermi wavevector and ε c the correlation energy per particle of the uniform electron gas. A = 1 4 belongs to the exchange-only ALDA, while the density-dependent term in A gives the correlation beyond the high-density limit.
The real-space representation of ALDA is a delta function which indicates the spatial locality of this kernel. ALDA gives reasonable accuracy for low-frequency, long-wavelength excitations, but is not the right choice for a general correction to RPA [27]. The ALDA kernel was applied to the ground state correlation energy of the uniform electron gas but makes an error of ∼ 0.5 eV [28]. This error is the same in magnitude but of opposite sign to the error that RPA makes for the same system.
The ALDA kernel can be made nonlocal, by applying a cutoff that makes the exchangecorrelation kernel cancel the Hartree kernel for q > k cut . The cutoff is introduced by the renormalized ALDA (rALDA) expression [29] f rALDA with the cutoff wavevector k cut = k F √ A . By construction the rALDAxc kernel keeps the correct q → 0 limit of ALDA, but improves the wrong q → ∞ behavior of ALDAxc. For inhomogeneous systems, more ingredients like the density gradient or the kinetic energy density give more flexibility for kernels, as for ground state density functional approximations. The kinetic energy density is one of the ingredients of the nonlocal energy optimized (NEO) exchange-only kernel [30].
The NEO kernel improves the ground state correlation energy and structural properties of real systems beyond RPA [30]. The NEO kernel is designed to satisfy further physical constraints beyond both ALDA or rALDA, and has the form where β = 1 4 c(1−z 2 ) . The one-electron limit is reached when the ingredient z = τ w τ equals 1, where τ w is the one-electron kinetic energy density [31] and τ the positive kinetic energy density constructed from the Kohn-Sham orbitals. In the uniform electron gas, z is zero.
When q → 0, the NEO kernel is properly independent of q for the uniform electron gas.
The parameter c has a key relevance to the current work. In the construction of NEO, c is designed to give a correction to the RPA correlation energy in the high-density limit. In other words, the standard c = 0.264 parameter in NEO provides a unique fit to the exact second-order correlation energy for the spin-unpolarized electron gas. The "second-order exchange" contribution to the second-order correlation energy of the uniform gas is the correction to direct RPA from wavefunction anti-symmetry that survives in the high-density limit. It can be evaluated from explicit expressions given by von Barth and Hedin for RPA [32] and by Langreth and Perdew [20] beyond RPA.
β can be chosen to satisfy another constraint relevant to the long-wavelength ( q → 0) limit. This choice is also made in the ALDA, rALDA, and in the CP07 dynamic exchangecorrelation kernel constructed by Constantin and Pitarke [33]. The compressibility sum rule, [34], as is an important requirement for frequency-dependent exchange-correlation kernels. Satisfying the compressibility sum rule, β becomes or c = 0.5 in the high-density limit.
Thus the energy optimized value of c = 0.264 in the high-density limit is different from the value c = 0.5 that yields the correct small-q kernel in the high-density limit. In the next section we will extensively discuss the impact of these physical constraints on the plasmon dispersion and provide a novel insight about the role of correlation effects. It will turn out that the difference between the ALDA and NEO kernels is important for the correlation energy, but not very important for the plasmon dispersion.
To set the scale of the problem, we plot in Fig. 1 three kernels as functions of wavevector q: the ALDAxc kernel (exact at small q), the NEO xc kernel with c = 0.264 (which we will argue later is more correct than ALDAxc at larger q), and minus the Hartree kernel (-4π/q 2 ).
Clearly, in the wavevector region 0 < q/k F < 1 that shapes the plasmon dispersion, the Hartree kernel (the only one present in RPA) dominates over the xc kernel. This dominant effect of the Hartree kernel explains the overall good performance of RPA for plasmon dispersion. For wavevectors q/k F > 1, which are important for the correlation energy, the xc kernels have a larger effect. An even better xc kernel might interpolate between ALDAxc at small q and NEO at larger q.

III. PLASMON DISPERSION WITH SPATIALLY NONLOCAL EXCHANGE-CORRELATION KERNELS IN NEARLY-FREE-ELECTRON METALS
The alkali metals Na and K are nearly-free-electron (NFE) systems and therefore realizations of the uniform electron gas. Correlation effects in alkali metals can be strong enough, especially for low electron densities, to impact electronic excitations. The high-resolution electron energy loss experiments indicate that the plasmon dispersion in heavy alkali metals such as Rb and Cs becomes negative, i.e., the plasmon frequency decreases with increasing wavevector q [9]. Since all the above experiments and calculations were performed for periodic crystals, it is difficult to decouple correlation and band structure effects [35] in the dispersion of the plasmon excitations.
Further calculations by Ku and Eguiluz in K crystal [12] indicate the relevance of band structure versus correlation and demonstrate that the exchange-correlation effects beyond RPA at the ALDA level have only a minor role in the dispersion of plasmons. A similar observation was made by Quong and Eguiluz [36] for Na crystal. Although crystal periodicity is important when modeling realistic conditions, the jellium model can offer an important way to separate the impact of many-body correlations and band structure. Within our work, we want to provide an additional theoretical and numerical insight only about the many-body correlation effects and explain why in general no beyond-RPA approximation for jellium can predict the correct plasmon dispersion and lifetime for Cs. The justification of our results is based only on exact physical constraints imposed on the construction of beyond-RPA approximations. By using the jellium model for alkali metals, we can build upon the conclusion that ALDA gives a minor improvement beyond RPA for the plasmon dispersion. With the nonlocal exchange-correlation kernels developed since the later 90's, we can make further conclusions about how these recent approximations compare to ALDA and RPA in terms of correlation effects.
We can make two groups of assessed approximations. The first group includes exchangeonly and exchange-correlation kernels based on the ALDA approximations. ALDAx and ALDAxc are both local kernels but differ in correlation contribution. rALDAx and rALDAxc are nonlocal. The second group consists of NEO exchange-only kernels [30] with the c parameter constructed by satisfying different physical constraints. The default version of the NEO kernel yields the exact correction to the RPA correlation energy of jellium in the high-density limit: where G x (K) refers to the kernel, according to the correspondence between kernels and local-field factors. K = q 2k F is a dimensionless wavevector, and W = Ω 2k 2 F is a dimensionless frequency. The explicit expression for the second-order exchange energy e 2X c comes from Langreth and Perdew [21] beyond RPA and uses the RPA correlation energy for the uniform electron gas given by von Barth and Hedin [32]. The c parameter that corresponds to the second-order correlation energy was found to be 0.264.
Alternatively, in the long-wavelength limit we can use the compressibility sum rule formulated as d 2 dn 2 [nε xc (n)] to determine " c " . Then the c parameter of NEO can be estimated from the compressibility sum rule of the ALDAxc expression. This fitting delivers a different c = (0.43 -0.47) at metallic densities considered here. While formally the NEO approximation remains an exchange-only kernel, this kind of fitting brings long-wavelength exchange-correlation effects into our NEO kernel [33].
Clearly " c " controls the correlation or screening within the NEO approximation starting from c→ ∞ in RPA. The impact of " c " as a screening parameter on ground state correlation energies was established by Bates et al in 2016 [30]. Changing " c " from its default 0.264 value was shown to yield different correlation energies in the uniform electron gas at r s = 4.
The exchange-only NEO kernel can be explicitly turned into an exchange-correlation approach by replacing " c " by an electron-density-dependent parameter. This approach was tested for jellium slab correlation energies at moderate densities, and resulted in improved integrated correlation energy [37]. For a given density, the density dependence can make " c " significantly smaller than its default value. In our analysis we also investigate the effect of a low " c " parameter on the plasmon dispersion of various NFE metals. At first, we discuss the plasmon dispersion up to the wavevector region where plasmons decay into single-particle excitations. We consider all the exchange-correlation kernels described above. Within the static approximation for the kernel f xc (q) , the plasmon frequency ω p (q) is found by solving the equation (q, ω) = 1 − (v c (q) + f xc (q)) χ 0 (q, ω) = 0, for λ = 1 [22], and the solutions are undamped outside the particle-hole contiunuum, i.e., for q < q c where ωp(qc) For the r s values considered here, q c /k F ≈ 1 where k cut is the cutoff wavevector for a kernel. Thus rALDA and ALDA kernels will yield the same plasmon dispersion.
Al is a metal with rather high density, and RPA becomes relatively exact in the highdensity limit [36]. Here, we model Al by jellium with r s = 2.07. The small wavevector behavior is demonstrated by the correct plasmon energy known from an EELS experiment [38]. All ALDA and rALDA kernels return the correct long-wavelength limit of f xc . The c = 0.47 NEO exchange kernel keeps the correct long-wavelength of f xc . Since the compressibility sum rule delivers direct information about the long-wavelength limit, among all the approximations NEO c = 0.47 has a direct impact on the curvature of the dispersion.
The fitting against this exact constraint designates that the plasmon dispersion must start out as horizontal at small wavevectors. NEO c = 0.47 exemplifies the best behavior in the long wavelength limit (q → 0) that any static kernel for jellium at r s = 2.07 can demonstrate.
Comparing the ALDA and rALDA kernels in the left panel of For jellium at r s =3.93 (as for Na), all ALDA and rALDA dispersion curves barely differ in the left panel of Fig. 3. As in Al, there is no significant change coming from the nonlocal kernels. According to the right subfigure, the NEO c = 0.264 and NEO c = 0.44 kernels differ more than they do in Al for a range of ∼ 0.5 for dimensionless wavevector Q. NEO c = 0.0037 leads to unphysically low plasmon frequencies.
Cs is the alkali metal with the lowest density [39]. We modeled it here as jellium with r s To summarize the role of the exact constraints, we compare all the exchange and exchange-correlation models described above. Figure  We have also extracted the plasmon dispersion for the real metals Na and Cs by calculating the electron energy loss spectrum using the GPAW [40][41][42][43] code. The results are presented in Figure 6. We used the projector-augmented wave (PAW) [44] pseudopotential provided with the GPAW code, an energy cutoff of 600 eV, and a 16×16×16 k-point mesh to sample the Brillouin zone. Our calculations including band structure effects confirm that the kernels within the jellium model. For both Na and Cs, the dispersion with RPA is also displayed with band structure effects obtained from the GPAW code, showing results close to experiment [9].
The momentum transfer (q) is in the [100] direction for bulk solids.

IV. THE DYNAMIC STRUCTURE FACTOR WITHIN AND BEYOND-RPA
The dynamic structure factor [21] or spectral function S(q, ω) shows the distribution of frequencies ω for density fluctuations of wavevector q in the ground state of the uniform electron gas. Although it arises from all density fluctuations of a given wavevector, the spectral function for small wavevector typically peaks around the frequency of a plasmon. The inverse frequency width of this peak, by the uncertainty principle, provides a lower bound on the decay time for such a fluctuation, while the dependence of the peak on q reflects the plasmon dispersion. Here we will investigate the effects of various model exchange-correlation kernels on the spectral function. The dynamic structure factor S (q, ω) is proportional to Imχ [45], the loss component of the dynamic density-density response function: This quantity has been investigated by Lewis and Berkelbach [46] using an equation-ofmotion coupled cluster singles and doubles formalism, which unlike our TDDFT methods allows for plasmon decay via multi-pair electron-hole decay channels. With our real static exchange-correlation kernels, the plasmon at small finite q does not decay.
In Figure 7, we analyze our approximations at three wavevectors: q = 0.1k F , q = 0.5k F , and q = k F . The latter wavevector is close to the one at which the plasmon decays into the single-pair electron-hole continuum at r s = 4. The first Figure 7 (a) compares RPA to the default NEO kernel with c = 0.264. The line shapes at q = 0.1k F and q = 0.5k F are very similar for both methods. At q = k F the line shape from NEO becomes broader than the one from RPA, and the plasmon peak from NEO is shifted to a lower frequency. Figure 7.b compares three NEO kernels at r s = 4 for the same three wavevectors.
The static structure factor S(q) is the integral over frequency of S(q, ω), divided by the electron number N. S(q) determines the well-known correlation energy of the uniform electron gas. We will now show that the range of q/k F (less than or about equal to 4) that distributes to the correlation energy is much greater than the range of q/k F (less than or about equal to 1) that contributes to the plasmon dispersion. In the smaller range (but not in the larger one), the ALDAxc kernel is nearly sufficient, while the default NEO kernel is nearly sufficient over the larger range. This completes our plasmon dispersion analysis with an explanation why a negative dispersion cannot exist in jellium at the density of Cs.
Exchange-correlation kernels can be applied to improve the ground state correlation energy of RPA through the adiabatic connection fluctuation dissipation theorem. This is the basis of the wavevector decomposition of the ground state exchange-correlation energy as known from Langreth and Perdew [21]: where λ is the coupling constant along the adiabatic connection path and S λ (q) is the static structure factor found by integrating S λ (q, ω) over frequency and dividing by the electron number. According to the expression given by Eq. 9, the exchange-correlation energy depends upon the dynamic structure factor. The exchange energy E x replaces S λ by S 0 , and the correlation energy is E c = E xc -E x . Figure 8 shows the wavevector decomposition of the correlation energy for all exchange and exchange-correlation kernels considered here. We plot this for r s = 4 and r s = 5.62.
The physical basis of our analysis is the exact exchange-correlation kernel f xc (q, ω) of the uniform electron gas. For the correlation energy, the static version of the kernel f xc (q, 0) can be applied to a good approximation [28]. Note that the frequency dependence at least qualitatively can also be ignored for ω ≈ ω p [47,48]. The ALDA exchange-correlation kernel f ALDA xc (q, 0) approaches the exact kernel for the uniform electron gas at q → 0. In the long wavelength limit, As we will see below, f ALDA xc breaks down for q/k F ≥ 1, where our constraint-based NEO kernels become less negative and more accurate. Therefore, as suggested by Fig. 1: The ALDA approximation becomes a lower bound to the static exact exchange-correlation kernel for the correlation energy of the uniform electron gas. Figure 8 visualizes the relation between ALDAxc, RPA and some other exchange-correlation kernels. The ALDAxc is shown in the left panel of Figure 8. The ALDAxc is very accurate for small wavevectors but starts to deviate from NEO at q = k F . At q = 2k F the ALDAxc yields a strong overestimation of the correction to RPA correlation energy. All static beyond-RPA kernels make the exchange-correlation energy of RPA less negative for any density including r s = 5.62. The correlation energies from NEO c = 0.264 and from the NEO kernel fitted against the compressibility sum-rule are close to each other but the unphysical NEO c = 0.0037 adds a much larger correction to the RPA correlation energy.
The constraints of Equations 10 and 11 control the plasmon dispersion and the correlation energy. The dynamic structure factor becomes a link that couples the physics in the correlation energy and plasmon dispersion. From the correlation energy, the exact uniform-electron gas-based kernel must be less negative than the ALDAxc kernel. The NEO c = 0.264 kernel is uniform electron gas-based only through its energy optimization to the high-density limit, and performs reasonably for both plasmon dispersion and correlation energy in jellium.

V. CONCLUSION
We have presented various model exchange-correlation kernels beyond-RPA for the plasmon dispersion within the jellium model for alkali metals. We have shown that the plasmon dispersion is strictly controlled by exact constraints. Additional physics beyond the ALDA kernel, such as nonlocality in space, can be unimportant for plasmon dispersion. In change, physical constraints such as the compressibility sum rule determine the plasmon dispersion with exchange-correlation kernels. Clearly none of our methods based on particle-hole RPA for jellium is able to predict the experimentally observed negative plasmon dispersion for the heavy alkali metal Cs which arises from band structure. The current exchange-correlation kernels do not have an explicit density dependence that could have a larger impact. For the exact exchange-correlation kernel, the ALDAxc is likely a lower bound (as suggested by Fig. 1). The ALDAxc is accurate for q k F < 1, the range of q that determines the plasmon dispersion, even though the ALDAxc kernel fails badly for q k F 1, a range that contributes significantly to the correlation energy.