The Relevance of Electronic Perturbations in the Warm Dense Electron Gas

Warm dense matter (WDM) has emerged as one of the frontiers of both experimental and theoretical physics and is challenging traditional concepts of plasma, atomic, and condensed-matter physics. While it has become common practice to model correlated electrons in WDM within the framework of Kohn-Sham density functional theory, quantitative benchmarks of exchange-correlation (XC) functionals under WDM conditions are yet incomplete. Here, we present the first assessment of common XC functionals against exact path-integral Monte Carlo calculations of the harmonically perturbed thermal electron gas. This system is directly related to the numerical modeling of X-Ray scattering experiments on warm dense samples. Our assessment yields the parameter space where common XC functionals are applicable. More importantly, we pinpoint where the tested XC functionals fail when perturbations on the electronic structure are imposed. We indicate the lack of XC functionals that take into account the needs of WDM physics in terms of perturbed electronic structures.

Warm dense matter (WDM) has emerged as one of the frontiers of both experimental and theoretical physics and is challenging traditional concepts of plasma, atomic, and condensed-matter physics. While it has become common practice to model correlated electrons in WDM within the framework of Kohn-Sham density functional theory, quantitative benchmarks of exchange-correlation (XC) functionals under WDM conditions are yet incomplete. Here, we present the first assessment of common XC functionals against exact path-integral Monte Carlo calculations of the harmonically perturbed thermal electron gas. This system is directly related to the numerical modeling of X-Ray scattering experiments on warm dense samples. Our assessment yields the parameter space where common XC functionals are applicable. More importantly, we pinpoint where the tested XC functionals fail when perturbations on the electronic structure are imposed. We indicate the lack of XC functionals that take into account the needs of WDM physics in terms of perturbed electronic structures.

I. INTRODUCTION
Understanding transient states in warm dense matter (WDM) is one of the grand challenges of plasma physics that is currently being tackled in a number of experimental facilities [1][2][3][4] . In these experiments, WDM is generated, for example, due to laser-induced shock compression [5][6][7] . At the foundational level, probing WDM facilitates a better understanding of astrophysical objects such as planetary interiors [8][9][10][11][12][13] and stars [14][15][16][17] . Furthermore, understanding WDM has great technological potential as it occurs during the fuel compression processes in inertial confinement fusion 18 . Finally, exploring novel materials properties by driving matter through the WDM regime is an active area of research [19][20][21][22] .
The interpretation of WDM experiments relies on a strong interplay with theory and simulation, because typical parameters like density and temperature cannot be inferred solely from the experimental data. The use of quantum Monte Carlo (QMC) techniques has recently paved the way for determining the exact properties of interacting electrons in the uniform electrons gas under WDM conditions [23][24][25][26][27] . However, this method is computationally expensive and does not take into account the coupling of the electrons to the ions explicitly. Kohn-Sham density functional theory (KS-DFT) 28 has therefore emerged as the standard method to model the electronic structure in WDM, because it includes both electron-electron and electron-ion interactions. Due to its balance of accuracy and computational efficiency 29 , KS-DFT enables calculating materials properties such as structural and electronic transport properties under the conditions relevant to the WDM community [30][31][32][33][34][35][36] .
Here, the central quantity is the KS potential -a meanfield potential that mimics the electron-electron interaction, in principle, exactly. The key ingredient to the KS potential is the exchange-correlation (XC) functional. In a) Electronic mail: a.cangi@hzdr.de practice, the XC functional needs to be approximated. From a historical perspective, the development of XC approximations has focused on the electronic structure of molecules and solids under ambient conditions 37 . These are commonly ranked in terms of increasing accuracy and computational cost on the so-called Jacob's ladder 38 . Starting at the lowest rung with the local density approximation (LDA) 28 , an array of XC approximations has been developed including, for instance, the generalized gradient approximation (GGA) 39,40 , the meta-GGA 41,42 , and hybrid functionals [43][44][45] . A key step most relevant for WDM modeling is the generalization of KS-DFT to finite temperature 46 . Based on this, several works have fleshed out the theoretical aspects of functional construction at finite temperature [47][48][49][50] . Most recently, these have led to the construction of XC functionals that have an explicit temperature dependence [51][52][53][54][55][56] .
Despite these efforts, functional development has not taken into account the needs of WDM physics to a large extent beyond the inclusion of finite-temperature effects in the electronic structure 57,58 . Some exact conditions imposed for chemistry and solid-state physics are not relevant to WDM. For example, in contrast to solid-state physics, the surface energy has no significance in WDM, because the generated samples have no well defined surface, but rather transition to a plasma state 59-61 . Conversely, there are properties that are highly important in WDM, but have little value for solid-state physics. The most prominent example are perturbations imposed on the electronic structure with respect to the wave number q. An accurate description of such perturbations across a large range of q is essential for WDM modeling. Accuracy is not only required in the long wavelength regime q < 2q F which is dominated by collective electronic excitations, but also at large wave numbers q > 2q F where so-called single-particle effects become important [62][63][64][65] . Achieving a high accuracy on a wide range of wave numbers is essential, e.g., for modeling X-Ray Thomson scattering (XRTS) experiments, which is a technique of paramount importance in WDM diagnostics 66 .
In this paper, we therefore analyze the accuracy of common XC functionals in the WDM regime when perturbations on the electronic structure are imposed across a large range of wave numbers q. To that end, we benchmark the results of KS-DFT calculations against pathintegral QMC data which are considered exact within the given error bars. Specifically, we employ four common XC functionals − the LDA in the Perdew-Zunger parametrization 67 , the GGAs PBE 40 and PBEsol 68 , and the meta-GGA SCAN 42 . We focus on the sensitivity of these XC functionals towards a perturbed electronic structure at a finite wave number q. Thereby, we uncover failures of these XC approximations. Our result strongly highlights the need for novel XC functionals that remain accurate under a perturbed electronic structure that is typically present in WDM. Furthermore, to our knowledge, this assessment of KS-DFT against QMC data with respect to electronic density perturbations in the WDM regime has not been performed yet. While recent developments introduced above have addressed the explicit temperature dependence of XC functionals, they do not account for the physics that stems from perturbations at finite q. In our analysis, we therefore separate the influence of finite temperature in the electronic states from the effect of perturbations in q. We follow the common practice of using the listed ground-state XC approximations and including only the implicit temperature dependence of the electronic structure in terms of a Fermi-Dirac occupation of the KS states.
The paper is organized as follows: we introduce the theoretical aspects and the simulation methods in Section II; we present our results on benchmarking common XC approximations in Section III; we address the issue of finite size effects in Section IV; we assess the performance of the aforementioned XC functionals in terms of the total energy in Section V; we provide conclusions and an outlook on future perspectives in Section VI. We provide methodological details on KS-DFT and pathintegral Monte Carlo (PIMC) in the Appendix.

