Assessing the accuracy of hybrid exchange-correlation functionals for the density response of warm dense electrons

The Kohn–Sham density functional theory (KS-DFT) is the most popular first-principles electronic structure method both at ambient and extreme conditions. This includes simulations of the warm dense matter (WDM),^1–4 where KS-DFT is widely used^5–7 for the interpretation of X-Ray Thomson scattering experiments.^8–10 WDM is characterized by high temperatures and densities and is generated using different techniques¹¹ at large research facilities like the European X-ray Free-Electron Laser (XFEL),¹² the National Ignition Facility (NIF),^13,14 the Linac coherent light source at SLAC,¹⁵ or the Sandia Z-machine.¹⁶ Both experimental and theoretical investigations of WDM are of high current interest due to its relevance for astrophysics, material science, and inertial confinement fusion.^14,17 Examples of astrophysical applications include white dwarfs¹⁸ and planetary interiors.^19–21 In addition, WDM is relevant for the discovery of novel materials at extreme conditions^19,22–24 and for hot-electron chemistry.²⁵ 

The rigorous understanding of the physics and chemistry of the WDM state is still in a stage of active development, since most widely used electronic structure calculations and diagnostic tools were developed for ambient conditions and require a consistent extension to high temperatures. For example, in the case of the KS-DFT, the key ingredient for an accurate simulation is the exchange-correlation (XC) functional. At ambient conditions, the quality of different XC functionals is tested by comparing them with the experimental data, e.g., with band gaps,²⁶ lattice constants,²⁷ and atomization energies.²⁸ 

In WDM simulations, it is a common practice to use ground-state XC functionals with thermal excitations included by the smearing of band occupation numbers according to the Fermi–Dirac distribution;^29,30 this is often referred to as zero-temperature approximation in the literature. Although there has been a certain progress in the development of finite-temperature XC functionals,^31–35 testing XC functionals at WDM conditions is challenging. First, there are significant uncertainties in the experimentally realized temperature and density.^9,10,36 Second, besides the very recent path integral Monte Carlo (PIMC) simulation results for hydrogen at certain WDM parameters,^37,38 the exact simulation data for complex materials like compressed solids at WDM conditions^24,39,40 that would enable reliable benchmarks of XC functionals are sparse.

This situation is highly unsatisfactory because the KS-DFT results are often used for the interpretation of experimental observations^6,24 and as an input to other methods and models, including machine-learning models,⁴¹ XC kernels for linear-response time-dependent density functional theory calculations,⁴² molecular dynamics with screened effective potentials,^43–49 and quantum hydrodynamics.^50–52 

Clearly, an extensive benchmark and analysis of the accuracy of KS-DFT calculations based on various XC functionals at high temperatures are needed. Recently, Moldabekov et al.^42,53,54 have reported the analysis of the local-density approximation (LDA),⁵⁵ PBE,⁵⁶ PBEsol,⁵⁷ AM05,⁵⁸ and SCAN⁵⁹ XC functionals at temperatures and densities relevant for WDM by considering an inhomogeneous electron gas in an external harmonic perturbation and comparing DFT results with exact PIMC reference data.^60–62 This has provided important insights to the performance of these functionals at WDM conditions. For example, in contrast to the ground state calculations, it was found that at high temperatures, SCAN performs worse than LDA, PBE, and PBEsol. Therefore, it is important to test the accuracy of XC functionals at WDM conditions. Furthermore, for the examples of the warm dense electron gas and hydrogen, it has been pointed out that even finite-temperature versions of the LDA and generalized gradient approximation (GGA) level functionals—created to extend the corresponding ground state versions to high temperatures—do not always lead to improved results at WDM conditions.⁴² Recently, Moldabekov et al.⁴² have analyzed the static XC kernel and density response function of the uniform electron gas (UEG) and hydrogen at WDM conditions using a finite-temperature LDA by Groth et al.³² It was found that for the UEG and warm dense hydrogen, the finite-temperature LDA provides less accurate results compared with the ground-state LDA. This seemingly surprising result was explained by a reduced validity (regarding wave number) of the long-wavelength approximation for the static XC kernel with an increase in temperature. It is clear that the explicit inclusion of the temperature dependence in an XC functional does not guarantee an improved description compared with its counterpart in the ground state. This highlights the importance of accurately testing and benchmarking XC functionals at WDM conditions.

