Density Functional Theory Perspective on the Nonlinear Response of Correlated Electrons across Temperature Regimes.

We explore a new formalism to study the nonlinear electronic density response based on Kohn-Sham density functional theory (KS-DFT) at partially and strongly quantum degenerate regimes. It is demonstrated that the KS-DFT calculations are able to accurately reproduce the available path integral Monte Carlo simulation results at temperatures relevant for warm dense matter research. The existing analytical results for the quadratic and cubic response functions are rigorously tested. It is demonstrated that the analytical results for the quadratic response function closely agree with the KS-DFT data. Furthermore, the performed analysis reveals that currently available analytical formulas for the cubic response function are not able to describe simulation results, neither qualitatively nor quantitatively, at small wavenumbers q < 2qF, with qF being the Fermi wavenumber. The results show that KS-DFT can be used to describe warm dense matter that is strongly perturbed by an external field with remarkable accuracy. Furthermore, it is demonstrated that KS-DFT constitutes a valuable tool to guide the development of the nonlinear response theory of correlated quantum electrons from ambient to extreme conditions. This opens up new avenues to study nonlinear effects in a gamut of different contexts at conditions that cannot be accessed with previously used path integral Monte Carlo methods.

Introduction

Quantum linear response theory (LRT) has been actively developed since the formulation of the foundations of quantum mechanics and has become one of the most fundamental theories for the computation of various properties.[1] At the same time, the ongoing development of technological and, along with it, experimental capabilities has resulted in the need for a theory that captures phenomena beyond the linear response regime. Specific examples include plasmonics,[2,3] optics,[4,5] and more recently warm dense matter (WDM)[6,7]—an extreme state that occurs in astrophysical objects[8,9] and that is also relevant for technological applications.[10−13]

However, in contrast to the LRT, the foundations of a quantum theory of the nonlinear response (NLRT) at finite wavenumbers is far from being established even for simple model systems such as a free electron gas.[14,15] In this regard, the lack of a reliable theoretical foundation makes the ab initio simulation an indispensable tool to guide the development of the NLRT. This was demonstrated recently for WDM by performing path integral quantum Monte Carlo (PIMC) simulations.[16,17] However, while these results are exact, the fermion sign problem[18,19] limits their application to moderate levels of quantum degeneracy. In contrast, the thermal Kohn–Sham density functional theory (KS-DFT) method[20] does not suffer from this limitation. Indeed, it has become standard practice to study the linear electronic response[21,22] based on the KS orbitals. In this work, we explore a new KS-DFT based approach to the nonlinear electronic response of arbitrary materials. First, this methodology allows us to compute higher-order (meaning quadratic, cubic, etc. with respect to the perturbation amplitude) response functions, with the only approximation being given by the choice of the exchange–correlation (XC) functional. In addition, we can straightforwardly estimate the validity range of LRT, which is highly important in its own right.

As a particular example, we apply this approach to the free electron gas[14,15]—the archetypical model system with general relevance for numerous applications in condensate matter physics and high-energy-density science. From a many-particle physics perspective, we note that it is imperative to first develop a NLRT for this general free electron gas model, before applying the NLRT to specific cases.

In this context, thermal KS-DFT[20] constitutes the method of choice because it allows calculations over large temperature ranges covering the strongly to partially degenerate regimes. Moreover, we note that the general nature of our present NLRT approach makes it directly useable for high-T DFT methods,[23,24] including orbital-free formulations.[25] Since the free electron gas model and the NLRT have important applications in WDM,[6,26] we start from relatively high temperatures relevant for laboratory astrophysics[27−29] as well as astrophysical models,[9] inertial confinement fusion,[10] and the synthesis of new materials at extreme conditions.[11−13] At these parameters, we can benchmark KS-DFT results against available PIMC results.[6,16,17] In addition, we consider lower temperatures down to the electronic ground state that are relevant for condensed matter physics.

It is convention to give the temperature T and density n0 of the free electron gas using the reduced temperature θ = T/TF (with TF being the Fermi temperature) and the mean inter-particle distance in a.u., rs = (4πn0/3)1/3. For example, a rather loose definition of the WDM regime corresponds to temperatures 0.1 ≲ θ ≲ 10 and densities 0.5 ≲ rs ≲ 10.

The prospect of the observation of nonlinear phenomena in WDM has triggered an active investigation of the nonlinear density response properties of the free electron gas by Dornheim et al.[16,17,30] using the ab initio PIMC method.[18,31] The focus of these PIMC studies was the static density response of WDM at temperatures θ ≥ 1. The main reason for not considering lower temperatures was the aforementioned fermion sign problem,[18] which results in an exponential increase in the computation time with decreasing temperature. Although there are other quantum Monte Carlo (QMC) methods that have different domains of applicability, such as the configuration PIMC approach,[32,33] the permutation blocking PIMC method,[34,35] and a phaseless auxiliary-field QMC technique,[36] there are always parameters at which QMC methods encounter significant difficulties. Approximately, the problematic domain for QMC methods corresponds to θ < 1 and rs ≳ 2.[33,36,37]