II. THEORY AND SIMULATION METHODS
We begin our assessment with the following Hamiltonian that imposes electronic perturbations on a warm dense uniform electron gas (UEG): whereĤ UEG denotes the standard Hamiltonian of the UEG [69][70][71] , N the number of electrons, N p the number of external harmonic perturbations, A their amplitude, and q i the wave vector of the perturbations. In this work, we consider both one harmonic perturbation N p = 1 and the combination of two harmonic perturbations N p = 2. Note that we work within atomic units throughout.
It has been shown that the electronic states described by Eq. (1) are generated in WDM experiments 72 . Furthermore, Eq. (1) has been used to rigorously examine various fundamental physical properties of the electronic structure in WDM, such as the local field correction (LFC) [73][74][75][76] and the non-linear response 72,77,78 that are used to describe XRTS signals. Moreover, it turns out that the UEG model provides an excellent description of many electronic properties in WDM, because the electron-ion coupling can be relatively weak 77,79,80 .
The Hamiltonian in Eq. (1) is a convenient device to control the degree of inhomogeneity and the wave number of the imposed perturbations on the UEG by tuning the parameters A and q, respectively. We analyze the performance of XC functionals with respect to electronic density perturbations by considering an amplitude range 0.02 ≤ A ≤ 5. This range covers regimes from weak density perturbations with δn/n 0 1 to strongly inhomogeneous electronic systems with δn/n 0 1, where δn = n − n 0 . Additionally, we consider perturbations imposed by tuning the wave number in the range 0.843 q F ≤ q i ≤ 5.9 q F . Thereby, we assess the performance of XC functionals at different length scales ranging from the long wavelength regime defined by collective behavior q < q F 62,63 to the scale defined by the singleparticle limit q q F 24,64 . This covers the entire range of wave numbers relevant to WDM generated in experiments.
WDM conditions prevail when we consider matter at solid density and a temperature T ∼ T F 81,82 , where T F = E F /k B denotes the Fermi temperature defined in terms of the Fermi energy E F and the Boltzmann constant k B . In our analysis, we consider the densities r s = 2 and r s = 6 and a temperature T = T F , where r s = a/a B defines the number density of electrons which is given as the ratio between the Wigner-Seitz radius a and the first Bohr radius a B . At T T F , the parameter r s also characterizes electronic non-ideality 83,84 . Therefore, it is used also as a coupling parameter. These parameters are typically encountered in experiments with laser-driven and shock-compressed solid targets [85][86][87] .
The KS-DFT calculations were performed with the GPAW code [88][89][90][91] , which is a real-space implementation of the projector augmented-wave method. Details of the KS-DFT simulation parameters such as the number of kpoints and number of bands are given in the Appendix.
We provide unassailable data for benchmarking our KS-DFT results by carrying out direct path-integral QMC calculations based on Eq. (1) without any restrictions on the nodal structure of the thermal density matrix. Therefore, the calculations are computationally expensive due to the fermion sign problem 92,93 , but exact within the given Monte Carlo error bars. The corresponding simulation details, such as the number of imaginary-time propagators, can be found in the Appendix VI.