On the other hand, a generalized KS-DFT⁶³ with orbital-dependent XC functionals has emerged as one of the most promising tools in the density-functional theory. The most often used orbital-dependent XC functionals^63,64 belong to the class of the hybrid functionals, which entail a certain weighted combination of the exact Hartree–Fock exchange contribution and of a standard XC functional like the LDA or PBE.^28,65 These orbital-dependent XC functionals naturally contain information about non-locality effects and reduce the self-interaction error.⁶⁶ Moreover, they provide a superior description of various material properties like lattice constants²⁷ and atomization energies.²⁸ 

Hybrid XC functionals—within generalized KS-DFT at finite temperature^31,35,67—can potentially overcome some problems of finite-temperature XC functionals since, in addition to the reduction of the self-interaction error, they contain extra information about thermal effects through orbitals in the exact exchange part. Therefore, in this work, we extend the previous benchmark and analysis of the LDA, GGA, and meta-GGA functionals for the harmonically perturbed electron gas at WDM conditions^53,54 to include a number of popular hybrid XC functionals. In our assessment, we consider the HSE06,^68,69 HSE03,⁷⁰ PBE0,⁷¹ PBE0-1/3,⁷² and B3LYP⁷³ hybrid functionals.

First, we consider the weakly perturbed electron gas and present results for the static hybrid XC kernel K_xc(q) of the uniform electron gas (UEG) at both ambient and WDM conditions using the recent framework by Moldabekov et al.⁴² In particular, K_xc(q) constitutes a wave-vector resolved measure for electronic XC effects and, therefore, is uniquely suited to benchmark the performance of XC functionals. In addition, the kernel constitutes a key input for a gamut of applications such as quantum hydrodynamics^50,51 or the interpretation of X-Ray Thomson-Scattering (XRTS) experiments.^74–77 

Second, we consider a strongly inhomogeneous, warm dense electron gas. The analysis is carried out by comparing with results for the PBE functional and exact PIMC reference data. PIMC allows us to unambiguously quantify the accuracy of the considered hybrid XC functionals, whereas the comparison with the PBE functional gives us insight into the physics behind the observed behavior of the hybrid functionals at different conditions.

The paper is organized as follows: we begin with providing the theoretical background in Sec. II, the details of our KS-DFT and PIMC calculations are given in Sec. III, the results of the calculations are presented in Sec. IV, and the paper is concluded in Sec. V with a summary of the main findings and a brief outlook on future research.

In Eq. (5), the implicit dependence on the electronic temperature is taken into account in $FxHF[ρσ]$ due to the smearing of the occupation numbers. Similarly, the usage of orbitals in the Fermi–Dirac distribution of the occupation numbers in the PBE0-1/3, HSE06, HSE03, and B3LYP allows for incorporating the implicit portion of thermal effects. Clearly, this is not a fully consistent way of generating a finite-temperature hybrid XC functional. Nevertheless, e.g., the HSE03 has been successfully used for the description of the experimental X-ray spectra of aluminum⁸⁰ and the reflectivity and dc electrical conductivity of liquid ammonia⁸¹ at high temperatures. Furthermore, it turns out that the straightforward use of the orbitals with smeared occupation numbers according to the Fermi–Dirac distribution in the PBE0, PBE0-1/3, HSE06, and HSE03 provides a significant improvement in the description of the electronic density response at high temperatures T ∼ T_F compared with PBE as reported in Sec. IV. Note that in the case of the HSE06 and HSE03, the Coulomb potential in $FxHF[ρσ]$ is substituted by a screened potential 1/r₁₂ → erfc(Ωr₁₂)/r₁₂⁷⁹ to avoid the computational difficulties due to the long-range character of the Coulomb potential. Consistently, the screened version for the exchange part of the PBE is also used. For example, Ω = 0.11 Bohr⁻¹ is used for the HSE03.⁶⁹ 

