Ab initio results for the plasmon dispersion and damping of the warm dense electron gas

Warm dense matter (WDM) is an exotic state on the border between condensed matter and dense plasmas. Important occurrences of WDM include dense astrophysical objects, matter in the core of our Earth, as well as matter produced in strong compression experiments. As of late, x-ray Thomson scattering has become an advanced tool to diagnose WDM. The interpretation of the data requires model input for the dynamic structure factor $S(q,\omega)$ and the plasmon dispersion $\omega(q)$. Recently the first \textit{ab initio} results for $S(q,\omega)$ of the homogeneous warm dense electron gas were obtained from path integral Monte Carlo simulations, [Dornheim \textit{et al.}, Phys. Rev. Lett. \textbf{121}, 255001 (2018)]. Here, we analyse the effects of correlations and finite temperature on the dynamic dielectric function and the plasmon dispersion. Our results for the plasmon dispersion and damping differ significantly from the random phase approximation and from earlier models of the correlated electron gas. Moreover, we show when commonly used weak damping approximations break down and how the method of complex zeros of the dielectric function can solve this problem for WDM conditions.


INTRODUCTION
There is growing interest in warm dense matter (WDM) [1,2] -an extreme state that occurs, e.g., in astrophysical objects [3][4][5][6] , in the core of our Earth [7] , in laser compression experiments [8] , and on the pathway towards inertial confinement fusion [9,10] . WDM is a complicated state due to an intricate interplay of many effects including quantum degeneracy and exchange of the electrons at finite temperatures, electronic and ionic correlations, including wide angle scattering, phase transitions, and partial ionization. The simultaneous occurrence of these effects makes the experimental and theoretic analysis of WDM extremely challenging. Among the experimental diagnostics, x-ray Thomson scattering (XRTS) has been established as an accurate and highly promising tool [11] . XRTS measures the dynamic structure factor of all the electrons in the system and yields information, among others, on the plasmon spectrum, density, temperature, and chemical composition. In Eq. (1), = 1∕ , is the density,̃ ( ) = 4 2 ∕ 2 is the Fourier transform of the Coulomb potential, and ( , ) denotes the dynamic dielectric function. In order to use XRTS as diagnostics, input from models for the dynamic structure factor, i.e., dynamic arXiv:2008.04605v1 [physics.plasm-ph] 11 Aug 2020 dielectric function ( , ) is required. Thus, the accuracy of the diagnostics of WDM crucially depends on the quality of the available models.
This concerns, in particular, the dispersion ( ) and damping ( ) of collective excitations. There has been extensive recent theoretical work on this subject which includes chemical models, sum rule models [12,13] , local field correction (LFC) theory [14,15] , time-dependent density functional theory (TDDFT) [16] , quantum kinetic theory [17,18] , non-equilibrium Green's functions [19][20][21] , and quantum hydrodynamics [22][23][24][25] . Standard assumptions that are used include the Chihara decomposition in multi-component plasmas (this incorporates the Born-Oppenheimer approximation and leads to the description of the response of the free electrons by an electron gas model) [26,27] , the decoupling of longitudinal and transverse modes, and the weak damping approximation, e.g. Ref. [28] . Due to the complexity of WDM, the influence of each of these approximations is often difficult to quantify and reliable predictive models are still missing.
It is, therefore, useful to disentangle these effects by using simpler but well-defined model systems. One such system is the uniform electron gas (UEG) at finite temperature (jellium). This is an example of an one-component system that constitutes an important test case for theory, as it allows one to focus on the treatment of quantum, correlation and finite temperature effects of the electrons and to benchmark models against first principle simulations. The plasmon dispersion of the UEG has been studied for many decades starting with the works of Bohm, Gross, Pines and Ferrell for metals [29][30][31] who developed a quantum mean field theory within the random phase approximation (RPA); for useful parametrizations see Ref [32] . Correlation effects were taken into account via local field corrections, e.g. in Refs. [33][34][35][36] . Recent applications to WDM include the analysis of experiments on beryllium and boron [14,37,38] .
However, the accuracy of these model results for the plasmon dispersion and damping, in particular under WDM conditions, is not known. Therefore, it is highly desirable to develop simulations that avoid approximations regarding correlation, quantum and finite temperature effects. The most accurate approach available is path integral Monte Carlo (PIMC) which, however, also faces fundamental difficulties such as the fermion sign problem [39][40][41] and the complicated access to frequency dependent observables. These problems were overcome by the combination of configuration PIMC (CPIMC) [42,43] , that is exact at strong degeneracy but exhibits a sign problem at weak degeneracy (low density), and of permutation blocking PIMC (PB-PIMC) [44][45][46] . Thus ab initio simulations of jellium now cover the entire range of WDM parameters [47,48] given by 0.1 ≲ Θ, ≲ 10, where Θ = ∕ and =̄ ∕ , with and denoting the Fermi energy and the Bohr radius. Finally, an accurate extrapolation to the thermodynamic limit was achieved in Ref. [49] and the connection to the ground state was realized in Ref. [50] , where also an accurate parametrization for the exchange-correlation free energy was reported. For an overview on the results and comparisons with earlier models and simulations, see Ref. [2,51] .
PIMC simulations can also be used to obtain dynamic quantities such as the dynamic structure factor ( , ) from an analytical continuation of the intermediate scattering function (density correlation function) ( , ) evaluated at imaginary times ∈ [0, ]. The scattering function is related to the dynamic structure factor ( , ) by a Laplace transform This is known to be an ill-posed problem that has occasionally been tackled using maximum entropy methods, see e.g. Ref. [52] and references therein. Recently it was found that a stochastic sampling of the LFCs ( , ) allows to very well reconstruct the imaginary time density response function and thus the dynamic structure factor, because additional exact constrains on the LFC makes the procedure very efficient and accurate [53] . Extensive further studies were reported in Ref. [54] . It was also noted in Ref. [53] , that in many cases the static LFC, ( ) = ( , = 0) is sufficiently accurate to recover the dynamic structure factor. For this purpose, a neural net representation for ( ) was constructed in Ref. [55] that is based on ab initio simulation data. Access to ( ) or even ( , ) allows for systematic extensions of the QMC based ab initio approach to other dynamic quantities such as the density response function ( , ) [24] , the dielectric function, and the dynamic conductivity [56] . For the plasmon dispersion, the key quantity is the dynamic dielectric function ( , ) where Π RPA is the polarization function in random phase approximation (RPA, Lindhard). Setting → 0, recovers the RPA dielectric function. On the other hand, Eq. (3) indicates that, with ab initio input for ( , ) also ab initio results for the dielectric function of correlated electrons under WDM conditions are becoming available.
The goal of this paper is to study the dynamic dielectric function in more detail with a focus on its zeroes because they determine the plasmon dispersion, ( ), and damping, ( ), of correlated electrons. In particular, 1. We present a detailed analysis of the wave number dispersion, ( ), at finite temperature. Starting with RPA, we review various analytical models and find that they exhibit significant deviations from the numerical result. We also present a novel analytical parametrization at finite temperature for ( ) in RPA.
2. We present results for the plasmon dispersion and damping that follow from the PIMC results for the local field correction ( , ) and compare the results to previous studies.
3. We carefully test the validity of the commonly used dispersion relation, Re ( , ) = 0. Since in XRTS experiments under WDM conditions the plasmon damping is not necessarily small, this relation has to be questioned. Therefore, we perform an analytical continuation of the retarded dielectric function to complex frequencies [17] . We present results for the RPA dielectric function and for Eq. (3) and observe significant deviations from the common approach based on the real part of . This has important implications for the correct interpretation of the XRTS measurements of warm dense matter.
The paper is organized as follows: In Sec. 2, we summarize the main ideas of our PIMC approach to the dielectric function that is based on the reconstruction of the dynamic local field correction. There, we also present a discussion of the longitudinal plasmon dispersion and of the analytical continuation. In Sec. 3, we present our ab initio simulation results for the local field correction, the dielectric function, and the plasmon dispersion. We conclude with a summary and outlook in Sec. 4.

