Momentum distribution function and short-range correlations of the warm dense electron gas -- ab initio quantum Monte Carlo results

In a classical plasma the momentum distribution, $n(k)$, decays exponentially, for large $k$, and the same is observed for an ideal Fermi gas. However, when quantum and correlation effects are relevant simultaneously, an algebraic decay, $n_\infty(k)\sim k^{-8}$ has been predicted. This is of relevance for cross sections and threshold processes in dense plasmas that depend on the number of energetic particles. Here we present extensive \textit{ab initio} results for the momentum distribution of the nonideal uniform electron gas at warm dense matter conditions. Our results are based on first principle fermionic path integral Monte Carlo (CPIMC) simulations and clearly confirm the $k^{-8}$ asymptotic. This asymptotic behavior is directly linked to short-range correlations which are analyzed via the on-top pair distribution function (on-top PDF), i.e. the PDF of electrons with opposite spin. We present extensive results for the density and temperature dependence of the on-top PDF and for the momentum distribution in the entire momentum range.

Here we consider another many-particle property -the momentum distribution function n(k) and how it is influenced by finite temperature and Coulomb interaction effects. It is well known that, for classical systems in thermodynamic equilibrium, n(k) is always of Maxwellian form regardless of the strength of the interaction. In contrast, in a quantum system the momentum and coordinate dependencies do not decouple which leads to fundamentally different behaviors of n(k) in ideal and nonideal quantum systems, and only for an ideal system the famil-iar Fermi distribution, n id (k) is being recovered (here we consider only Fermi systems). However, in a non-ideal Fermi system, the momentum distribution decays much slower with k, exhibiting a power law asymptotic. The importance of a power law asymptotic has been pointed out by Starostin and co-workers [33][34][35] and many others, e.g. [36], because an increased number of particles in high-momentum states could have a significant effect on scattering and reaction cross sections, in particular on fusion reaction rates in dense plasmas [37][38][39]. The main goal of the present paper is, therefore, to present accurate theoretical results for the tail of the momentum distribution function. Before outlining our goals in more detail, we briefly recall the main available theoretical results on the large-k asymptotic of the momentum distribution function.
It was first demonstrated by Wigner [40] how to incorporate quantum uncertainty between coordinate and momentum into n(k). Following the development of perturbation theory for the electron gas in the 1950's, e.g. [41], [42], Daniel and Vosko [43] calculated the momentum distribution for an interacting electron gas. They used the approximation due to Gell-Mann and Brueckner for the correlation energy [44] which corresponds to the random rhase Approximation (RPA). For the ground state, T = 0 K, they derived an analytical expression for the large-k asymptotic of the momentum distribution, i.e. they found an algebraic decay, in striking contrast to the exponential asymptotic of an ideal classical or quantum system. Galitskii and Yakimets [45] used Matsubara Green functions and the Kadanoff-Baym relation [46] between the energy distribution in equilibrium, f EQ (ω) [which is always a Fermi or Bose distribution], and the spectral function A(k, ω), Correlation effects enter only via the spectral function A, which is given by A id (k, ω) = 2πδ[ ω − E(k)], for an ideal gas. Ref. [45] computed the leading correction to the ideal spectral function and confirmed the asymptotic, Eq. (1). For a systematic improvement of this result higher order selfenergies have been computed, e.g. by Kraeft et al. Ref. [47], and we also refer to the text books Refs. [46,48,49]. The exact limiting behavior in the asymptotic (1) was found independently by Kimball [50] via a shortrange ansatz to the two-electron wave function, and by Yasuhara and Kawazoe [51] who analyzed the largemomentum behavior of the ladder terms in Goldstone perturbation theory. An important result of Yasuhara et al. is the proof [51] that, at T = 0 K, the asymptotic can be expressed via the on-top pair distribution function (on-top PDF), i.e. the PDF of a particle pair with different spin projections at zero distance, g ↑↓ (r = 0), lim k→∞ n(k) = 4 9