On the other hand, the parameter range corresponding to densities rs ≳ 2 and temperatures 0.01 ≲θ < 1 is highly important for numerous applications. Recently, it was shown that the static nonlinear density response functions of the electron gas can be used for the construction of advanced kinetic energy functionals required for orbital-free density functional theory (OF-DFT)-based simulations[38,39] with applications at ambient[40] and extreme conditions.[41,42] Additionally, nonlinear density response functions can extend quantum fluid models (quantum hydrodynamics and time-dependent OF-DFT) beyond the weak perturbation regime.[43−48] Moreover, static nonlinear density response functions are needed for the systematic improvement of effective pair interaction models for WDM[49−52] and liquid metals.[53−55] Finally, Moldabekov and co-workers[56,57] have recently suggested to deliberately probe the nonlinear regime in X-ray Thomson scattering experiments[58] as an improved method for the inference of plasma parameters such as the electronic temperature. However, these applications remain in their infancy since the NLRT of correlated electrons is significantly less developed compared to the LRT.[37] One of the reasons is that the derivations in the NLRT are much more mathematically involved.[59−61] In fact, the NLRT is burdened with easy-to-make-mistake mathematical tasks and poorly converging integrals. Therefore, the ab initio calculation of the NLRT properties across parameter ranges is required not only to describe certain phenomena but also to guide and test new theoretical developments.

The key goal of this work is to demonstrate the high value of KS-DFT to study the nonlinear density response across temperature regimes as an alternative to much more expensive—for certain parameters even prohibitively expensive—QMC simulations. This is achieved by developing and testing the KS-DFT based methodology for the analysis and investigation of the higher order static density response functions. Therefore, first of all, we show that KS-DFT can be effectively used to compute static nonlinear density response properties of correlated electrons at low temperatures (θ ≲ 1) and is able to reproduce available PIMC results at θ = 1. Second, we provide an analysis of the available theoretical results for the diagonal parts of quadratic and cubic response functions by combining the KS-DFT simulation of correlated electrons, the KS-DFT calculations with the XC functional set to zero, and recently developed machine learning (ML) representation of the local field correction (LFC) of the free electron gas.[62,63] This confirms the high accuracy of the analytical results for the quadratic response function and reveals the significant deficiency of the available analytical results for the cubic response function. Finally, we are able to show the change in the characteristic features of the NLRT functions on the way from moderately to strongly degenerate regimes.

The paper is organized as follows: in Section , we provide the theoretical background of the studied NLRT characteristics; in Section , we give the description of the performed simulations; the new results are presented and discussed in Section ; and the paper is concluded by summarizing the main findings and providing an outlook over future investigations in Section .

Theory

First, we briefly discuss the state-of-the-art theory of the static nonlinear density response functions. Along with that, we establish the terminology used throughout the paper. In general, the definition of the NLRT functions follow from the perturbative expansion of the density n(r) around its unperturbed value n0(r).[16,59,60] Specifically, we consider the response of the uniform electron gas to an external harmonic perturbation,[6,64]V(r) = 2Acos(q·r), with amplitude A and wavenumber q. In this case, the Fourier expansion of the density distribution reads[16]where we have introduced the density perturbation components in Fourier space ⟨ρ̂⟩. The latter quantity is essentially the density perturbation in k space induced by an external perturbation with amplitude A and wavenumber q.

From eq , we see that ⟨ρ̂⟩ has non-zero components at multiples of the perturbing field wavenumber, that is, at k = ηq with η being an integer number. We refer to ⟨ρ̂η⟩ at η = 1, η = 2, and η = 3 as density perturbations at the first, second, and third harmonics, respectively. Next, using the density response ⟨ρ̂⟩, we arrive at the following definitions of the density response functions[16,64,65]where χ(1) (q) is the linear response function, χ(1,cubic) (q) is the cubic response function at the first harmonic, χ(2) (q) is the quadratic response function, and χ(3) (q) is the cubic response function at the third harmonic. Evidently, eqs –4 are given by expansions in terms of the perturbation amplitude A and are accurate up to the third order. While the density response is only given by a single term at the wavenumber of the original perturbation within LRT, the consideration of nonlinear effects leads to a richer picture including the excitation of higher-order harmonics.

In the ideal Fermi gas approximation, the linear response function χ0(1) is given by the (temperature-dependent) Lindhard function.[66] On the same level of description, Mikhailov expressed the ideal quadratic response function χ0(2) (q) and ideal cubic response function at the third harmonic χ0(3) (q) in terms of the Lindhard function[65,67]