Path integral Monte Carlo
The basic idea of the standard path integral Monte Carlo method [57] is to stochastically evaluate the thermal density matrix in coordinate space, with = ( 1 , … , ) containing the coordinates of all particles and = 1∕ B being the ususal inverse temperature. As a direct evaluation of ( , , ) is not possible, one performs a Trotter decomposition [58] , and the final result for the partition function is given as the sum over all closed paths of particle coordinates in the imaginary time ∈ [0, ], see Refs. [41,55] for details. We note that this formulation in the imaginary time is particularly convenient in the context of the present work, as it allows for a straightforward computation of imaginary-time correlation functions, such as the density autocorrelation function All PIMC data presented in this work have been obtained without any nodal restrictions [59] in Eq. (4). Therefore, the simulations are computationally demanding due to the fermion sign problem (see Ref. [41] for a review article), but exact within the given error bars. Moreover, we use a canonical adaption [60] of the worm algorithm introduced by Boninsegni et al. [61,62] .

Stochastic sampling of the dynamic LFC
The numerical inversion of Eq. (2) is a notoriously hard problem [63] . Solutions for ( , ) are, in general, not unique as the information contained in the PIMC data for ( , ) does not fully determine the DSF. To overcome this obstacle, Dornheim and co-workers [53,54,64] have introduced a stochastic sampling scheme for the dynamic local field correction, which automatically satisfies a number of exact constraints on ( , ) and, in this way, sufficiently constraints the space of possible solutions for ( , ).
The basic workflow of this method is as follows [53,54] : 1) Generate a random trial solution for Im ( , ) that already incorporates a number of well-known exact relations. 2) Use the Kramers-Kronig relations [65] to compute the corresponding real part, Re ( , ). (2) and measure the deviation to the PIMC data for ( , ) for all -points. Only those ( , ) which lead to an imaginary-time density-density correlation function in agreement to the PIMC data constitute valid solutions.