A. Single Harmonic Perturbation
Let us begin our analysis by considering the case of a single harmonic perturbation, i.e., N p = 1 in Eq. (1). First, we investigate the electronic density distribution at a density (coupling) parameter r s = 2 for a range of perturbation amplitudes and wave numbers. Then, we consider a stronger coupling regime at r s = 6. These parameters are of particular interest in WDM experiments.

Metallic Density, rs = 2
The density distribution along the direction of the density perturbation is shown for two amplitudes A = 0.1 and A = 1 and increasing wave number in Fig. 1. The top panel displays the results for q 1 = q min = 0.843 q F , the middle panel for q 1 = 2q min , and the bottom panel for q 1 = 3q min . On the scale of the total density, all considered XC functionals are in overall agreement with each other and with the QMC data.
Let us now consider the actual deviation of the KS-DFT data from the reference QMC data in closer detail. To that end, we consider the relative density deviation ∆n/max{δn} between the KS-DFT data and the reference QMC data, where max{δn} is the maximum deviation of the QMC data from the mean density. We use max{δn} for the analysis of the KS-DFT results because the physically important quantity in the case of a weak perturbation δn/n 0 1 is the deviation of the density from n 0 rather than the total density n itself. Indeed, δn defines the density response of the system which is a cornerstone of linear response theory describing all related physical properties of electrons in equilibrium. More specifically, when the external field defined via Eq. (1) is weak, the density perturbation δn = 2Aχ(q) cos(q · r) defines the static response function χ(q). Therefore, ∆n/max{δn} = ∆χ/χ QMC measures the error in the response function. On the other hand, to keep the present analysis general, we provide the values of max{δn} along with ∆n/max{δn}. This allows a simple conversion of data from ∆n/n 0 or ∆n/n.
In Fig. 2, we show the relative difference between the KS-DFT data and the reference QMC data at r s = 2 for the amplitudes A = 0.02, 0.1, 0.5, 1 (from left to right) and for the wave numbers q 1 = q min , 2q min , 3q min (from top to bottom). Additionally, the corresponding largest absolute value of the discrepancy is given in Table I. Our assessment of the results presented in Fig. 2 and Table I can be summarized as follows. At A = 0.02 and q 1 = q min , the results computed with the LDA, PBE, and PBEsol functionals have about the same level of accuracy. They are consistent with the QMC reference densities to within 7 %. With a maximum deviation of 5.66 % from the QMC data, the SCAN functional yields slightly more accurate results. Note that here the statistical uncertainty of the QMC results is of the same order as the relative difference between the KS-DFT and QMC data (see the top left corner in Fig. 2). Upon increasing the perturbation strength up to A = 1.0, the relative error of the KS-DFT obtained from the LDA and GGA functionals remains less than ≤ 3 %, while the use of SCAN provides a remarkable accuracy better than 1 %. Therefore, the KS-DFT calculations using SCAN are virtually exact at q 1 = 0.843 q F for both weak and strong perturbations. Contrarily, the accuracy provided by the SCAN functional is not maintained with an increase of the perturbation in terms of its wave number to q 1 = 2q min and q 1 = 3q min . Across the perturbation amplitdes 0.02 ≤ A ≤ 0.5 (δn/n 0 ≤ 0.244), the relative differences for calculations using SCAN are about 4 % and 8 %. The LDA and PBEsol functionals provide an accuracy of about 2 % (5 %) at q 1 = 2q min (q 1 = 3q min ).
The PBE functional provides a comparable accuracy to LDA and PBEsol at q 1 = 2q min , but becomes less accurate at q 1 = 3q min and A = 0.5 reaching an error of 7.5 %. In the strong-perturbation regime (A = 1.0 and max{δn} > 1.3 n 0 ), all functionals LDA, PBE, PBEsol TABLE I. The performance of common XC functionals in terms of the relative density deviation ∆n/max{δn}100 %. A single harmonic perturbation at a fixed density rs = 2 (metallic density) and varying perturbation amplitude 0.02 ≤ A ≤ 1 and wave number qmin ≤ q ≤ 3qmin is considered, where qmin = 0.843 qF . The largest absolute values of the deviation are listed in this table. functionals now provide about the same level of accuracy when compared to the QMC data. In summary, overall the LDA and PBEsol show a more robust performance compared to SCAN and PBE for all perturbation amplitudes and relatively small wave numbers at the typical mass density of metals.
Next in Fig. 3, we assess the performance of the XC functionals in the limit of large perturbation wave numbers by setting q 1 = 7q min = 5.9 q F , but keeping the mass density metallic (r s = 2). We consider the per-turbation amplitudes A = 1 (with max{δn} = 0.28 n 0 ) and A = 5 (with max{δn} = 1.268 n 0 ). At A = 1, we find that for LDA, PBEsol, and SCAN, the relative difference between the KS-DFT data and the QMC data is less than 2.5 %, while the PBE functional yields an error of 5 %. When we increase the perturbation amplitude to A = 5, the KS-DFT results from all considered functionals are in agreement with each other, but exhibit a strong disagreement with the QMC data of up to about 12 %. In this extreme regime, the deviation of the KS- DFT data from the QMC results is significant even on the scale of the total density n. This is illustrated in Fig. 4. In summary, all tested XC functionals fail dramatically in yielding an accurate electronic density for q q F and max{δn}/n 0 > 1 at metallic density. We argue that a large number of other XC functionals that are derived from LDA, GGA, and meta-GGA classes are afflicted with the same limitation.  Next, we investigate strongly correlated electronic systems where r s = 6. Such low densities can be realized, for example, in evaporation experiments 94 . From a theoretical perspective, these conditions are particularly challenging due to the substantial impact of electronic XC effects 95,96 on physical observables like the electrical conductivity or the density.
Again, we consider weak as well as strong perturbation amplitudes corresponding to a density enhancement in the range 0.034 n 0 ≤ δn ≤ 2.044 n 0 . The electronic density distribution along the direction of the perturbation is shown in Fig. 5 for q 1 = q min , 2q min , and q 1 = 3q min (from top to bottom) at A = 0.05 (left) and A = 1.0 (right). The agreement between the KS-DFT data and QMC data is excellent when q 1 = q min . With increasing wave number, the accuracy decreases. For example, the errors are larger than in the previous case where a metallic density is considered (see Fig. 1).
We further delineate the relative difference between the KS-DFT data and the reference QMC results in Fig. 6. There, the relative difference is illustrated for the ampli- tudes A = 0.01, 0.05 and A = 0.1 (from left to right) and the wave number q 1 = q min , 2q min , and q 1 = 3q min (top to bottom) at r s = 6. The largest absolute values of the differences are listed in Table II. The main conclusions of our assessment in this parameters range are as follows. At q 1 = q min and all perturbation amplitudes, the SCAN functional yields the most reliable results with an accuracy better than 1.64 %. The other tested functionals provide an accuracy better than 3.13 %. In contrast to that, when we increase the wave number to q 1 = 2q min , SCAN performs much worse with a maximum deviation of 12.58 % at A = 0.01 and 9.97 % at A = 0.05. Both the LDA and PBEsol functional yield an accuracy better than 3.5 %. Similarly, PBE yields an error of 4.02 % and 5.08 % at A = 0.05 and A = 0.1, respectively. However, irrespective of the wave number, the SCAN functional provides the most accurate results with the maximum deviation of 2.63 % at a stronger perturbation amplitude A = 0.1. Finally, increasing the wave number of the perturbation to q 1 = 3q min renders KS-DFT data less accurate with a maximal deviation in the range of about 8 % and 12 % (for 0.01 ≤ A ≤ 0.1) when LDA and PBE are used. The PBEsol functional is in better agreement with the QMC results showing a maximum deviation in the range between 5.63 % and 9.45 %. The SCAN functional provides very low accuracy when the perturbation is weak(A = 0.01), with a deviation of almost 30 %. This improves in the case of strong perturbation A = 0.1, where the maximum deviation is 6.2 %. In summary, the overall performance of the considered XC functionals is worse when strongly coupled electronic systems are considered. We observe failures of the functionals, in particular, when the wave number of the perturbation increases.

