Density Response of the Warm Dense Electron Gas beyond Linear Response Theory: Excitation of Harmonics

In a recent Letter, Dornheim et al. [PRL 125, 085001 (2020)] have investigated the nonlinear density response of the uniform electron gas in the warm dense matter regime. More specifically, they have studied the cubic response function at the first harmonic, which cannot be neglected in many situations of experimental relevance. In this work, we go one step further and study the full spectrum of excitations at the higher harmonics of the original perturbation based on extensive new ab initio path integral Monte Carlo (PIMC) simulations. We find that the dominant contribution to the density response beyond linear response theory is given by the quadratic response function at the second harmonic in the moderately nonlinear regime. Furthermore, we show that the nonlinear density response is highly sensitive to exchange-correlation effects, which makes it a potentially valuable new tool of diagnostics. To this end, we present a new theoretical description of the nonlinear electronic density response based on the recent effective static approximation to the local field correction [PRL 125, 235001 (2020)], which accurately reproduces our PIMC data with negligible computational cost.


I. INTRODUCTION
Linear response theory (LRT), i.e., the assumption of a linear response of a system of interest to a sufficiently small perturbation, is ubiquitous throughout physics and related disciplines. Prominent examples include the density and current response of many-body systems to external fields, e.g. [1], the computation of optical absorption or plasmon spectra [2,3], the probing of superfluidity [4,5], the description of electron-phonon interactions [6], and dynamical simulations based on density functional theory [7]. From a practical point of view, the linear response only depends on the equilibrium properties of the unperturbed system, which often makes an accurate theoretical description feasible in the first place. This also allows one to characterize a system in terms of its response, which is of fundamental importance for diagnostics.
Yet, while it is clear that LRT is only accurate for small perturbations, often the validity of this assumption is not checked. On the one hand, this is quite understandable, as the general nonlinear response is substantially more complicated than LRT and, consequently, much poorer understood. On the other hand, it is well known that nonlinear effects play an important role for many applications, e.g., in the optical excitation and ionization of atoms [8,9], in nonlinear optics [10], in nuclear physics [11], in laser plasmas [12], in the inverse bremsstrahlung heating of plasmas [13,14], nonlinear plasmons in magnetized plasmas [15][16][17], the excitation spectrum of graphene [18][19][20], the optical prop- * t.dornheim@hzdr.de erties of large molecules [21], or intersubband transitions in quantum-well semiconductor heterostructures [22].
The present work is devoted to the nonlinear density response of warm dense matter (WDM), an exotic state with extreme temperatures (T ∼ 10 4 − 10 8 K) and densities in the vicinity of solid states [23][24][25]. These conditions occur naturally in astrophysical objects such as giant planet interiors [26,27] and brown dwarfs [28,29]. In addition, the compression path of a fuel capsule towards inertial confinement fusion has to traverse the WDM regime [30]. Finally, we mention the potential of hot electrons as a catalyst for chemical reactions [31], and the recent discovery of new materials such as lonsdaleite [32] and nanodiamonds [33] at extreme conditions. Due to this gamut of applications, WDM has emerged as one of the most active frontiers in plasma physics and material science, and WDM is nowadays routinely realized in experiments at large research facilities around the globe; see the recent review article by Falk [34] for an overview of different experimental techniques.
Unfortunately, the theoretical description of WDM is notoriously difficult due to the intriguingly intricate interplay of Coulomb correlations, thermal excitations, and quantum effects. The correlation strength rules out perturbative expansions such as Green functions [35], a nonvanishing temperature precludes the applications of the well-stocked arsenal of ground-state methods, and quantum effects make semi-classical methods like molecular dynamics insufficient. Therefore, computationally expensive quantum Monte Carlo methods [36,37] have emerged as the method of choice to determine the basic properties of WDM [38][39][40][41][42][43][44].
Like theoretical concepts in many other fields, current WDM theory is based, to a large extent, on LRT. For example, the theory of X-ray Thomson scattering [45,46], the most capable diagnostics in WDM experiments, presupposes a purely linear relation between perturbation and response. Other examples include the construction of electronically screened effective potentials [47][48][49], the estimation of energy-relaxation rates [50][51][52], and the description of energy-loss properties like the stopping power [53,54]. In this context, the central property is often given by the linear density response function, χ(q, ω) (with q and ω being the wave number and frequency of the perturbation), of the electrons in the system. These can be approximated by the uniform electron gas (UEG) [1], one of the most fundamental model systems in physics and quantum chemistry [36,55]. Consequently, many recent works [56][57][58][59][60] have been devoted to the estimation of χ, which has culminated in the complete characterization [58] of the static local field correction (LFC) G(q) containing the full wave-number resolved information about electronic exchange-correlation effects (see also Eq. (5) below). Furthermore, Dornheim et al. [61,62] have introduced the effective static approximation to the LFC, which is available as a simple analytical representation [62] and allows for the computation of highly accurate electronic properties within LRT with negligible computation cost.
In a recent Letter [63], three of us carried out extensive ab initio path integral Monte Carlo (PIMC) simulations of the harmonically perturbed electron gas in the WDM regime without any assumptions about the response being linear. This allowed us to check the validity range of LRT, and we found that nonlinear effects cannot be neglected in many situations of experimental relevance. In addition, we used the exact PIMC data for the density response to extract the cubic density response function at the wave number of the perturbation itself.
In the present work, we substantially extend these efforts by considering the full excitation spectrum over the entire range of wave numbers. In particular, we find strong nonlinear excitations at the integer higher harmonics of the original perturbation. The quadratic response at the second harmonic constitutes the largest nonlinear contribution for moderate perturbations. Furthermore, we find that the nonlinear response functions strongly depend on electronic exchange-correlation effects, thus making them a potentially valuable new tool of diagnostics. To this end, we introduce a new theory of the nonlinear density response based on available representations of the LFC [58,62], which is capable to accurately reproduce our PIMC data with negligible computational cost.
The paper is organized as follows: In Sec. II, we introduce the relevant theoretical background, including the PIMC method (II A), the simulation idea and setup (II B), the extended theory of the nonlinear density response functions (II C), and the specific way how the latter can be estimated from our PIMC data (II D). In Sec. III, we present our extensive new simulation results covering a metallic density (III A), the temperature dependence (III B), and a strongly coupled case (III C). The paper is concluded by a brief summary and outlook in Sec. IV.

II. THEORY
We assume Hartree atomic units throughout this work. From a theoretical perspective, the WDM regime is typically defined by two characteristic parameters that are both of the order of one: (i) the density parameter (Wigner-Seitz radius) r s = r/a B , where r and a B are the average inter-particle distance and first Bohr radius, and (ii) the reduced temperature θ = k B T /E F , with E F being the usual Fermi energy [1]. A third principle parameter is given by the spin-polarization ξ = (N ↑ − N ↓ )/N . Throughout this work, we restrict ourselves to the case of a purely unpolarized system, i.e., N ↑ = N ↓ = N/2 and ξ = 0.