Next, on the mean-field level, usually called RPA, the results for χ(2) (q) and χ(3) (q) can be obtained by taking into account screening on the level of the linear response and dropping quadratic or higher order corrections to screening[16]and

Finally, some electronic correlation effects beyond the mean-field level can be taken into account using a LFC G(k) in the denominator[16]andSimilarly, the screened eqs –10 take into account electronic XC on the basis of the LRT. However, contrary to the case of the LFC in the linear response function, the insertion of the LFC here cannot give an exact result as further terms are missing. Equations and 10 are easy-to-compute solutions for the case of a harmonically perturbed electron gas since the static LFC at the parameters of interest is readily available from a ML representation, which is based on QMC simulation results.[62,63]

Finally, we note that there are no satisfactory analytical results for the ideal cubic response function at the first harmonic χ0(1,cubic) (q). Nevertheless, there is a formal relation between χ0(1,cubic) and the cubic response function of correlated electrons, which follows from perturbative analysis based on the Green function method[16]

Note that the mean-field result for the cubic response function at the first harmonic follows from eq by setting G = 0

In this work, we use the KS-DFT method to compute the set of density response functions defined by eqs –4 and subsequently verify the quality of the KS-DFT results by comparing with PIMC results at θ = 1. Then, we use KS-DFT results to analyze the analytical approximations given by eqs –10 in the wide range of parameters inaccessible for QMC methods. This allows us to unambiguously assess the importance of the negligible higher order (nonlinear) screening and LFC effects.

Simulation Details

The computational workflow consists of four main steps: First, the thermal KS-DFT simulations[20] of the free electron gas perturbed by an external field V(r) = 2Acos(q·r) are performed for different A and q values; second, the wave functions from KS-DFT simulations are used to compute the total density distribution along the direction of the wave vector q; third, the density perturbation components in k space ⟨ρ̂⟩ are computed using eq ; and finally, the density response functions are found by fitting data for ⟨ρ̂⟩ using eqs –4.

To begin with, we consider a strongly correlated electron gas with rs = 6 at θ = 1 and rs = 5 at θ = 0.01. At rs = 6, we compare the results with the available finite temperature PIMC data for the linear and nonlinear density response functions.[16] At rs = 5, we compare with diffusion quantum Monte Carlo (DMC) results for the linear density response function computed by Moroni, Ceperley, and Senatore.[64] Furthermore, we investigate a metallic density with rs = 2 at three different values of the degeneracy parameter, θ = 1, θ = 0.5, and θ = 0.01. In this case, we also benchmark results against PIMC data at θ = 1 and compare with the linear density response function from the DMC simulations at θ = 0.01.

The KS-DFT simulations of the free electron gas were performed using the GPAW code,[68−71] which is a real-space implementation of the projector augmented-wave method. The number of particles in the main simulation box varied in the range from 14 to 66. Accordingly, the main cubic cell size is defined by rs and N as L = rs(4/3πN)1/3. The direction of the perturbation is set to be along the z-axis. Due to periodic boundary conditions, the value of the perturbation wavenumber of the external harmonic field is defined by the reciprocal lattice vectors of the main simulation cell q = η × 2π/L, with η being a positive integer number. We used a Monkhorst-Pack[72] sampling of the Brillouin zone with a k-point grid of Nk × Nk × Nk total points, with Nk = 12 at rs = 2 and Nk = 8 at rs = 6 and rs = 5. The calculations were performed using a plane-wave basis where the cutoff energy has been converged to 800 eV at θ = 1 and rs = 2 and to 440 eV at the rest of the rs and θ values. The number of orbitals is set to Nb = 500 at rs = 2 and θ = 1 with the smallest occupation number fmin ≲ 10–7. We set Nb = 240 at rs = 6 and θ = 1 and Nb = 140 at rs = 2 and θ = 0.5 (fmin ≲ 10–6). At θ = 0.01, we set Nb = 70 for N = 66 particles and Nb = 2 N for N = 20 and N = 14 particles (with fmin = 0).

At rs = 2, the perturbation amplitudes are set in the range from A = 0.01 to A = 0.1 with a step of ΔA = 0.03 (here, A is in Hartree atomic units). At rs = 6 and rs = 5, the perturbation amplitudes are in the range from A = 0.002 to A = 0.017 with the step ΔA = 0.005. These values of the perturbation amplitudes used for the calculation of the density response functions were found empirically guided by the requirement Δn/n0 ≪ 1 and by testing the validity of eqs –4. Examples of the dependence of the density perturbation ⟨ρ̂⟩ on A and the application of eqs –4 are illustrated in the Appendix.