It allows us to investigate both the weakly and strongly inhomogeneous electron gas by tuning the external perturbation. As a measure of the density inhomogeneity, one can use the maximum value of the density perturbation $δnmax/n0=maxδn(r)/n0$⁠, where δn(r) = n(r) − n₀ is the density deviation from the UEG with density n₀.

The XC kernel constitutes the key input to linear-response time-dependent density functional theory. Clearly, the XC kernels from the hybrid XC functionals are also interesting in their own right since they provide insight into the XC features at different wavenumbers. We stress that the static XC kernels obtained from hybrid XC functionals for the UEG had not been reported prior to our work—neither for the ground state nor for WDM.

We use the available exact quantum Monte Carlo data for G(q) of the UEG to analyze various XC kernels.^{60,76,82,87–89}

Beyond the linear response theory, accurate data for the XC kernel (LFC) in combination with ideal density response functions allows one to compute higher order quadratic and cubic density response functions.^62,90,91

Commonly used parameters for the description of the electrons in the WDM regime are the coupling parameter r_s and the reduced temperature θ = T/T_F, where the former is the mean inter-particle distance in Hartree atomic units and the latter is the ratio of the temperature to the Fermi temperature.⁹² The WDM state is defined by the regime where electrons are non-ideal r_s ≳ 1, and thermal excitations and quantum degeneracy are simultaneously important, i.e., θ ∼ 1. Often in experiments, the WDM state is generated by shock-compression and (isochoric) heating of solid-state samples.^24,39,40 In this regard, understanding the WDM state around metallic densities is particularly important. Thus, we consider the characteristic metallic density of r_s = 2 and the partially degenerate case with θ = 1. This corresponds to an electron temperature T ≃ 12.5 eV. In the case of a weak perturbation δn_max/n₀ ≪ 1, we additionally consider the ground state by setting θ = 0.01. For this case, we analyze the quality of the considered hybrid XC functionals in the UEG limit in terms of quantum Monte Carlo data for the density response function by Moroni et al.⁸² Furthermore, the analysis of the ground state allows us to gauge the impact of thermal excitation effects on the performance of the hybrid XC functionals at T ≃ T_F.

For the perturbation in Eq. (6), we consider A = 0.01, A = 0.1, A = 0.5, and A = 1.0 (in Hartree). These parameters correspond to various degrees of the density inhomogeneity from δn_max/n₀ ≲ 10⁻² (at A = 0.01) to δn_max/n₀ ≳ 1 (at A = 0.5 and A = 1.0). Thus, A = 0.01 allows one to probe the properties of the system in the UEG limit. We particularly extract the static density response function and the static XC kernel following the approach introduced in Ref. 42 and compare them with the PIMC data.

The wavenumber of the perturbation q is directed along the z-axis. Note that we drop the vector notation in the following. The values of the perturbation wavenumber are defined by the length of the simulation cell as q = 2πj/L, where j ≥ 1 is a positive integer number. The smallest value of q corresponding to j = 1 is denoted as q_min. We consider q in the range from q_min to 4q_min and up to 5q_min at A = 0.01 for the static density response and XC kernel calculations.

The KS-DFT calculations were performed using the ABINIT package.^93–98 The calculations at θ = 1 were performed with the cubic main cell length L = 7.7703 Bohr, N_b = 200 orbitals, and the maximal kinetic energy cutoff 13 Ha. The used cell length corresponds to N = L³ × n₀ = 14 electrons, where at r_s = 2, we have n₀ ≃ 0.02984 Bohr⁻³. The smallest perturbation wavenumber for 14 particles is q_min[N = 14] ≃ 0.842 68 q_F, where q_F is the Fermi wavenumber. Periodic boundary conditions and a k-point grid of size 8 × 8 × 8 were used. Self-consistent-field cycles for the solution of the Kohn–Sham equations were converged with absolute differences of the total energy δE < 10⁻⁷ Ha. A large number of bands N_b ≫ N are needed due to significant thermal smearing of the occupation numbers. At the considered parameters, N_b = 200 corresponds to the smallest occupation number f_min < 10⁻⁴, which is sufficiently small to acquire convergence.⁵³ In prior works, it was shown that, in the case of the harmonically perturbed electron gas, the results for the response functions with N = 14 particles are not affected by finite size effects.^{42,53,99–102} In this work, keeping other parameters equivalent, the convergence is cross-checked by comparing the data for N = 14 with the results computed with N = 20 and N_b = 300 at q = 2q_min (with A = 0.01). Additionally, we used three different particle numbers in the simulation cell for the ground state calculations: N = 38, 54, and 66; with a k-point grid of 8 × 8 × 8. The agreement and consistency of the results for all considered hybrid XC functionals computed with different particle numbers reassure the convergent character of the calculations. Finally, for both θ = 1 and θ = 0.01, we performed KS-DFT calculations with the XC functional set to zero yielding the static density response function of the UEG in the random phase approximation (RPA) (see Sec. II B) and compared it with the analytical RPA (in N → ∞ limit). The perfect agreement of our KS-DFT data with the analytical RPA [shown in Figs. 1(a) and 2(a)] substantiates our confidence in the convergence of our calculations.