LFCs and the dynamic dielectric function
Having obtained ab initio results for the dynamic local field correction, it is straightforward to obtain the dynamic dielectric function via Eq. (3). This can be rewritten in terms of the polarization function Π where In the mean field limit, → 0, and we recover the RPA dielectric function As already mentioned, an important approximation is obtained by replacing ( , ) with its static limit, ( , ) → ( , 0) = ( ) in Eq. (3) or Eq. (7). This is still a dynamic dielectric function which will be denoted SLFC ( , ), while the full dynamic result will be called DLFC ( , ). Comparing results for the dynamic structure factor revealed, that the static approximation is accurate for ≲ 4, for all wave numbers [53] .

Longitudinal plasmon dispersion
The existence of longitudinal collective plasma oscillation follows from Maxwell's equations which, after Fourier transform, reduce to a wave equation for the Fourier components of the electric field strength [28] which are three coupled equations for the complex cartesian components (q, ), , = , , with being the dielectric tensor of the medium. Non-trivial solutions for the field strength exist if the determinant of the matrix in the braces vanishes. In an isotropic medium, such as a plasma, the dielectric tensor has only two non-vanishing components -the longitudinal, wherê ( ) is the plasmon frequency for wavenumber which is, in general, a complex function. This means the solution of Maxwell's equation in the plasma, following a (possibly random) excitation with wave number q have a time-dependence that is governed by the solution of Eq. (10), ( , ) ∼ − ̂ ( ) . In thermodynamic equilibrium, this solution has to vanish in the long time limit, which requires Im̂ ( ) < 0.
In case of weak damping, this solution can be approximated by the roots of the real part, and the plasmon damping ( ) follows in perturbation theory [17] ,

Weak damping approximation
Let us now turn to the dispersion of collective plasmon oscillations which, for weak damping, is derived from Re RPA = 0. In a plasma in 3D in the limit → 0, the dispersion starts at = . For finite , there appear corrections that involve even powers of . There is a large body of work for classical and quantum plasmas at zero and finite temperature in the absence of correlation effects. Let us summarize some common approximations that include the terms of order 2 and 4 (a derivation of the full RPA results at finite temperature is presented in the appendix): 1. The first result for the 2 term of a classical plasma is due to Bohm and Gross [29] , where 2 th = 3 ∕ is the thermal velocity.
2. In a degenerate quantum plasma at = 0, this dispersion is replaced by [28] 2 ( ) 2 = 1 + 3. The first account of the 4 corrections to the dispersion for an ideal Fermi gas at = 0 is due to Bohm and Pines [66] , and the result was subsequently improved by Ferrell [31] who reported , and the subscript "0" indicates the average with the ground state Fermi function. Note that the term with (Δ 2 0 ) 2 [i.e. with 16 2 ] is not present in Ref. [66] and most other works. 4. An extension of the ideal Fermi gas parametrization to finite temperatures was reported by Arista and Brandt [32] and, more recently, by Thiele et al. [38] 2 ( ) 2 = 1 + who then replaced, in the denominator, ( ) → . For finite temperatures on the order of = 1, the authors of Ref. [38] neglected the ⟨ 4 ⟩ term and proposed the following parametrization: with the degeneracy parameter = Λ 3 = ℎ 3 (2 ) −3∕2 .
5. The finite temperature RPA dispersion (16) can be further improved if the terms ( ) in the denominator are not replaced by but, instead, the full result is used iteratively: where the velocity moments are computed with the finite temperature Fermi function. This gives the most accurate result for the 4 coefficient (see Appendix). Evaluating the Fermi integrals we find the following parametrization of the dispersion, where, for the coefficients, very accurate analytical expressions are presented in the Appendix, The result (19) is shown in Fig. 1 by line "I" and is the most accurate RPA result in the weak damping approximation for small . The result that follows when the term (Δ 2 ) 2 is replaced by 4 is shown by line "II". Eq. (17) leads to the line labelled "III". It follows from neglecting (Δ 2 ) 2 . The approximation that neglects, in Eq. (18), all terms of order 4 , is shown by the line "IV". Note that all of these approximations do not take into account that the dispersion relation Re = 0 has solutions only for a finite range of wavenumbers.