The XC functional in our KS-DFT simulations is the local density approximation (LDA) in the Perdew–Zunger parametrization.[73] Recently, it was demonstrated for θ = 1 that commonly used GGA functionals such as PBE,[74] PBEsol,[75] AM05,[76] and the meta-GGA functional SCAN[77] are not able to provide a superior description compared to LDA in the case of the free electron gas perturbed by an external field with fixed wavenumber q when Δn/n0 < 1.[78,79] We do not aim to further study this problem in this work. Therefore, we do not consider other types of XC functionals beyond LDA.

In addition to the LDA-based calculations, we performed simulations with zero XC functional [to which we refer to as DFT (RPA)]. This allows us to obtain exact results for the density response on the mean-field level. The value of this type of KS-DFT calculations allows us to assess the accuracy of the corresponding theoretical mean-field expressions given in eqs and 8. Furthermore, once the analytical results have been verified, the KS-DFT calculations on the mean-field level can be combined with the LFC to compute a highly accurate response function. We demonstrate that this is the case for the quadratic response function and the cubic response function at the first harmonic using eqs and 11.

Results

Strongly Correlated Hot Electrons

We start the discussion of our simulation results by considering the strongly correlated electron gas with the density parameter rs = 6 and at the reduced temperature θ = 1. This corresponds to WDM generated in evaporation experiments.[80] At these parameters, we can benchmark the KS-DFT calculations against previous PIMC calculations.[6,16]

Linear Density Response in the WDM Regime

First, we verify that our KS-DFT calculations provide accurate data for the linear static density response function χ(1) (q), which allows us to systematically analyze and exclude the possibility of finite size effects.[81] We start this analysis by comparing the linear static density response function computed using KS-DFT with the exact analytical results on the mean-field level and with χ(1) of the correlated electron gas computed using the exact data for the LFC.[62]

In Figure , we present the KS-DFT results computed using N = 14 electrons. In this case, the cell size is L = 12.335 Å, and the accessible values of the wavenumbers are multiples of qmin ≃ 0.8427qF.

Figure 1

Linear static density response function at rs = 6 and θ = 1.

From Figure , first of all, we see that the linear density response function computed using KS-DFT with zero XC functional, χ̃RPA(1), accurately reproduces the exact random phase approximation (RPA) result for the static linear response function,where v(q) = 4π/q2 and χ0(1) (q) is the Lindhard function. This shows that finite size effects in our KS-DFT calculations with as few as 14 electrons is negligible in the considered case. For completeness, we note that this is consistent with previous findings of PIMC simulations at similar conditions.[62]

Second, we combine the density response function computed using the KS-DFT with zero XC functional, χ̃RPA(1), with the LFC to find the linear density response function of correlated electronswhere G(q) is computed using the ML representation of the LFC.[62]

From Figure , we see that χ̃LFC(1) (q) [labeled as DFT(RPA) + LFC] is in excellent agreement with the exact value computed using the Lindhard functionwhere after the second equality, we used eq to express χLFC(1) (q) in terms of χRPA(1) (q).

Note that eq follows from eq after the substitution χ̃RPA(1) → χRPA(1) and χ̃LFC(1) → χLFC(1).

Now, after verifying that our calculations are not affected by the finite size effect, we compare the LDA-based KS-DFT calculations of the linear density response function, χ̃LDA(1) (q), with the exact result χLFC(1) (q). From Figure , we see that χ̃LDA(1) (q) is in good agreement with χLFC(1) (q) at q < 2qF (with qF being the Fermi wavenumber) and exhibits significant disagreements at q > 2qF. To understand this finding, we recall that the LDA corresponds to the long wavelength approximation of the LFC with GLDA = γk2, where γ is defined by the compressibility sum rule.[82] This approximation is applicable at q ≲ qF and increasingly deviates from the exact result with the increase in the wavenumber beyond 2qF.[78] Note that all χ̃LDA(1) (q), χ̃LFC(1) (q), and χLFC(1) (q) tend to χRPA(1) (q) in the limit of large wavenumbers since the screening factor dominates over XC effects in this limit. This can be seen from eqs and 15, where the LFC is suppressed by the factor q–2.

The insights that we have gained considering χ̃RPA(1), χ̃(1) (q), and χ̃LDA(q) will help us understand the KS-DFT results for the higher-order nonlinear density response functions discussed in the next subsection. Furthermore, we use a tilde over a symbol to differentiate response functions calculated using KS-DFT from the theoretical definitions.

Nonlinear Density Response in the WDM Regime

Our new results for the nonlinear density response functions at rs = 6 and θ = 1 are presented in Figure . In particular, the left panel shows the results for the quadratic response function (defined by eq ), the middle panel presents data for the cubic response function at the first harmonic (the cubic term in eq ), and the right panel shows the results for the cubic response function at the third harmonic (defined by eq ).