A. Path Integral Monte Carlo
We carry out path-integral Monte-Carlo (PIMC) simulations [4,64,65] of N electrons in the canonical ensemble, i.e., the volume V = L 3 (with L being the box length), inverse temperature β = 1/k B T , and particle number N are constant.
The basic idea of the PIMC method is to stochastically sample the thermal density matrix where R = (r 1 , . . . , r N ) T contains the coordinates of all N electrons. In the end, the partition function is re-cast into the integral over P density matrices at a P times higher temperature, and the resulting high-dimensional integrals are evaluated using the Metropolis Monte-Carlo method [66] to avoid the curse of dimensionality of standard quadrature methods; see Ref. [4] for an extensive review article on the PIMC method. An additional obstacle regarding the PIMC simulations of electrons is the requirement of anti-symmetry under particle exchange, which renders the partition function into a sum over both positive and negative contributions. The resulting cancellation of terms with opposite sign is the origin of the notorious fermion sign problem [67,68] that leads to an exponential increase in computation time both with increasing the system size N and with decreasing temperature T [68].
For this reason, one often employs approximate methods like the restricted path-integral Monte-Carlo method [69] that are not afflicted with the sign problem. Unfortunately, this advantage comes at the cost of an in practice uncontrolled systematic error [43,70,71]. Therefore, we here avoid such nodal restrictions and carry out the computationally expensive, yet exact standard PIMC simulations and deal with the sign problem by increasing the computation time, where this is feasible.
For completeness, we mention that all PIMC results presented in this work have been obtained using a canonical adaption [72] of the worm algorithm by Boninsegni et al. [73,74].