B. Double Harmonic Perturbation
Furthermore, we assess the accuracy of the considered XC functionals when more complex perturbations are applied. To that end, we consider a double harmonic perturbation N p = 2 with q 1 = q min and q 2 = 2q min . This allows us to check if the observed poor performance of the XC functionals manifests itself when the perturbations of different wave numbers are superimposed. Again, we consider both metallic densities (r s = 2) and strongly coupled systems (r s = 6). When r s = 2, we set the perturbation amplitude to A = 0.1 resulting in a density perturbation of the order of 0.1 n 0 . When r s = 6, we set A = 0.01 leading to a similar density perturbation of the order of 0.1 n 0 .
The resulting electronic density distributions are shown in Fig. 7. At r s = 2, the KS-DFT results show good agreement with the QMC data when gauged on the scale of the total density. At r s = 6, the SCAN functional deviates from the QMC data significantly, while the results obtained from the other XC functionals are indistinguishable from the QMC data on this scale. We take a closer look at the performance of the various XC functionals in Fig. 8, where we use the relative deviation in the density ∆n/max{δn} between the KS-DFT data and the reference QMC results. The corresponding largest absolute values of the differences are given in Table III. The results in Fig. 8 and Table III show that LDA, PBE, and PBEsol provide an accuracy better than 2.3 % for both r s = 2 and r s = 6. Contrarily, the densities computed using SCAN have a maximum deviation of 4.31 % and 12.05 % at r s = 2 and r s = 6, respectively. These numbers are similar to the deviations observed in the case of a single harmonic perturbation at the wave number 2q min . This provides a strong indication about the general applicability of the present findings to other systems, as any external potential can be expressed as a superposition of harmonic perturbations in reciprocal space. Relative deviation in the density ∆n/max{δn}100 % between the KS-DFT data and the reference QMC data for q = qmin (top) and q = 3qmin (bottom) at rs = 2 and A = 0.1.

