We assess the accuracy of common hybrid exchange-correlation (XC) functionals (PBE0, PBE0-1/3, HSE06, HSE03, and B3LYP) within the Kohn–Sham density functional theory for the harmonically perturbed electron gas at parameters relevant for the challenging conditions of the warm dense matter. Generated by laser-induced compression and heating in the laboratory, the warm dense matter is a state of matter that also occurs in white dwarfs and planetary interiors. We consider both weak and strong degrees of density inhomogeneity induced by the external field at various wavenumbers. We perform an error analysis by comparing with the exact quantum Monte Carlo results. In the case of a weak perturbation, we report the static linear density response function and the static XC kernel at a metallic density for both the degenerate ground-state limit and for partial degeneracy at the electronic Fermi temperature. Overall, we observe an improvement in the density response when the PBE0, PBE0-1/3, HSE06, and HSE03 functionals are used, compared with the previously reported results for the PBE, PBEsol, local-density approximation, and AM05 functionals; B3LYP, on the other hand, does not perform well for the considered system. Additionally, the PBE0, PBE0-1/3, HSE06, and HSE03 functionals are more accurate for the density response properties than SCAN in the regime of partial degeneracy.

## I. INTRODUCTION

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^{11} at large research facilities like the European X-ray Free-Electron Laser (XFEL),^{12} the National Ignition Facility (NIF),^{13,14} the Linac coherent light source at SLAC,^{15} or the Sandia Z-machine.^{16} 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^{18} 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.^{25}

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,^{26} lattice constants,^{27} and atomization energies.^{28}

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,^{41} XC kernels for linear-response time-dependent density functional theory calculations,^{42} 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),^{55} PBE,^{56} PBEsol,^{57} AM05,^{58} and SCAN^{59} 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.^{42} Recently, Moldabekov *et al.*^{42} 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.*^{32} 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^{63} 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.^{66} Moreover, they provide a superior description of various material properties like lattice constants^{27} and atomization energies.^{28}

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,^{70} PBE0,^{71} PBE0-1/3,^{72} and B3LYP^{73} 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.*^{42} 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.

## II. THEORETICAL BACKGROUNDS

### A. Hybrid XC functionals for systems at finite temperature

*adiabatic connection*formula.

^{28}At finite temperature, the thermal adiabatic connection formula for the exchange-correlation functional, in terms of the potential energy of exchange-correlation, $Uxc\lambda $,

^{78}is given by

^{31}

^{,}

*λ*is the electron–electron coupling parameter (interaction strength). In the

*λ*= 0 limit, $Uxc\lambda =0$ reduces to the finite temperature Hartree–Fock exchange, with $Uxc\lambda =0=FXHF$ given by

^{35,79}

*f*

_{iσ}being the orbital occupation number defined by the Fermi–Dirac distribution.

*λ*= 1 limit $Uxc\lambda =1$ is approximated by a standard XC density functional, $ExcDF[n]$, such as LDA or GGA

^{28,65}(where DF is an abbreviation for density functional). In their essence, hybrid XC functionals are based on the approximation of the $Uxc\lambda $ in the form of a certain (usually polynomial) interpolation between the limits $Uxc\lambda =0=FxHF[\rho \sigma ]$ and $Uxc\lambda =1=ExcDF[n]$.

^{28,65,71–73}For example, employing

^{65}

^{,}

*n*= 4 and $ExDF$ being the exchange part of $ExcDF$ and implementing the PBE as the DF, one arrives at a finite temperature version of the PBE0 hybrid XC functional,

^{71}

In Eq. (5), the implicit dependence on the electronic temperature is taken into account in $FxHF[\rho \sigma ]$ 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^{80} and the reflectivity and dc electrical conductivity of liquid ammonia^{81} 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[\rho \sigma ]$ is substituted by a screened potential 1/*r*_{12} → erfc(Ω*r*_{12})/*r*_{12}^{79} 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^{−1} is used for the HSE03.^{69}

### B. Density response functions and hybrid XC kernel

^{62,82,83}

^{3,30,84}and

*N*is the number of electrons. Note that we restrict ourselves to the unpolarized case with an equal number of spin-up and spin-down electrons,

*N*

^{↑}=

*N*

^{↓}=

*N*/2, throughout this work. In addition to $H\u0302UEG$, the Hamiltonian in Eq. (6) has an extra harmonic perturbation defined by the amplitude

*A*and wavelength

*λ*= 2

*π*/|