B. Simulation idea
Following Refs. [75][76][77][78][79][80], we simulate a uniform electron gas that is subject to an external static harmonic perturbation. The corresponding Hamiltonian is given byĤ whereĤ UEG is the standard (unperturbed) UEG Hamiltonian [1,36,55]. The second term on the RHS. of Eq. (2) is the external potential corresponding to a single cosine of wave vector q = 2π/L(n x , n y , n z ) T (with n i ∈ Z, and L being the length of the simulation box) and the perturbation amplitude A. As A corresponds to an energy in this notation, it is given in Hartree throughout this work. Moreover, we restrict ourselves to a strictly static perturbation as PIMC simulations of a time-dependent system are severely hampered by an additional dynamic sign problem [81].
To measure the density response, we compute the induced density where . . . q,A indicates the parameters of the perturbation in Eq. (2), and k corresponds to the wave vector at which we measure the response of the system. In particular, within linear response theory it holds with χ (1) (q) being the standard linear response function [1]. The latter is conveniently expressed as where v(q) = 4π/q 2 is the Fourier transform of the Coulomb potential (we use atomic units throughout), χ 0 (q) is the density response function of the ideal (noninteracting) Fermi gas that is known from the literature [1], and G(q) is the static local field correction containing the full wave-number resolved information about exchange-correlation effects. In particular, highly accurate data for G(q) have recently become available as a neural-net representation that was constructed on the basis of extensive PIMC simulations [58,61].
In addition, we mention that the full wave-number dependence of χ (1) (q) can be computed from a single simulation of the unperturbed UEG via the imaginary-time version of the fluctuation-dissipation theorem, which states that with F (q, τ ) being the intermediate scattering function [45] evaluated at an imaginary time argument τ ∈ [0, β], see Refs. [58,60,[83][84][85][86][87][88][89] for different applications of this quantity. Finally, we mention the density of the perturbed electron gas in coordinate space, which, in LRT, is given by with n 0 being the average value of the density.

C. Theory of nonlinear density response
The nonlinear density response of plasmas and condensed matter systems in and out of equilibrium has been studied in some detail theoretically and by simulations, e.g., [90][91][92][93]. Here we concentrate on the nonlinear response functions in thermodynamic equilibrium. The general definitions of the second and third order response functions are given by an expansion of the induced density, n ind (r) = n(r) − n 0 , according to a weak perturbation assumption [94,95] n ind (r) = dr χ(r, r )V (r ) + dr dr Y (r, r , r )V (r )V (r ) + dr dr dr Z (r, r , r , r )V (r )V (r )V (r ) where V is the perturbing potential, whereas Y and Z are the second order and third order response functions, respectively. Taking into account that, for a homogeneous system, (8) can be rewritten in Fourier space as (Ω denotes the volume) where wave-vector notation k is used to avoid confusion with the wave-vector q of the external harmonic field V (r) = 2A cos q · r.
By substituting into Eq. (9) the Fourier representation of the external potential (A = A/Ω), we arrive at The first line contains the linear response χ (1) at the incoming wave vector q. The second line collects the quadratic response Y at the second harmonic 2q. The third line gives the cubic response at the first harmonic χ (1,cubic) where the following notation for brevity was introduced Finally, the third and fourth lines contain the cubic response function Z at the third harmonic 3q. We will drop all negative wave vectors from now on as the behavior is symmetric for positive and negative q. It should also be noted that only integer multiples of the original wave vector appear as possible higher order excitations. Now, by performing an inverse Fourier transform of Eq. (11), we obtain, for the induced density in real space, where, for η = 1, η = 2, and η = 3: and, using Eqs. (11) and (12) we identify Equations (17) and (19) show the connection between the general second order, Y , and third order, Z , response functions with the response functions of the system at the second and third harmonics [χ (2) (q) and χ (3) (q), respectively]. Furthermore, Eq. (18) shows the connection of the cubic response at the first harmonic, χ (1,cubic) (q), to the third order response function Z .
By definition, the total density in coordinate space n(r) can as well be written as the Fourier series over all wave numbers q, Yet, as only the harmonics of the perturbation wave vector q exhibit a nonzero response, the density is simply
Mikhailov [19,106] expressed the ideal response functions at the second harmonic, χ 0 (q), without invoking Y 0 (k, q) and Z 0 (k, k , q), but by directly expanding the induced density at harmonics. The recursion relations derived by Mikhailov [19,106] are and [19] Such a recursion formula does not exist for the ideal cubic response at the first harmonic χ (1,cubic) (q). In fact, the terms making up the ideal cubic response at the first harmonic diverge.

Random Phase Approximation
Taking into account screening in linear approximation, the RPA result for the response function from the Green's functions based consideration has the following form [94] At the second harmonic, χ RPA (q) is found using the relation (17) and Eq. (23) To illustrate the level of approximation in Eqs. (23) & (24), it is instructive to provide an alternative derivation of Eq. (24). To begin with, we recall that, in contrast to the ideal response function, the response function in RPA takes into account inter-electronic interactions in a mean field approximation. Therefore, introducing a total potential as a sum of the external potential and the Hartree potential due to the induced electronic density, Φ tot = Φ ext +Φ ind , the induced density in RPA at the second harmonic is obtained in terms of the ideal response functions  (14)].
On the other hand, the induced density in RPA is expressed as the response to the external field using the quadratic response function in RPA From Eqs. (26) and (25), we find To solve Eq. (27), we combine the Poisson equation for the induced potential at the second harmonic with Eq. (26) and find a relation between Φ ind (2q) and Φ ext (q) where it was taken into account that there is no contribution from the external potential to the total potential at the second harmonic. Secondly, we neglect the contribution from the cubic term, i.e. set χ (1,cubic) 0 (2q) = 0 in Eqs. (25) & (27), and, thus, approximate the total potential at the first harmonic using linear response theory, which reads in RPA Using Eqs. (28) and (29), we can solve Eq. (27) with respect to χ RPA (q) to find Eq. (24). For further discussion of the QMC data, it is important to emphasize that, in χ (2) RPA defined by Eq. (24), (i) the screening is taken into account only in linear approximation and (ii) the contribution due to the cubic response function is neglected.
Next, by following the steps from Eq. (25) to Eq. (27), but for the induced density at the third harmonic, and making use of the screening in linear approximation, we find for the response function at the third harmonic in RPA . (30) Note that, for a weak perturbation, the neglected contribution to the screening due to the cubic response, χ (1,cubic) 0 , in Eq. (24) is a higher order correction compared to the quadratic response. In contrast, for χ has the same order contribution as χ (3) 0 . Therefore, it can be expected that Eq. (24) performs better than Eq. (30) when applied for the description of the QMC data at the corresponding harmonics.
The corresponding cubic response in RPA at the first harmonic is given by It is known that there are further terms that are technically contributing to the cubic response at the level of RPA. These terms are made up from non-diagonal quadratic response functions entirely [107,108].

Going beyond the RPA description by using LFC
In Eqs. (24) and (30), the electronic interactions effect is included on the basis of the linear response functions in the denominator. Therefore, in analogy to the usual practice in linear response theory, we can go beyond RPA by introducing a local field correction (LFC), G(q). In this way, we arrive at the following equations for the response functions at the second and third harmonics with higher order electronic exchange-correlations effect included by using LFC and Equations (33) and (34) provide an improved description of the response functions on the second and third harmonics with electronic exchange-correlations effect taken into account on the level of the linear response theory. We still give the respective result for the cubic response at the first harmonic The temperature and density dependent static LFC has been obtained from QMC simulations recently [56][57][58][59][60]. Furthermore, Dornheim et al. [61,62] have introduced the new concept of the effective static approximation to the LFC, which is available as a simple analytical representation [62].

D. PIMC approach to nonlinear density response functions
To extract the nonlinear density response functions of different harmonics from our PIMC simulation, we evaluate Eq. (3) given above for different perturbation amplitudes A and subsequently perform fits to these data.
In particular, the density response in reciprocal space for the same wave number as the perturbation is fitted to [77,78] where the coefficients χ (1) (q) and χ (1,cubic) (q) are the free parameters. Note that the determination of χ (1) (q) from the PIMC data via Eq. (36) is redundant as it is already known from previous simulations of the unperturbed UEG via Eq. (6). However, these independent benchmark data can be compared to our new results for χ (1) (q) and, thus, constitute a valuable consistency check of our approach, see also Ref. [75]. Furthermore, it is important to note that Eq. (36) only holds up to a maximum perturbation strength, beyond which contributions with a higher order of A start to substantially contribute; see Sec. III for a hands-on discussion of this point. The density response of the second harmonic is, in first order in A, given by where χ (2) (q) is the free parameter that is obtained from a fit to the A-dependence for a fixed perturbation wave number q. Similarly, the density response of the third harmonic is obtained by fitting the PIMC data for the response at k = 3q to Let us start our investigation by considering the response of the electron gas at r s = 2 and θ = 1. These conditions are located at the center of the WDM regime and are realized experimentally for example in experiments with aluminum, e.g., Ref. [109]. In Fig. 1, we show the density response of the system in dependence on the perturbation amplitude A for the same wave number as the perturbation, i.e., for the first harmonic ρ q q,A . The green crosses show the raw PIMC data which have been obtained by evaluating Eq. (3). The solid red line shows the prediction by LRT and has been obtained from a previous simulation of the unperturbed UEG [58] by evaluating Eq. (6). Evidently, LRT is in excellent agreement with the PIMC data for small A only and the disagreement rapidly increases with A. This can be seen particularly well in the bottom panel of Fig. 1 showing the relative deviation to the PIMC data. For completeness, we mention that a similar analysis has been presented in Ref. [63], although for a different wave number.
The dotted blue curve has been obtained by fitting Eq. (36) to the PIMC data within the interval A ∈ [0, 0.225], see also the vertical dashed grey line. Evidently, this curve captures the emerging deviations between the green crosses and the red line and remains accurate for substantially larger values of A. In particular, it remains accurate even beyond the fitting interval, which constitutes a strong empirical confirmation of the functional form in Eq. (36). In addition, the dashed black line shows the LRT prediction using the first coefficient from the fit and perfectly agrees with the independent red line. Still, we note that the blue curve, too, eventually becomes inaccurate, as the higher-order terms in A that are neglected in Eq. (36) start to significantly contribute.
A possibly more intuitive illustration of the impact of the harmonic perturbation on the system is presented in Fig. 2 where we show the density of the system along the direction of the perturbation. More specifically, the left column corresponds to A = 0.02, which falls well into the LRT regime, cf. Fig. 1. Consequently, the red curve that has been obtained from Eq. (7) is in excellent agreement to the green crosses depicting the PIMC data. The relative deviation of LRT to the latter is shown in the bottom panel of the same figure, and we find perfect agreement within the given Monte-Carlo error bars that are of the order of ∆n/n 0.1%. Overall, we find that the density modulation attains an amplitude of almost 5% of the unperturbed density.
Let us next turn to the right column of Fig. 2 showing the same quantity for a larger perturbation amplitude, A = 0.2. In this case, the density modulation is substantially larger and the amplitudes of n(x) are of the order of 50% with respect to the unperturbed density. Thus, the system is indeed to a large degree shaped by the presence of the external potential. Consequently, the red LRT curve only gives a qualitative description of the density and systematically underestimates (overestimates) the height of the peaks (depth of the minima) of n. This can again be seen most clearly in the deviation plot in the bottom panel, where we find differences between LRT and the PIMC data exceeding 5%.
As a first step towards an improved description of the density profile n(x), we go beyond LRT by including the cubic contribution χ (1,cubic) to the first harmonic. Indeed, Fig. 1 clearly indicates that this provides a fully adequate description of ρ q q,A , and the dash-dotted black curve in Fig. 2 corresponds to Eq. (20) truncated after η = 1. Let us list a few observations: i) the impact of χ (1,cubic) on n(x) is relatively small despite the large density modulation amplitude and the inaccurate description of LRT; ii) from the top panel, we see that χ (1,cubic) leads to a somewhat improved agreement with the PIMC data around the minima, but increased deviations around the peaks. This becomes more clear by looking at ∆n/n 0 in the bottom panel. The LRT curve exhibits a deviation profile that oscillates twice as fast in space as the original density modulation, but the amplitude of this deviation is not constant. In fact, this deviation profile is a combination of a) the insufficient description of ρ q q,A by LRT (cf. Fig. 1) and b) the omission of the second harmonic that does, by definition, oscillate twice as fast in space as the first one. By including χ (1,cubic) in the description of ρ q q,A , we have removed effect a), and the residual deviation profile (black squares in the bottom panel of Fig. 2) does indeed have a constant amplitude and corresponds to the contribution of the second harmonic.
Evaluating Eq. (20) up to η = 2 gives the dashed blue curve, which is in striking agreement with the green crosses, over the entire x-range. This is again confirmed by the deviation plot in the bottom, where the corresponding blue diamonds fluctuate around zero deviation with a relative accuracy of ∆n ∼ 0.1%.
A more systematic investigation of the emerging impact of higher harmonics upon increasing the perturbation amplitude is presented in Fig. 3, where we show the full wave-number dependence of ρ k q,A for different values of A. The left panel shows results for q ≈ 1.69q F , i.e., the same perturbation wave number as in Figs. 1 and 2. First and foremost, we empirically confirm that a non-zero value for the density response can indeed only be found at the integer harmonics of the original perturbation wave vector q. For the smallest depicted value of the perturbation amplitude, A = 0.02 (black squares), we only find a signal at q itself, which agrees with the prediction known from LRT. Increasing the perturbation strength by a factor of five (A = 0.1) leads to the results given by the red circles. For the first harmonic, the response remains almost unchanged and varies only by ∼ 1% (cf. the bottom panel of Fig. 1). At the same time, we find a significant response at the second harmonic that cannot be neglected and substantially contributes to observables like the density profile in coordinate space n(r) that was discussed above.
A further increase of A by a factor of two is depicted by the green crosses. In this case, the response at the first harmonic clearly deviates from the prediction by LRT, and we have already seen that the cubic contribution χ (1,cubic) (q) is required for an adequate description. Further, we observe an even more pronounced deviation from LRT at the second harmonic, which was shown to be responsible for the bulk of the deviation between LRT and the PIMC data for n(r) in Fig. 2. At the same time, no significant response is observed for higher harmonics (η > 2) at these conditions. This changes for the largest depicted perturbation strength shown in Fig. 3, i.e., A = 0.7 (blue diamonds). In this case, the signal at the first harmonic is substantially reduced compared to the other data points, and even the cubic response function χ (1,cubic) (q) cannot provide a reasonable description of the density response at these conditions, see Fig. 1. Moreover, we find a large signal at the second harmonic, as it is expected, and a significant contribution for the third harmonic, k = 3q.
The corresponding density profile in coordinate space is depicted in Fig. 4, with the green crosses again being the exact PIMC data. Evidently, the external perturbation dominates the behaviour of the system under these conditions and we find an almost shell-like structure with two pronounced peaks at the positions of the minima of the harmonic potential exceeding 2.5 times the original unperturbed density n 0 . These shells are separated by two deep minima where the cosinusoidal potential has its maxima.
Let us next consider the solid red line showing the prediction by LRT, i.e., Eq. (7). Evidently, the deviation to the PIMC data is dramatic and almost attains 80% in terms of the unperturbed density n 0 . In particular, LRT predicts an unphysical negative density around the minima of the density profile. Furthermore, we find that the deviation profile (bottom panel) is quickly oscillating and exhibits a complicated behaviour, which means that the systematic deficiencies of LRT are both quantitative and qualitative. The next step towards an improved description of n(x) is given by the inclusion of the cubic response function χ (1,cubic) (q) into ρ q q,A , i.e., ρ q q,A ≈ χ (1) (q)A + χ (1,cubic) (q)A 3 . The resulting evaluation of Eq. (20) truncated at η = 1 is depicted as the dash-double-dotted yellow curve and leads to a qualitative improvement of the description of the exact PIMC data, even though the deviations are still systematic and partly attain ∼ 60% of the unperturbed density n 0 . This is somewhat expected as we know from Fig. 1 that even the first harmonic is not properly described by χ (1) (q) and χ (1,cubic) (q) alone, and more terms in the expansion would be needed. This is remedied by the dash-dotted black curve, where we again truncate Eq. (20) after η = 1, but include the full PIMC expectation value for ρ q q,A . On the one hand, this does not substantially improve the agreement with the exact PIMC data and there remain systematic deviation of ∼ 40%. On the other hand, it makes the deviation profile shown in the bottom panel more uniform. In particular, it exhibits oscillations with twice the wave number q and a nearly constant amplitude.
Naturally, these deviations are caused by the second harmonic, and including the η = 2 term with the exact PIMC expectation value for ρ 2q q,A in Eq. (20) leads to the dotted purple curve. The inclusion of the second harmonic leads to a drastic improvement, and the purple curve is in qualitative agreement with the green crosses everywhere. In particular, we note that the unphysical negative density points no longer appear. Still, the examination of the deviation profile of this curve shown in the bottom panel reveals that there remain systematic errors exceeding 5%. Further, the deviation has a rapidly oscillating form with a constant amplitude, and is readily identified as the contribution of the third harmonic, k = 3q.
In other words, it is necessary to evaluate Eq. (20) up to η = 3 for such a large perturbation amplitude, and the results are shown as the dashed blue curve. Evidently, it is in excellent agreement with the PIMC data over the entire x-range, and the deviation plot shows that the difference between the two is fluctuating around zero within the given statistical uncertainty.
Let us next consider the right panel of Fig. 3, where we show the k-dependence of the density response ρ k q,A for a smaller wave number, q ≈ 0.84q F . In particular, this is the smallest wave number that is accessible within a PIMC simulation for N = 14, see also Sec. II B above.
For the lowest depicted value of the perturbation amplitude (A = 0.05, black squares), the signal at the first harmonic k = q cannot be distinguished from LRT [63]. In addition, we find a small yet significant signal for k = 2q, which indicates that the second harmonic is the dominant contribution to nonlinear effects in the density response of electrons in the WDM regime. This is an important finding that is discussed extensively in the context of Fig. 11 below.
For A = 0.15 (red circles), LRT still relatively accu-rately describes the signal at the first harmonic, and we find a deviation to the PIMC data of less than 1%. At the same time, the signal at k = 2q is significantly increased, which further corroborates the previous observation about the respective importance of the different harmonics to the nonlinear density response. Upon further increasing the perturbation strength to A = 0.3 (green crosses), the systematic error of LRT in the description of ρ q q,A increases to ∼ 3%, while the second harmonic continues to have a higher impact. Moreover, we observe a small but significant density response at the third harmonic, that was absent for A 0.15.
The blue diamonds show ρ k q,A for A = 0.7, which was the largest perturbation strength considered for q ≈ 1.69q F shown in the left panel of Fig. 3. For the smaller q, we observe a large systematic error of LRT for the first harmonic, which, however, is very accurately described by the cubic response function χ (1,cubic) (q). This is a direct consequence of screening effects, which make nonlinear effects less important for small wave numbers, see Ref. [63] for the first description of this finding. In the present work, we find similar effects for the second and third harmonics as well, see below. Further, we find a large signal for k = 2q and a smaller, but noticeable signal for k = 3q.
Finally, we consider the case of strong perturbation strength (A = 1.5), which is depicted by the yellow triangles. At these parameters, the density response at the first harmonic is substantially reduced compared to LRT, and the signal at k = 2q is of the same order of magnitude. Further, the response at the third harmonic is dras- tically increased compared to A = 0.7, and we even find a significant response for the fourth harmonic, k = 4q.
Let us proceed with a more rigorous quantification of the nonlinear electronic density response by obtaining the respective generalized response functions defined in Eqs. (36)- (38). For the case of the cubic response function of the first harmonic χ (1,cubic) (q), this is demonstrated in Fig. 5 for three different wave numbers. More specifically, we show PIMC results for the dependence of ρ q q,A on the perturbation strength for q ≈ 0.84q F (green crosses), q ≈ 1.69q F (blue diamonds), and q ≈ 2.53q F (red circles), and the corresponding curves show the respective fits according to Eq. (36). The bottom panel shows the deviation between fits and PIMC data and the vertical lines indicate the respective maximum values of A up to which data have been included in the fitting procedure.
Let us briefly touch upon the A-dependence of the Determination of the cubic response function χ (1,cubic) (q). Top panel: density response of the UEG for N = 14, rs = 2, and θ = 1 for q ≈ 0.84qF (green crosses) q ≈ 1.69qF (blue diamonds), and q ≈ 2.53qF (red circles). The corresponding curves of the same colour have been obtained by fitting Eq. (36) to the PIMC data. Bottom panel: relative deviation of the fits to the PIMC data. The vertical lines show the respective maximum value of A that was included in the fits.
Monte-Carlo error bars in the deviation plot. For weak perturbations A, the absolute values of the density response ρ q q,A are small, which means that the denominator of the relative deviations shown in the bottom of Fig. 5 is small, too. At the same time, the absolute values of the Monte-Carlo error bars are nearly independent of A, which means that the relative accuracy of our PIMC data for Eq. (3) decreases for weak perturbations, and, eventually, the density response cannot be resolved within the given statistical uncertainty.
Regarding the A-dependence of ρ q q,A itself, we find that the blue diamonds exhibit the largest density response at the first harmonic, whereas the signal is weakest for the largest wave number, q = 2.53q F . While not being trivial, this is somewhat expected as the LRT function χ (1) (q), too, exhibits a maximum response around q ∼ 1.5q F , cf. Wave-number dependence of the cubic response function of the first harmonic χ (1,cubic) (q) for the warm dense electron gas at rs = 2 and θ = 1.
mains accurate for larger perturbation amplitudes A for the green crosses compared to the other two cases. This can be explained by the comparatively smaller value of χ (1,cubic) in this case, see Fig. 6.
Repeating the fitting procedure demonstrated in Fig. 5 for different values of the perturbation wave numbers gives us access to the full q-dependence of χ (1,cubic) (q), and the results are shown in Fig. 6. More specifically, the green data points have been obtained from our PIMC simulations and the stars and diamonds show results for N = 14 and N = 20 particles. While being available at different discrete q-points [41,110], we note that no dependence on the system size can be resolved within the given confidence intervals [63]. In addition, the grey stars, too, have been obtained from PIMC data, but for N = 14 ideal (noninteracting) fermions. In this case, there are no screening effects and χ (1,cubic) 0 (q) attains a finite value in the limit of q → 0.
Unfortunately, no simple recursion relation is available for the ideal function χ (1,cubic) 0 (q) [see Sec. II above], so that we cannot evaluate the RPA and LFC expressions for the cubic response at the first harmonic given in Eqs. (31) and (35) for all wave numbers. Yet, we can still use the three PIMC data points, and the results are depicted by the blue (RPA) and red (LFC) data points in Fig. 6. Evidently, the RPA curve seems to qualitatively reproduce the correct behaviour, but severely underestimates the true magnitude of the cubic response around its maximum. In contrast, the LFC curve nicely agrees with the PIMC data points for the interacting systems for all three wave numbers.
Let us next repeat this analysis for the quadratic response function of the second harmonic χ (2) (q), that has been shown to constitute the dominant nonlinear effect above. To this end, we show the A-dependence of the FIG. 7. Determination of the quadratic response function χ (2) (q). Top panel: density response of the second harmonic ρ2q q,A of the UEG for N = 14, rs = 2, and θ = 1 for q ≈ 0.84qF (green crosses) q ≈ 1.69qF (blue diamonds), and q ≈ 2.53qF (red circles). The corresponding curves of the same colour have been obtained by fitting Eq. (37) to the PIMC data. Bottom panel: relative deviation of the fits to the PIMC data. The vertical lines show the respective maximum value of A that was included in the fits. density response of the second harmonic ρ 2q q,A in Fig. 7 for the same wave-numbers as in Fig. 5.
Firstly, we observe the same absolute ordering of the three data sets with q as for the first harmonic, with q = 1.69q F (q = 0.84q F ) exhibiting the largest (weakest) signal for the density response. Further, we see that the parabolic expansion given in Eq. (37) is in excellent agreement with the PIMC data points for small q, which gives us direct access to the corresponding quadratic density response function.
The results for the full wave-number dependence of χ (2) (q) are shown in Fig. 8 with the same key as χ (1,cubic) (q) shown in Fig. 6 above. There is excellent agreement between the exact quadratic response function of the ideal Fermi gas [yellow curve, cf. Eq. (21)] and our PIMC data for N = 14 noninteracting fermions (grey stars) over the entire depicted q-range. In addition, we find that the quadratic response function of the UEG χ (2) (q) exhibits a qualitatively similar behaviour to χ (1,cubic) (q) shown above. More specifically, χ (2) (q) van-  ishes for small wave-numbers due to the perfect screening in the UEG, and, too, vanishes in the limit of large q, albeit slower than for q → 0.
The dashed blue RPA curve [Eq. (24)] is in qualitative agreement with our PIMC data for the UEG (green stars and diamonds) and correctly reproduces both of these limits. Yet, there appear systematic deviations exceeding 20% for intermediate wave-numbers q ∼ 1.5q F (i.e., in the vicinity of the maximum), which resembles the known deficiencies of RPA within LRT, see, e.g., Refs [58,111,112].
Finally, the dotted red curve has been obtained by including the neural-net representation of the static LFC G(q) given by Dornheim et al. [58] [Eq. (33)] and is in excellent agreement with the PIMC data for all q.
The last generalized response function to be considered in this work is the cubic response function of the third harmonic χ (3) (q), that can be obtained by fitting Eq. (38) to PIMC data for ρ 3q q,A . This is shown in Fig. 9, where the A-dependence of the density response at k = 3q is shown for the same three wave numbers as for the other response functions investigated above. In this case, we observe a comparatively large increase of ρ 3q q,A with A at q = 1.69q F that is approximately four times as large as for the other two wave-numbers, which are similar in magnitude.
Moreover, we note that the absolute value of the density response at k = 3q is an order of magnitude smaller than for k = 2q for the same values of the perturbation Top panel: density response of the third harmonic ρ3q q,A of the UEG for N = 14, rs = 2, and θ = 1 for q ≈ 0.84qF (green crosses) q ≈ 1.69qF (blue diamonds), and q ≈ 2.53qF (red circles). The corresponding curves of the same colour have been obtained by fitting Eq. (38) to the PIMC data. Bottom panel: relative deviation of the fits to the PIMC data. The vertical lines show the respective maximum value of A that was included in the fits.
amplitude A. This, in turn, means that the relative error bars in ρ 3q q,A are larger than for the other harmonics shown above, which makes the determination of χ (3) (q) more challenging.
The corresponding wave-number dependence of this function is shown in Fig. 10. Let us first focus on the noninteracting case, with the grey stars depicting PIMC results for N = 14 ideal fermions and the solid yellow curve the exact ideal response function χ (3) 0 (q) that we obtain from Eq. (22). Remarkably, even without the electronic Coulomb repulsion and the associated screening effects, the density response function of the third harmonic exhibits a complicated, non-monotonous behaviour, with a maximum in magnitude around the Fermi wave number. This is confirmed by the corresponding ideal PIMC data point at q = 0.84q F , which nicely agrees with the exact result.
For the interacting UEG, this interesting behaviour at small q is masked by screening effects, and the generalized response function exhibits similar trends as χ (1,cubic) (q)  and χ (2) (q), that is, it vanishes both in the limits of small and large q. As seen before, the RPA of the cubic response at the third harmonic underestimates the response, but the inclusion of LFCs gives a nice agreement between QMC simulations and theory.
Let us conclude this section with a systematic investigation of the respective contribution of the different generalized response functions to the total nonlinear density response. In Fig. 11, we show the relative contribution of χ (2) (q) (blue), χ (1,cubic) (q) (green), and χ (3) (q) (black) compared to the density response at the first harmonic predicted by LRT, Eq. (4). More specifically, the data points have been obtained from fits to PIMC data, and the dotted (dashed) lines are theoretical results including an LFC (RPA only).
The top panel of Fig. 11 corresponds to a weak perturbation, A = 0.02. At these conditions, LRT is relatively accurate (see, e.g., the left panel of Fig. 2) and the response due to the cubic response functions for both the first and third harmonic is smaller than 0.1% compared to LRT. By far the largest nonlinear effect can be observed at the second harmonic with a maximum signal of almost 1% in terms of ρ LRT . In particular, the response for k = 2q is more than an order of magnitude larger than the other two depicted curves.
The bottom panel of Fig. 11 corresponds to a stronger perturbation, A = 0.2. This is a particularly interesting example, as the density response for the first three harmonics is fully described by the different response func- tions at these conditions. In contrast, additional terms of the respective expansions in powers of A would be needed for an accurate description of ρ k q,A at much larger values of A. Let us first consider the density response of the second harmonic which, again, constitutes the largest nonlinear effect, with a maximum contribution of almost 10% around q ∼ 1.5q F compared to Eq. (4). Further, we note that the relative impact of ρ 2q q,A only diminishes in the small wave-number limit for q 0.2q F , and disappears even more slowly for large values of q, with ρ 2q q,A /ρ LRT exceeding 1% even at q = 5q F .
In addition, the cubic contribution to the density response at the first harmonic ∆ρ q q,A = χ (1,cubic) (q)A 3 is substantially more important compared to A = 0.02, and we find a maximum impact of approximately 4% around Density profile of the harmonically perturbed electron gas for N = 14 and rs = 2, with q ≈ 0.84qF and q = 2π/L(1, 0, 0) T for A = 0.1. Top: PIMC results for the density n(x) along the direction of the perturbation for θ = 1 (blue diamonds), θ = 2 (red circles), and θ = 4 (green crosses). The solid curves correspond to LRT [cf. Eq. (7)]. Bottom: relative deviation of LRT to the PIMC data. q = 1.5q F . Indeed, this is the same order of magnitude as the signal of the second harmonic, which means that both χ (2) (q) and χ (1,cubic) (q) have to be taken into account for an accurate description of the system, cf. the right panel of Fig. 2 above showing the radial density n(r) for this case. At the same time, we find that the impact of the cubic contribution to the first harmonic rapidly vanishes in both the limits of large and small wave-numbers and is only relevant for intermediate q.
Finally, the black points and curves correspond to the contribution of the third harmonic, which is an order of magnitude smaller than the green curve. Therefore, it appears to be mostly negligible for practical purposes.