IV. FINITE SIZE EFFECTS
The results presented so far were computed for 14 electrons in a simulation cell. Naturally, it is important to inquire whether the large deviations of the KS-DFT results from the QMC data observed at 2q min and 3q min Relative deviation in the density ∆n/max{δn}100 % between the KS-DFT data and the reference QMC data for q = qmin (top) and q = 3qmin (bottom) at rs = 2 and A = 0.1.
compared to the case with q 1 = q min are indeed genuine results or an artifact due to finite size effects. To exclude this possibility, we consider N = 20 and N = 34 electrons in the simulation cell at r s = 2 and A = 0.1. The density distribution along the direction of the density perturbation for the case with 20 electrons is shown in Fig. 9. The top panel displays the results for q 1 = q min = 0.843 q F and the bottom panel for q 1 = 3q min . We observe that on the scale of the total density, the KS-DFT results are in agreement with the QMC data. The corresponding relative difference between the KS-DFT results and the reference QMC data is displayed in Fig. 10, where the top panel is for q 1 = q min and the bottom panel is for 3q min . Fig. 10 shows that at q 1 = q min the relative difference between KS-DFT results obtained using the SCAN functional and the QMC data is within the given uncertainty range, whereas it increases to up to 8% with an increase in the wave number of the perturbation to 3q min . The other considered XC functionals lead to slightly worse agreement at q 1 = q min , with about 2% maximum relative difference, and to much better results at q 1 = 3q min , with a maximum relative difference of about 4%.
Next, we further increase the number of electrons in the simulation cell to 34. The density distribution for 34 particles along the direction of the density perturbation is shown in Fig. 11 at q 1 = q min (top panel) and at q 1 = 3q min (bottom panel). Again, on the scale of the total density we observe agreement of KS-DFT results with the QMC data. The corresponding relative differences, shown in Fig. 12, confirm this assessment. At q 1 = q min (top panel) the KS-DFT results agree with the QMC data within the given numerical uncertainty. At q 1 = 3q min (the bottom panel in Fig. 12), KS-DFT results obtained with the SCAN functional exhibit a relative deviation from QMC data of about 8%. The other XC functionals yield a maximum relative difference of about 4%. Note the statistical uncertainty of the QMC data at N = 34 is noticeably larger due to the infamous fermion sign problem. See Appendix for more details and Ref. 92 for a topical review of this issue.
This analysis of the densities and relative differences between KS-DFT results and QMC data for N = 20 and N = 34 electrons clearly shows that our conclusions at N = 14 electrons are not affected by finite size effects.

V. ENERGY PERTURBATION
Finally, we analyze the correlation of the revealed inaccuracies in the densities with the total energy. To that end, we focus on the change in the total energy compared to the unperturbed case. This is quantified by the reduced perturbation in the total energy where E(A) is the total energy of electrons in the presence of the external perturbation with the amplitude A, E(A = 0) is the total energy of the unperturbed system, and E QMC (A = 0) is the total energy of the unperturbed system from QMC simulations. The results for δ E at r s = 2 and at r s = 6 for different amplitudes and wave numbers are presented in Fig. 13