Figure 2

Nonlinear static density response functions at rs = 6 and θ = 1. Left: quadratic response function at the second harmonic. Middle: cubic response function at the first harmonic. Right: cubic response function at the third harmonic. We note that analytical results for the cubic response function at the first harmonic are presently not available.

First of all, we observe that the LDA XC functional-based calculations are generally in good agreement with the PIMC results at q < 2qF and overestimate the considered nonlinear density response functions at q > 2qF. The reason for this behavior of the LDA-based calculations is the inaccuracy of the LFC incorporated in the LDA at q > 2qF as it has been discussed in Section above.

Next, from the left panel of Figure , we see that KS-DFT calculations with XC set to zero [i.e. DFT(RPA)] are in excellent agreement with the theoretical RPA curve for the quadratic response function. This confirms the high quality of the analytical result eq in the WDM regime. In the case of the cubic response function at the third harmonic, as shown in the rightmost panel of Figure , we observe that eq is accurate at q > 2qF but overestimates the response at q < 2qF. Note that the LDA-based KS-DFT results, the KS-DFT calculations without XC, DFT(RPA), and the PIMC results all have positive sign at q < qF, while the theoretical curves fail to capture the change in the sign of the cubic response function at the third harmonic with decrease in wavenumber.

Let us next combine the KS-DFT data for the quadratic response computed with zero XC functional, χ̃RPA(2) (q), with the LFC. For that, we express χLFC(2) (q) via χRPA(2) (q) using eqs and 9. Then, we perform substitutions χ̃LFC(2) (q) → χLFC(2) (q) and χ̃RPA(2) (q) → χRPA(2) (q). As the result, we have the following relation

where χ̃0(1) (q) can be extracted from χ̃RPA(1) (q) as

Comparing the results calculated using eq with the PIMC data, we conclude that the relationship (eq ) is fulfilled with high accuracy.

Similarly, we derive the connection between χ̃LFC(1,cubic) (q) and χ̃RPA(1,cubic) (q) using eqs and 12 and replacing χLFC(1,cubic) (q) → χ̃LFC(1,cubic) (q) and χRPA(1,cubic) (q) → χ̃RPA(1,cubic) (q). As the result, we find

Using χ̃RPA(1,cubic) (q) obtained from KS-DFT simulations with zero XC and the LFC computed using the ML representation,[62] we have found that eq reproduces the PIMC results in the entire range of the wavenumbers as it can be seen from the middle panel of Figure . It is only the availability of χ̃RPA(1,cubic) (q) that allows us to estimate the cubic response at the first harmonic with PIMC accuracy as no analytical theory for χ0(1,cubic) currently exists.

To further explore the combination of the KS-DFT calculations with zero XC functional and the ML representation of the LFC, we next analyze the quality of the theoretical result (eq ) for the cubic response at the third harmonic. Using eqs and 10 and replacing χLFC(3) (q) by χ̃LFC(3) (q) and χRPA(3) (q) by χ̃RPA(3) (q), we arrive at the following relation between χRPA(3) (q) computed using the KS-DFT calculations with zero XC functional and the LFC

The comparison of χ̃LFC(3) (q) with the PIMC results is presented in the right panel of Figure . From this figure, we see that χ̃LFC(3) (q) significantly deviates from the PIMC data at q < 2qF. This means that the relation (eq ) does not provide an adequate description of the correlated electron gas. This is expected since we have already demonstrated above that the RPA result (eq ) is inadequate at q < 2qF. Therefore, the description of the screening on the mean-field level must first be improved to describe the actual system.

Strongly Correlated and Strongly Degenerate Electrons

Next, we investigate the strongly degenerate case with rs = 5 and θ = 0.01. In this regime, we are able to verify our KS-DFT calculations by comparing with the accurate DMC calculations of the linear static density response function, χ(1), by Moroni, Ceperley, and Senatore.[64]

Additionally, we further assess possible finite size effects at low temperature by comparing the simulation results for N = 14 particles to the results computed using N = 20, N = 38, and N = 66 particles. In this case, the cell size is L = 10.28 Å (for N = 14), L = 11.577 Å (for N = 20), L = 14.34 Å (for N = 38), and L = 17.236 Å (for N = 66). Correspondingly, the accessible values of the wavenumbers are multiples of qmin ≃ 0.8427qF (for N = 14), qmin ≃ 0.74822qF (for N = 20), qmin ≃ 0.604qF (for N = 38), and qmin ≃ 0.50qF (for N = 66).

Linear Density Response in the Limit of Strong Degeneracy