9π
2/3 r s where k F denotes the Fermi momentum, and the coupling (Brueckner) parameter r s =r/a B is the ratio of the mean interparticle distance,r = [3/(4πn)] 1/3 , to the Bohr radius [52]. A more general derivation has been presented by Hofmann et al. [53] who have shown that Eq. (3) holds also for finite temperature. An extension of the results of Yasuhara et al. and Kimball to arbitrary spin polarizations of the electron gas was performed by Rajagopal et al. in Ref. [54] who derived the next order in the asymptotic which becomes dominant in the case of a ferromagnetic electron gas because the on-top PDF vanishes: Aside from dense plasmas, the tail of the momentum distribution is also relevant for the electron gas in metals, e.g. [55], as well as cold fermionic atoms [56,57]. In the latter case, however, the short-range character of the pair interaction leads to a modified large-momentum asymptotic, n(k) ∼ k −4 , instead of (1). A second approach to the high-momentum tail is based on quantum Monte Carlo simulations. Here one can either directly compute the asymptotic of n(k) or determine it from the Fourier transform of the density matrix. While the former requires to extend the simulations to very large momenta and to resolve the occupations over many orders of magnitude, the latter way is potentially more efficient. Here one calculates the on-top PDF (which is called "contact" in the cold atomic gas community). In addition to its use in Eq. (3), we mention that an accurate description of g ↑↓ (0) is interesting in its own right, and is important for many other applications, like the description of the static local field correction [58][59][60][61][62].
Accurate QMC results for n(k) of the UEG in the ground state were obtained in Refs. [55,63], whereas the on-top PDF was studied in multiple QMC-based works [55,[64][65][66], most recently by Spink and coworkers [67]. At finite temperatures, the momentum distribution n(k) has been investigated by Militzer et al. [68,69] who carried out restricted path integral Monte Carlo (RPIMC) simulations and recently by Filinov et al. [70] based on a version of fermionic PIMC that is formulated in phase space. Furthermore, the only comprehensive data set for g(0) in this regime was presented by Brown et al. [71], again on the basis of RPIMC simulations.
Note that fermionic PIMC in coordinate space is limited to moderate degeneracy [12,72], due to the notorious fermion sign problem, see Ref. [73] for an accessible topical discussion. On the other hand, RPIMC has been shown to exhibit significant systematic errors of the thermodynamic quantities, for example the error for the exchange-correlation energy reaches 10% at r s = 1 and Θ = 0.25 [74]. In addition, RPIMC is substantially hampered by an additional sampling problem (reference point freezing [75]) at high densities, r s 1.
Therefore, it is of high interest to perform alternative simulations that can access the momentum distribution of the uniform electron gas at high degeneracy without any systematic errors. In this context, a suitable approach is given by the recently developed configuration PIMC (CPIMC) method that is formulated in Fock space (Slater determinant space) and is highly efficient at high to moderate quantum degeneracy [74,76]. In particular, CPIMC simulations were the basis for the first ab initio thermodynamic results for the warm dense UEG [74]. In combination with the likewise novel permutation blocking PIMC [77][78][79] scheme, it was possible to avoid the fermion sign problem and to obtain ab initio thermodynamic results for the UEG at warm dense matter conditions [2,80]. In addition, also ab initio results for the static density response [81] have been obtained with CPIMC.
The goal of this paper is to utilize CPIMC to obtain ab initio data for the momentum distribution of the uniform electron gas at finite temperature and high density corresponding to r s 0.7. To access stronger coupling, we also employ a recently developed approximate methodrestricted CPIMC [82] as well as direct fermionic propagator PIMC simulations in coordinate space -an extension of permutation blocking PIMC [77]. In particular, i: we verify that the high-momentum asymptotic does obey a k −8 behavior, and that it is solely determined by the on-top PDF; ii: we present detailed CPIMC results for g ↑↓ (0) and analyze its temperature and density dependence; iii: investigate the momentum distribution function in the vicinity of the Fermi momentum and for small momenta; iv: investigate the momentum range of the onset of the large-momentum asymptotic.
This paper is organized as follows: In Sec. II we present a brief overview on earlier theoretical work pertaining to the uniform electron gas, together with the main predictions. This is followed by an introduction into our quantum Monte Carlo simulations in Sec. II B and by a presentation of the numerical results in Sec. III.

A. On-top pair distribution
Since the high-momentum tail of the momentum distribution function can be expressed in terms of the on-top pair distributions, cf. Eq. (3), we start by considering the pair distribution of electrons with spin projections σ 1 and σ 2 [83], ] is a fermionic field operator annihilating [creating] an electron in spin state |r 1 σ 1 . Note that the two-particle density in the numerator is normalized to the single-particle spin densities, n σ (r) = Ψ † σ (r)Ψ σ (r) , in the denominator. Thus, in the absence of correlations and exchange effects, g σ1σ2 (r 1 , r 2 ) ≡ 1. For electrons there exist four spin combinations. Assuming a homogeneous paramagnetic system, we have g ↑↑ (r 1 , r 2 ) ≡ g ↓↓ (r 1 , r 2 ) and g ↑↓ (r 1 , r 2 ) ≡ g ↓↑ (r 1 , r 2 ).
The total pair distribution function follows from the spin-resolved functions (5) according to where the normalization assures that, in the absence of exchange and correlation effects, g ≡ 1. In a spatially homogeneous system, such as the UEG, the PDFs depend only on the distance of the pair, g σ1σ2 (r 1 , r 2 ) = g σ1σ2 (|r 2 − r 1 |). Of particular importance is the case of zero separation. Then, the Pauli principle leads to g ↑↑ (0) ≡ g ↓↓ (0) ≡ 0. On the other hand, the probability of finding two electrons with different spins "on top of each other" yields the on-top PDF, g ↑↓ (0), which is related to total PDF in the paramagnetic case by [cf. Eq. (6)] which is a fundamental property for the characterization of short-range correlations. While in a non-interacting system (r s → 0), g ↑↓ id (0) = 1, Coulomb repulsion leads to a reduction of this value. Thus for the UEG a monotonic reduction with r s is expected which will directly influence, via Eq. (3), the tail of the momentum distribution.
There exist a variety of analytical parametrizations of the on-top PDF. The ground state on-top PDF of correlated electrons was investigated in Ref. [84] by using the Overhauser screened Coulomb potential in the radial two-particle Schrödinger equation. The results were parametrized for r s ≤ 10 according to where A = 0.0207, B = 0.08193, C = −0.01277, D = 0.001859 and E = 0.7524. These results will be called "Overhauser model" and used for comparison below.
On the other hand, the high-temperature asymptotic of the on-top PDF of a classical non-degenerate electron gas where χ = nΛ 3 1, and Θ = k B T /E F 1, is also known. Here n is the density depending on the mean inter-particle distance,r ∼ n −1/3 , and Λ is the thermal DeBroglie wavelength, Λ 2 = h 2 /(2πmk B T ). A quantummechanical expansion was given in Ref. [53], where the result depends on the order the high-temperature limit, T → ∞, and the classical limit, → 0, are taken. The reason is the existence of a third length scale [48,85], the Bjerrum length, l B = βe 2 , where β = (k B T ) −1 , giving rise to a second dimensionless parameter, the classical coupling parameter, Γ = βe 2 /r = l B /r.
In the case Γ χ 1/3 (i.e. l B Λ), the result is [53] where the behavior is still dominated by the ideal Fermi gas properties with deviations scaling like Γχ −1/3 , or, On the other hand, in the case χ 1/3 Γ (i.e., Λ l B ), which corresponds to classical plasmas at moderate temperatures, the on-top PDF becomes [53] This value is exponentially small due to the moderate Coulomb repulsion and is not influenced by quantum effects. Nevertheless, quantum effects (finite Λ) show up in the algebraic momentum tail, according to Eq. (3), but only on length scales much smaller than Λ or, correspondingly, at momenta strongly exceeding Λ −1 . The latter case is out of the range of WDM and not relevant for the present analysis.
Finally, there exists a more recent parametrization of the ground state on-top-PDF that is based on QMC simulations [67]: which will be used for comparison below. For an overview about different models of g(0) for the ground state, the reader is referred to the paper by Takada [62]. With explicit results for the on-top PDF and, using Eq. (3), the large-k asymptotics of the momentum distribution function can be reconstructed.
For finite temperature one can relate the PDF to an effective quantum pair potential, g ↑↓ (r) = e −βV Q (r) , an idea that was put forward by Kelbg [86] and further developed, among others, by Deutsch, Ebeling, and Filinov and co-workers, cf. Refs. [87][88][89][90] and references therein. We will return to this issue in Sec. III B 2.