We report results for a total of 155 KS-DFT calculations for the HSE06, HSE03, PBE0, PBE0-1/3, and B3LYP hybrid XC functionals with different choices of parameters (amplitude A and wavenumber q) in the density perturbation.

A fraction of the PIMC data (with q = q_min, 2q_min, and 3q_min) used in this work has been reported in the analysis of the ground state LDA, PBE, PBEsol, AM05, and SCAN functionals in Ref. 53. It was shown there that the quality of the KS-DFT calculations based on these XC functionals has the tendency to deteriorate with an increasing wavenumber of the perturbation. Here, we test the hybrid XC functionals in the regime where the functionals in the previous analysis have failed. We have, therefore, generated new PIMC data for q = 4q_min at A = 0.1, 0.5, and 1.0.

Since a detailed introduction to the direct PIMC method¹⁰³ and the associated fermion sign problem¹⁰⁴ have been presented elsewhere, we here restrict ourselves to a brief summary of the main simulation details. We have used the extended ensemble approach introduced in Ref. 105, which is a canonical adaption of the worm algorithm by Boninsegni et al.^106,107 Moreover, we use P = 200 imaginary-time propagators from a primitive factorization scheme, which is sufficient to ensure convergence; see the supplementary material of Ref. 61 for a corresponding analysis.

We start with the weakly perturbed electron gas where |δn(r)|/n₀ ≪ 1 and consider both the ground state and the parameters relevant for WDM. We analyze the results both in terms of the static density response function χ(q) and the static XC kernel K_xc(q). In the next step, we consider strongly perturbed systems with |δn(r)|/n₀ ≳ 1 and provide the analysis of the accuracy of the electron density calculations based on the considered hybrid XC functionals.

First, the analysis of the T → 0 limit helps to establish a baseline for further analyzing the performance of the hybrid XC kernels at high temperatures. Second, for the construction of the LDA and GGA (e.g., PBE) XC functionals, the reproduction of the LFC of the UEG at q < 2q_F was one of the main criteria for an adequate description of metals.⁵⁶ Thus, it is interesting to see to what degree the considered hybrid XC functionals can reproduce the LFC of the UEG.

The results for the static density response function and the static XC kernel are shown in Figs. 1(a) and 1(b), respectively. In Fig. 1, we compare the KS-DFT results with the diffusion quantum Monte Carlo (DMC) data by Moroni et al. (black squares)⁸² computed for T → 0. In addition, we show the results from the machine-learning representation (ML-PIMC) by Dornheim et al. of extensive PIMC data at finite temperatures⁶⁰ and the DMC data⁸² for T = 0. In addition, the results based on the hybrid XC functionals, in Fig. 1, we also show results computed using PBE. The latter comparison is meaningful because PBE0, PBE0-1/3, HSE06, and HSE03 are based on PBE.

Before discussing our results on the hybrid XC functionals, we validate the absence of any significant finite size effects in Fig. 1(a) by comparing the analytical RPA with the static density response from the KS-DFT calculations where the XC functional is set to zero (labeled NullXC⁴²). The latter is equivalent to the RPA and does indeed agree with the analytical RPA as expected.