Beyond weak damping
If, however, the damping is not small, it is straightforward to improve approximations (11) and (12) by extending the Taylor expansion of the complex dispersion relation (10) to terms that are second order in | |∕ . The next order result has the form of two coupled equations (for details, see Appendix B) This approximation improves the accuracy of the dispersion in the case of moderate damping, as we will demonstrate in Sec. 3.2.
If the damping is large, the above Taylor expansion will fail, and we have to return to the full condition (10). We rewrite it explicitly in terms of real and imaginary parts where we introduced the complex frequencŷ . [Note that we defined = −Im̂ , which is positive in an isotropic equilibrium plasma, as shown above and, moreover also in an isotropic plasma out of equilibrium [67] .] The pair of equations (21), if solved for several points in the complex frequency plane, yields the dispersion and damping simultaneously, cf. Sec. 3.2. The relations (21) imply that the retarded dielectric function has been analytically continued into the complex frequency plane which has occasionally been done for the mean field approximation [68][69][70] . Here we will report results for the RPA dielectric function, Eq. (8) and for the correlated dielectric function, using the static LFC.

Analytic continuation of the dynamic dielectric function
Taking Eq. (10) or Eqs. (21) seriously requires the knowledge of the retarded dielectric function in the complex plane. Therefore, it needs to be analytically continued from the real axis into the complex plane. The analytic continuation of the whole dielectric function is based on the complex continuation of the finite-temperature Lindhard polarization function. For real frequencies, the retarded and advanced polarization functions are given by Starting from the spectral function we obtain the analytic continuation to the complex frequency plain using Cauchy's integral formula yielding Π in the upper and Π in the lower half-plain. As discussed above, cf. Eq. (21), collective modes follow from zeroes of the retarded dielectric function, which may appear in the lower half-plane, where the retarded polarization function is given bỹ In order to evaluate the spectral functionΠ at complex frequencies, the integration in Eq. (23) has to be performed first, yielding: Using this result, the replacement → =̂ can be made, and the analytic continuation of Π to the lower half-plane can be carried out. The analytic continuation of the real and imaginary part of the RPA dielectric function follows from this easily. These results will be used in our numerical analysis of the plasmon dispersion in Sec. 3.2.
To extend this result beyond the weak coupling case, we have to perform an analytic continuation of the correlated dielectric function. This problem is solved in the following way. Our starting point is the -formally exact -relation (3), which includes correlation effects via our ab initio data for the dynamic local field correction ( , ). Using the static approximation ( , 0), this expression is easily continued (we denote the analytic continuation by a "hat") because it only involves the analytic continuation of the Lindhard function.

Ab initio dynamic dielectric function
Using Eq. (3), the dynamic dielectric function is directly expressed by the local field correction to which we have access in our ab initio simulations. Thus, it is straightforward to directly compare the RPA dielectric function to correlated results that use either the static or dynamic LFC. A first typical result for the dielectric function of the correlated electron gas is shown in Fig. 2 , for the cases of = 2 (bottom row) and = 8 (top row) and = 1. In the left (right) panel, we show the real (imaginary) part of the dielectric function for Im  two wave numbers. At large frequencies, ≳ , the correlated results are in close agreement with the RPA. However, below deviations occur that increase with . The peak of the imaginary part narrows and shifts to much lower frequencies. Due to the Kramers-Kronig relations, the same trend is observed for the real part. The statistical uncertainty of the reconstruction of ( , ) leads to an uncertainty in the region of the peak of Im that is indicated by the red band. Interestingly, the static approximation is very close to the full dynamic results, at the present parameters. In the left column of Fig. 2 , we also show the imaginary part of the inverse dielectric function, -Im −1 , which is proportional to the dynamic structure factor, cf. Eq. (1). For the case of = 8 at the lower wave number (top left figure), its peak is close to the larger zero of the real part of , whereas for = 2 no zeroes of the real part of exist at these -values.
The wavenumber dependence of the dielectric function for the case of = 2 is explored more in detail in Fig. 3 , for smaller than in Fig. 2 . Here we include only the result using the SLFC in addition to the RPA, because the difference to the full dynamic result is very small. The existence of zeroes of Re on the real frequency axis sensitively depends on the wave number: for small wavenumbers, Re has two zeroes, but with increasing , the zeroes vanish as in the examples of wave numbers = 0.4 and 0.5 in Fig. 3 . During this transition, the value of the imaginary part of at the upper zero of the real part and thus the plasmon damping ( ) [Eq. (12)] increase drastically. Hence, the peak of -Im −1 [and, with it, the peak of ( , ), cf. Eq. (1)] broadens strongly. However, this transition is clearly beyond the validity of the weak damping approximation for the plasmon, (11) and (12), and requires to consider improved approximations that were discussed in Sec. 2.4 and 2.4.3. We will study these effects on the plasmon dispersion below.