In Figure , we present results for the static linear density response function. Evidently, the results for χ̃RPA(1) computed using different numbers of particles accurately reproduces the exact mean-field level result χRPA(1), eq . Furthermore, at all considered numbers of particles, the combination of χ̃RPA(1) with the LFC using eq allows us to reproduce the exact result given by eq . Therefore, the reduction of the number of particles from N = 66 to N = 38, then to N = 20, and further to N = 14 does not lead to a deterioration of the quality of the data for χ̃RPA(1). This confirms the remarkable convergence of the KS-DFT simulations for as few as N = 14 particles.

Figure 3

Comparison of the DMC data by Moroni, Ceperley, and Senatore[64] with DFT results for the linear response function at rs = 5 and θ = 0.01.

To get a picture about the quality of the LDA-based calculations, we compare χ̃LDA(1) with the DMC results by Moroni et al.[64] and with χLFC(1) computed using the ML representation of the LFC by Dornheim et al.[62] Despite the fact that the LDA is designed to describe only the long wavelength limit of the LFC, we observe that in the strongly degenerate case, the LDA-based KS-DFT calculations provide high-quality results for the linear density response function with a level of accuracy similar to the ground-state QMC calculations.

Nonlinear Density Response in the Limit of Strong Degeneracy

After successfully testing the accuracy of our KS-DFT simulations on the linear density response function, we analyze results for the higher order density response functions presented in Figure . In the left panel, we see that χ̃RPA(2) is in excellent agreement with χRPA(2) and that χ̃LFC(2) also reproduces χLFC(2). This confirms the correctness of the analytical results for the quadratic response function given by eqs and 9 in the limit of strong degeneracy. The LDA-based data χ̃LDA(2) provides an adequate description of the quadratic response function and captures the effect of the stronger response when XC effects are included compared to the mean-field level results χRPA(2) and χ̃RPA(2). Certain quantitative disagreements between χ̃LDA(2) and χLFC(2) can be understood by noting that accurate data for LFC (beyond LDA) is needed to correctly describe the quadratic response. In fact, Dornheim et al.[57] recently pointed out that the quadratic response is directly related to three-body correlations, which explains this sensitivity to XC effects.

Figure 4

NLRT functions at rs = 5 and θ = 0.01. Left: quadratic response function at the second harmonic. Middle: cubic response function at the first harmonic. Right: cubic response function at the third harmonic. We note that analytical results for the cubic response function at the first harmonic are presently not available.

The middle panel of Figure presents data for the cubic response at the first harmonic. Compared to the partially degenerate case with θ = 1, the results exhibit a much sharper peak and much stronger response at 1.5qF < q < 2qF. At these wavenumbers, the difference between χLFC(1,cubic) and χ̃LDA(1,cubic) are most likely a direct consequence of the fact that the cubic response depends on the fourth power of the LFC, meaning that any deviations from the correct response in the first order gets amplified when a higher-order response is considered.

The right panel of Figure shows the cubic response at the third harmonic. In this case, the theoretical result for the response at the mean-field level given by eq , χRPA(3), fails to reproduce the exact data χ̃RPA(3). This extends the conclusion that eq fails to correctly describe the screened response from the WDM regime considered earlier to the case of strong degeneracy. As a consequence, being built upon eq , the LFC result χLFC(3) defined by eq also does not provide the correct description of the cubic response of the correlated electron gas at the third harmonic. This makes the analysis based on the comparison of χ̃LFC(3) and χ̃LDA(3) less meaningful. Since the LDA is an approximation to the true XC effects, χ̃LDA(3) cannot be considered to be the exact result. Nevertheless, it provides the correct quantitative outcome. Particularly, we see from the right panel of Figure that XC effects lead to a stronger response of the system compared to χRPA(3). We stress that χ̃RPA(3) is still the exact ab initio result for the cubic response on the mean-field level. Therefore, it can be used to verify theoretical derivations. Once a more accurate theoretical result for χRPA(3) that includes nonlinear screening effects that are negligible in eq is derived, the correct way to include the LFC should directly follow.

Free Electron Gas at Metallic Density

As a particularly important regime from the point of view of applications, we next consider rs = 2, which is a characteristic metallic density. In this case, we investigate three different values of the degeneracy parameter, namely, θ = 1, θ = 0.5, and θ = 0.01. For θ = 1, we have performed series of calculations with N = 14, N = 20, and N = 34. At θ = 0.5, we have considered N = 14 and N = 34 particles. At θ = 0.01, we have performed simulations with N = 14, N = 20, N = 38, and N = 66 particles. In agreement with the calculations in the strongly correlated case, there is no noticeable finite size effect for rs = 2 at these numbers of particles.

Linear Density Response