In Fig. 1(a), we observe that the results for χ(q) based on PBE0, PBE0-1/3, HSE06, HSE03, and PBE are closely grouped and are in rather good agreement with the reference DMC data. In contrast, the B3LYP results for χ(q) significantly deviate from the reference ML-PIMC data at q > q_F. This is expected since B3LYP does not attain the UEG limit by design. This was elucidated by Paier et al. for the unperturbed UEG.²⁷ Here, we further substantiate it for finite perturbation wavenumbers in the range 0.6 ≤ q ≤ 3q_F as it can be seen in Fig. 1(b) for the static XC kernel, where a strong deviation of the B3LYP results from the DMC and ML-PIMC data is clearly visible. Furthermore, from Fig. 1(b), we see that the B3LYP data do not obey the long wave-length expansion of the LFC Eq. (12) with quadratic q dependence of the LFC at q ≲ q_F; this can be seen particularly well in the inset that shows a magnified segment.

In contrast to B3LYP, the LFCs of PBE0, PBE0-1/3, HSE06, and HSE03 reproduce the quadratic dependence on the wavenumber at q ≲ 1.5q_F in accordance with Eq. (12), as shown in Fig. 1(b). We also observe in Fig. 1(b) that the PBE results have a quadratic dependence on q at all considered wave numbers since it was designed in a way that the gradient corrections to the exchange and to the correlation parts exactly cancel each other in the limit of the UEG.⁵⁶ 

At large wavenumbers, we observe from Fig. 1(b) that the PBE0, PBE0-1/3, HSE06, and HSE03 results deviate from the DMC data at q ≳ 2q_F. Nevertheless, at 2q_F ≤ q ≤ 3q_F, the results from the considered hybrid XC functionals exhibit smaller absolute deviations from the DMC data compared with the PBE results. An interesting feature of the LFC based on the PBE0, PBE0-1/3, HSE06, and HSE03 is that it has a maximum somewhere around 2q_F and a monotonically decreasing trend with increasing q at q > 2q_F. A similar behavior—with a somewhat worse agreement with the DMC data—was reported by Armiento and Mattsson for the AM05⁵⁸ functional, which is a semi-local GGA. Another related observation is that AM05 provides as accurate a description as the PBE0 and HSE06 hybrid XC functionals for lattice constants, bulk moduli, and surfaces of solids.¹⁰⁸ The latter property, i.e., the electronic behavior at surfaces, corresponds to a domain of strongly inhomogeneous electron densities characterized by large wave numbers q ≫ q_F.¹⁰⁹ Thus, it is an interesting observation that the AM05 XC kernel has a similar behavior as the XC kernel computed using the HSE06, HSE03, PBE0, and PBE0-1/3 functionals at q > 2q_F. We do not further extend the analysis of the observed similarities between the AM05 XC kernel and the hybrid XC kernels since it would deviate from the main focus of this work, but we leave such investigation as a possible direction for future research.

In Fig. 1(b), one can see that at q ≲ 1.5q_F, the PBE0, PBE0-1/3, HSE06, HSE03, and PBE functionals provide static XC kernels in excellent agreement with the DMC data. Additionally, the PBE0 and PBE0-1/3 results are much closer to the DMC data at q ≃ 2q_F compared with HSE06, HSE06, and PBE. At q ≳ 2q_F, PBE0 provides a much better description of the static XC kernel of the UEG compared with the other hybrid XC functionals.

In addition, we find that increasing the contribution of HF exchange from 1/4 (as in PBE0) to 1/3 (as in PBE0-1/3) lowers the resulting values of K_xc(q) at q > 2q_F. Keeping the coefficient 1/4 in front of the HF exchange fixed in Eq. (5), the same effect is achieved by using a screened potential instead of the Coulomb potential as illustrated by the static XC kernel based on the HSE06 and the HSE03. On the other hand, the behavior of the K_xc(q) at q < 2q_F is not sensitive to the screening of the Coulomb potential or the variation of the HF exchange contribution. One could adjust the screening parameter or the weight of the HF exchange contribution in order to accurately reproduce the exact data for the K_xc(q) in a wider range of wavenumbers.