Perturbation results for the RPA plasmon dispersion
As we discussed in Sec. 2.4 the approximation (11) applies only for sufficiently weak damping. There we also discussed a straightforward way to relax this restriction. To this end we have improved the weak damping approximation, Eqs. (11) and (12) by extending the Taylor expansion of the dielectric function to third order in | |∕ . The result is given by Eqs. (20), for details see Appendix B. The results of the weak damping approximation and of the second order Taylor expansion [using the iterative procedure, cf. Appendix] are compared in Fig. 4 for a moderate temperature, = 1, and two densities, = 1 and = 4. The figure confirms that the weak damping approximation, Eqs. (11) and (12), agrees well with the next order of the expansion, Eqs. (20), for small and moderate wavenumbers. Deviations start growing for ∕ ≳ 0.25, for = 1, and ∕ ≳ 0.5, for = 4. Notice that, at these wavenumbers, the ratio | ( )|∕ ( ) is below 0.1, but nevertheless the expansion already breaks down. Since it is not clear up to what -values the higher order expansion, Eqs. (20), will yield an improvement, we also include,  The results of the analytic continuation, for the parameters of Fig. 4 , are shown by the green lines [the numerical results will be discussed below, in Sec. 3.3]. The comparison confirms that the higher order expansion is more accurate than the lowest order one and is valid up to larger wavenumbers than the latter, i.e. up to ∕ ≈ 0.3, for = 1, and up to ∕ ≈ 0.6, for = 4. The main conclusion from this comparison is that, even though the weak damping condition is well fulfilled [even at these wavenumbers, | ( )|∕ ( ) ≲ 0.2], the approximation Re = 0 -which is the commonly used condition for plasma oscillations -fails badly. Not only are the values for the plasmon frequency wrong and, even more so, for the damping (except for small wavenumbers), the weak damping approximation also makes grossly incorrect predictions for the wavenumber range where plasmons exist. We will return to this issue in Sec. 3.4.

Analytic continuation of the dielectric function
Due to the convergence problems of the weak damping expansion, it is important to improve the computation of the plasmon dispersion and damping by avoiding any weak damping ansatz and Taylor expansion. This is indeed possible, by resorting to the solution of the complex dispersion relation, Eq. (21), i.e. by performing the analytic continuation of the dielectric function, as we discussed in Sec. 2.4. We have carried out this procedure for both, the RPA and the static LFC (SLFC) approximation, Eq. (26). The result of the analytic continuation, for the same parameters and two wave numbers of Fig. 3 , is plotted in the top row of Fig. 5 . This figure shows iso-lines of the zeroes of the real and imaginary parts of the dielectric function in the complex frequency space (full and dotted lines, respectively). Zeroes of the real part on the real axis are clearly visible for = 0.4 , but no real zeroes survive for = 0.5 in agreement with Fig. 3 . At the same time, even for = 0.5 , complex zeroeŝ ( ), i.e. crossings of the full and dotted lines, exist. In the left panel ( = 0.4 ) the plasmon is located at̂ ∕ ≈ 1.36 − 0.09 which is close the solution of the weak damping approximation. In the right part ( = 0.5 ) the plasmon parameters arê ∕ ≈ 1.52 − 0.2. There exist further crossings of the real and imaginary parts located at higher imaginary parts of the frequency. These excitations are strongly damped and therefore suppressed.
A second case of the analytic continuation, for the larger coupling of = 10 and = 1, is presented in the bottom row of Fig. 5 . As in the previous case, collective mode solutions in the complex plane extend to much larger wave numbers than the weak damping approximation predicts. Again the static LFC data are at lower frequency and higher damping than the RPA plasmon. This is particularly clearly visible at larger wavenumbers (right column). This figure also shows another peculiarity of the correlated dielectric function: for small wave numbers, the real part of has only a single branch Re (̂ ) = 0. This is not an artifact of the analytic continuation but a correlation effect that is visible also in the behavior of Re as a function of real frequencies. A further analysis of this effect is given in Ref. [56] .