FIG. 13.
Relative deviation in the total energy from UEG (A = 0) due to external perturbation with amplitude A (see Hamiltonian (1)) at rs = 2. From top to bottom: q1 = qmin, 2qmin, and q1 = 3qmin. and Fig. 14, respectively. In these figures we also plot δ E ∼ A 2 using the data point at A = 0.1 (A = 0.01) for r s = 2 (r s = 6). This dependence is known to be exact in the limit of small A in the ground state at T T F and is referred to as the stiffness theorem 97 . At finite temperature, the stiffness theorem applies to the total free energy. Figs. 13 and 14 can be understood as an empirical indication of validity of stiffens theorem for the total energy at finite temperature. The deviation of the simulation data from the A 2 dependence is related to the emergence of the higher-order non-linear effects.
In Fig. 13, in addition to KS-DFT results within LDA, PBE, PBEsol, and SCAN, we also show results computed using null XC functional, i.e., the electron-electron interaction approximated only using the mean Hartree field. First of all, in Fig. 13 we clearly see that neglecting the XC energy leads to less accurate results compared to the QMC data at amplitudes A = 0.1, A = 0.5 and A = 1.0. At A = 0.02, the change in the energy due to the exter- nal field is comparable with statistical uncertainty. This prevents us from quantifying the accuracy of the KS-DFT results for a reduced perturbation in the energy at that perturbation strength. At A = 0.1, all KS-DFT data points, except the null XC functional case, are in agreement with the QMC data within the given statistical uncertainty. Furthermore, in Table IV, we provide the relative difference of the reduced perturbation in the total energy, ∆ E, which is defined by From Table IV we infer that in the case of null XC functional, ∆ E is about 15% (9%) and 13% (7%) for A = 0.5 and A = 1.0 at q = q min (q = 2q min ). This large difference mainly stems from an inaccurate density due to neglecting the XC functional (see Appendix). This clearly demonstrates the importance of including the XC IV. The performance of common XC functionals in terms of the relative deviation in the total energy as defined by Eq. (3). A single harmonic perturbation at a fixed density rs = 2 and perturbation amplitudes A = 0.5 and A = 1.0 and wave numbers qmin ≤ q ≤ 3qmin is considered, where qmin = 0.843 qF . The presented data has about 1% uncertainty. functional in the considered KS-DFT calculations. Also, as shown in Table IV, the LDA, PBE, and PBEsol functionals show similar performance. The fact that SCAN is more accurate for ∆ E at q 1 = q min and less accurate at q 1 = 2q min correlates with the previous assessment of the densities.
In the case of r s = 6 ( Fig. 14), the dependence of δ E on A is similar to that at r s = 2. At A = 0.01, the KS-DFT results computed using different XC functionals is in agreement with the QMC data within the given statistical uncertainty. At A = 0.05 and A = 0.1, the SCAN results tend to overestimate δ E with increasing wave number from q 1 = q min to q 1 = 2q min and q 1 = 3q min . The LDA, PBE, and PBEsol results are very similar across all considered A and q 1 values.
In Table V we list ∆ E as defined by Eq. (3) in order to better delineate the differences. Similar to r s = 2, also at r s = 6 the SCAN results are most accurate at q 1 = q min . Contrarily, at q 1 = 2q min and q 1 = 3q min , the largest deviation from the QMC data is exhibited by the SCAN results. This reflects the fact that there are more inaccuracies in the density when SCAN is used over the other XC functionals. The KS-DFT results for δ E computed using LDA, PBE, and PBEsol have similar values at q 1 = q min and q 1 = 2q min . At q 1 = 3q min , the PBE results are less accurate by a few percent than the LDA and PBEsol results. Overall, by comparing data in Table II and in Table V, we observe that the variations in ∆ E reflect those in ∆n/max{δn}.