B. Configuration PIMC (CPIMC) approach to g(0)
and n(k) of the warm dense electron gas

Idea of CPIMC simulations
CPIMC was first formulated in Ref. [91] and applied to the UEG in Refs. [74,76,92]. For a detailed description of the CPIMC formalism we refer to the overview articles [2,93] and to the recent developments [82]. Here we only summarize the main idea. The thermodynamic expectation value of an arbitrary operatorÂ is determined by the density operatorρ and its normalization -the partition function Z, where we use the canonical ensemble,ρ Since the Hamiltonian involves only one-and two-body operators, its expectation value can be described via the reduced one-and two-particle density matrices, d ij and d ijkl , see the definitions (16) and (17). Here the sums are over arbitrary complete sets of single-particle states which below will be specified to momentum eigenstates. Quantum Monte Carlo estimators for these quantities are obtained through differentiation of the partition function [93,Eq. (5.88)] with respect to the single-particle matrix element and the two-particle matrix element respectively. The resulting expressions depend on the order and choice of the indices (i, j) and (i, j, k, l), respectively.
Let us now present explicit expressions for the oneparticle and two-particle density matrices in CPIMC. Configuration PIMC is path integral Monte Carlo formulated in Fock space [91], i.e. in the space of N -particle Slater determinants, |{n} = |{n 1 , n 2 , . . . } , constructed from the single-particle orbitals |i where n i is the associated occupation number.
In CPIMC the canonical partition function (13) is written as a Dyson series in imaginary time, for details see Ref. [82]. A configuration C determining a MC state is given by a set of initially occupied orbitals {n}, along with a set of K changes κ i to this set, called kinks at their respective times t i , 1 ≤ i ≤ K, Due to the Slater-Condon rules for fermionic 2-particle operators, each interaction matrix element yields either a 2-particle term, corresponding to κ = (i, j), or a 4particle term, κ = (i, j, k, l). Thus the kinks are given by either two or four orbital indices, respectively. The kink matrix element q i,i−1 (κ i ) represent the off-diagonal matrix elements with respect to the possible choices of 2or 4-tuples κ i . The final result for the partition function is [82] where paths with K = 1 violate the periodicity and have to be excluded. Configurations can be sampled from the partition function with the weight function which allows one to rewrite thermodynamic expectation values (14) as An example configuration (path) is illustrated in Fig. 1.
With three particles present, horizontal solid lines represent diagonal matrix elements, as given by the exponential factor in the partition function (19) and the occupation number state at a given time-interval is specified by the set of all these lines in this interval. On the other hand, the vertical solid lines represent interaction terms, where the occupation changes according to the specified kink κ i , weighted by the respective kink matrix element q i,i−1 (κ i ). Due to the periodicity of the expectation values (14), the kinks must add to yield the initial occupation vector at 0 < t < t 1 again: The three kinks 1, 3, 5, at times t1, t3, t5, each involve four orbitals: κ1 = (1, 4; 3, 5), κ3 = (1, 3; 0, 2), κ5 = (4, 5; 1, 3) respectively. The two kinks 2 and 4, at t2 and t4, involve two orbitals, each: κ2 = (0; 1) and κ4 = (2; 4). The fourth Slater determinant |n (4) exists between the imaginary "times" t3 and t4 and contains three occupied orbitals { 0, 2, 5 }.
This representation of the partition function can now be applied to the observables of interest. For the oneparticle density matrix we obtain, for i = j, For the uniform electron gas, the off-diagonal matrix elements vanish in a momentum basis, whereas the diagonal ones yield the momentum distribution, as will be discussed in Sec. II B 2 Let us now turn to the CPIMC estimator for the twoparticle density matrix. Here we have to distinguish several cases of index combinations [94,Eq. 3.14]. If i < j, k < l are pairwise distinct The term under the sum (without the Kronecker-delta) will be abbreviated as the weight of the kink κ ν , In the case of i = k, but with all other indices being different, Finally, if i = k and j = l, but i = j, the matrix elements are given by The expectation value of this estimator is given by the weighted sum over all possible configurations C, Due to the large single-particle basis sizes that have to be used in the CPIMC simulations, the variances of these estimators may be very large for some transitions [i.e. combinations of indices (i,j) or (i,j,k,l)]. However, special cases can be used to derive the estimators needed to measure short-range properties of the system: The momentum distribution and the on-top PDF.