Ab intio plasmon dispersion of the correlated electron gas
So far, we have considered the dispersion of plasma oscillations ( ) and their damping ( ) focusing on the RPA. Having the correlated dielectric function available, cf. Sec. 3.1, we can now extend this analysis to the correlated electron gas at finite temperature. Due to the close agreement of the static and dynamic approximations, we will only use SLFC to evaluate the correlated plasmon dispersion. Figure 6 summarizes the results for the plasmon dispersion and damping over a broad range of densities and temperatures corresponding to 2 ≤ ≤ 10 and 0.5 ≤ ≤ 2. We include the results from the weak damping approximation and from the analytical continuation of both, RPA and SLFC , and also the analytical approximations (18) and (16) for the RPA dispersion. The main observations on the complex behavior are as follows.
1. The general trend is that of steeper slopes of the dispersion with increased temperature and reduced slope with decreasing density. This overall trend is expected from the prefactor of the 2 term in the RPA dispersion which is given by the thermal velocity, in the classical case, and the Fermi velocity, at strong degeneracy, for details see Sec. 2.4.1. This trend is reproduced well by all computational approaches.
2. The inclusion of correlation effects via the SLFC shifts the plasmon to lower values, compared to the RPA, and reduces the slope (smaller prefactors of 2 and 4 terms). Further, as expected, the SLFC results are systematically more strongly damped than the RPA data because they include, in addition to Landau damping also correlation induced damping effects.
3. All approaches agree very well for small wave numbers where local field effects are negligible and the RPA collective modes exist and are accurate. Thus, for sake of efficiency and accuracy, the analytical low-series expansions I to IV (green lines) are to be preferred in this regime, with "I" corresponding to the parametrization (19), being the most accurate one.
4. The perturbative approaches that are based on the weak damping expansion and use Eqs. 5. However, the wave number range for which this is the case cannot be easily estimated beforehand, see the results in section 3.2 and Fig. 4 . In the presented range, it is found that one could use expansion IV (Eq. (18), no 4 terms) for an estimate of the plasmon location for higher temperatures. Near the one-particle continuum (at low T), "I" [Eq. (18)] is to be preferred. 6. In accordance with the complex dispersion relation, Eq. (10), the method of analytical continuation (black lines) is the most accurate approach to the plasmon dispersion. It also yields solutions for intermediate wave numbers, where the perturbative approaches fail. At the same time for sufficiently large wavenumbers, the complex approach to the dispersion also ceases to find collective modes. Interestingly, the termination is observed close to the point where the dispersion ( ) enters the pair continuum.
7. We also observe striking deviations of the location of the maximum of -Im −1 from the plasmon dispersion before the line = 1, which gives a rough limit to the dominance of collective effects in the dielectric function, and also before the pair continuum is reached. We will discuss the significance of the fact below, in Sec. 3.5.
We conclude that the perturbation approach (weak damping expansion, orange lines) and complex theory (black lines) give good agreement for the damping, for RPA and SLFC both, in the low-q range where the former is valid. On the other hand, the FWHM of the plasmon peak of -Im −1 on the real axis is in good agreement with the damping derived from the complex method. This can obviously only be the case for parameters for which a plasmon as such exists and shows the limits of interpreting the peak of -Im −1 and its FWHM as "plasmon" or "plasmon damping". Naturally, as soon as the peak width becomes of the order of the mode frequency, as is observed e.g. in the case = 2, = 10, in the bottom right part of Fig. 6 , the interpretation as a collective excitation is not appropriate. We discusse this issue in more detail in the next section.

