{'fraction_non_alphanumeric': 0.03956642776757685, 'fraction_numerical': 0.02330746905170882, 'mean_word_length': 4.48789208516171, 'pattern_counts': {'":': 0, '<': 13, '<?xml version=': 0, '>': 11, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Benchmarking Exchange-Correlation Functionals in the Spin-Polarized Inhomogeneous Electron Gas under Warm Dense Conditions

Warm dense matter is a highly active research area both at the frontier and interface of material science and plasma physics. We assess the performance of commonly used exchange-correlation (XC) approximation (LDA, PBE, PBEsol, and AM05) in the spin-polarized inhomogeneous electron gas under warm dense conditions based on exact path-integral quantum Monte-Carlo calculations. This extends our recent analysis on the relevance of inhomogeneities in the spin-unpolarized warm dense electron gas [Z.~Moldabekov et al., J. Chem. Phys. 155, 124116 (2021)]. We demonstrate that the predictive accuracy of these XC functionals deteriorates with (1) a decrease in density (corresponding to an increase in the inter-electronic correlation strength) and (2) an increase of the characteristic wave number of the density perturbation. We provide recommendations for the applicability of the considered XC functionals at conditions typical for warm dense matter. Furthermore, we hint at future possibilities for constructing more accurate XC functionals under these conditions.

Warm dense matter is a highly active research area both at the frontier and interface of material science and plasma physics. We assess the performance of commonly used exchange-correlation (XC) approximation (LDA, PBE, PBEsol, and AM05) in the spin-polarized inhomogeneous electron gas under warm dense conditions based on exact path-integral quantum Monte-Carlo calculations. This extends our recent analysis on the relevance of inhomogeneities in the spin-unpolarized warm dense electron gas [Z. Moldabekov et al., J. Chem. Phys. 155, 124116 (2021)]. We demonstrate that the predictive accuracy of these XC functionals deteriorates with (1) a decrease in density (corresponding to an increase in the inter-electronic correlation strength) and (2) an increase of the characteristic wave number of the density perturbation. We provide recommendations for the applicability of the considered XC functionals at conditions typical for warm dense matter. Furthermore, we hint at future possibilities for constructing more accurate XC functionals under these conditions.

I. INTRODUCTION
Matter under extreme conditions [1,2] is an active area of research both in terms of theory and experiment. Interest in this research area is fueled by its importance for astrophysics [3][4][5][6], controlled fusion [7,8], and discovery of novel materials [9][10][11]. In high-energy-density plasma science, extreme conditions are created by various laser-and shock-driven compression techniques [12]. Highly compressed matter at temperatures around the electronic Fermi temperature is referred to as warm dense matter (WDM). Loosely speaking, WDM is a transient state between solids and dense plasmas. As such, understanding the physics in WDM has emerged as a new interdisciplinary challenge for condensed matter and plasma physics [2,13].
Commonly the Born-Oppenheimer approximation is applied which couples the ion dynamics to the electrons through a potential energy surface in the electronic Schrödinger equation. Thus, obtaining an accurate electronic structure is of great significance for WDM modeling. Over the years, a wide range of theoretical methods has been developed for dealing with the electronic structure problem at typical WDM parameters. These include path-integral Quantum Monte-Carlo (QMC) [14][15][16], restricted QMC [17][18][19], Kohn-Sham density functional theory (KS-DFT) [20,21], and orbital-free density functional theory [22].
KS-DFT is often the method of choice due to its balance of reasonable accuracy and affordable computational cost. The temperature generalization of KS-DFT was originally performed by Mermin [23]. Several formal aspects of functional construction at finite temperature have been investigated more recently [24][25][26][27]. The fundamental ingredient to KS-DFT is the exchange-correlation * z.moldabekov@hzdr.de (XC) functional which needs to be approximated in practice. Hence, the predictive capability of KS-DFT relies on the accuracy of the XC functional. Therefore, developing and assessing XC functionals for WDM applications is of pivotal importance.
From a theoretical perspective, the relevance of both electronic correlations and thermal effects strongly impedes our ability to generate accurate data for response properties of electrons in WDM [2,13]. Obtaining highly accurate, first-principles data for both the static [28][29][30][31][32][33][34][35] and dynamic response properties [36][37][38][39][40] of electrons has become possible only very recently by aid of the pathintegral QMC method. This, in turn, has allowed us to present the first assessment of commonly used XC functionals in the inhomogeneous electron gas for the relevant WDM parameter range [41].
In the present work we extend our prior analysis by considering the spin-polarized inhomogeneous electron gas under WDM conditions. To this end, we have performed path-integral QMC calculations of the spinpolarized electron gas under the impact of an external harmonic perturbation [42]. The generated QMC dataset, primarily the electronic densities, at various wave numbers q and amplitudes A of the external perturbation serves as a reference point to assess the accuracy of several XC functionals.
Understanding the accuracy of existing XC functionals for spin-polarized electrons is important for a number of reasons [43]. First of all, the impact of the spinpolarization is crucial for modeling WDM in an external magnetic field. Specifically, non-quantizing strong magnetic fields with an amplitude B ∼ 10 T...10 4 T are generated in experiments related to inertial confinement fusion [44,45], where non-quantizing means that the characteristic quantum kinetic energy T 2 F + T 2 is dominant over the electron cyclotron energy ω c [46] with T denoting the electronic temperature and T F the Fermi temperature [47,48]. The latter condition defines the range of non-quantizing magnetic fields as B/B 0 18.4/r 2 s (T /T F ) 2 + 1 where B 0 2.25 × 10 5 T. Furthermore, electromagnetic pulses applied to WDM can induce spin-polarized states due the effect of the spinponderomotive force [49]. Moreover, the spin-resolved fluid description of quantum plasma dynamics in an external fields requires an XC functional that is capable of describing strongly, non-ideal electrons [49][50][51][52][53][54]. Finally, spin-polarized systems are ubiquitous in quantum chemistry and, consequently, benchmarking XC functionals is relevant for modeling hot-electron phenomena [55].
While a wide range of XC functionals is available in the literature, we select the most relevant approximations used for thermal KS-DFT calculations of WDM. We focus on the basic and most commonly used functionals, including the local density approximation (LDA) in the Perdew-Zunger parametrization [56] and the generalized gradient approximations PBE [57] and PBEsol [58]. Additionally, we consider the Armiento-Mattsson functional (AM05) [59] which is a semi-local GGA based on the notion of interpolating between different model systems, namely the uniform electron gas (UEG) and the Airy gas. It was demonstrated for solids [60] that AM05 is comparable in its accuracy to the hybrid functionals PBE0 [61,62] and HSE06 [63]. This motivated the use of AM05 for the calculation of the equation of state and electronic structures at WDM parameters [64][65][66].
While a number of XC functionals with an explicit temperature dependence has become available recently [15,[67][68][69][70][71], we purposely do not include them in this benchmark study. We are interested in assessing solely the impact of an inhomogeneous electronic structure due to external perturbations at finite q on the accuracy of groundstate XC functionals. Thus, the temperature dependence is included in our assessment implicitly in terms of a Fermi-Dirac occupation of the KS states. An assessment of explicitly temperature-dependent XC approximations shall be the focus of future studies.
Our paper is organized as follows: we begin with providing the theoretical background and technical details of our calculations in Sec. II; we present our benchmark study in Sec. III; finally, we conclude with a summary of our main findings and an outlook on future research in Sec. IV.

II. THEORY AND SIMULATION METHODS
Our assessment of XC functionals is based on the Hamiltonian of a harmonically perturbed, interacting electron gas [72][73][74][75][76] whereĤ UEG denotes the standard Hamiltonian of the UEG [16,47,77], N the number of electrons, and N p the number of harmonic perturbations. The strength of the perturbation is controlled by the amplitude A i and its wavelength by the wave vectors q i . This Hamiltonian has been used in prior work to investigate fundamental physical properties of both the electron gas and electron liquids, e.g., in terms of the non-linear response [29,30,33,78], the local field correction (LFC) [72,73,75,76], plasmons, and other types of excitations [79]. Employing the Hamiltonian in Eq. (1) allows us to asses the performance of different XC functionals in a wide range of conditions ranging from the weakly perturbed UEG where δn/n 0 1 to the strongly inhomogeneous electron gas with regions of large electronic localization where δn/n 0 1 with δn = n − n 0 . While any combination of external perturbations is possible, we limit ourselves to an external perturbation with either one harmonic perturbation (N p = 1) and a combination of two harmonic perturbations (N p = 2). Also, we consider harmonics of the perturbation with the same amplitudes, i.e, A i = A throughout this work.
From a physical perspective, we consider typical WDM conditions which correspond to densities r s = 2 and r s = 6, where r s = (4/3πn) −1/3 is the mean interelectronic distance in Hartree atomic units. The temperature is set equal to the Fermi temperature of an unpolarized electron gas [15]. Accordingly, at r s = 6 and r s = 2 we have T 1.4 eV and T 12.5 eV, respectively. These parameters are achieved in solid targets by isochoric heating and laser-induced shock compression experiments [80][81][82].
The KS-DFT calculations are carried out using the GPAW code [83][84][85][86]. GPAW is a real-space implementation of the projector augmented-wave method which allows us to fix the spin polarization as an external parameter. A Monkhorst-Pack [87] sampling of the Brillouin zone was used with a k -point grid of 8 × 8 × 8. The calculations are performed using a plane-wave basis where the cutoff energy has been converged to 540 eV. The number of orbitals is set to 250 at r s = 2 and to 100 at r s = 6. At the considered temperatures, the smallest occupation number at r s = 2 is less than 10 −8 and at r s = 6 the smallest occupation number is less than 10 −5 . It is common practice in thermal KS-DFT calculations to converge the occupation threshold to this value [88]. Moreover, this criterion was used successfully in our prior study on the spin-unpolarized electron gas where the same Hamiltonian was used [41].
We generate accurate QMC benchmark data by carrying out direct path-integral Monte Carlo calculations without any nodal restrictions [42] based on the Hamiltonian in Eq. (1). Therefore, the simulations are extremely costly, but are exact within the given error bars. The compute time is up to O 10 5 CPUh for the most challenging parameter sets. The direct path-integral QMC method is afflicted with the notorious fermion sign problem [89,90], which leads to an exponential increase in the compute time with increasing the system size N or decreasing the temperature T . This problem is further exacerbated in the spin-polarized case, which results in a more pronounced impact of fermionic exchange effects and, therefore, a more severe sign problem. Yet, we stress that our results are not afflicted with any nodal errors, which have been shown to be considerable in the WDM regime [91,92]. In addition, we note that we have employed a canonical adaption [93] of the worm algorithm by Boninsegni et al. [94,95]. Finally, we have utilized a straightforward primitive factorization [96,97] of the density matrix using P = 200 high-temperature factors. The convergence with P has been checked carefully.
We choose N = 14 electrons within a simulation cell where its size is defined by 4/3πnL 3 = 1. Accordingly, the smallest wave number of the external perturbation is q min = 2π/L 0.843q F , where q F = (3π 2 n 0 ) 1/3 and n 0 denotes the mean density when there is no external perturbation. Note that our calculations do not suffer from finite size effects at the considered parameter range. This has been demonstrated previously in the assessment of the spin-unpolarized electron gas [41] . Additionally, it was shown in previous QMC studies of the warm dense electron gas [28,30,32,33] that the electronic density response is well converged with respect to the number of electrons for as few as N = 14 electrons.

III. RESULTS
Our assessment of XC functionals centers on the WDM application domain, i.e., partially degenerate plasmas. To better quantify these conditions, we consider the ratio of the Coulomb interaction potential between two elec- trons separated by the mean distance (1/r s ) with the Fermi energy (E F ∼ r −2 s ). As this is directly proportional to the density parameter r s , we employ r s as the characteristic coupling parameter of electrons [48] for the partially degenerate plasma state considered here. For ensuring a comprehensive analysis, we perform calculations for both r s = 2 and r s = 6 which are conditions representative for WDM, where we characterize r s = 2 as moderate coupled and r s = 6 as strongly coupled regimes. At r s = 2, the kinetic energy of the system is dominant over the XC energy. Contrarily, at r s = 6, the XC energy prevails over kinetic energy [98][99][100].
A. Moderate coupling, rs = 2 We start with r s = 2 and a single harmonic perturbation. This choice of parameters allows us to probe the performance of the XC functionals at different values of the perturbation wave number and amplitude. In principle, any other external perturbation can be represented as a sum of harmonic perturbations. Therefore, the conclusions obtained for the single harmonic perturbation are generally valid and should be transferable. As we show below, this claim is validated by considering the combination of two harmonic perturbations with different wave numbers.
In Fig. 1, we show the total electron density along the direction of the perturbation with an amplitude A = 0.1 and at three different values of the wave numbers q min , 2q min , and 3q min . From the top panel we deduce that the combination of A = 0.1 and q 1 = q min leads to small changes in the electron density δn/n 0 < 1. Here, the KS-DFT results are overall in good agreement with the QMC data. The same is observed for both q 1 = 2q min and q 1 = 3q min , as shown in the bottom panel of the same figure.
Let us next consider a substantially stronger external perturbation. As shown in Fig. 1, we find that an increasing perturbation amplitude (A = 0.5) results in a stronger density change with δn/n 0 1. Also here the agreement between KS-DFT and QMC data is quite excellent on the scale of the total density for all considered wave numbers. We note that the spin-polarization of the electrons leads to a weaker perturbation at a given amplitude A compared to the spin-unpolarized case considered in Ref. [41]. The physics behind this observation is discussed in detail in the Appendix.
To gain deeper insight, we consider the density perturbation ∆ n which is the physically more important quantity here. It defines the (linear) static response function in the case of a weak external field since δn = 2Aχ(q) cos(q · r). Therefore, we further quantify differences in the density perturbation by considering the Note that at A = 0.1 (A = 0.5), the uncertainty is ±0.39% (±1.5%) at q1 = qmin, ±0.55% (±0.5%) at q1 = 2qmin, and ±0.78% (±0.3%) at q1 = 3qmin.
where max{δn} is the maximum deviation of the QMC data from the mean density. This quantity provides a meaningful point of reference to gauge the quality of the DFT data. In Figs. 3 and 4 we illustrate the relative density devi-  Table I. At q 1 = q min and q 1 = 2q min , the maximum deviation of the KS-DFT data from the QMC data is about 2% and all considered XC functionals provide a similar quality within the given statistical uncertainty of the QMC data. When we increase the perturbing wave number to 3q min , we observe larger errors in the KS-DFT data, where the calculations based on the PBE functional exhibit the largest errors of about 6%. In contrast to that, the AM05 functional performs well with a maximum deviation of about 3%. The cal-culations based on LDA and PBEsol functionals yield a comparable accuracy in ∆ n with maximum deviations of 5% and 4%, respectively. Increasing the perturbation strength to A = 0.5, as shown in Fig. 4, generally leads to less accurate KS-DFT results. The PBE performs worst for the two largest considered wave numbers. At A = 0.5 and q 1 = 3q min , the error in PBE is close to 9%, while the LDA yields an error of about 4%. Within the statistical uncertainty of the QMC data, the overall quality of the KS-DFT calculations based on LDA, PBE, and AM05 at A = 0.5 is similar to that at A = 0.1. Nevertheless, we can single out LDA as the XC functional providing the best quality on average for q 1 = q min , q 1 = 2q min , and q 1 = 3q min . Additionally, in both cases δn/n 0 < 1 and δn/n 0 1, we conclude that all considered XC functionals provide an accurate description of the total density at perturbation wave numbers q ≤ 2q min 1.7q F . At larger wave numbers 3q min the considered XC functionals start becoming less accurate.

B. Strong coupling, rs = 6
Next, we investigate the strongly coupled regime at r s = 6. We again consider a single harmonic perturbation with different amplitudes and wave numbers. The strongly coupled regime turns out to be more challenging for the XC functionals under examination.
We begin with considering δn/n 0 1 which is achieved by a weak perturbation with A = 0.01. The total electron density along the perturbation direction is shown in Fig. 5 for three different values of the perturbation wave number, q min , 2q min , and 3q min . While at small wave numbers, such as q 1 = q min in the top panel, all XC functionals perform well, significant deviations from the QMC reference data become evident for most XC functionals at larger wave numbers, e.g., at 2q min and 3q min (bottom panel).
In Fig. 6, we further examine the disagreement between the KS-DFT and QMC calculations using the relative density deviation. We infer that at q 1 = q min all considered XC functionals have a maximum deviation of about 2% in the central region around the maximum of the density and deviate up to about 3% around the density minima. At q 1 = 2q min , the AM05 functional has an error abour 5% around the density minima and maxima. Contrarily, at q 1 = 2q min other XC functionals provide a rather accurate description, the PBE data yielding a maximum error of about 2% and both the LDA and PBEsol being virtually exact. When the perturbing wave number increases to 3q min , the LDA, PBE, and PBEsol functionals suffer from large errors as shown in the bottom panel of Fig. 6. Specifically, LDA and PBE have an error of about 15% and PBEsol of up to 10%. In contrast to that and to the case of small wave numbers, the AM05 functional provides relatively accurate results with a maximum error of about 2.5% at large wave numbers. The high accuracy of the KS-DFT calculations based on the considered XC functionals in the case of a weak perturbation is due to connection of the XC energy density with the long wavelength of XC kernel of the static response function via the compressibilty sum-rule [101]. For the same reason, the accuracy of the KS-DFT results becomes worse with an increasing wave number of the perturbation. We discuss this point in more detail in Sec. IV.
It is noteworthy that the AM05 functionals performs better for perturbations at large wavenumbers (compare the bottom panel of Fig. 6 with the top and middle). This can be understood heuristically by the fact that AM05 is based on an two electronic reference systems, namely the UEG and the Airy electron gas [102]. It interpolates between the slowly varying limit in jellium bulk, q → 0, and the limit far from the jellium surface with a characteristic wave number q/q F 1 representing the density changes. The latter is the reason why the AM05 functional performs better at 3q min . Note that this assessment is valid at r s = 2 and r s = 6 if the perturbation is weak δn/n 0 < 1, but not when the density perturbation is strong δn/n 0 > 1. We extend this discussion towards the end in Sec. IV.
Next, we consider strong perturbations A = 0.1 which cause strong density changes δn/n 0 > 1 in Fig. 7. It illustrates the total density at the perturbation wave numbers 2q min and 3q min . Similar to weak perturbations, the KS-DFT results differ significantly from the QMC data. We analyze these differences further by computing the the relative density deviations shown in Fig. 8. Here, the PBE functional performs worst with maximum errors of about 8% and 15% maximum deviations at 2q min and 3q min , respectively. The AM05 and PBEsol functionals exhibit both comparable errors of up to 6% at 2q min and about 13% at 3q min . The LDA performs best with a maximum error of about 4% at 2q min , but still with large errors of about 10% at 3q min .
At the smallest considered wave number of the perturbation, q min , the KS-DFT data is relatively accurate with a maximum error of about 3% and 4%. The maximum values of the relative density deviation are listed in Table II.
We conclude our assessment of the strong-coupling regime by summarizing that LDA, PBE, and PBEsol provide an accurate description of the total density for relatively small wave numbers q ≤ 2q min 1.7q F and small density changes δn/n 0 < 1. We point out that AM05 does not show such a clear trend with a change in the perturbing wave number when δn/n 0 < 1. We provide a heuristic explanation in the Sec. IV. However, when the perturbation is strong (δn/n 0 > 1), all considered XC functionals yield significant errors with respects to the QMC reference data at both q = 2q min 1.7q F and q = 3q min 2.53q F .

C. Combined effect of two harmonic perturbations
In general, any perturbation can be represented as a linear combination of harmonic perturbations. To provide supporting evidence that our analysis is not an artefact of using a single harmonic perturbation, we consider the combination of two harmonic perturbations with the wave numbers q 1 = q min and q 2 = 2q min . Note that TABLE II. The performance of common XC functionals in terms of the relative deviation in the density ∆ n [%]. A single harmonic perturbation at a fixed density rs = 6, perturbation amplitudes A = 0.01 and A = 0.1, and wave numbers qmin ≤ q ≤ 3qmin is considered, where qmin = 0.843 qF . The largest absolute values of the error are listed. Note that at A = 0.1 (A = 0.5), the uncertainty in the QMC data is ±0.35% (±0.25%) at q1 = qmin, ±0.18% (±0.13%) at q1 = 2qmin, and ±0.25% (±0.11%) at q1 = 3qmin.  we choose the same amplitudes, namely A = 0.1 in the moderately coupled regime (r s = 2) and A = 0.01 in the strongly coupled regime (r s = 6). These correspond to weakly perturbed states for which the analysis based on the single harmonic perturbation showed failure of AM05 at r s = 6 and 2q min (see middle panel in Fig. 6). Now, we can test whether the same holds when the perturbation is the combination of harmonics with different wave numbers. The resulting electron density and the relative density deviation are shown in Figs. 9 and Fig. 10. From these we infer that LDA, PBE, and PBEsol provide an accurate description of the total density for both r s = 2 and r s = 6 when the perturbation is weak. The same is valid for AM05 at r s = 2. However, AM05 becomes less accurate for strong coupling (r s = 6). We note that maximum errors in the relative density deviation follow a similar trend as observed for a single harmonic perturbation with the same values of the perturbation amplitudes (cf . Figs 3 and 6). Therefore, we confirm that our conclusions based on a single harmonic perturbation are generally applicable.

IV. CONCLUSIONS AND OUTLOOK
We have benchmarked the accuracy of KS-DFT calculations based on the LDA, PBE, PBEsol, and AM05 XC functionals against exact QMC data for a spin-polarized, partially degenerate, inhomogeneous electron gas in the regime of both weak and strong perturbations. Furthermore, we have considered both moderate and strong coupling regimes which are relevant for current WDM appli-cations. We stress that a comprehensive assessment of the vast amount of the available XC functionals has not been attempted in this work. Instead, we focused on the LDA functional and it is generalizations, which are most often used for WDM simulations. We therefore provide access to our QMC data to facilitate assessments of any existing or newly developed XC functional by the scientific community [103]. With the present work devoted to the spin-polarized case and our earlier work to the spin-unpolarized case [41], we conclude assessing common, ground-state XC functionals within the harmonically perturbed, warm dense electron gas.
Our analysis has revealed (1) a parameter range where the KS-DFT method can be used to obtain QMC-level accuracy and (2) conditions when the considered XC functionals become inaccurate and eventually fail to yield an accurate electron density.
We summarize the main conclusions regarding the applicability of the considered XC functionals as follows: • For moderate coupling (r s = 2) and perturbing wave numbers q 2.53q F , LDA, PBEsol, and AM05 provide an accurate electron density for both weak (δn/n 0 < 1) and strong (δn/n 0 > 1) density perturbations.
• For strong coupling (r s = 6) and strong density perturbations (δn/n 0 > 1), all considered XC functionals provide an accurate electron density at small perturbing wave numbers q < q F . However, all of them become unreliable for larger wave numbers, q > q F . Some computational aspects regarding the choice of XC functional are discussed in Appendix B. While this discussion is rather technical, we believe it can be helpful for users of KS-DFT to save compute time and resources.
Apart from delivering useful benchmarks, our assessment advances our understanding on the expected accuracy of KS-DFT when used for spin-polarized simulations at parameters relevant for WDM applications. Clearly, the considered XC functionals struggle to provide accurate results for perturbations with large wave numbers. These conditions are, however, relevant to contemporary WDM experiments where the impact of the external driving fields and of shock formation causes such perturbations in the electronic structure.
We can rationalize the observed errors at large wave numbers by considering the static linear response function.
It is defined as χ(q) −1 = χ 0 (q) −1 − 4π/q 2 + f xc (q) , where χ 0 denotes the response function of the UEG and f xc (q) the XC kernel. Using the LDA is equivalent to working within the long wavelength limit of the XC kernel f xc (q 0), which is connected to the XC energy via compressibility sum rule [101]. The same applies to PBE and PBEsol, for which exchange and correlation gradient corrections are designed to cancel out each other in the limit of the UEG [104]. Therefore, our analysis points out a significant deficiency of common XC functionals in terms of the wave number q. When perturbations at values larger than the Fermi wave number are present, common XC functionals will not be reliable.
An interesting trend is observed for AM05 at δn/n 0 < 1. The accuracy of AM05 first drops with increasing wave number from q min 0.84q F to 2q min . However, with a further increase in the wave number to 3q min , its accuracy improves. We can shed light into this behavior by looking more closely how AM05 is constructed. The exchange energy in AM05 is an interpolation of two subsystems: the LDA exchange energy and the exchange energy of the Airy gas. Similarly, the correlation energy of AM05 is interpolated between LDA correlation and a scaled LDA correlation such that it reproduces the jellium surface energy. As an integral part of the AM05 functional, the Airy gas is a model system which describes the electrons on surfaces. It is characterized by a decaying electron density with increasing distance from the surface. The connection between the UEG and Airy gas is established using the characteristic wave number of the density variation q ∼ |∇n|/n. Essentially, AM05 reproduces the strongly inhomogeneous Airy gas at large wave numbers s = q /(2q F ) > 1. Therefore, in our assessment, AM05 performs better than LDA, PBE, and PBEsol, for weak density perturbations δn/n 0 < 1 with 3q min > 2q F . This is no longer the case when the perturbation is strong δn/n 0 > 1. This might be related to the fact that only the exchange energy of the Airy gas was used in the construction of AM05, while its correlation energy is based on a scaled LDA due to the lack of exact QMC data for the Airy gas. Therefore, this result hints at the need for accurate QMC data for the Airy gas. Furthermore, our findings reveal that the interpolation scheme used for AM05 struggles with density perturbations at intermediate wave numbers q F < q < 2q F when δn/n 0 < 1. At δn/n 0 > 1, the AM05 is less accurate than LDA. We note that these conclusions are valid also for unpolarized electrons. An in-depth analysis of the AM05 parametrization shall be presented elsewhere.
Our presented findings suggest that interpolating between reference systems in the fashion of the AM05 functional might be a promising route to XC functionals accurate across a large range of wave numbers needed for WDM applications.
Finally, we note that our conclusions are relevant for further developing the non-linear response theory of the electron gas and liquids at WDM conditions and beyond [30,33,78]. Indeed, various available QMC methods are restricted with respect to the parameter range due to the infamous fermion sign problem [16,89,90,FIG. 11.
Comparing the electron density in the spinpolarized electrons gas with the spin-unpolarized case. The density is shown along the perturbation direction for A = 0.1 and rs = 2 at q1 = qmin (top) and q1 = 3qmin (bottom). 105,106]. KS-DFT allows us to study non-linear response phenomena at temperatures T < T F and densities 2 r s 6 that are difficult to reach with QMC methods.

ACKNOWLEDGMENTS
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 A: Spin-polarized versus spin-unpolarized electron gas While the spin-unpolarized electron gas has been assessed elsewhere [41], we provide a concise comparison with the spin-polarized case for completeness. To this end, we contrast the electron densities at r s = 2 and a temperature T = T F = 12.5 eV for a single harmonic perturbation with A = 0.1 in Fig. 11, The top panel corresponds to a perturbing wave number q 1 = q min and the bottom panel to q 1 = 3q min . From the top panel (q 1 = q min ) we infer that the density perturbation induced in the spin-polarized electron gas is weaker than that in the spin-unpolarized case. This can be understood by recalling that the spin-polarized state has a larger Fermi energy. Therefore, the external field has a smaller amplitude on the scale of the Fermi energy which sets the energy scale for the system. Contrarily at q 1 = 3q min (bottom panel), the induced density perturbation is virtually identical for both the spin-polarized and spin-unpolarized cases. At large wave numbers, the response is dominated by single-particle effects. Consequently, both spin polarizations of the electron gas exhibit a similar response to an external field. For completeness, we note that the external field considered here does not interact with electronic spin.

Appendix B: Computational scaling
The amount of CPU hours required to obtain converged results for different parameters at r s = 2 and r s = 6 is shown in Figs. 12 and 13.
At r s = 2, the calculations with PBEsol turned out most expensive, while the other XC functionals were similar in their computational cost (except AM05 at A = 0.1 and 2q min ). In contrast to this, at r s = 6, the calculations with PBEsol were overall much cheaper. Note that, when δn/n 0 > 1, the calculations based on AM05 required substantially more compute time than those for LDA, PBE, and PBEsol.

DATA AVAILABILITY
The data supporting the findings of this study are available on the Rossendorf Data Repository (RO-DARE) [103].
FIG. 12. Comparison of CPU hours required to obtain converged results using different XC functionals at rs = 2 for weakly (top) and strongly (bottom) perturbed spin-polarized electrons.