Momentum distribution with CPIMC
The momentum distribution function is given by the diagonal part of Eq. (16), if a plane wave basis is being used, cf. Sec. II B 3. For i = j, we obtain where the contribution of each time-slice is weighted by the length of horizontal paths.

On-top pair distribution function with CPIMC
The definition (5) of the spin-resolved PDF requires the two-particle density matrix in coordinate representation which is obtained from the two-particle density matrix, Eq. (17), in momentum representation, i.e. using plane wave orbitals, To shorten the notation, the wave vector k will be represented by an index i ↔ k i of the corresponding single-particle basis eigenvalue. The field operators in a position-spin basis are related to the creation and annihilation operators in a momentum-spin basis |i := |k i s i byΨ The on-top PDF, Eq. (8), follows from the PDF, Eq. (5), for different spin projections, σ 2 = σ 1 , With the basis transformation (31) of the field operators, a straightforward calculation yields the CPIMC estimator for the on-top PDF (for details see Appendix A),

III. SIMULATION RESULTS
We have performed extensive CPIMC simulations with N = 54 particles. Due to the fermion sign problem, these simulations are restricted to small coupling parameters, r s 0.7. To extend the range of parameters, we also performed simulations with N = 14 particles. As shown before, important structural properties, such as the static structure factor [16,95] and the pair distribution function only weakly depend on the particle number. A quantitative analysis of the N -dependence of the results will be performed for the tail of the momentum distribution in Sec. III B. The CPIMC results are complemented by restricted CPIMC simulations [82]. To access larger values of the coupling parameter, we also include fermionic PIMC simulation results in coordinate space for the ontop PDF.
A. Momentum distribution

Overview
Let us start by analyzing the general trends of the momentum distribution when either the temperature or the coupling strength are varied. In Fig. 2 we present CPIMC data for N = 54 particles showing the entire momentum range for moderate coupling, r s = 0.5, and three temperatures and indicating that the occupation of high-momentum states is coupled in a non-trivial way to occupation of lower momentum states. Interestingly, an increase of temperature not only leads to the familiar broadening of n(k) around the Fermi edge and depletion below it, but may also lead to a lower population of the tail (see below). The most striking observation is the strong deviation, in the tail region, from the exponential decay in case of an ideal Fermi gas. Our simulations clearly confirm the correlation-induced enhanced population of high-momentum states with the asymptotic, n(k) ∼ k −8 . Let us now turn to the dependence on the coupling parameter. To this end, we present, in Fig. 3, the momentum distribution for a fixed temperature, Θ = 2, and two values of r s and also compare to the ideal Fermi gas. For large momenta, k 6k F , we observe an increase of the population when r s grows. However, for intermediate momenta, k F k 6k F , the ideal distribution is significantly above the correlated distributions. Finally, below the Fermi momentum, the correlated distributions are again above the ideal momentum distribution.
This behavior seems counter intuitive, and we analyze it more in detail in the next section. CPIMC results with N = 54 particles are compared to the ideal Fermi-Dirac distribution n id (full black line). For momenta below approximately 6kF the correlated distributions are indistinguishable from n id . For comparison, the ground state distributions, as given by Ref. [84], are shown by the dashed lines of the same color as the finite temperature result.