Comparing the complex plasmon dispersion to the peak of the dynamic structure factor
In Sec. 2.2 we discussed how to obtain ab initio PIMC data for the dynamic structure factor of the correlated electron gas. Even though this function is of prime importance for comparison with Thomson scattering data of warm dense matter, the physical interpretation of the results -even for jellium -is not trivial. The reason is that ( , ) is due to a number of different processes, including single-particle (particle-hole excitations) and collective plasma oscillations, e.g. [11,38] . From the relation to the dielectric function, Eq. (1), it is clear that sharp peaks of ( , ) at small wave numbers most likely correspond to a collective mode. However, in the general case, a broad peak of is made from a mix of collective and single particle effect and thus does not directly yield information about the plasmon dispersion. Since correlation effects lead to an additional broadening of the peaks, this correspondence between plasmons and ( , ) becomes even more difficult. Having ab initio results for both, ( , ) and ( , ), available, we are able to answer this question rigorously, for the first time.
Before analyzing the numerical results, let us recall the commonly used criteria for a separation of collective and single-particle excitations.
A. At zero temperature, the boundary is given by the pair continuum that is shown in Fig. 6 by the grey shaded area. Of course, at finite temperature, the pair continuum has merely qualitative relevance.
B. The second criterion for collective excitations is the existence of zeroes of the dispersion relation. If no solutions exist at a given wave number, the excitations are expected to be entirely of single-particle nature. At the same time we have seen, that the existence condition of solutions strongly differs for the weak damping approximation and for the solution of the complex dispersion relation, where obviously the latter constitutes the most accurate result.
C. Finally, a third qualitative separation of the two types of processes is often performed on the basis of the so-called scattering parameter [11] = 1 , where is the screening length: for ≳ 1 ( < 1) collective (single-particle) processes dominate the response of the plasma to the radiation of wavenumber , as was already noted by Bohm and Gross [29] . The reason is that, for wavelengths exceeding the screening length all particles inside of a sphere of radius will be excited simultaneously ("collectively"). Thus, Eq. (27) predicts the existence of a critical wave number, max = 1∕ beyond which no collective excitations exist.
The existence of a maximum wavenumber for collective modes is in full qualitative agreement with the trends predicted by the other two criteria. It is, therefore, interesting to perform also a quantitative comparison of the criteria A-C over the relevant range of densities and temperatures. Such an analysis has already been performed in Fig. 6 , and more detailed results will be given below in Fig. 7 . For the evaluation of the scattering parameter, Eq. (27), we use the stating long wavelength limit of the RPA polarization, as is described in Appendix C.
Let us now discuss the data displayed in Figs. 6 and 7 with respect to the distinction between single particle and collective effects. First, the wave number max following from = 1 is shown in all plots, cf. the vertical black dashed lines. Second, the pair continuum is shown by the grey shaded area. Complex zeroes of the dispersion relation have been plotted for all wavenumbers for which solutions of the corresponding dispersion relations exist. Remarkably, even for finite temperature, complex zeroes cease to exist in the vicinity of the grey areas (pair continuum). The estimation using = 1 agrees excellently with the other two criteria for low temperatures and small , whereas for higher temperatures, it predicts a significantly too small wavenumber max . This is particularly striking at small , e.g. the bottom left part of Fig. 6 . Interestingly, for larger , the criterion = 1 works better again. Of course, there is only a single strict criterion which is the existence of solutions of the complex dispersion relation (21) that follows from the analytical continuation of the dielectric function.
Thus, Figure 6 allows one to verify the reliability of the other criteria and approximations. In particular, by comparing with the complex dispersion relation for RPA , we conclude that use of the RPA static screening length in Eq. (27) is most accurate for a particular -value, at a given temperature: For Θ = 0.5, for ≈ 2, for Θ = 1.0, for ≈ 6, and for Θ = 2.0, for ≈ 10. For smaller (larger) values of , at the same temperature, the screening length is overestimated (underestimated). Note that this behavior of the screening length is not due to the neglect of correlations but reflects a deficiency of the oversimplified criterion (14).
After having discussed the wave number range where collective electronic plasma oscillations exist under WDM conditions, let us now turn to the dynamic structure factor and analyze how reliably it captures plasmons. To this end, we compare in Figs. 6 & 7 the results for the plasmon dispersion that follow from the complex zeroes of the dielectric function (black lines) to the data for the imaginary part of the inverse of the dielectric function (red lines, -Im −1 ) and to data for the dynamic structure factor ( , ) (blue lines) where, for the latter two cases, the position and FWHM of the peak are being used. This comparison is again performed for two sets of approximations -the RPA, RPA , and the static local field correction, SLFC , respectively. The peak positions of ( , ) and -Im −1 ( , ) agree well with the plasmon dispersion obtained from the complex dielectric function for small wave numbers, both, for RPA and the correlated result following from the static LFC. However, when the wave number increases and reaches the wider vicinity of any of the criterias A-C, the peak of -Im −1 ( , ) and even more so the peak of the dynamic structure factor quickly fall below the plasmon dispersion. The effect is particularly strong for the correlated (SLFC) result when and thus the coupling increases. In particular, for ≳ 6 the peak positions of -Im −1 ( , ) and ( , ) even decrease with in a finite wavenumber range (for = 10 this range is approximately between and 2 ). This was first observed in Ref. [53] where it was interpreted as indication of a negative plasmon dispersion. As the width of the peak of ( , ) in that range of wavenumbers is of the same order as the frequency ( ) of the peak, one would rather expect an overdamped oscillation. Even though this effect should be observable in XRTS experiments and is an exciting correlation effect in WDM, it is, however, mainly due to single-particle excitations. This conclusion can now be clearly made, based on our analysis of the complex plasmon dispersion that reveals a strictly monotonic increase of ( ) where the solution ceases to exist at ≈ 1.4 ( ≈ 1.5) for SLFC (for RPA), in the case of = 4. For = 10, the critical wavenumbers are correspondingly ≈ 1.43 for both, SLFC and RPA.
In general, the value of the FWHM of ( , ) and -Im −1 ( , ) is only equal to the damping of the plasmons for very small wavenumbers. As for the location of the plasmon that was discussed above, deviations between plasmon damping and FWHM start to appear early in such a way that the FWHM can be larger or smaller than the damping of the plasmon. If the FWHM is larger than the plasmon damping from the complex theory, we take that as an indicator for (correlated) single particle effects to contribute to the signal. When, for large wavenumbers, the FWHM is lower than the plasmon damping, we find an overdamped state with mostly single particle excitations contributing to the spectrum.
Thus, the procedure to extract plasmon dispersions from an XRTS experiment should never be based on the dynamic structure factor alone, except for small wavenumbers and narrow peaks. In general, one should use, in addition the spectral function which is the imaginary part of the inverse dielectric function as well as the dielectric function itself and solve the dispersion relation. Since Im −1 naturally contains all possible excitations, not just plasmons, to extract the collective excitations one has to analyze the dielectric function.