VI. CONCLUSIONS AND OUTLOOK
We benchmarked the performance of KS-DFT based on the LDA, PBE, PBEsol, and SCAN XC functionals against exact QMC data in the WDM regime. Our assessment revealed a set of conditions for the successful simulation of WDM with KS-DFT at QMC level accuracy. Our comparative analysis unambiguously demonstrates when the considered XC functionals fail to correctly describe the electronic density.
We found that the KS-DFT results are sufficiently accurate for small wave numbers of the density perturbation, q < q F . In particular, using the SCAN functional yields an excellent agreement with the QMC reference data. However, with increasing wave number q > q F (as tested for q = 1.686q F and q = 2.529q F ), the SCAN functional performs much worse than LDA, PBE, and PBEsol. This is somewhat surprising, because the SCAN functional is ostensibly designed to be superior over LDA and GGA functionals. In contrast to that, LDA and PBEsol show a robust performance with an accuracy better than 4 % at q ≤ 1.686 q F for both r s = 2 and r s = 6.
As a key finding of our assessment, we highlight that the overall performance of the considered XC functionals deteriorates upon increasing the wave number of the density perturbation. At r s = 2 and q = 2.529 q F , LDA and PBEsol still provide an accuracy with an error of less than 6 % in δn = n − n 0 . At r s = 6, the same is valid for PBEsol in the regime of weak perturbations, δn/n 0 1. Other XC functionals essentially fail at r s = 6, δn/n 0 < 1, and q = 2.529 q F . At the largest considered wave number, q = 5.9 q F , all considered XC functionals yield errors of less than 6 % if δn/n 0 1 for metallic densities (r s = 2.0). Finally, we highlight the failure of the considered XC functionals in the regime of strong perturbations, δn/n 0 > 1, and large wave-numbers, q = 5.9 q F , where they exhibit a maximum deviation of about 10 %. We also showed that the reported errors in the density correlate with the errors in the total energy.
Based on the performed analysis we can formulate the following general recommendations for using XC functionals within the typical WDM regime at temperatures T ≤ T F and when the total density is of interest: When characteristic wave numbers q < q F are of interest, the SCAN functional is the most reliable choice, it provides accuracy at the level of QMC. For a wider range of wave numbers q ≤ 5.9 q F (at r s = 2) and q ≤ 1.686 q F (at r s = 6), the LDA and PBEsol functionals should be used if δn/n 0 < 1, because they provide consistent results with a relative error not exceeding a few percent.
One can understand the performance of the consid-ered XC functionals by recalling that, in the case of a weak perturbation, the density deviation from the mean value is defined by the static response function χ(q) since δn(r) = 2Aχ(q) cos q · r. Within linear response theory of the UEG, XC effects are included in terms of the local field correction G(q) via the relation χ(q) −1 = χ 0 (q) −1 − 4π q 2 (1 − G(q)), where χ 0 (q) is the response function of ideal electron gas, i.e., the non-interacting UEG. In this case, the XC functional is K xc (q) = − 4π q 2 G(q). In the UEG, the LDA functional reduces to the long wavelength limit q 2q F of K xc . Clearly, this behavior can be attributed to the failure of the LDA at perturbation wave numbers q > 2q F . The same is valid for PBE and PBEsol, where gradient corrections to exchange and correlation exactly cancel each other in the limit of the UEG. On the other hand, when the density perturbation is strong, the characteristic wave number of the density perturbation defined by the local gradient of the density is larger than the wave number of the external harmonic perturbation. In this regime, regions with a strong density localization are formed. These lead to an increase of the kinetic energy and reduce the importance the wave-number dependence in XC functionals. This may explain the better performance of the considered XC functionals when the perturbation amplitude is increased. Additionally, we note that the SCAN functional is designed to yield PBE-like results in the limit of the UEG. However, it is not trivial to explain the bad performance of SCAN compared to PBE in the case of a weak perturbation and q > 2q F due to the large number of exact conditions enforced. A more in-depth analysis is needed to explain this behavior.
The presented data along with our assessment constitute an indispensable guide on the choice of the XC functional for KS-DFT calculations when an inhomogeneous electronic structure of WDM is investigated. This is of paramount importance for the diagnostics of XRTS experiments. We highlight the importance of this application by pointing out that KS-DFT results are used to extract electronic parameters like temperature and density from experimental observations. Besides that, our findings advance our understanding on how well KS-DFT is capable of capturing the electronic structure under WDM conditions. They also point to the parameter space where XC functionals ought to be improved for their use in the WDM application domain. We highlight the need for XC functionals that are accurate when perturbed electronic states are present. This goes beyond the inclusion of explicit temperature effects in the XC free energy. This particular outcome of our assessment is valuable for DFT developers.
In the light of the vast amount of available XC functionals, a comprehensive assessment was beyond the scope of this work. An extensive comparison of available XC functionals shall be presented elsewhere. Instead, in this work, we focus on a representative set of XC functionals − the basic LDA and its common generalizations. Furthermore, we have set up a numerical workflow for benchmarking XC functionals under perturbed electronic structures in WDM. The presented KS-DFT data, input scripts, and QMC data will be made accessible online 98 . In doing so, we provide tools to test any existing or newly developed XC functional under perturbed electronic states in WDM. We believe that our work is a valuable contribution that facilitates both the rigorous assessment and construction of XC functionals adapted to the needs of WDM.

ACKNOWLEDGMENTS
We acknowledge helpful feedback from M. Bussmann. ZM gratefully acknowledges stimulating discussions with Timothy Callow. This work was funded by the Center for Advanced Systems Understanding (CASUS) which is financed by the German Federal Ministry of Education and Research (BMBF) and by the Saxon Ministry for Science, Art, and Tourism (SMWK) with tax funds on the basis of the 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 shp00026, and on the Bull Cluster at the Center for Information Services and High Performance Computing (ZIH) at Technische Universität Dresden.