Our new results for the static density response function and the static XC kernel of the UEG at θ = 1 are presented in Figs. 2(a) and 2(b), respectively. Similar to the ground state, we observe that the B3LYP results for the static density response function [Fig. 2(a)] and the static XC kernel [Fig. 2(b)] disagree significantly with the reference ML-PIMC data at q ≳ q_F. Also, the B3LYP data for the static XC kernel do not obey a quadratic dependence on the wavenumber at q ≲ q_F and approaches the PBE result from above with a decrease in q and crosses the PBE data at q ≃ 0.84q_F. The PBE data show a quadratic dependence on q in the entire considered range of wavenumbers and agree well with the reference data (ML-PIMC) at q ≲ 1.5q_F. The PBE0, PBE0-1/3, HSE06, and HSE03 results are in close agreement with the reference data (ML-PIMC) at q ≲ 2.5q_F (with the PBE0 being the most accurate). Therefore, these functionals lead to a better description of the density response compared with PBE. At q ≳ 2.5q_F, the PBE0, PBE0-1/3, HSE06, and HSE03 functionals produce static XC kernels that decrease with increasing q and eventually become negative at q ≳ 3.5q_F. A similar behavior has been reported for the AM05 functional at T ≃ T_F.⁴² In general, the trends in the q-dependence and main features of χ(q) and K_xc(q) of the considered hybrid XC functionals are similar to those in the ground state discussed above. This indicates that the partial occupation of the orbitals at high temperatures does not lead to fundamentally new features in the hybrid XC functional results.

With that, we conclude our consideration of the weak perturbation limit. We continue further with the case of the strong density inhomogeneity.

To generate a strong density perturbation, we apply the external harmonic perturbation as defined in Eq. (6) with amplitudes A = 0.1, 0.5, and 1.0 and wavenumbers in the range from q = q_min ≃ 0.84 q_F up to q = 4q_min. The amplitude of the density inhomogeneity is qualitatively illustrated in Fig. 3, where we plot the largest values of the density accumulation (δn_⊕) and the density depletion (δn_⊖) in the simulation cell due to the external harmonic field with different wavenumbers. The gray area indicates density perturbations less than 0.1n₀. We see that at the considered amplitudes of the perturbation, we always have |δn_⊕|≥ 0.1n₀ and |δn_⊖|≥ 0.1n₀. These density perturbation magnitudes are well beyond the linear response regime considered in Sec. IV A. At A = 0.5, we have |δn_⊕|∼ n₀ and |δn_⊖|∼ 0.5n₀. Furthermore, at A = 1.0, the electron density values almost vanish in the density depletion regions with δn_⊖ → −n₀ at q = q_min and q = 2q_min, where a physical lower bound δn_⊖ > − n₀ is in place (indicated by the dashed horizontal line in Fig. 3). Correspondingly, most of the electrons are localized in the density accumulation regions with δn_⊕ ≈ 2n₀.

An important point for our analysis is that at fixed values of A, the density perturbation magnitude can have a non-monotonic dependence on q due to a stronger response of electrons around q = 2q_min ≈ 1.5q_F as it is clearly illustrated in Fig. 3 for A = 0.1 and A = 0.5. A stronger density response at the wavelength corresponding to the mean inter-particle distance d ≃ 2r_s, i.e., at q ≃ 2π/d ≃ 2q_F, is a manifestation of the fact that it is energetically more favorable for particles to be localized at the distance r ≃ d from each other. This was demonstrated and used for the explanation of the roton feature in the excitation spectrum of the UEG by Dornheim et al.,.¹¹⁰ Contrarily, at fixed values of q, the increase in A always leads to a monotonic increase in the density perturbation amplitude. Therefore, to analyze the error of the hybrid XC functionals in the KS-DFT calculations, it is more transparent to compare the deviations from the exact PIMC data for various A values at a fixed q, rather than vice versa.