SUMMARY AND DISCUSSION
In this paper we have extended the ab initio quantum Monte Carlo approach to the dynamic structure factor, that was presented in Ref. [53] , to the dynamic dielectric function. In particular, the ab initio results for the dynamic and static local field correction have allowed us to obtain unbiased accurate results for the dispersion of longitudinal electronic plasma oscillations that fully include correlation and finite temperature effects. To explore the effect of exchange and correlations, we performed a detailed comparison of SLFC results to the random phase approximation (RPA).
To have reference results, we revisited the commonly used RPA dispersion relation, RPA ( , ) = 0, and compared different analytical parametrizations of ( ). The most accurate one was found to be the ground state result of Ferrell, Eq. (15), which we extended to arbitrary finite temperatures, and for which we presented an accurate analytical parametrization given by Eq. (19).
In addition to the analytical parametrizations for the RPA plasmon dispersion, we explored several numerical methods of obtaining the plasmon dispersion and damping: perturbation theory (with respect to the damping) for the zeroes of the dielectric function on the real frequency axis, higher order perturbation theory results, and the analytical continuation of the retarded dielectric function into the lower frequency half plane. Finally we also performed a direct determination of the location and width of the plasmon from the peaks of the plasmon spectral function, -Im −1 .
The true plasmon dispersion and damping are determined by the complex dispersion relation, requiring one to solve simultaneously the two equations, Re =Im = 0, with the complex solution̂ ( ) = ( ) − ( ). The spectral function reflects all aspects of the longitudinal excitations, not just plasmons, and the position, FWHM and shape of its peak can be directly compared to the XRTS measurements.
However, the interpretation of the experimental signal in terms of collective or single-particle excitations, requires an independent evaluation of the plasmon dispersion. Here all of the presented methods give identical results for very small wavenumbers, deep in the collective regime. In this case, the perturbation theory has clear advantages as it does not suffer much from numerical noise as is the case for the complex method or the spectral functions. For mixed regimes of collective and single-particle excitations, we observed that a separation based on the scattering parameter can be very inaccurate. Instead, the analytic continuation method is clearly the best in extracting the plasmon excitations and its wave number range.
The present first ab initio results for the plasmon dispersion, including correlation effects, will be particularly important for WDM experiments at elevated temperatures, ∼ , i.e. in the plasma phase. The biggest correlation effects should occur for comparatively low densities where ∼ 4 … 10, and wavenumbers of the order of 1 … 2 . A suitable candidate could be hydrogen.

APPENDIX APPENDIX A: ANALYTICAL PARAMETRIZATION OF THE RPA PLASMON DISPERSION RELATION
In this Appendix we present details on the analytical plasmon dispersion for a uniform electron gas at finite temperature that were discussed in Sec. 3.2. We consider the regime with ≫ ℏ 2 ∕(2 ) and ≫ . In this limit, the expansion of the real part of the dielectric function in RPA has the following form [32] : wherẽ = ( ∕ )( ∕ ). In Eq. (1), the velocity moments can be expressed in terms of the Fermi integrals of order , , as where = ∕ and is the Fermi velocity. From the condition of the plasmon resonance, in the weak damping approximation, Re ( , ) = 0, and Eq. (1), the dispersion relation follows: 2 ( ) 2 = 1 + which contains the plasmon dispersion also on the r.h.s., so we proceed by iteration. The first iteration is obtained by substituting ( ) = into the right hand side of Eq. (3) [32,38] : The parametrization (19) with the coefficients 2 ( , ) and 4 ( , ) agrees with the exact numerical results with a precision better than 3% in the entire range of , , and and also reproduces the exact analytical limits at ≪ 1 and ≫ 1.

APPENDIX B: PLASMON DISPERSION RELATION FOR MODERATE DAMPING
Here we extend the weak damping result for the plasmon dispersion and damping, Eqs. (11) and (12), to stronger damping. To this end, we extend the Taylor expansion of the complex dispersion relation to terms of order ( ∕ ) 3 : 2 2! Re ′′ ( ) + (− ) 3 3! Re ′′′ ( ) There are several ways to solve this equation.

Iterative solution
The first is an iterative solution, for a fixed .
iteration 1: The first approximation for the frequency, 1 ( ), follows from the lowest order to the real part of Eq. (9): Re ( 1 ) = 0, and coincides with approximation (11).
iteration 2: Inserting this result into the imaginary part of Eq. (9) yields [we denote ∕ by a prime] 1 = ( 1 ): which is the previous result (12) that is improved in the following.
and so on. While this iterative approach significantly improves the result for the dispersion and damping, it has the disadvantage that the existence condition for solutions ( ) is determined by the first iteration, i.e. by the weak damping dispersion relation.

Selfconsistent solution
The iterative procedure can be avoided by solving Eq. (9) directly, simultaneously for ( ) and ( ). This leads to the system of two coupled equations for the real and imaginary parts of Eq. Im ′′ ( ) .
If the damping is moderate, one can restrict the second equation to the first two terms on the right because the last two terms are of third order in the damping. The simplest approximation for the screening length is that for an ideal Fermi gas which, in the ground state, is given by the Thomas-Fermi length

APPENDIX C: SCATTERING PARAMETER AND SCREENING PARAMETER
For finite temperature this result is generalized by computing the static long wavelength limit of the RPA polarization function Π, Eq. (24). The result is shown in Fig. A1. For zero temperature, we recover the analytical result (13) whereas, with increasing temperature, the screening length (the screening parameter = −1 ) increases (decreases). For very high temperatures, the result is identical to the Debye screening. Thus, the maximum wavenumber for collective modes, max = predicted by criterion C can be directly read off from Fig. A1.