Interaction-induced enhanced population of low-momentum states
Let us now investigate in more detail the behavior of the momentum distribution in the range from k = 0 to momenta on the order of several k F . To focus on correlation effects we plot, in Fig. 4, the difference of the correlated distribution and the Fermi distribution for the case of r s = 0.5. Clearly, we observe an enhanced population of low-momentum states, k 1.5k F , compared to the Fermi function. The effect is biggest at the lowest temperature and decreases monotonically with Θ. On the other hand, it is clear that, upon further reduction of Θ, this effect will decrease again and vanish in the ground state. The reason is that, at T = 0 K, all lowmomentum states are completely occupied, and, due to the Paui principle, correlations can only enhance the population of unoccupied states, at k > k F . The same analysis is performed, for a fixed temperature but different coupling parameters, in Fig. 5. Here we observe a monotonic trend: with increasing r s , the difference of the populations increases with respect to the ideal case.
This interaction-induced enhanced population of low-k state has been reported before, e.g. based on restricted PIMC simulations, by Militzer and Pollock [68], and on thermodynamic Green functions by Kraeft et al. [47]. The origin of this effect is interaction-induced lowering of the energy eigenvalues, E(k) < E id (k) [68]. Here, the interacting energy contains, in addition, an exchange and a correlation contribution, The behavior reported here is dominated by the exchange contribution, i.e. by the Hartree-Fock selfenergy (the Hartree term vanishes due to homogeneity and charge neutrality) which is negative, The negative Hartree-Fock selfenergy shift is largest at small momenta and decreases monotonically with k. As a consequence, the system tends to increase the population of low-momentum states.
An interesting consequence of this population increase is that the mean kinetic energy of the correlated electron gas may be lower than that of the ideal electron gas at the same temperature [47,68]. Our simulations clearly confirm this prediction. This effect is illustrated in the lower panels of Figs. 4 and 5 where we plot the k-resolved difference of kinetic energy densities. For the parameters shown in theses figures, the excess kinetic energy (compared to the ideal UEG) concentrated in low-momentum states (positive difference) is smaller than the kinetic energy reduction (negative difference) at larger momenta. This is evident from the areas under the curves in the lower panels of Figs. 4 and 5. As a result the total kinetic energy difference of the interacting system compared to the ideal system is negative for a broad range of parameters. The corresponding kinetic energies for the interact-ing and ideal systems are presented in the appendix, in tables II and III, for 54 and 14 particles, respectively. Our argument, so far, was based on the negative sign of the Hartree-Fock selfenergy. However, for a complete picture we also need to consider the energy shift due to correlations, ∆E c . In contrast to the Hartree-Fock shift, the correlation corrections to the energy dispersion are typically positive, but smaller, as was shown for the Born approximation (Montroll-Ward approximation), in Ref. [47]. However, this result applies only for weak coupling. For stronger coupling, in particular, r s 1, at least T-matrix selfenergies would be required. An alternative are QMC simulations, as presented in Ref. [68], which allow one to map out the range of density and temperature parameters where the difference of correlated and ideal kinetic energies changes sign.
The present CPIMC simulations are not directly applicable to the range r s 1. However, we can take advantage of the accurate parametrization of the exchangecorrelation free energy f xc of Groth et al. [80] that is based on a combination of CPIMC, PB-PIMC and ground state QMC results. In particular, the exchangecorrelation contribution to the kinetic energy is obtained by evaluating [26] K xc = −f xc − θ ∂f xc ∂θ rs − r s ∂f xc ∂r s θ , and the corresponding results are depicted in Fig. 6. The line where the kinetic energy difference changes sign is in good agreement with the results of Ref. [68], for r s 1, but we find significant deviations at smaller r s and lower temperatures.
It is interesting to compare the parameter values where the kinetic energy difference changes sign to the occupation of the zero-momentum state, n(0), relative to the ideal distribution, n id (0). For most temperatures considered, the interacting zero-momentum state n(0) has a larger population than the corresponding ideal state. Only for the lowest temperatures, Θ ∈ { 1/16, 1/8 }, we observe the opposite behavior.

High-momentum asymptotics of n(k)
In Figs. 7 and 8 we present data for low to moderate temperatures focusing on momenta beyond the Fermi edge. We directly compare the CPIMC data to the asymptopic behavior where a k −8 tail is expected, with the coefficient determined by the on-top PDF g(0), cf. Eq. (3), where g(0) is taken from the same CPIMC simulation. As can be seen in these figures, the CPIMC data clearly exhibit the expected algebraic decay, for sufficiently large k. To make a quantitative comparison, we also plot, in the lower panels, the relative difference between CPIMC data, n(k), and the asymptotic, n ∞ (k), according to The results for δ ∞ clearly confirm that our ab initio CPIMC data approach the asymptotic. Moreover, we can estimate the momentum range where the asymptotic behavior dominates. For low temperatures of E F /16, the asymptotic is reached at about 6k F , cf. Fig. 7. With increasing temperature, the asymptotic is approached only at larger momenta, e.g. for Θ = 2, around 11k F , cf. Fig. 8. A systematic analysis of the onset of the asymptotic will be given in Sec. III C.
In these figures we also included ground state data for the momentum distribution (green lines) which allows us to analyze finite temperature effects. In all figures we observe that the finite temperature distribution, n(k; Θ), intersects the ground state function, n(k; 0), coming from above, before it reaches the asymptotic. In the range of the algebraic tail the finite temperature function is always below the ground state result, for the same r s and k, in agreement with Fig. 2. This behavior is, at first sight, counter intuitive because one expects that finite temperature effects increase the population of high momentum states. As we will show in Sec. III B 2 this temperature dependence is, in fact, non-monotonic and is due to a competition between Coulomb repulsion and exchange effects.
Finally, we note that our simulations reveal that the k −8 asymptotic is observed independently of the particle number, in agreement with the predictions of Refs. [53,97]. We will return to the question of the particle number dependence in Sec. III B 1. After analyzing CPIMC data for the large momentum tail of the distribution function we now concentrate on the coefficient in front of the asymptotic k −8 term. According to Eq. (3), this coefficient is entirely determined by the on-top PDF g(0) which is directly accessible in quantum Monte Carlo simulations. For PIMC in coordinate space, the straightforoward way is to analyze the r-dependence of the PDF and subsequently extrapolate to r = 0. Typical results are shown in Fig. 9 for direct fermionic (labeled "PIMC") and restricted ("RPIMC") PIMC simulations. In contrast, in CPIMC a direct estimator for the on-top PDF is available, cf. Eq. 33, and the results are included in Fig. 9 with the red symbols. These results depend on the size of the single-particle basis and the corresponding cut-off energy E max (top xaxis). Overall, for a sufficiently large basis, very good agreement of the two independent fermionic simulations -PIMC and CPIMC -is observed for the parameter combinations where both are feasible.
This gives additional support for our CPIMC data, in particular for its use at low temperatures, where CPIMC provides the only ab initio approach. In fact, CPIMC data for g(0) were already used for comparisons above. In this section we investigate the density and temperature dependence of g(0). But first we explore how sensitive this value depends on the number of particles in the simulation cell.