In Figure , we present results for the linear density response function at rs = 2 and compare the KS-DFT data with PIMC results and with χLFC(1) at θ = 1 in the left panel. We observe that the LDA-based results χ̃LDA(1) are in good agreement with both the PIMC data and χLFC(1). At θ = 0.5, too, we find good agreement of χ̃LDA(1) with χLFC(1) as it can be seen from the middle panel of Figure . In the limit of θ → 0, we compare χ̃LDA(1) with the DMC data by Moroni et al.[64] as well as with χLFC(1). Evidently, the LDA-based KS-DFT simulations provide an accurate description of the linear density response function in the strongly degenerate case too. We note that KS-DFT is more accurate in particular for q > 2qF compared to the previously considered case of rs = 6 due to the reduced impact of electronic XC effects at the higher density.

Figure 5

Linear static density response function for rs = 2 at θ = 1 (the left panel), at θ = 0.5 (the middle panel), and at θ = 0.01 (the right panel).

The comparison of χ̃RPA(1) with χRPA(1) that is also presented in Figure confirms the high accuracy of the KS-DFT results for the description of the linear response function on the mean-field level across temperature regimes. As a consequence, χ̃LFC(1) is in an excellent agreement with χLFC(1) over the entire wavenumber range.

Quadratic Density Response

The results for the quadratic response function at rs = 2 are shown in Figure . The quadratic response χ̃LDA(2) closely reproduces both the PIMC data and χLFC(2) at θ = 1 as it is demonstrated in the left panel of Figure . The agreement between χ̃LDA(2) and χLFC(2) somewhat deteriorates with the decrease in the temperature from θ = 1 to θ = 0.5 and further to θ = 0.01. This is shown in the middle and right panels of Figure . Nevertheless, the LDA, which is an XC functional purely based on the uniform electron gas model, provides an overall impressively accurate description of the quadratic response function at all considered wavenumbers of the perturbation.

Figure 6

Quadratic static density response function for rs = 2 at θ = 1 (the left panel), at θ = 0.5 (the middle panel), and at θ = 0.01 (the right panel).

Similarly, the discussed cases of the strongly coupled electrons, χ̃RPA(2) and χRPA(2) are in close agreement with each other at θ = 1, θ = 0.5, and θ = 0.01. This is a clear illustration of the high accuracy of the theoretical result eq for the mean-field description. Consequently, we find almost the same result using χ̃LFC(2) and χLFC(2).

From comparing amplitudes of the quadratic response function in Figure at different temperatures, one can see that the response of the system at the second harmonic becomes stronger upon decreasing the temperature of the electrons. For example, the decrease in the temperature from the partially degenerate case (θ = 1) to the strongly degenerate case (here represented by θ = 0.01) leads to an increase in the maximum value of the quadratic response function by a factor of 2.5. For comparison, the amplitude of the linear response function increases about two times with the decrease in the temperature from θ = 1 to θ = 0.01 at the same conditions.

Cubic Density Response at the First Harmonic

Let us next investigate the cubic response function at the first harmonic, which is shown in Figure . The comparison with the PIMC data at θ = 1 is provided in the left panel of Figure . From this comparison, it is evident that the LDA-based KS-DFT calculations are able to provide a very accurate description of the cubic response at rs = 2 and θ = 1 as they are in good agreement to both the PIMC data and χ̃LFC(1,cubic). This is an indication that the relation (eq ) holds since χ̃LFC(1,cubic) is computed using the KS-DFT data χ̃RPA(1,cubic) and the LFC according to eq .

Figure 7

Cubic static density response function at the first harmonic for rs = 2 at θ = 1 (the left panel), at θ = 0.5 (the middle panel), and at θ = 0.01 (the right panel).

Decreasing the electronic temperature leads to a substantial increase in the amplitude of the cubic response function at the first harmonic. This is visible from the comparison of the amplitudes of the cubic response at θ = 1 (left), at θ = 0.5 (middle), and at θ = 0.01 (right) in Figure . For example, the maximum of the cubic response function at the first harmonic at θ = 0.01 is about 16 times larger than at θ = 1. The consequences of this remarkable behavior are discussed in more detail in Section .

Cubic Density Response at the Third Harmonic

Finally, we present the results for the cubic response at the third harmonic in Figure , and the left panel depicts the comparison with PIMC data at θ = 1. Evidently, the LDA-based calculations χ̃LDA(3) are in good agreement with the latter. Next, χRPA(3) overestimates the strength of response compare to χ̃RPA(3) at q < qF. In contrast to the response at stronger coupling (rs = 6 and θ = 1), at rs = 2, we do not see the change in the sign of the cubic response at the third harmonic upon increasing the wavenumber from small q < qF to large q > qF values. However, this behavior is restored with a decrease in the temperature to θ = 0.5 as it is clearly visible in the middle panel of Figure . Importantly, such behavior at θ = 0.5 is not captured by the analytical results given in eqs and 10. At q < 1.5qF, this leads to a significant disagreement between the analytical result in the mean-field level approximation, eq , and the exact result computed using zero XC in the KS-DFT simulations. However, in the limit of strong degeneracy depicted in the right panel of Figure , eq again provides the qualitatively correct result by reproducing the change in the sign of the corresponding response function but remains in quantitative disagreement with χ̃RPA(3) in the vicinity of the positive maximum and the negative minimum. These disagreements between eq and the DFT (RPA) results could be caused by the approximation of the screening on the linear response level in eq as it was shown by Dornheim et al.[16]