B. Temperature dependence
Let us next turn our attention to the dependence of nonlinear effects in the static density response on the temperature θ. In Fig. 12, we show the density profile in coordinate space, n(r), along the direction of the perturbation, for N = 14 electrons at r s = 2 for q ≈ 0.84q F and a relatively small value of the perturbation amplitude, A = 0.1. The blue diamonds have been obtained for θ = 1 and exhibit a deviation from the unperturbed density n 0 of almost 20% at the maxima and minima. The solid blue curve shows the prediction from LRT [cf. Eq. (7)] and qualitatively captures the correct behaviour. This can again be seen particularly well in the bottom panel where we show the relative deviation between the PIMC data and LRT. Evidently, LRT results in an oscillating deviation profile that is mainly due to the density response at the second harmonic ρ 2q q,A (see the previous section III A), with a maximum deviation below 0.8%.
Let us next focus on the red circles corresponding to twice the temperature, θ = 2. Firstly, we note that the overall density response is reduced compared to θ = 1, as it is expected, see also the corresponding discussion in Refs. [63,111]. In particular, an increasing temperature leads to an increased internal energy scale of the system, which, in turn, means that the relative impact of an external potential is reduced if the perturbation amplitude A is kept constant. The green crosses correspond to a further increase of the temperature by a factor of two, and the total density response is again reduced compared to the other data sets.
In addition, a reduced density response means that nonlinear effects become less important and we find maximum deviations between LRT and the PIMC data of 0.4% (0.2%) for θ = 2 (θ = 4). At the same time, we note that the nonlinear response is significant for all three depicted temperatures and is mainly given by the signal at the second harmonic, k = 2q.
Let us next consider Fig. 13, which shows the same information as Fig. 12, but for a significantly increased perturbation strength A = 0.7. In this case, the density response is strongly nonlinear for all three depicted temperatures, although the density profile n(r) is substantially different for all three cases. For θ = 1, the external perturbation constitutes the dominant effect and fully shapes the physical behaviour of the system. Most particles are, on average, located around the center of the simulation cell where the cosinusoidal potential has its minimum, whereas the density is reduced by roughly 90% for small and large x. The solid blue curve has been computed by truncating Eq. (20) after the second harmonic (η = 2) using the exact PIMC expectation values for the coefficients ρ q q,A and ρ 2q q,A . In this case, we find that the agreement to the PIMC data is only very qualitative, and there appear rapidly oscillating deviations (bottom panel) with an amplitude exceeding 3%. Naturally, these deviations constitute the signal of the third harmonic, k = 3q, see also the discussion of Fig. 14 below.
For θ = 2, the density at the maxima of the external perturbation is substantially larger compared to the lower temperature, and the impact of the third harmonic on the density profile n(r) does not exceed 2%. Finally at the highest temperature, we find oscillations in the total density n of approximately 50%, which is substantially smaller compared to the other cases. Consequently, the first two harmonics already provide a rather accurate description, and the systematic error is smaller than 1% everywhere.
A more systematic investigation of the respective strength of the density response at different harmonics is shown in Fig. 14, where we show the full k-dependence of ρ k q,A for the two examples investigated in Figs. 12 and 13. The top panel corresponds to the weaker perturbation, A = 0.1. We note that we again only observe a nonzero signal for the integer harmonics of q, as it is expected. For the first harmonic (k = q), we see the expected order of the signals for the three depicted temperatures, with the response at θ = 4 being less than half the response at θ = 1.
For the second harmonic (k = 2q), the relative strength of the signals for different temperatures is even FIG. 14. Full wave-number dependence of the density response ρ k q,A of the UEG for N = 14, rs = 2, and θ = 1 with q ≈ 0.84qF and A = 0.1 (top) and A = 0.7 (bottom). Shown are PIMC results for Eq. (3) for θ = 1 (blue diamonds), θ = 2 (red circles), and θ = 4 (green crosses). The vertical dotted lines indicate the location of the first four harmonics. more different, with the response for θ = 4 (θ = 2) being approximately 25% (50%) of the response for θ = 1.
Finally, we mention that the signal at the third harmonic (k = 3q) vanishes within the given level of accuracy for all θ.
The bottom panel of Fig. 14 shows the same results for A = 0.7 and we find similar trends. In fact, the ratio of density response for θ = 2 and θ = 4 for both the first and second harmonic compared to the respective response at θ = 1 is remarkably unchanged, although the total response at k = 2q has substantially increased. Moreover, we find a nonzero signal at k = 3q for all three temperatures, whereas the response at the fourth harmonic vanishes even for θ = 1.
The next step towards a more complete description of the temperature dependence of the nonlinear density response is given by the determination of the different generalized response functions. To this end, we show the dependence of the density response at the second harmonic ρ 2q q,A on the perturbation amplitude A in Fig. 15 for the same conditions as in the previous figures in this section. The points depict PIMC data for the three different temperatures, and the curves parabolic fits to these data according to Eq. (37). The bottom panel shows the relative deviation between PIMC and the fits, and the vertical bars show the respective maximum value of A up to which points have been included into the fitting procedure. Evidently, Eq. (37) is capable to accurately reproduce the PIMC data for A 0.4 in all three cases. Going back to the top panel, we see that both the PIMC results and the parabolic fits are systematically FIG. 16. Wave-number dependence of the quadratic response function of the second harmonic χ (2) (q) for the warm dense electron gas at rs = 2 and θ = 1 (blue diamonds), θ = 2 (red circles), and θ = 4 (green crosses). The dashed and solid lines show theoretical results within RPA and using the LFC from Ref. [58], cf. Eqs. (24) and (32). ordered with θ, as it is expected.
Repeating this analysis for different values of the perturbation wave-number q gives us access to the full wavenumber dependence of the quadratic response function of the second harmonic χ (2) (q), and the results are shown in Fig. 16 for different relevant temperatures. The symbols have been obtained on the basis of fitting Eq. (37) to PIMC data, and the solid [dashed] curves correspond to the theory with a LFC [within RPA], cf. Eq. (33) [Eq. (24)].
First and foremost, we note that χ (2) exhibits a substantial dependence on θ, and the maximum value of this function at θ = 0.01 (black) is larger by an order of magnitude compared to θ = 4 (green). Moreover, the shape of the curves strongly depends on θ as well: while the curves for all temperatures are relatively similar both in the limit of small and large wave-numbers q, there emerges a sharp peak around q ∼ 1.3q F with decreasing θ.
Let us next analyse the level of accuracy of both RPA and LFC, which, too, strongly depends on θ. For θ = 4, the system is only weakly coupled and both RPA and LFC are in very good agreement to the PIMC data (green crosses) over the entire depicted wave-number range. Upon decreasing the temperature to θ = 2, RPA systematically underestimates the quadratic response function in the vicinity of its maximum, whereas the LFC still nicely reproduces the PIMC results (red circles) everywhere. For θ = 1, the deviation between RPA and LFC becomes substantial and partly exceeds 20%, whereas the LFC gives both a qualitatively and quantitatively sufficient description. At the same time, we like to re-iterate the point that Eq. (33) does only include screening effects to a linear order, and, thus does not give an exact description even with the exact static limit of the LFC G(q) taken from Ref. [58]. Indeed, the systematic inaccuracy of the approximate description of screening effects in Eq. (33) increases with the effective coupling strength of the system. With the density parameter r s being constant, the effective coupling increases towards small temperature and small deviations between the PIMC data points and the theoretical prediction become noticeable for θ 1.
Finally, the black curves show results for θ = 0.01, where the system is in the ground state. In this case, the systematic deviation between LFC and RPA is even larger compared to θ = 1, as the effective coupling strength of the system is further increased. Unfortunately, PIMC simulations are not feasible at these conditions due to the fermion sign problem [68], and groundstate QMC simulations are not available for this quantity.
As a reference, we show the analogous temperature dependence of the LRT response function χ (1) (q) in Fig. 17. Overall, χ (1) (q) exhibits a qualitatively similar behaviour as χ (2) (q), with a screening-induced decay at small wave-numbers, and a maximum of its modulus around q ∼ 1.5q F . Still, the decrease of the density response towards high temperature is substantially less Determination of the cubic response function χ (1,cubic) (q). Top panel: density response of the first harmonic ρq q,A of the UEG for N = 14, rs = 2, and q ≈ 0.84qF for θ = 1 (blue diamonds) θ = 2 (red circles), and θ = 4 (green crosses). The corresponding curves of the same color have been obtained by fitting Eq. (36) to the PIMC data. Bottom panel: relative deviation of the fits to the PIMC data. The vertical lines show the respective maximum value of A that was included in the fits.
pronounced for χ (1) (q), which makes the experimental probing of the second harmonic a potentially valuable method of diagnostic for the temperature [45,46]. The next generalized response function of interest is given by χ (1,cubic) (q), the determination of which is demonstrated in Fig. 18 for the usual three values of θ. The data points are PIMC results for ρ q q,A and the different curves have been fitted to the latter via Eq. (36). Evidently, this functional form is again well suited to represent the PIMC data, and deviations can only be spotted with the naked eye at large perturbation amplitudes and are most pronounced for θ = 1. In the bottom panel, we show the corresponding relative deviations ∆ρ/ρ, with the vertical lines representing the maximum values of the perturbation amplitude A that have been included into the fitting procedure. In particular, the fit function remains accurate for large values of A beyond the respective fitting interval for all three temperatures, which is a strong indication for the high quality of our FIG. 19. Wave-number dependence of the cubic response function of the first harmonic χ (1,cubic) (q) for the warm dense electron gas at rs = 2 and θ = 1 (blue diamonds), θ = 2 (red circles), and θ = 4 (green crosses).