Particle Number Dependence
We have performed extensive CPIMC simulations for g(0) for a broad range of particle numbers, from N = 14 to N = 66. Two typical examples are shown, for Θ = 0.0625, in Fig. 10, and for Θ = 2, in Fig. 11. In these figures we use the case N = 14 as the reference for comparison because, for this number, the widest range of parameters is feasible, although, naturally, simulations with larger N are more accurate. All figures confirm that finite size effects are very small in g(0) and to not exceed 2%, even for N = 14. Regarding simulations with the two approximate CPIMC variants that were discussed above [82], the analysis reveals that RCPIMC+ is reliable for intermediate temperatures, 0.1 Θ 0.5. Even at lower temperatures, cf. Fig. 10, we observe that RCPIMC+ data points for N = 54 are close to CPIMC simulations for N = 54 particles (and more accurate than CPIMC for N = 14) and, therefore, can be well used for larger r s -values, where CPIMC is not possible, due to the sign problem. At the same time, RCPIMC [82] turns out to be not sufficiently accurate for computing g(0) and is not being used in this paper.

Temperature Dependence
We now turn to the temperature dependence of the ontop PDF. In Figs. 12 and 13, we plot g(0) from CPIMC data over a broad range of temperatures for r s = 0.2, and 0.2 ≤ r s ≤ 0.7, respectively. The figures display an interesting non-monotonic behavior: the on-top PDF increases, both towards low and high temperatures. This is easy to understand: At very low temperatures, the system approaches an almost ideal Fermi gas for which g(0) would be exactly 0.5. The (weak) Coulomb repulsion gives rise to an additional depletion of zero distance pair states. This is confirmed by the lower absolute values of g(0) when r s is increased from r s = 0.2 to 0.4 and 0.7.
On the other hand, for increasing temperature, in the range where the electron gas is dominated by classical behavior (Θ > 1), both, exchange and Coulomb repulsion effects are suppressed, as compared to thermal motion, and the probability that two particles approach each other closely, tends to unity, as it would be in a non-interacting classical gas. A non-trivial ques-tion is the position of the minimum. It appears around Θ = k B T /E F ∼ 0.63, with a depth of 0.42, for r s = 0.2, around Θ = 0.63, with a depth of 0.365, for r s = 0.4, and around Θ = 0.63, with a depth of 0.296, for r s = 0.7.
This minimum can be understood as due to the balance of two opposite trends: depletion of g(0), due to Coulomb repulsion and increase of g(0), due to quantum delocalization effects. At high temperatures and low densities, the PDF can be expressed in binary collision (ladder) approximation where V is the Coulomb potential, which reproduces the behavior right of the minimum. At small interparticle distances, r Λ, however, quantum effects have to be taken into account in the pair interaction. Averaging over the finite spatial extension of electrons leads to the replacement of the Coulomb potential by the Kelbg potential (quantum pair potential) [86,98,99], whereΛ = Λ. Note that V K has the asymptotic V K (0; β) = e 2 Λ(β) ∼ T 1/2 which removes the Coulomb singularity at zero separation. While this potential has the correct derivative, dV K (0)/dr = − e 2 Λ 2 , its value at r = 0 is accurate only at weak coupling. At the same time, this potential can be extended to arbitrary coupling by retaining the same analytical form, but correcting the standard thermal DeBroglie wavelength Λ (referring to an ideal gas) to the wave length of interacting particles, which gives rise to the so-called improved Kelbg potential [88,89], At low temperature the effective wavelength of the electrons increases, γ(β) ∼ T −1/2 , which ensures that g ↑↓ IK (0) = e −βV IK (0) is finite. Accurate values for the function γ in a two-component plasma and for different spin projections were presented in Refs. [88,89] from a fit to PIMC data. In similar manner, the present ab initio QMC results for the on-top PDF can be used to compute an effective DeBroglie wavelength of the warm dense uniform electron gas, and the concept of an effective quantum pair potential allows for a simple physical interpretation of some of its thermodynamic properties.
As we already saw for the example of three densities, the location of the minimum changes with the coupling strength r s . This effect is analyzed systematically in Fig. 14. We observe an increase of the minimum position, Θ min , with r s (full squares, left axis). The reason is that, with increasing coupling, the interaction strength increases, as is seen by the increasing depth of the minimum (open symbols, right axis). Therefore, the monotonic increase of g(r) with temperature sets in already at a higher temperature, when r s is increased. In addition to CPIMC data which are restricted to r s 1 we also included an analytical fit ("ESA" [58]) that agrees well with CPIMC and extends the data to r s = 8. More information on this approximation is given in the discussion of Fig. 16.