APPENDIX: KS-DFT SIMULATION PARAMETERS
All KS-DFT calculations were performed with the GPAW code [88][89][90][91] . A k -point grid of N k × N k × N k with N k = 12 at r s = 2 and N k = 8 at r s = 6 with a Monkhorst-Pack sampling of the Brillouin zone (kpoints) was used. At T = T F , 180 orbitals (with the smallest occupation number of about 10 −4 ) were used for a total of N = 14 electrons. The grid spacing was set to 0.15Å for 0.02 ≤ A ≤ 1 and r s = 2, 0.05Å for A = 5 and r s = 2, 0.25Å for r s = 6 and 0.05 ≤ A ≤ 0.1. At A = 5, 240 bands were used. Convergence criteria used for the self-consistency cycle: the energy change (last 3 iterations) must be less than 0.5 meV per valence electron, the change in integrated absolute value of density change must be less than 0.0001 electrons per valence electron, and the integrated value of the square of the residuals of the Kohn-Sham equations should be less than 4 × 10 −8 eV 2 per valence electron (see GPAW documentation).
In Fig. 15 we compare results computed using 180 bands and 280 bands using LDA and SCAN. From Fig.  15 we see that there is no notable difference in the results for densities.
In Fig. 16, we present results for the total density computed using null XC functional in comparison with the data obtained using QMC simulation and KS-DFT simulation with LDA XC functional. From Fig. 16 we see that it is essential to use XC functional to correctly compute the total density using KS-DFT method.

APPENDIX: PIMC SIMULATION DETAILS
The basic equation behind the PIMC method 99-101 is the canonical partition function × dR R| e −βĤ |π σ ↑π σ ↓ R , where R = (r 1 , . . . , r N ) T contains the coordinates of all N particles, andπ σ ↑ (π σ ↓ ) is the permutation operator corresponding to a particular element σ ↑ (σ ↓ ) from the permutation group S N . In addition, sgn(σ ↑ , σ ↓ ) is the sign function, which is positive (negative) for an even (odd) number of pair exchanges 102 . Unfortunately, the matrix elements of the density operatorρ = e −βĤ cannot be evaluated in a straightforward way as the kinetic (K) and potential (V ) contributions to the full Hamiltonian do not commute, R| e −βĤ |R = R| e −βK e −βV |R .
FIG. 16. The electronic density distribution along the perturbation direction for two different amplitudes A at wave number qmin for rs = 2 and T = TF . The comparison of the data computed using LDA XC functional and null XC functional with QMC data shows the importance of the XC functional for the calculation of the electronic density.
As a practical workaround, one uses the exact semi-group property ofρ, which allows one to re-write Eq. (4) as where we have defined = β/P . The comparison between Eqs. (4) and (6) reveals that the partition function has been transformed into a high-dimensional integration over P density matrices that have to be evaluated at P times the original temperature T . For a sufficiently large P , one can introduce a suitable hightemperature approximation like the simple primitive factorization e − Ĥ ≈ e − V e − K . In fact, the associated factorization error decays as P −2103,104 , and the convergence in the limit of P → ∞ is ensured by the well-known Trotter formula 105 . For completeness, we mention that higher-order factorizations ofρ are frequently employed in PIMC simulations 104,106,107 , although we do not find them necessary for the present study. In practice, we use P = 200 primitive high-temperature factors such that any factorization error is substantially below the given level of statistical uncertainty; see the Supplementary Material of Ref. 108 for corresponding convergence plots. The high dimensionality of Eq. (6) requires a stochastic evaluation, which can be done efficiently using methods that are based on the celebrated Metropolis algorithm 109 . Specifically, we employ a canonical adaption 110 of the worm algorithm by Boninsegni and co-workers 111,112 . Additional care has to be taken due to the antisymmetric nature of the fermionic density matrix, which can result in positive and negative contributions to Eq. (6). The corresponding cancellation of positive and negative terms is the origin of the notorious fermion sign problem 92,113 , which leads to an exponential increase in compute time with increasing the system size N or decreasing the temperature T .
In practice, this bottleneck is often avoided by imposing the fixed-node approximation 114 . While this allows one to formally remove the sign problem, this advantage comes at the cost of an uncontrolled approximation 115 . Obviously, this is a disadvantage that severely limits the value of such data to benchmark other approximations like thermal DFT. Therefore, we do not impose any nodal restrictions in our simulations, and our PIMC results are exact within the given level of uncertainty. This is possible for the present parameters by using modern supercomputer clusters, and we have spend O(10 5 ) CPUh for the most costly data points with N = 34 at r s = 2 and θ = 1.

DATA AVAILABILITY
The data supporting the findings of this study are available on the Rossendorf Data Repository (RO-DARE) 98 .