**q**|.

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 $\delta nmax/n0=max\delta n(r)/n0$, where *δn*(**r**) = *n*(**r**) − *n*_{0} is the density deviation from the UEG with density *n*_{0}.

*δn*

_{max}/

*n*

_{0}≳ 1 and sufficiently small perturbation wavelength, we access the regime where standard XC functionals fail due to the strong degree of density inhomogeneity and self-interaction.

^{53}In the limit

*δn*

_{max}/

*n*

_{0}≪ 1, the static density response function is defined only by the properties of the UEG and is independent of

*A*. In this case, for the static density response function,

*χ*(

**q**), we have

*χ*(

**q**,

*ω*) can be represented in terms of the ideal response function

*χ*

_{0}(

**q**,

*ω*) and the XC kernel,

*χ*

_{0}(

*q*,

*ω*) is the temperature-dependent Lindhard function.

*χ*(

**q**) =

*χ*(

**q**,

*ω*= 0), and

*χ*

_{0}(

**q**) =

*χ*

_{0}(

**q**,

*ω*= 0),

*χ*

_{RPA}(

**q**) is the density response function in the random phase approximation (RPA) defined by setting

*K*

_{xc}(

**q**) ≡ 0. In the static case,

*ω*= 0, from Eq. (8) we get

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.

^{30,85}The relation between the static XC kernel and the static LFC reads

^{86}

*n*

_{0}=

*N*/

*V*denoting the average number density of the electrons in the UEG or the mean value of the conduction electron density in metals.

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}

## III. SIMULATION SETUP

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.^{92} 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*_{0} ≪ 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.*^{82} 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*_{0} ≲ 10^{−2} (at *A* = 0.01) to *δn*_{max}/*n*_{0} ≳ 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 4*q*_{min} and up to 5*q*_{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*^{3} × *n*_{0} = 14 electrons, where at *r*_{s} = 2, we have *n*_{0} ≃ 0.02984 Bohr^{−3}. 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^{−7} 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^{−4}, which is sufficiently small to acquire convergence.^{53} 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* = 2*q*_{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}, 2*q*_{min}, and 3*q*_{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* = 4*q*_{min} at *A* = 0.1, 0.5, and 1.0.

Since a detailed introduction to the direct PIMC method^{103} and the associated fermion sign problem^{104} 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.

## IV. RESULTS

We start with the weakly perturbed electron gas where |*δn*(**r**)|/*n*_{0} ≪ 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*_{0} ≳ 1 and provide the analysis of the accuracy of the electron density calculations based on the considered hybrid XC functionals.

### A. Weakly perturbed electron gas

#### 1. Ground state results

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* < 2*q*_{F} was one of the main criteria for an adequate description of metals.^{56} 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)^{82} 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^{60} and the DMC data^{82} 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^{42}). 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.^{27} Here, we further substantiate it for finite perturbation wavenumbers in the range 0.6 ≤ *q* ≤ 3*q*_{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.5*q*_{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.^{56}

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* ≳ 2*q*_{F}. Nevertheless, at 2*q*_{F} ≤ *q* ≤ 3*q*_{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 2*q*_{F} and a monotonically decreasing trend with increasing *q* at *q* > 2*q*_{F}. A similar behavior—with a somewhat worse agreement with the DMC data—was reported by Armiento and Mattsson for the AM05^{58} 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.^{108} 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}.^{109} 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* > 2*q*_{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.5*q*_{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* ≃ 2*q*_{F} compared with HSE06, HSE06, and PBE. At *q* ≳ 2*q*_{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* > 2*q*_{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* < 2*q*_{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.

#### 2. Results at a high temperature *T* ≃ *T*_{F}

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.84*q*_{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.5*q*_{F}. The PBE0, PBE0-1/3, HSE06, and HSE03 results are in close agreement with the reference data (ML-PIMC) at *q* ≲ 2.5*q*_{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.5*q*_{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.5*q*_{F}. A similar behavior has been reported for the AM05 functional at *T* ≃ *T*_{F}.^{42} 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.

### B. Strongly inhomogeneous electron gas

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* = 4*q*_{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.1*n*_{0}. We see that at the considered amplitudes of the perturbation, we always have |*δn*_{⊕}|≥ 0.1*n*_{0} and |*δn*_{⊖}|≥ 0.1*n*_{0}. These density perturbation magnitudes are well beyond the linear response regime considered in Sec. IV A. At *A* = 0.5, we have |*δn*_{⊕}|∼ *n*_{0} and |*δn*_{⊖}|∼ 0.5*n*_{0}. Furthermore, at *A* = 1.0, the electron density values almost vanish in the density depletion regions with *δn*_{⊖} → −*n*_{0} at *q* = *q*_{min} and *q* = 2*q*_{min}, where a physical lower bound *δn*_{⊖} > − *n*_{0} 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*_{⊕} ≈ 2*n*_{0}.

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* = 2*q*_{min} ≈ 1.5*q*_{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* ≃ 2*r*_{s}, i.e., at *q* ≃ 2*π*/*d* ≃ 2*q*_{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.*,.^{110} 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.

*δn*

_{QMC}} is the maximum deviation of the QMC data from