approach.
The full wave-number dependence of the thus determined cubic response function of the first harmonic χ (1,cubic) (q) is shown in Fig. 19. First and foremost, we find that this function exhibits an even more pronounced dependence on the temperature compared to both χ (1) (q) and χ (2) (q), which further highlights the potential utility of the nonlinear density response as a method of diagnostics for WDM. More specifically, the magnitude of the maximum at θ = 4 is more than order of magnitude smaller compared to its analogue at θ = 1. Unfortunately, no RPA and LFC data are available in this case, as the ideal cubic response function is not known analytically and, hence, would require extensive PIMC simulations of the harmonically perturbed ideal Fermi gas.
The final response function considered in this work is given by the cubic response function of the third harmonic, χ (3) (q). The determination of this function on the basis of our PIMC data for ρ 3q q,A is shown in Fig. 20, where the curves have been fitted to the latter via Eq. (38). Since the observed signal for the third harmonic is comparably quite small for all three temperatures (cf. Fig. 14), the relative statistical uncertainty of the PIMC data points is large, especially at small A. Yet, Eq. (38) only constitutes the lowest order in A of ρ 3q q,A , which makes the determination of the coefficient χ (3) (q) difficult. This is particularly true for θ = 4, where we can only give a qualitative estimation of the order of magnitude of this function.
The full wave-number dependence of Determination of the cubic response function χ (3) (q). Top panel: density response of the third harmonic ρ3q q,A of the UEG for N = 14, rs = 2, and q ≈ 0.84qF for θ = 1 (blue diamonds) θ = 2 (red circles), and θ = 4 (green crosses). The corresponding curves of the same colour have been obtained by fitting Eq. (37) to the PIMC data. Bottom panel: relative deviation of the fits from the PIMC data. The vertical lines show the respective maximum value of A that was included in the fits.
in Fig. 21, again for a gamut of different values of θ. For the highest temperature (θ = 4, green curve), the signal is two orders of magnitude smaller than in the groundstate (θ = 0.01, black curve) and can hardly be resolved with the bare eye. Upon decreasing the temperature to θ = 2 (red curve), the signal is increased by a factor of five and can be resolved with the PIMC procedure described in the discussion of Fig. 20. Moreover, the difference between RPA [Eq. (30)] and the LFC-based description [Eq. (34)] is quite small, and the PIMC data points are in excellent agreement with the latter.
Decreasing the temperature by an additional factor of two leads to the blue curves and data points, which display approximately thrice the magnitude in the maximum of the density response around q ∼ 1.5q F . In addition, there appears a much larger difference between RPA and LFC, and-although to a much smaller degreebetween the LFC curve and the PIMC data; see also the discussion of Fig. 6 above.
Lastly, the black curves correspond to the ground state (θ = 0.01), where no PIMC simulations are available due to the fermion sign problem [68]. Evidently, RPA underestimates the depth of the minimum of χ (3) (q) around q ≈ 1.5q F by almost 50% and, thus, only gives a very qualitative description of the density response at the third harmonic of the original perturbation. In addition, we find a complicated, nontrivial behaviour as χ (3) (q) changes its sign around the Fermi wave number. We expect this to be a real physical effect, as a quite similar behaviour has been detected on the basis of our PIMC data at a lower density, r s = 6, shown in Fig. 22 below.
C. Strong coupling: rs = 6 Let us next investigate the nonlinear density response of the UEG at stronger coupling, r s = 6 and θ = 1. Such conditions can be realized experimentally for example in evaporation experiments [113][114][115][116] and constitute a challenging benchmark for models and theories due to the pronounced impact of electronic exchange-correlation effects [61,115].
Since the determination of the different generalized response functions from ab initio PIMC data for ρ k q,A is completely analogous to the case of r s = 2 that was discussed in the previous sections, here we restrict ourselves to a concise discussion of the final results summarized in Fig. 22.
The top left panel corresponds to the usual LRT func-tion χ (1) (q), where the green stars have been obtained from a single simulation of the unperturbed UEG by evaluating Eq. (6) for N = 34 unpolarized electrons. These data are in good agreement with the dotted red curve that has been computed from the static LFC presented in Ref. [58], whereas RPA (dashed blue) substantially underestimates the magnitude of the minimum of χ (1) (q) around q = 2q F by more than 25%. This is of course expected and a direct consequence of the increased importance of electronic correlation effects at low density [112]. The solid yellow curve depicts results for the ideal Fermi gas and becomes increasingly inaccurate towards small wave-numbers, when screening effects start to manifest in the UEG [117]. The top right panel shows the wave-number dependence of the quadratic response function of the second harmonic, χ (2) (q). We find that this function exhibits a significantly faster decay towards zero both in the smalland large-q limits, similar to our findings for r s = 2 above. In addition, the electronic LFC G(q) has a substantially larger impact on χ (2) (q) compared to χ (1) (q), and the RPA underestimates the peak at q ≈ 1.75q F by 50%. In contrast, the LFC curve is in good agreement with the PIMC data points for both N = 14 (green diamonds) and N = 34 (green stars) over the entire depicted wave-number range, and the deviation due to the approximate treatment of screening effects [cf. Sec. II] is small. In addition, the PIMC data points for the two system sizes exhibit a smooth progression, which means that no finite-size effects can be resolved within the given level of accuracy. This is in good agreement with the recent investigation of χ (1,cubic) (q) presented in Ref. [63].
Let us next focus on the bottom left panel showing the cubic density response function of the first harmonic, χ (1,cubic) (q) at the same conditions. In particular, this function is sharply peaked around q ≈ 2q F and quickly vanishes both in the limits of small and large wave numbers. We again mention that the evaluation of the RPA and LFC expressions given in Eqs. (31) and (35) is not possible as χ (1,cubic) 0 (q) cannot be readily evaluated. Finally, we investigate the cubic density response function of the third harmonic, χ (3) (q) shown in the bottom right panel of Fig. 22. While the associated density response at the third harmonic is relatively small and unimportant compared to the signals at k = q and k = 2q, this function still deserves close attention as it exhibits interesting, nontrivial behaviour. More specifically, the PIMC data points are positive for small wave numbers and χ (2) (q) changes its sign around the Fermi wave number, q = q F . This trend is further substantiated in Fig. 23, where we show the dependence of the actual density response at the third harmonic ρ 3q q,A on the perturbation amplitude A. The green crosses show results for the smallest wave number accessible for N = 34 electrons, q ≈ 0.62q F , and we can clearly resolve a positive density response that monotonically increases with A. In contrast, the blue diamonds and red circles corresponding to two larger wave FIG. 22. Wave-number dependence of the linear response function χ (1) (q) (top left), the quadratic response function of the second harmonic χ (2) (q) (top right), the cubic response function of the first harmonic χ (1,cubic) (q) (bottom left), and the cubic response function of the third harmonic χ (3) (q) (bottom right). The green stars and diamonds correspond to PIMC data points for N = 34 and N = 14, the solid yellow curve to the ideal Fermi gas, the dashed blue curves to RPA, and the dotted red curves have been obtained using the LFC from Ref. [58]. All data have been obtained for rs = 6 and θ = 1.
numbers (q ≈ 1.25q F and q ≈ 2.51q F , respectively) both exhibit a negative density response. As usual, the different curves have been obtained by fitting Eq. (38) to the PIMC data, and the bottom panel of Fig. 23 shows the relative deviation. Interestingly, the fit function remains accurate for larger A at both q = 0.62q F and q = 1.25q F than for q = 2.51q F , even though the response is substantially larger at the intermediate wave number (see also the vertical lines in the bottom panel showing the maximum value of A that has been included into the fitting procedure).
Let us get back to the full wave-number dependence of χ (3) (q) shown in the bottom right panel of Fig. 22. Evidently, neglecting electronic exchange-correlation effects within the RPA (dashed blue curve) results in a very poor description of the cubic response function, and the real depth of the minimum is underestimated by a factor of two. Including the static LFC gives the dotted red curve that is in very good agreement with the PIMC data points for q q F . Yet, the LFC result, Eq. (34), is not capable to describe the sign change in χ (3) (q) around the Fermi wave number and, thus, stays negative over the entire q-range. This, in turn, means that the approximate linear treatment of screening effects [cf. Sec. II above] leads to a description where this sign change is missed.
Finally, the dotted black curve has also been obtained by evaluating Eq. (34) including the LFC, but for θ = 0.01, i.e, in the ground state. While this leads to a curve with a substantially deeper minimum around q ∼ 1.5q F than at the Fermi temperature, it does exhibit the sign change in χ (3) (q) that we find in our PIMC data at the higher temperature.