Density Dependence
Let us now discuss the density dependence of the ontop PDF. As we have seen above, with increasing coupling strength, r s , the value of g ↑↓ (0) decreases, due to the increased interparticle repulsion. This connection can be qualitatively understood from Eq. (38) if it is used with an effective potential that includes many-body effects beyond the pair interaction. This monotonic decrease with r s is confirmed by our simulations for all temperatures. As an illustration, we show in Figs. 15 and 16 the behavior for Θ = 0.0625 and Θ = 1, respectively.
At low temperature and weak coupling, the temperature dependence of g(0) is very weak, cf. Fig. 15, in agreement with Fig. 13. At θ = 1, finite temperature effects increase the particle repulsion due to stronger localization of electrons, and g(0) falls slightly below the ground state value, cf. Fig. 16. This confirms the non-monotonic temperature dependence of g(0) discussed above, since this temperature is in the vicinity of the minimum of g(0).
Let us now discuss the consequences of this density and temperature dependence of g(0) for the high-momentum asymptotics of n(k). According to Eq. (3), the number of electrons occupying large-k states is proportional to n(k; r s , Θ) ∝ r 2 s · g(0; r s , Θ) where we made the dependence on the coupling parameter explicit. Taking into account that k F ∝ n 1/3 ∼ r −1 s , the absolute value of the asymptotic occupation number, at a given k and fixed Θ, scales as n(k) ∝ r −6 s g(0; r s , Θ)·k −8 . On the other hand, considering the occupation number as a function of the momentum normalized to the Fermi momentum, κ = k/k F , the density dependence becomes s(r s , Θ) = 9 2 α 8 r 2 s · g(0; r s , Θ) , α := 4 9π 1 3 .
Given the monotonic decrease of g(0) with r s , the function n(κ) may exhibit non-monotonic behavior as a function of r s , including a maximum at an intermediate r svalue. This is clearly seen in Fig. 17 for the temperatures Θ = 2, 4.
As expected, at all temperatures, the coefficient s(r s ) increases monotonically, for small r s , starting from zero. The decrease, governed by the monotonic decrease of g(0) sets in only at large r s where CPIMC simulations are not possible any more. On the other hand, an extensive set of restricted PIMC data [71] for g(r) is available, for 1 ≤ r s ≤ 40, which has recently been used by Dornheim et al. [58] to construct an analytical parametrization of g(0; r s , θ). The results are denoted as ESA because they constitute an important ingredient to the effective static approximation for the static local field correction that was presented in Ref. [58].
An example is shown in the lower part of Fig. 17 for two temperatures, Θ = 2 and Θ = 4. The maximum of s is observed around r s = 4, for Θ = 2 and r s ≈ 5, for Θ = 4. We have performed a systematic parameter scan on the basis of the analytical fit (ESA) over a broad range of temperatures. The results are collected table I. These results show that the maximum of s(r s ) is generally located in the range 3.5 r s 6.0. Interestingly r max s the r s -value where the maximum is located -exhibits a non-monotonic temperature dependence. The reason is the non-monotonic temperature dependence of g(0) that was discussed in detail in Sec. III B 2. Finally, the comparison with the ab initio results contained in Fig. 17 suggests that the ESA fit can be further improved using our CPIMC and FP-PIMC data.

C. Onset of the large-k asymptotic of n(k)
Let us now find an approximate value of the momentum k ∞ where the k −8 -asymptotic starts to dominate the behavior of the distribution function. In particular, we are interested to understand how this value depends on density and temperature. First, we observe that the significant broadening of the low-momentum part of the distribution that is observed when the temperature is increased pushes the value k ∞ to larger momenta.  by the Fermi-Dirac distribution function n id (k): This approach is demonstrated in Fig. 18, and the results are presented for a broad range of densities, in the range of r s = 0.2 . . . 1.6, and temperatures Θ ≤ 4, in Fig. 19. For this procedure, to obtain the asymptotic n ∞ we used the value of g(0) that was computed in CPIMC simulations. This figure shows that, with an increase of correlations (increase of r s ) the onset of the asymptotic is shifted to lower momenta, even though the dependence is weak. The figure also shows that an algebraic tail of the momentum distribution exists also in a weakly quantum degenerate plasma with Θ > 1. With increasing temperature, the onset of this asymptotic is pushed to larger momenta with k ∞ /k F increasing slightly faster than Θ 0.5 .

A. Summary
In this paper we have performed an analysis of the momentum distribution function of the correlated warm dense electron gas using recently developed ab initio quantum Monte Carlo methods. We have presented extensive data obtained with CPIMC, for small r s . This was complemented with new fermionic propagator PIMC data, for r s 1, so the entire density range hase been covered. Our CPIMC results for the momentum distribution of the warm uniform electron gas achieve an unprecedented accuracy -the asymptotic is resolved up to the eleventh digit for momenta up to approximately 15k F , cf. Figs. 2 and 3. For all parameters the existence of the 1/k 8 asymptotic is confirmed. Moreover, based on accurate data for the on-top PDF the absolute value of n(k) in the asymptotic is obtained.
While the value of the on-top PDF decreases monotonically with r s , it exhibits an interesting non-monotonic temperature dependence with a minimum around Θ = 0.656, e.g. Fig. 13. This was explained by a competition of Coulomb correlations and exchange effects. We also investigated the density and temperature dependence of the momentum where the algebraic decay begins to dominate the tail of the momentum distribution.
In addition to the large-momentum tail we also investigated the occupation of low-momentum states in the warm dense electron gas. An interesting observation is that Coulomb interaction may lead to an enhanced occupation of low-momentum states (compared to the ideal case), which is mostly due to exchange effects, cf. Fig. 5. Together with an enhanced population of highmomentum states this leads to a depopulation of intermediate momenta in the range k F k 3k F . This non-trivial re-distribution of electrons may give rise to a counter-intuitive interaction-induced decrease of the kinetic energy of the finite temperature electron gas. This confirms earlier results [47,68] and, at the same time, complements them extensive new and more accurate data in a broad range of parameters.