*n*

_{0}. From a physical perspective, Eq. (13) defines an error in the density response computed using KS-DFT in comparison with the exact density response from the PIMC calculations. It becomes obvious by considering a weak perturbation where

*δn*(

*r*) = 2

*A*cos(

*qr*)

*χ*(

*q*) and Eq. (13) reduces to the error evaluation of the linear density response function.

^{111}

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*_{0} and the density error $\Delta n\u0303$ profiles when the harmonic perturbation with *A* = 0.1 and *q* = 3*q*_{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 $\Delta n\u0303$ in Fig. 4(b), where $|\Delta n\u0303|\u22721%$ 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 $\Delta n\u0303$. 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 $\Delta n\u0303$ in the density accumulation region, $\Delta n\u0303\u2295$, and the largest magnitude of the $\Delta n\u0303$ in the density depletion region, $\Delta n\u0303\u2296$, 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 $\Delta n\u0303\u2295$ and $\Delta n\u0303\u2296$ at *q* = *q*_{min} ≃ 0.84 *q*_{F} and different values of *A*. We observe that the B3LYP data have the smallest error with $\Delta n\u0303\u2295<1%$ and $\Delta n\u0303\u2296<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* = 2*q*_{min}, 3*q*_{min}, and 4*q*_{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.1*n*_{0} and |*δn*_{⊖}|∼ 0.1*n*_{0} (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*_{0} ≪ 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 $\Delta n\u0303\u2295\u223c\u22122%$ and $\Delta n\u0303\u2296\u223c1%$. As shown in Fig. 6, at *q* = 2*q*_{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* = 3*q*_{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 $\Delta n\u0303\u2295$ and $\Delta n\u0303\u2296$. In contrast, the PBE0, PBE0-1/3, HSE06, and HSE03 results are highly accurate with deviations of $|\Delta n\u0303\u2295|\u22641%$ and $|\Delta n\u0303\u2296|\u22641%$. From Fig. 7, one can see that at *q* = 3*q*_{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* = 3*q*_{min}, all considered XC functionals show similar accuracy with $|\Delta n\u0303\u2295|\u22722%$ and $|\Delta n\u0303\u2296|\u22722%$.

Finally, in Fig. 8, we show the results for *q* = 4*q*_{min} ≃ 3.371 *q*_{F}. Similar to the case with *q* = 3*q*_{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 $|\Delta n\u0303\u2295|$ and/or $|\Delta n\u0303\u2296|$. At last, for *A* = 1.0 and *q* = 4*q*_{min}, we observe similar accuracy for all considered XC functionals with $|\Delta n\u0303\u2295|\u22724%$ and $|\Delta n\u0303\u2296|\u22724%$.

## V. CONCLUSIONS AND OUTLOOK

Before drawing general conclusions, we recall that in previous work, Moldabekov *et al.*^{53} 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* ≲ 3*q*_{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* < 3*q*_{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.^{42} At *q* > 3*q*_{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*^{2}/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,^{42} in particular.

## ACKNOWLEDGMENTS

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.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**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).

## DATA AVAILABILITY

The data that support the findings of this study are available from the Rossendorf Data Repository (RODARE).^{114}

## REFERENCES

*Frontiers and Challenges in Warm Dense Matter*

*Ab initio*simulation of warm dense matter

*Quantum Theory of the Electron Liquid, Masters Series in Physics and Astronomy*

*Ab initio*exchange–correlation free energy of the uniform electron gas at warm dense matter conditions

*Ab initio*path integral Monte Carlo simulations

*ab initio*path integral Monte Carlo simulations

*ab initio*path integral Monte Carlo study and machine learning representation

*Ab initio*path integral Monte Carlo approach to the momentum distribution of the uniform electron gas at finite temperature without fixed nodes

*Ab initio*path integral Monte Carlo simulation of the uniform electron gas in the high energy density regime