Before we dive into the details of the error analysis of the hybrid XC functionals, in Fig. 4, we illustrate the density perturbation δn(r) = n(r) − n₀ and the density error $Δñ$ profiles when the harmonic perturbation with A = 0.1 and q = 3q_min is applied. From Fig. 4(a), we observe that the PBE data exhibit a certain deviation from the PIMC data at the maximum and minimum values of δn(r). The B3LYP results show a significant disagreement with the PIMC data in both the density accumulation (with δn(r) > 1) and density depletion regions (with δn(r) < 1). The results for the other considered hybrid XC functionals are in a close agreement with the exact PIMC data. This is further reaffirmed by $Δñ$ in Fig. 4(b), where $|Δñ|≲1%$ for the PBE0, PBE0-1/3, HSE06, and HSE03 results. Additionally, from Fig. 4(b), we clearly see that the PBE (the B3LYP) data have an error of up to ±5% (±15%) in the relative density deviation $Δñ$⁠. Obviously, in this case, the PBE0, PBE0-1/3, HSE06, and HSE03 hybrid XC functionals provide a significant improvement due to mixing in a fraction of HF exchange.

To perform a general analysis of the errors in the density response, we consider the largest value of the $Δñ$ in the density accumulation region, $Δñ⊕$⁠, and the largest magnitude of the $Δñ$ in the density depletion region, $Δñ⊖$⁠, for different A and q.

In Fig. 5, we show the histogram representation of the largest value of the error in the relative density deviation as defined by $Δñ⊕$ and $Δñ⊖$ at q = q_min ≃ 0.84 q_F and different values of A. We observe that the B3LYP data have the smallest error with $Δñ⊕<1%$ and $Δñ⊖<1%$ compared with other KS-DFT data. However, it appears to be the crossing point between the B3LYP density response and the exact density response as the function of the parameters A and q, where one can say that the good performance of B3LYP is due to chance. We have already observed such a crossing point when considering the linear density response function in the case of the weak perturbation in Sec. IV A. Indeed, the B3LYP data significantly deviate from the reference PIMC data at other considered q > q_min values at A = 0.1, 0.5, and 1.0 (see Figs. 6–8). In fact, at A = 0.1 and q = 2q_min, 3q_min, and 4q_min, the B3LYP data have the largest errors in both density accumulation and depletion regions among all considered XC functionals. This is consistent with the conclusions drawn from the linear density response function in Sec. IV A. The reason is that at A = 0.1, with |δn_⊕|∼ 0.1n₀ and |δn_⊖|∼ 0.1n₀ (see Fig. 3), the linear density response function has the largest contribution to the total density response within non-linear response theory, which includes quadratic and cubic density response functions.^62,91 This reasoning is not applicable at A = 0.5 and A = 1.0 since density deviation from the UEG case is well beyond of the regular perturbation expansion, which is valid for |δn(r)|/n₀ ≪ 1.

From Fig. 5, one can see that at q = q_min ≃ 0.84 q_F and A = 0.1, 0.5, and 1.0, PBE0, PBE0-1/3, HSE06, and HSE03 have a similar performance as PBE, with $Δñ⊕∼−2%$ and $Δñ⊖∼1%$⁠. As shown in Fig. 6, at q = 2q_min ≃ 1.685 q_F, the PBE0, PBE0-1/3, HSE06, and HSE03 results are a bit more accurate than PBE with the improvement of the accuracy approximately between 1% and 2% at A = 0.5 and A = 1.0 in the density accumulation region.

At q = 3q_min ≃ 2.528 q_F and A = 0.1, we observe the largest improvement in the accuracy of the density response using PBE0, PBE0-1/3, HSE06, and HSE03 compared with PBE (see Fig. 7). At these parameters, the PBE data deviate significantly with about 5% in $Δñ⊕$ and $Δñ⊖$⁠. In contrast, the PBE0, PBE0-1/3, HSE06, and HSE03 results are highly accurate with deviations of $|Δñ⊕|≤1%$ and $|Δñ⊖|≤1%$⁠. From Fig. 7, one can see that at q = 3q_min and A = 0.5, PBE0, PBE0-1/3, HSE06, and HSE03 provide an improvement of a few percent in accuracy compared with PBE. At A = 1.0 and q = 3q_min, all considered XC functionals show similar accuracy with $|Δñ⊕|≲2%$ and $|Δñ⊖|≲2%$⁠.