IV. SUMMARY AND DISCUSSION
In this work, we have presented an in-depth analysis of nonlinear effects in the density response of the warm Determination of the cubic response function χ (3) (q). Top panel: density response of the third harmonic ρ3q q,A of the UEG for N = 34, rs = 6, and θ = 1 for q ≈ 0.62qF (green crosses) q ≈ 1.25qF (blue diamonds), and q ≈ 2.51qF (red circles). The corresponding curves of the same colour have been obtained by fitting Eq. (38) to the PIMC data. Bottom panel: relative deviation of the fits to the PIMC data. The vertical lines show the respective maximum value of A that was included in the fits. dense electron gas, substantially extending the first investigation of Ref. [63]. We have obtained extensive new PIMC results carrying out simulations of a harmonically perturbed electron gas, which has allowed us to obtain the full spectrum of excitations at the integer harmonics of the original perturbation. First and foremost, we have found that the dominant nonlinear contribution is given by the quadratic response at the second harmonic, for weak to moderate values of the perturbation amplitude A. The second potentially important nonlinear term is given by the cubic response at the first harmonic studied in Ref. [63], whereas the cubic response at the third harmonic is practically negligible in most realistic situations. In addition, we have found that the nonlinear response functions more strongly depend on system parameters such as the temperature compared to the usual LRT, which makes them a potentially valuable new tool of diagnostics for WDM experiments. In addition, the nonlinear response is strongly shaped by electronic exchangecorrelation effects, which makes a reliable theory based on ab initio PIMC data indispensable. Furthermore, we have extended previous works on the nonlinear density response of warm dense electrons and presented new expressions for the response functions in terms of the wellknown LFC that is readily available based on previous PIMC simulations [58,61,62]. In particular, this approach is capable to reproduce the benchmark PIMC data for the nonlinear density response with negligible computational cost.
While the present analysis was based on simulations for a purely static density perturbation, even stronger nonlinear effects are to be expected for dynamic excitations. For example, the excitation of warm dense matter by a strong laser field will give rise to a time-dependent density response that will involve harmonics. If proper resonance conditions are fulfilled, as e.g. in the case of high harmonics generation in rare gases [9], or laser plasmas [118], the excitation of even higher harmonic orders should be expected.
We are convinced that the findings of our paper will open up many new avenues of future research in multiple directions. For example, a realistic perturbation in a WDM experiment will be given by a superposition of many distinct harmonic perturbations. While this has no profound consequence within LRT, a nonlinear treatment of the density response will give rise to mode-coupling between different perturbations, which will make the excitation spectrum more interesting and complicated. One can imagine two-color x-ray beam experiments to intentionally use the mode coupling for diagnostics. A second topic for future research is given by the investigation of the dynamic density response, which can be done on the basis of our new theoretical expressions and a suitable local field correction such as the ESA presented in Refs. [61,62]. In addition to these basic questions, we stress that an accurate theory of nonlinear effects can be directly used for the improved computation of different material properties such as the electronic stopping power, screened potentials, or energy relaxation rates.