B. Outlook
Part of our results for the on-top PDF were obtained with help of the recent extended static approximation (ESA) [58]. Its advantage is that it allows for relatively easy parameter scans in a broad range of densities and temperatures. Therefore, an important task is to further improve this approximation with the present highquality data for g(0). The present simulations concentrated on the range of r s 10 which is of relevance for warm dense matter. At the same time the jellium model is also of interest for the strongly correlated electron liquid, e.g. [101,102]. It will, therefore, be interesting to extend this analysis to larger r s -values, which should be straightforward based on an analysis of the on-top PDF.
Finally, the momentum distribution function is of crucial importance for realistic two-component plasmas for which extensive restricted PIMC simulations, e.g. [11,24] and fermionic PIMC simulations, e.g. [12,103] have been performed. Therefore, an extension of the present analysis of the on-top PDF two two-component QMC simulations if of high interest.
This will also be the basis for the application of the present results to estimate the effect of power law tails in n(k) in fusion rates, e.g. [37][38][39], and other inelastic processes, that involve the impact of energetic particles. An example for the latter are electron impact excitation and ionization rates of atoms in a dense plasma. Such effects were predicted for various chemical reactions in Ref. [35] based on a approximate treatment of collision rates and phenomenological Lorentzian-type broadening of the electron spectral function in Eq. (2). However, such approximations are known to violate energy conservation, e.g. [104]. The present approach to n(k) makes such approximations obsolete and, moreover, eliminates the multiple integrations over the energy variables in Ref. [35], substantially simplifying the expressions for the rates.
Finally, the relevance of algebraic tails of n(k) for nuclear fusion rates in dense plasmas was discussed by many authors, e.g. [33,35,36,105], but the agreement with experimental data remains open. The results of the present work are applicable to many fusion reactions of fermionic particles, such as the proton-proton or 3 He-3 He fusion reactions in the sun or supernova stars that were considered e.g. in Refs. [38,105]. For quantitative comparisons the present simulations should be extended to multicomponent electron-ion plasmas and include screening effects of the ion-ion interactions, e.g. [39], which does not pose a principal problem.

ACKNOWLEDGMENTS
This work has been supported by the Deutsche Forschungsgemeinschaft via project BO1366-15/1. TD acknowledges financial support by the Center for Advanced Systems Understanding (CASUS) which is financed by the German Federal Ministry of Education and Research (BMBF) and by the Saxon Ministry for Science, Art, and Tourism (SMWK) with tax funds on the basis of the budget approved by the Saxon State Parliament.
We gratefully acknowledge CPU-time at the Norddeutscher Verbund für Hoch-und Höchstleistungsrechnen (HLRN) via grant shp00026 and on a Bull Cluster at the Center for Information Services and High Performance Computing (ZIH) at Technische Universität Dresden.
The expectation value, Eq. (A4), and the spin density, Eq. (A6), contain products of plane wave single-particle orbitals (30) for which |ϕ i (r)| 2 = 1 V , and, due to momentum conservation, With the definition (5)  The estimator can be read off the expression in the braces, g ↑↓ 0 (C) = 1 2N σ1 (C)N σ2 (C) ijkl δ si,s l δ sj ,s k (1 − δ si,sj )d ijkl (C) where the sum can be rearranged as 1 2 ijkl g ijkl = k =i<j =l k<l δ si,s l δ sj ,s k − δ sj ,s l δ si,s k (1 − δ si,sj )d ijkl − i<j (1 − δ si,sj )d ijij , using the symmetry properties of the two-particle density matrix. The first sum is over the off-diagonal matrix elements, where the conditions of Eq. (25) apply. The latter sum is diagonal in creation and annihilation operators, and the conditions of Eq. (27) are met. Finally, we obtain, g ↑↓ 0 (C) = 1 β K ν=1 k =i<j =l k<l δ sj ν ,s lν δ si ν ,s kν − δ si ν ,s lν δ sj ν ,s kν (1 − δ siν ,sjν )w(κ ν ) The first sum extends over all kinks κ ν := (i ν , j ν , k ν , l ν ) with the proper ordering of the indices ensured by the Kronecker deltas. The second sum extends over all occupation numbers of occupied orbitals, i < j, with opposite spin projections, at all imaginary time intervals weighted by the relative extension of the time slice in imaginary time.
Appendix B: Modification of the kinetic energy by interaction effects The influence of Coulomb interaction on the kinetic energy of the warm dense uniform electron gas was studied in the main text in Sec. III A 2, see in particular Figs. 4 and 5. In this Appendix we provide tables with extensive benchmark data for the kinetic energy of the UEG compared to the kinetic energy of the ideal system, based on ab initio CPIMC simulations, for temperatures 0.0625 ≤ Θ ≤ 4 and r s ≤ 2.
TABLE II. Comparison of the kinetic energy of an ideal system with the one of an interacting system. CPIMC results using 54 particles. The errors in the 4th column are lower than 10 −5 .