Finally, in Fig. 8, we show the results for q = 4q_min ≃ 3.371 q_F. Similar to the case with q = 3q_min ≃ 2.528 q_F, PBE0, PBE0-1/3, HSE06, and HSE03 yield significantly higher accuracy than PBE for the density response at A = 0.1 and A = 0.5, with up to two times smaller errors in $|Δñ⊕|$ and/or $|Δñ⊖|$⁠. At last, for A = 1.0 and q = 4q_min, we observe similar accuracy for all considered XC functionals with $|Δñ⊕|≲4%$ and $|Δñ⊖|≲4%$⁠.

Before drawing general conclusions, we recall that in previous work, Moldabekov et al.⁵³ had assessed the accuracy of both the LDA and GGA functionals for the density response in the warm dense electron gas. They had shown that the PBEsol and LDA functionals were approximately as accurate as PBE for the density response in the partially degenerate case with T ≃ T_F for both weak and strong perturbations. The SCAN functional—being highly accurate in the ground state—had deviated significantly from the reference PIMC data at q > q_F.^42,53

In this work, the benchmark against the exact PIMC data at q ≲ 3q_F for both weak and strong perturbations demonstrates that the PBE0, PBE0-1/3, HSE06, and HSE03 hybrid XC functionals yield a more accurate description of the density response compared with PBE. In contrast, B3LYP is significantly less accurate for the perturbed electron gas compared with the PBE or the other hybrid functionals.

Furthermore, we presented new results for the static XC kernel computed with the PBE0, PBE0-1/3, HSE06, HSE03, and B3LYP functionals for both the ground state and partially degenerate case with T ≃ T_F. Among the considered XC functionals, it has been revealed that PBE0 yields the closest agreement with the exact PIMC and DMC reference data. Particularly, at q < 3q_F, the PBE0 static XC kernel is practically exact. For the partially degenerate regime, this is a significant improvement over LDA, PBEsol, AM05, and SCAN static XC kernels.⁴² At q > 3q_F and T ≃ T_F, none of the aforementioned LDA, GGA, meta-GGA, and hybrid XC functionals are able to accurately describe the static LFC of the UEG, G(q) = (q²/4π)K_xc(q). This is not a critical point for the density response function since K_xc(q) ≪ 1 and K_xc(q) do not introduce a significant correction into χ(q) as evident from Eq. (8). Experimentally, large wave numbers are probed in WDM and dense plasmas using the XRTS signal collected at large scattering angles.^75,112,113

The static XC kernel, being a crucial quantity in linear-response time-dependent density functional theory, enables a unique insight into the performance of the various XC functionals. Thus, in combination with the exact quantum Monte Carlo data for the UEG,^{60,76,82,87–89} it also can help with the construction of improved XC functionals for WDM applications, where the XRTS signal from large q values is crucial for the diagnostics of system parameters in experiments.

We are convinced that our present findings will be useful for future DFT calculations of WDM, in general, and for the many applications that utilize the XC kernel,⁴² in particular.

This work was funded by the Center for Advanced Systems Understanding (CASUS), which was financed by Germany’s Federal Ministry of Education and Research (BMBF) and by the Saxon state government out of the State budget approved by the Saxon State Parliament. We gratefully acknowledge computation time at the Norddeutscher Verbund für Hoch-und Höchstleistungsrechnen (HLRN) under Grant No. shp00026, and on the Bull Cluster at the Center for Information Services and High Performance Computing (ZIH) at Technische Universität Dresden. The authors also thank Henrik Schulz and Jens Lasch for providing very helpful support on high performance computing at HZDR.

The authors have no conflicts to disclose.

Zhandos A. Moldabekov: Conceptualization (lead); Data curation (lead); Formal analysis (lead); Investigation (lead); Methodology (equal); Visualization (lead); Writing – original draft (lead). Mani Lokamani: Data curation (equal); Writing – review & editing (equal). Jan Vorberger: Methodology (equal); Writing – review & editing (equal). Attila Cangi: Conceptualization (equal); Methodology (equal); Writing – review & editing (equal). Tobias Dornheim: Conceptualization (equal); Data curation (equal); Funding acquisition (equal); Methodology (equal); Project administration (equal); Writing – review & editing (equal).

The data that support the findings of this study are available from the Rossendorf Data Repository (RODARE).¹¹⁴