Figure 8

Cubic static density response function at the third harmonic for rs = 2 at θ = 1 (the left panel), at θ = 0.5 (the middle panel), and at θ = 0.01 (the right panel).

Interestingly, despite the poor performance of eq , the agreement between χ̃LDA(3) and χ̃LFC(3) is rather good at all considered temperature regimes. This is explained by the fact that XC effects are less pronounced compared to the above considered strongly correlated cases.

Conclusions and Outlook

In this work, we have explored a new methodology for the study of the nonlinear electronic density response based on KS-DFT. As a particular example, we have investigated the free electron gas model and demonstrated that the KS-DFT method is an effective and valuable tool for the investigation of various nonlinear electronic density response functions across temperature regimes. This conclusion is important for parameters where QMC methods experience significant difficulties or fail to converge at all due to the fermion sign problem. This approximately corresponds to θ < 1 and rs ≳ 2.[33,36,37] A particularly effective method to gauge and guide the development of new theoretical approaches is given by the KS-DFT simulation with zero XC functional. This is due to the fact that theoretical models of correlated electrons are built upon mean-field approximations, in combination with the electronic LFC. Therefore, an accurate check of the quality of the analytical results for the nonlinear density response functions in the mean-field approximation is required to build a reliable theory. The presented DFT-based methodology provides such a tool in a wide range of parameters.

As a demonstration of the KS-DFT method-based analysis of the theoretical results, we have considered the quadratic and cubic response functions at different values of the density and degeneracy parameters. First of all, we have confirmed the validity of the analytical results for the quadratic response function (eqs and 9) for partially to strongly degenerate electrons. This confirms and complements the earlier PIMC-based analysis at θ ≥ 1.[16] Second, it has been shown that the analytical results for the cubic response function at the third harmonic (Eqs and 10) are quantitatively inaccurate at q ≲ 1.5 qF for all considered values of θ. Moreover, eqs and 10 are qualitatively inadequate at θ = 0.5.

The application of the KS-DFT method to study the cubic response at the first harmonic (as defined in eq ) has allowed us to observe a change of its characteristics with the decrease in the temperature to θ < 1. We have revealed that the decrease in the temperature from the partially degenerate regime with θ ∼ 1 to the strongly degenerate regime (θ ≪ 1) leads to a significant increase in the maximum of the cubic response at the first harmonic. Let us demonstrate the implication of this finding for electrons at a metallic density, rs = 2. Using eqs –4 and the data presented in Figures –8, one can deduce that, at rs = 2, both the quadratic response and the cubic response at the third harmonic remain inferior to the linear density response function if the perturbation amplitude A < 1 at all considered θ values and wavenumbers of the perturbation. In contrast, the cubic response function at the first harmonic becomes dominant over the linear response function if A > 0.8 at θ = 0.5 and if A > 0.36 at θ = 0.01. The applicability of the perturbative analysis requires the smallness of the higher order correction compared to the first order term in eq . Therefore, at θ ≪ 1, not the quadratic response but the cubic response at the first harmonic leads to the strongest restriction on the applicability of the nonlinear density response theory of free electrons with respect to the perturbation amplitude.

Another important finding is that the LDA XC functional provides a remarkably accurate description of the linear and nonlinear density response at metallic density, rs = 2, across the entire considered temperature range. This strongly indicates that our new KS-DFT-based approach constitutes a reliable tool for the investigation of the nonlinear density response both at ambient conditions and in the WDM regime. At stronger coupling parameters, rs = 6 and rs = 5, the LDA XC functional-based KS-DFT calculations of the correlated electron gas provide qualitatively correct behavior but ought to be considered with caution from a quantitative point of view.

We conclude this study by pointing out that the present methodology is very general and can be directly applied to arbitrary materials; the only approximation is given the choice of the XC-functional. In particular, simulations including ions are much more problematic and computationally expensive for the QMC methods compared to the considered case of a free electron gas model. Therefore, the KS-DFT method is particularly valuable for the multi-component systems. This proof of concept of the capability of KS-DFT for the estimation of the NLRT of an electron gas is thus a pivotal first step before extending our considerations to real materials.