Ab initio path integral Monte Carlo simulations of warm dense two-component systems without fixed nodes: structural properties

We present extensive new \emph{ab initio} path integral Monte Carlo (PIMC) results for a variety of structural properties of warm dense hydrogen and beryllium. To deal with the fermion sign problem -- an exponential computational bottleneck due to the antisymmetry of the electronic thermal density matrix -- we employ the recently proposed [\textit{J.~Chem.~Phys.}~\textbf{157}, 094112 (2022); \textbf{159}, 164113 (2023)] $\xi$-extrapolation method and find excellent agreement with exact direct PIMC reference data where available. This opens up the intriguing possibility to study a gamut of properties of light elements and potentially material mixtures over a substantial part of the warm dense matter regime, with direct relevance for astrophysics, material science, and inertial confinement fusion research.


I. INTRODUCTION
Warm dense matter (WDM) is an extreme state that is characterized by the simultaneous presence of high temperatures, pressures, and densities [1].In nature, WDM occurs in a host of astrophysical objects, including giant planet interiors [2], brown dwarfs [3], and the outer layer of neutron stars [4].In addition, WDM plays an important role for technological applications such as hotelectron chemistry [5], and the discovery and synthesis of exotic materials [6][7][8].A particularly important example is given by inertial confinement fusion (ICF) [9], where the fuel capsule (typically a mixture of the hydrogen isotopes deuterium and tritium) has to traverse the WDM regime in a controlled way to reach ignition [10].In fact, the recently reported [11] net energy gain from an ICF experiment has opened up the intriguing possibility to further develop ICF into an economically viable and sustainable option for the production of clean energy [12].
Yet, the rigorous theoretical description of WDM constitutes a notoriously difficult challenge [1,[13][14][15], as it must capture the complex interplay of Coulomb coupling and nonideality, quantum degeneracy and diffraction, as well as strong thermal excitations out of the ground state.From a theoretical perspective, these conditions are conveniently characterized in terms of the Wigner-Seitz radius (also known as density parameter or quantum coupling parameter ) r s = (3/4πn) 1/3 , where n = N /Ω is the electronic number density, and the degeneracy temperature Θ = k B T /E F , with E F the usual Fermi energy [16].In the WDM regime, one has r s ∼ Θ ∼ 1 [17], which means that there are no small parameters for an expansion [1].
In this situation, thermal density functional theory (DFT) [18] has emerged as the de-facto work horse of WDM theory as it combines a generally manageable computation cost with an often acceptable level of accuracy.
At the same time, it is important to note that DFT crucially relies on external input i) to provide the required exchange-correlation (XC) functional, ii) as a benchmark for its inherent approximations.At ambient conditions where the electrons are in their ground state, the availability of highly accurate ab initio quantum Monte Carlo (QMC) results for the uniform electron gas (UEG) [19][20][21][22] was of paramount importance to facilitate the arguably unrivaled success of DFT with respect to the description of real materials [23].Yet, it is easy to see that thermal DFT simulations of WDM require as an input a thermal XC-functional that depends on both density and the temperature [1,24,25], the construction of which, in turn, must be based on thermal QMC simulations at WDM conditions [26].
In this context, the ab initio path integral Monte Carlo (PIMC) method [27][28][29] constitutes a key concept.On the one hand, PIMC is, in principle, capable of providing an exact solution to a given quantum many-body problem of interest for arbitrary temperatures; this has given insights into important phenomena such as superfluidity [27,30,31] and Bose-Einstein condensation [32,33].On the other hand, PIMC simulations of Fermi systems such as the electrons in WDM are afflicted with the notorious fermion sign problem (FSP) [34][35][36].It leads to an exponential increase in the required compute time with increasing the system size N or decreasing the temperature T , and has prevented PIMC simulations over a substantial part of the WDM regime.
As a consequence, a number of strategies to attenuate the FSP in PIMC simulations have been discussed in the literature [13,37,39,45,[68][69][70][71][72][73][74][75][76][77][78].Note that complementary QMC methods like configuration PIMC [13,38,44,[79][80][81] and density matrix QMC [41,42,82] lie beyond the scope of the present work and have been covered elsewhere.Three decades ago, Ceperley [72] suggested a reformulation of the PIMC method, where the sampling of any negative contributions can be avoided by prohibiting paths that cross the nodal surfaces of the thermal density matrix.This fixed-node approximation, also known as restricted PIMC (RPIMC) in the literature, is both formally exact and sign-problem free.Unfortunately, the exact nodal surface of an interacting many-fermion system is generally a-priori unknown, and one has to rely on de-facto uncontrolled approximations in practice.For the warm dense UEG, Schoof et al. [38] have reported systematic errors of ∼ 10% in the XC-energy at high density (r s = 1), whereas Dornheim et al. [83] have found better agreement for the momentum distribution for r s ≥ 2. In addition, the nodal restrictions on which RPIMC is based break the usual imaginary-time translation invariance and, thus, prevent the straightforward estimation of imaginary-time correlation functions.In spite of these shortcomings, the RPIMC method has been successfully applied by Militzer and coworkers [73,74,84,85] to a variety of WDM systems, and their results form the basis for an extensive equation-of-state database [86].
A different line of thought has been pursued by Takahashi and Imada [28], who have suggested to use inherently anti-symmetric imaginary-time propagators, i.e., determinants.This leads to the blocking (grouping together) of positive and negative contributions to the partition function into a single term, thereby alleviating the sign problem.Similar considerations have been used by Filinov et al. for PIMC simulations of the UEG [87], warm dense hydrogen [88], and exotic quark-gluon plasmas [89].In practice, this strategy works well if the number of imaginary-time propagators P is small enough so that the thermal wavelength of a single imaginary-time slice λ ϵ = √ 2πϵ, with ϵ = β/P , is comparable to the average interparticle distance r.The key point is thus to combine the determinant with a high-order factorization of the density matrix [90], which allows for suffi-cient accuracy even for very small P .This is the basic idea of the permutation blocking PIMC (PB-PIMC) method [13,39,45,69,91,92], which has very recently been applied to the simulation of warm dense hydrogen by Filinov and Bonitz [93].While being arguably less uncontrolled than the RPIMC method, the PB-PIMC idea has three main bottlenecks: i) even though the FSP is significantly attenuated, the method still scales exponentially with β and N ; ii) the convergence with P is often difficult to ensure, especially at low T where this would be most important; iii) it is difficult to resolve imaginarytime properties due to the necessarily small number of imaginary-time slices.
A third route has been recently suggested by Xiong and Xiong [75,76] based on the path integral molecular dynamics simulation of fictitious identical particles that are defined by the continuous spin-variable ξ ∈ [−1, 1] [cf.Eq. ( 2) below].To be more specific, they have proposed to carry out simulations in the sign-problem free domain of ξ ≥ 0, and to subsequently extrapolate to the fermionic limit of ξ = −1.This ξ-extrapolation method has subsequently been adapted to PIMC simulations of the warm dense UEG by Dornheim et al. [77,78], who have found that it works remarkably well for weak to moderate levels of quantum degeneracy.Although such an extrapolation is empirical, it combines a number of strong advantages in practice.First and foremost, only effects due to quantum statistics have to be extrapolated.These are known to be local in the case of fermions [94], which means that it is possible to verify the applicability (or lack thereof) of the ξ-extrapolation method for relatively small systems (e.g.N ∼ 4, ..., 14) before applying it to larger numbers of particles where no direct check is possible [77].Second, the method has no sign-problem, which allows one to simulate very large numbers of electrons [78].Finally, it gives one access to all observables that can be computed with direct PIMC, including the ITCF.This has very recently allowed Dornheim et al. [95] to compare extensive PIMC simulations of warm dense Be to XRTS experiments carried out at the National Ignition Facility (NIF) [96], resulting in an unprecedented consistency between theory and experiment without the need for any empirical parameters.
In the present work, we report a detailed study of the application of the ξ-extrapolation method to warm dense two-component plasmas, focusing on various structural characteristics of hydrogen and beryllium.The paper is organized as follows.In Sec.II, we introduce the relevant theoretical background, including a definition of the all-electron Hamiltonian governing the entire twocomponent system (II A), a brief introduction to the ab initio PIMC method (II B) and a subsequent overview of the ξ-extrapolation method (II C).In Sec.III, we present our extensive new simulation results, starting with the fermion sign problem (III A); interestingly, it is substantially more severe for two-component plasmas than for the UEG at the same conditions, reflecting the more complex physics in a real WDM system.In Secs.III B and III C, we present a detailed analysis of hydrogen at the electronic Fermi temperature (Θ = 1) at a metallic (r s = 2, ρ = 0.34 g/cc, T = 12.53 eV) and solid (r s = 3.23, ρ = 0.08 g/cc, T = 4.80 eV) density.In Sec.III D, we consider the substantially more complex case of strongly compressed (r s = 0.93, ρ = 7.49 g/cc) Be at T = 100 eV (Θ = 1.73), which is relevant e.g. for experiments at the NIF [95][96][97].Remarkably, the ξ-extrapolation method is even capable of reproducing the correct interplay of XC-effects with double occupation of the atomic K-shell as they manifest in observables such as the spin-resolved pair correlation function.Finally, we consider the spatially resolved electronic density in the external potential of a fixed ion snapshot in Sec.III E. The paper is con-cluded by a summary and outlook in Sec.IV.

II. THEORY A. Hamiltonian
Throughout this work, we restrict ourselves to the fully unpolarized case where N ↑ = N ↓ = N /2 with N being the total number of electrons in the system.Moreover, we consider effectively charge neutral systems where N = ZN atom with Z being the nuclear charge and N atom the total number of nuclei.The corresponding Hamiltonian governing the behaviour of the thus defined twocomponent plasma then reads where r and Î denote electron and nucleus position operators, and we assume Hartree atomic units (i.e., m e = 1) throughout this work.The pair interaction is given by the usual Ewald potential, where we follow the convention given by Fraser et al. [98].

B. Path integral Monte Carlo
The ab initio PIMC method constitutes one of the most successful tools in statistical physics and quantum chemistry.Since more detailed introductions to PIMC [27][28][29] and its algorithmic implementation [83,99] have been presented in the literature, we restrict ourselves to a brief summary of the main ideas.As a starting point, we consider the canonical (i.e., number of particles N , volume Ω, and inverse temperature β = 1/k B T are fixed) partition function in coordinate representation, where R contains the coordinates of all electrons and nuclei.The double sum over all possible permutations of coordinates of the spin-up and spin-down electrons is required to properly realize the fermionic antisymmetry of the thermal density matrix, where the exponent N pp corresponds to the number of pair permutations that is required to realize a particular combination of permutation elements σ N↑ and σ N↓ .In general, the cases ξ = 1, ξ = 0, and ξ = −1 correspond to Bose-, Boltzmann-, and Fermi-statistics, with the latter applying to the electrons studied in the present work.We note that we treat the nuclei as distinguishable Boltzmann particles (i.e., boltzmannons), which is appropriate for the conditions studied in the present work.Further note that this does still take into account nuclear quantum effects due to the finite extension of the nucleonic paths, which are however small (∼ 0.1%) in the warm dense matter regime even for hydrogen.
The problem with Eq. ( 2) concerns the evaluation of the matrix elements of the density operator ρ = e −β Ĥ that is not directly possible in practice, because the potential and kinetic contributions to the full Hamiltonian Ĥ = V + K do not commute.The usual workaround is to evoke the exact semi-group property of the density operator combined with a Trotter decomposition [100] of the latter, leading to the evaluation of P density matrices at a P -times higher temperature.In the case of electron-ion systems, an additional obstacle is given by the diverging Coulomb attraction between electrons and ions at short distances; this prevents a straightforward application of the Trotter formula, which only holds for potentials that are bounded from below.We over-come this problem by using the well-known pair approximation [27,[101][102][103] that is based on the exact solution of the isolated electron-ion two-body problem, and utilize the implementation presented in Ref. [103].In essence, the quantum many-body problem defined by Eq. ( 2) has then been mapped onto an ensemble of quasiclassical ring polymers with P segments [104]; these are the eponymous paths of the PIMC method.Alternatively, one can say that each quantum particle is now represented by an entire path of length τ = β along the imaginary time t = −i ̵ hτ with P discrete imaginarytime steps.Therefore, PIMC gives one straightforward access to a number of imaginary-time correlation functions [66,99,105,106], with the density-density correlator F ab (q, τ ) = ⟨n a (q, τ )n b (−q, 0)⟩ between species a and b being a particularly relevant example for WDM research [58][59][60][107][108][109].The PIMC formulation becomes exact in the limit of P → ∞, and the convergence with P has been carefully checked.We find that P = 200 proves sufficient for all cases studied here.
The basic idea of the PIMC method is to randomly generate all possible path configurations X (where the meta variable X contains both the coordinates of electrons and nuclei) using an implementation of the celebrated Metropolis algorithm [110].We note that this also requires the sampling of different permutation topologies, which can be accomplished efficiently via the worm algorithm that was introduced by Boninsegni et al. [99,111].In practice, we use the extended ensemble algorithm from Ref. [83] which is implemented in the ISHTAR code [112].

C. The ξ-extrapolation method
In the case of fermions, the factor of (−1) Npp leads to a cancellation of positive and negative terms, which is the root cause of the notorious fermion sign problem [34][35][36]; it results in an exponential increase in the compute time with N and β.In other words, the signal-to-noise ratio of the direct PIMC method vanishes for large systems, and at low temperatures.As a partial remedy of this bottleneck, Xiong and Xiong [75] have suggested to carry out path integral molecular dynamics simulations of fictitious identical particles where ξ ∈ [−1, 1] is treated as a continuous variable.It is straightforward to carry out highly accurate simulations in the sign-problem free domain of ξ ∈ [0, 1]; one can then extrapolate to the fermionic limit of ξ = −1 using the empirical relation ( For completeness, we note that an alternative extrapolation method has been introduced very recently [76], but the possibility of extending it to observables beyond the total energy remains a subject for dedicated future works.Subsequently, Dornheim et al. [77] have adapted the ξ-extrapolation approach for PIMC simulations of the warm dense UEG, and found that is works remarkably Snapshot of an all-electron PIMC simulation of Natom = 25 Be atoms at rs = 0.93 (ρ = 7.49 g/cc) and Θ = 1.73 (T = 100 eV).The green orbs depict the approximately pointlike nuclei, and the red-blue paths represent the quantum degenerate (i.e., smeared out) electrons.
well for weak to moderate quantum degeneracy.A particular practical advantage of this method is that it can be rigorously verified for small systems [N = O(1) − O (10)], where comparison with exact direct PIMC calculations is feasible.Then, a breakdown of the ξ-extrapolation [i.e., Eq. ( 3)] for larger systems at the same conditions is rendered highly unlikely due to the well-known local nature of quantum statistics effects for fermions [94].Other effects, such as the interplay of long-range Coulomb coupling with the quantum delocalization of the electrons is still fully taken into account by the PIMC simulations in the sign-problem free domain of ξ ∈ [0, 1] for larger systems.This is illustrated in Fig. 1, where we show a snapshot of a PIMC simulation of N atom = 25 Be atoms at r s = 0.93 (ρ = 7.49 g/cc) and Θ = 1.73 (T = 100 eV).The green orbs represent the Be nuclei, which are indistinguishable from classical point particles at these conditions, even though this assumption is not hardwired into our set-up, see Sec.II B above.The red-blue paths represent the quantum degenerate electrons, with their extension being proportional to the thermal wavelength λ β = √ 2πβ.Thus, effects such as quantum diffraction are covered on all length scales, which has recently allowed four of us to present unprecedented simulations of the warm dense UEG with up to N = 1000 electrons [78].Furthermore, the ξ-extrapolation method has been used in Ref. [95] to carry out extensive PIMC simulations of strongly compressed Be (ρ = 7.5−30 g/cc), resulting in excellent agreement with XRTS measurements taken at the NIF [96].Here, we present a much more detailed tech- Dashed black: exponential fit to the data at rs = 3.23.The UEG data for rs = 2 are partially taken from Ref. [114].
nical analysis of this approach, and apply it to hydrogen and beryllium for a set of different conditions.

III. RESULTS
All PIMC results that are presented in this work are freely available in an online repository, see Ref. [113].

A. Fermion sign problem
The FSP constitutes the main bottleneck for PIMC simulations of WDM systems.It is a direct consequence of the cancellation of positive and negative contributions to the partition function Eq. ( 2) due to the antisymmetry of the thermal density matrix under the exchange of particles.This cancellation is conveniently characterized by the average sign S that is inversely proportional to the relative uncertainty of a given observable [34].In Fig. 2, we show the dependence of S on the number of electrons N based on PIMC calculations.The blue diamonds and yellow squares have been obtained for the UEG at Θ = 1 and r s = 2 and r s = 3.23, respectively.For an ideal Fermi gas, the degree of quantum degeneracy would be exclusively a function of Θ and the results independent of the density.For the UEG, r s serves as a quantum coupling parameter [16].The stronger coupling at the solid density thus leads to an effective separation of the electrons within the PIMC simulation, leading to a decreased probability for the formation of permutation cycles [115].In other words, a sparser electron gas is less quantum degenerate than a dense electron gas, resulting in the well-known monotonic increase of S with r s for the UEG.In addition, we find the expected exponential decrease of the sign with the number of electrons that is well reproduced by the exponential fit to the yellow squares, and which is ultimately responsible for the exponentially vanishing signal-to-noise ratio for direct PIMC simulations of fermions.In practice, simulations are generally feasible for S ∼ O (10 −1 ).
Let us next consider the red circles, which show PIMC results for hydrogen at r s = 2 (and Θ = 1).Evidently, the average sign is systematically lower compared to the UEG at the same conditions.This is a direct consequence of the presence of the protons, leading to local inhomogeneities.To be more specific, electrons tend, on average, to cluster around the protons; this even holds at conditions where the majority of electrons can be considered as effectively unbound.The reduced average distance between the electrons then leads to an increased sampling of permutation cycles and, therefore, a reduced average sign compared to the UEG.
Finally, the green crosses show results for hydrogen at r s = 3.23.Interestingly, we find the opposite trend compared to the UEG: the average sign decreases with increasing r s .In other words, quantum degeneracy effects are even somewhat more important at lower density compared to the high-density case.This is a consequence of the more complex physics in real two-component plasmas compared to the UEG, leading to two competing trends for hydrogen.On the one hand, increasing the density compresses the electronic component, making exchange effects more important.On the other hand, reducing the density leads to the formation of H 2 molecules, where the two electrons are in very close proximity.In fact, the PIMC sampling by itself cannot directly distinguish between a spin-polarized or spin-unpolarized H 2 molecule.As we will see below, the correct predominance of the unpolarized case is exclusively a consequence of the fermionic antisymmetry that is realized by the cancellation of positive and negative contributions in PIMC, i.e., a subsequent re-weighting of the actually sampled set of configurations.In practice, this crucial effect is nicely captured by the ξ-extrapolation method, making it reliable beyond the comparably simple UEG model system.

B. Hydrogen: metallic density
In Fig. 3, we consider hydrogen at r s = 2 (ρ = 0.34 g/cc) and Θ = 1 (T = 12.53 eV).These conditions are at the heart of the WDM regime [1,13] and are realized, for example, on the compression path of a fuel capsule in an ICF experiment at the NIF [116].From a physical perspective, hydrogen is expected to be strongly ionized in this regime [85,93,117], which means that the electrons can be expected to behave similarly to the UEG at the The ξ-dependence of the electron-electron SSF See(q) for three wavenumbers and N = 14, rs = 2, Θ = 1.Symbols: PIMC data for q = 1.53 Å−1 (black stars), q = 3.06 Å−1 (blue diamonds), and q = 4, 59 Å−1 (green crosses).Solid red lines: fits according to Eq. ( 3) based on PIMC data in the FSP free domain of ξ ∈ [0, 1].same conditions.Given the excellent performance of the ξ-extrapolation method for the warm dense UEG [77,78], this regime constitutes a logical starting point for the present investigation.
In panel a), we show the electron-electron static structure factor (SSF) S ee (q) as a function of the wave number q for N = 14.The green crosses correspond to direct PIMC results for the fermionic limit of ξ = −1; in this case, the simulations are challenging due to the FSP, but still feasible, and we find an average sign [34] of S = 0.08466 (14).The red circles have been computed by fitting Eq. ( 3) to the sign-problem free domain of ξ ∈ [0, 1] (shaded grey area) and are in excellent agreement with the exact results for all q.Interestingly, the SSF exhibits a more pronounced structure in the bosonic limit of ξ = 1, which is likely a consequence of the effective attraction of two identical bosons reported in earlier works [118,119].
In Fig. 4, we show the ξ-dependence of S ee (q) for three selected wavenumbers.Specifically, the symbols show our exact direct PIMC results, which are available for all ξ at these parameters, and the solid red lines fits via Eq.( 3) that have been obtained exclusively based on input data from the sign-problem free domain of ξ ∈ [0, 1].The ξdependence is most pronounced for the smallest depicted q-value, and Eq.(3) nicely reproduces the PIMC results over the entire ξ-range in all cases, as it is expected.
In Fig. 3b), we analyze the thermal structure factor [120], defined by the imaginary-time density-density correlation function evaluated at its τ = β/2 minimum, i.e., F ee (q, β/2).From a theoretical perspective, the ITCF can be expected to depend even more strongly on quantum effects in general, and quantum statistics in particular, and thus constitutes a potentially more challenging case compared to the static S ee (q) = F ee (q, 0).Nevertheless, we observe that the ξ-extrapolation method is capable of giving excellent results for the thermal structure factor over the entire wavenumber range, just as in the previously investigated case of the UEG [77].
Let us next focus more explicitly on the effect of the nuclei (i.e., protons).To this end, we show the electronproton SSF S ep (q) and proton-proton SSF S pp (q) in panels c) and d).First and foremost, we find the same excellent agreement between the exact, direct PIMC results and the ξ-extrapolation results in both cases for all q.Second, the impact of quantum statistics is comparably reduced in particular for S pp (q), which is very similar to S ee (q) in the fermionic limit, but qualitatively differs substantially in the bosonic limit of ξ = 1.Such direct insights into the importance of quantum degeneracy effects on different observables are a nice side effect of the ξ-extrapolation method.
To focus more closely on such effects, we investigate the spin-resolved electron-electron pair correlation functions (PCFs) g ud (r) and g uu (r) in Fig. 5a) and b), respectively.In the spin-offdiagonal case, hardly any effects of quantum statistics can be resolved.While somewhat expected, this is still different from S pp (q), for which quantum statistics prove to be important even though the protons themselves are effectively distinguishable, see Sec.II B above.In stark contrast, the behaviour of the spin-diagonal PCF is predominantly shaped by quantum statistics, in particular for small separations r → 0. As FIG. 6. PIMC results for the electronic static structure factor See(q) at rs = 2 and Θ = 1 for N = 14 (diamonds), N = 32 (crosses) and N = 100 (circles).Left: impact of quantum statistics, with green, blue and red symbols showing results for ξ = −1 (extrapolated), ξ = 0, and ξ = 1; the yellow squares show direct PIMC results for N = 14 and ξ = −1.Right: comparing PIMC results for hydrogen to UEG results at the same parameters computed via the ESA [121,122].
mentioned above, bosons tend to cluster around each other, resulting in a large contact probability g ξ=1 uu (0), and a nontrivial maximum around r = 0.4 Å.For distinguishable Boltzmann quantum particles, there is still a finite contact probability of g ξ=0 uu (0) ≈ 0.5, but the bosonic maximum for small but finite r disappears.In the fermionic limit, the contact probability completely vanishes due to the Pauli exclusion principle.It is striking that the ξ-extrapolation method very accurately captures these stark qualitative differences for all r, and the red circles nicely recover the exact direct PIMC results that are shown by the green crosses.For completeness, we note that the on-top PCF g(0) constitutes a very important property even for the UEG [121][122][123][124].The applicability of the ξ-extrapolation method thus opens up the possibility to study g ee (0) and to reduce finitesize effects [123] and the possible impact of nodal errors in the case of restricted PIMC [37] by simulating larger system sizes as demonstrated in the recent Ref. [78].
In Fig. 5c) we show the electron-proton PCF g ep (r) [rescaled by a factor of r 2 ], which contains information about the electronic localization around the protons.In an atomic system, one would expect a maximum around the Bohr radius of r = 0.529 Å [93], see the dotted blue curve showing the probability density of an electron around an isolated proton at the hydrogen ground state.For this observable, quantum statistics effects are mostly restricted around intermediate distances, while the contact range between an electron and a proton is mostly unaffected.Finally, we show the proton-proton PCF g pp (r) in Fig. 5d).It is largely featureless, except for a pronounced exchange-correlation hole around r = 0. Interestingly, fermionic exchange effects play a more prominent role compared to g ep (r), and they are perfectly captured by the ξ-extrapolation at these conditions.
Let us conclude our analysis of the metallic density case with an investigation of finite-size effects.To this end, we show the electronic SSF S ee (q) in Fig. 6 for N = 14 (diamonds), N = 32 (crosses), and N = 100 (circles) hydrogen atoms.The red, blue, and green symbols in the left panel correspond to ξ = 1, ξ = 0, and ξ = −1 (extrapolated), respectively; the yellow squares show the corresponding direct PIMC results for ξ = −1 which are available only for N = 14.Apparently, no dependence on the system size can be resolved for the fermionic limit of prime interest for the present work.This is consistent both with previous findings for the UEG model [13,40,[125][126][127][128], and also with the wellknown principle of electronic nearsightedness [94].Similarly, we cannot resolve any clear dependence on N in the case of ξ = 0 despite the smaller error bars.Interestingly, this situation somewhat changes for the bosonic case of ξ = 1, where e.g. the results for N = 14 exhibit small differences compared to the other data.Heuristically, this can be attributed to a breakdown of nearsightedness for bosons, which, in the extreme case, are known for exhibiting off-diagonal long-range order e.g. in the case of superfluidity [129].It is important to note that the ξ-extrapolation still gives the correct fermionic limit despite the larger finite-size effects in the ξ > 0 data, as it becomes evident from the excellent agreement between the green diamonds and yellow squares.
The right panel of Fig. 6 shows a comparison of our new PIMC results for hydrogen with the UEG model at the same conditions; the latter has been computed from the effective static approximation (ESA) [121,122], which is known to be quasi-exact in this regime.For large q, S ee (q) approaches the single-particle limit for both hydrogen and the UEG, and they agree for q ≳ 3 Å−1 .In the long wavelength limit of q → 0, S ee (q) is described by a 0.8 0.9 parabola vanishing for q = 0; this is a well-known consequence of perfect screening in the UEG [13,16,130].In contrast, S ee (q) attains a finite value for real electronion systems that is governed by the compressibility sumrule [131].This can be explained by considering the definition of S ee (q) as the normalization of the dynamic structure factor, For the UEG, S ee (q, ω) consists of a single (collective) plasmon peak for small q.For hydrogen, this free electron gas feature is complemented by a) a contribution due to transitions between bound and free states [61] and b), more importantly in this context, a quasi-elastic feature due to effectively bound electrons and the screening cloud of free electrons [132,133].While the plasmon weight vanishes for small q, this trend does not hold for the other contributions in the case of hydrogen, leading to a finite value of S ee (0).

C. Hydrogen: solid density
Let us next consider hydrogen at solid density, i.e., r s = 3.23 (ρ = 0.08 g/cc) at Θ = 1 (T = 4.80 eV).Such conditions can be realized e.g. in experiments with hydrogen jets [134] which can be optically heated and subsequently be probed with XRTS [135].From a physical perspective, such conditions may give rise to interesting effects such as a non-monotonic dispersion relation of the dynamic structure factor at intermediate wave numbers [136,137] resembling the roton feature known e.g. from ultracold helium [138][139][140][141]. From a technical perspective, lower densities are, generally, more challenging for theoretical methods due to the increased impact of XC-effects and the larger degree of inhomogeneity [142,143].In the case of PIMC, the same trend holds w.r.t the FSP as it has been explained during the discussion of Fig. 2 above.
In Fig. 7, we show extensive new PIMC results for the ξ-extrapolation of various structural properties at Θ = 1.Panel a) corresponds to S ee (q) and is flat to a high degree.This is mainly a consequence of the electronic localization around protons.At the same time, we find that the ξ-extrapolation works with very high accuracy and reproduces the exact direct PIMC results for ξ = −1 for all q within the given level of accuracy.The same holds for the thermal structure factor F ee (q, β/2), electron-proton SSF S ep (q), and proton-proton SSF S pp (q) shown in panels b), c), and d), respectively.
In Fig. 8a), we show the spin-offdiagonal PCF g ud (r).First, we find a somewhat more pronounced impact of quantum statistics in the solid density case compared to r s = 2 that has been investigated in Fig. 3 above.Second, g ud (r) exhibits an interesting and nontrivial structure with a significant maximum around r = 0.5 Å.It is a sign of the formation of H − ions and the incipient formation of molecules in the system.Third, the ξ-extrapolation method again works well for all r.The corresponding spin-diagonal PCF g uu (r) is shown in Fig. 8b) and exhibits a strikingly different behaviour.In this case, the impact of quantum statistics is approximately 100%, and we find a pronounced exchange-correlation hill around r = 0.5 Å for ξ = 1.Nevertheless, the ξ-extrapolation method accurately captures the correct XC-hole in the fermionic limit.We note that it holds g ud (r) = g uu (r) for ξ = 0. From a physical perspective, the bosonic effective attraction leads to a clustering of spin-aligned electrons and, in this way, to the formation of bosonic molecules.This can be seen particularly well in Fig. 8d), where we show the proton-proton PCF g pp (r).It exhibits a pronounced peak around the molecular distance of r = 0.74 Å [53] that vanishes with decreasing ξ.It is entirely absent for ξ = −1, which is fully captured by the ξ-extrapolation method.Finally, panel Fig. 8c) shows the electron-proton PCF g ep (r).While it is qualitatively similar to the case of r s = 2 shown in Fig. 5 above, it exhibits a nearly flat progression for a B ≲ r ≲ 2a B ; this feature is indicative of a somewhat lower degree of electronic delocalization, as it is expected.b) thermal electron-electron structure factor Fee(q, β/2), c) electron-proton SSF Sep(q), d) proton-proton SSF Spp(q), e) spin up-down PCF g ud (r), f) spin up-up PCF guu(r).Green crosses: direct (exact) PIMC results for ξ = −1; red circles: ξ-extrapolated results [Eq.( 3)]; grey area: FSP free domain of ξ ∈ [0, 1].The ξ-dependence of the electron-electron SSF See(q) for Natom = 4, rs = 0.93, and Θ = 1.73.Symbols: PIMC data for q = 3.14 Å−1 (black stars), q = 6.29 Å−1 (blue diamonds), and q = 14.06 Å−1 (green crosses).Solid red lines: fits according to Eq. ( 3) based on PIMC data in the FSP free domain of ξ ∈ [0, 1].

D. Strongly compressed beryllium
As a final example, we consider compressed Be at r s = 0.93 (ρ = 7.5 g/cc) and Θ = 1.73 (T = 100 eV).Such conditions have been realized in experiments at the NIF [96,97], and have recently been studied with the ξ-extrapolation method in Ref. [95].In Fig. 9, we show an analysis of the familiar set of structural properties computed for N atom = 4 Be atoms (i.e., N = 16 electrons) with the usual color code.We observe the same qualitative trends as reported for hydrogen in the previous sections and restrict ourselves to a discussion of the key differences here.First and foremost, we find excellent agreement between the ξ-extrapolation and the exact direct PIMC simulation for ξ = −1 for all considered observables.This can be discerned more directly for the electron-electron SSF in Fig. 10, where we explicitly show the ξ-dependence for three selected wavenumbers q.Interestingly, S ee (q) exhibits a slightly though definite non-monotonic behaviour with a shallow minimum around intermediate q.It is a consequence of the interplay between the plasmon weight that decreases with q, and the increasing weight of the elastic feature in the long wavelength limit.In contrast, F ee (q, β/2), S eI (q), and S II (q) exhibit the same qualitative trends as reported previously for hydrogen.The spin-offdiagonal PCF g ud (r) shown in Fig. 9e) exhibits a substantial increase towards r → 0; this is a direct consequence of the presence of ions with a fully occupied K-shell at these conditions, which is more pronounced in the bosonic case.This trend is completely absent for the spin-diagonal pendant g uu (r), where g(r) = 0 holds by definition.At the same time, we point out the remarkable impact of quantum statistics for Be at these conditions with g ξ=1 uu (0) ≈ 4 for bosons.We also note that we find an average sign of S ≈ 0.11 for N atom = 4 Be atoms in our direct PIMC calculations.
Finally, we study the dependence of the electronelectron SSF on the system size in Fig. 11.Specifically, the diamonds and crosses show our PIMC results for N atom = 4 and N atom = 10, and the red, blue and green symbols in the top panel correspond to ξ = 1, ξ = 0, and ξ = −1.Overall, we find the same qualitative trends as observed for hydrogen in Fig. 6 above, namely a substantial dependence on the number of particles in the bosonic case which is absent for boltzmannons and fermions.This conclusion is further substantiated by the yellow squares that show the exact, direct fermionic PIMC simulations for N atom = 4 and nicely agree with the ξ-extrapolated results for both system sizes.
In the bottom panel of Fig. 11, we compare our PIMC results for Be with the corresponding electronic SSF for the UEG at the same conditions; it has been computed within the aforementioned ESA [121,122] and is depicted by the solid blue curve.Evidently, the Be system does not even qualitatively resemble the UEG model despite the comparably high temperature.More specifically, the localized K-shell electrons, as well as the loosely local- ized electronic screening cloud give rise to a pronounced elastic feature in the dynamic structure factor S ee (q, ω) in this regime [61,95,96].This, in turn, leads to an increasing normalization S ee (q) [cf.Eq. ( 4)] for small q despite the vanishing plasmon.

E. Strongly compressed beryllium: snapshot
Let us conclude our investigation with a brief study of the ξ-extrapolation method applied to the electronic problem in a fixed external ion snapshot potential.A similar set-up has recently allowed Böhme et al. [103,117] to present the first exact results for the electronic density response and the associated exchange-correlation kernel of warm dense hydrogen.In addition, such calculations provide direct information about the impact of the ions on the density response [144][145][146], and about the formation of bound states and molecules [53].In Fig. 12, we show PIMC results for the electronic density computed for a snapshot of N atom = 4 Be atoms at the same conditions studied in Sec.III D. These results nicely illustrate a high degree of electronic localization around the nuclei, and a substantially reduced density in the interstitial region.From the perspective of the ξ-extrapolation method that constitutes the focus of the present work, one might expect that such a spatially resolved observable constitutes a most challenging benchmark case due to the comparably larger impact of quantum statistics effects in the vicinity of the nuclei.
In Fig. 13, we show the corresponding electronic density in the y-z-plane, and the x-positions have been chosen so that panels a) and b) approximately contain two ions, whereas panels c) and d) show the interstitial region without such nuclei.Further, the left column shows exact, direct fermionic PIMC results, whereas the right column has been computed from the ξ-extrapolation via Eq.(3) using as input PIMC results from the signproblem free domain of ξ ∈ [0, 1].Most importantly, we find no significant differences between the direct and the extrapolated results in both cases.This can be seen more clearly in Fig. 14, where we show scan lines over the 2D densities.Let us first consider the left column that has been obtained for the layer that includes two ions, and the top and bottom panels correspond to the dotted and dashed lines in Fig. 13.As it is expected, we find a comparably larger localization in the bosonic limit around the nuclei.Moreover, the ξ-extrapolation works very well, although there appear small differences around the maximum of the density in the top panel.This, however, is not indicative of a statistically significant systematic underestimation of the true fermionic density in this region, and does not occur for the other nuclei (i.e., see the bottom panel).The right column shows the same analysis for the layer that is located in the interstitial region.In this case, the bosonic density is comparably reduced to the other data sets; this is required to balance the increased bosonic localization around the nuclei.Most importantly, the ξ-extrapolation fully captures the correct fermionic limit everywhere.

IV. SUMMARY AND DISCUSSION
In this work, we have presented a detailed investigation of a variety of structural properties for warm dense hydrogen and beryllium.As the first test case, we have considered hydrogen at r s = 2 and Θ = 1, where most of the electrons are assumed to be free [53,85,93].As a consequence, the system behaves qualitatively similar to the free UEG, where it is known that the ξ-extrapolation works well [77,78]; the same does indeed hold for H at these conditions.The second test case concerns hydrogen at solid density, r s = 3.23 and Θ = 1.These conditions can be created in experiments with hydrogen jets [134,135], and might give rise to a nontrivial rotontype feature at sufficiently high temperature [136,137].From a theoretical perspective, they constitute a more challenging example due to the comparably larger impact of quantum statistics that is reflected by the somewhat more severe fermion sign problem.In practice, we find that simulations based on Bose-Einstein statistics result in the formation of molecules, which are completely absent for a proper fermionic treatment of the electrons.At the same time, the ξ-extrapolation method nicely captures these qualitative differences and works equally well for all the investigated properties.The third physical system studied in this work concerns compressed Be at r s = 0.93 and Θ = 1.73; these conditions are relevant for experiments at the NIF [96,97] and have very recently been studied with the ξ-extrapolation method by Dornheim et al. [95].Unsurprisingly, the complex interplay between the more strongly charged Be nuclei and the electrons gives rise to a richer physics that is reflected by an even more pronounced impact of quantum statistics compared to hydrogen despite the larger value of Θ.This includes the partial double-occupation of the K-shell, which strongly depends on the spin-polarization of the involved electrons, giving rise to substantial differences in the spin-resolved pair correlation functions.Nevertheless, the ξ-extrapolation method works exceptionally well despite these complexities.Finally, we have used PIMC to solve the electronic problem in the external potential of a fixed configuration of N atom = 4 Be nuclei.From a technical perspective, analyzing the corresponding electronic density distribution might constitute the most challenging benchmark case that we have considered in this work as it allows us to spatially resolve the impact of quantum statistics, whereas potential systematic errors might be averaged out in aggregated properties such as the SSFs and PCFs considered before.In practice, the ξ-extrapolation method works remarkably well both in the vicinity of the nuclei, as well as in the interstitial region.All PIMC results that have been presented here are freely available online [113] and can be used to benchmark new methods and approximations.We are convinced that these findings open up a variety of potential projects for future works.While the ξ-extrapolation method has worked exceptionally well in all presently studied cases, previous findings for the warm dense UEG model [77,78] indicate that this idea breaks down for Θ < 1.The particular limits will likely strongly depend both on the density and the elemental composition.A detailed study of the applicability range of the method for various light elements and their mixtures thus constitutes an important task for dedicated future works.
A particular strength of the ξ-extrapolation method is given by its polynomial scaling with the system size, which combines a number of advantages.First, this ac- cess to comparably large systems allows one to study finite-size effects, which are small for the structural properties studied here, but may play a more important role for the computation of an equation-of-state [85,93].The latter are of direct importance for a number of applications including laser fusion [10], and are thus of considerable interest in their own right.Second, simulating large systems allows one to compute properties such as the ITCF in the limit of small q.This is important e.g. to describe XRTS experiments in a forward-scattering geometry [95], and might be needed to fully capture dynamic long-range effects in other properties [147].
An additional direction for future research is the dedicated study of density response properties.In this regard, we note the excellent performance of the ξ-extrapolation method w.r.t. the ITCF F ee (q, τ ), which can be used as input for the imaginary-time version of the fluctuationdissipation theorem [60,148] In this way, one can get direct access to the speciesresolved static density response function χ ab (q, 0) and, in this way, a set of exchange-correlation properties such as the electron-ion local field correction of different systems.Moreover, the corresponding estimation of higherorder imaginary-time correlation functions can be used as input for similar relations [66] that give one access to a variety of nonlinear response properties [63][64][65]149].
Finally, we note that the impact of quantum statistics is comparably large for Coulomb interacting systems such as in warm dense matter, but may be substantially smaller for more short-range repulsive interactions [119].This opens up the intriguing possibility to study ultracold fermionic atoms such as 3 He [138,140,150] with unprecedented accuracy.Given the above, returning to two-component systems, the ξ−extrapolation method can also improve our rather poor understanding of liquid 3 He-4 He mixtures at low temperatures.Binary bosonfermion mixtures have traditionally gathered strong interest owing to the possibility of double superfluidity and the unique correlation-driven interplay between different quantum statistics [151][152][153].Nevertheless, ab initio simulations of such systems are essentially non-existent.To our knowledge, finite temperature isotopic helium mixtures have been so far investigated within the fixed-node approximation [154][155][156][157] or with the PIMC method but neglecting fermionic exchange effects [158].The quantities of interest have been limited to the superfluid fraction, the component kinetic energies, the momentum distribution functions and the pair correlation functions.Thus, the ξ−extrapolation technique, being sign problem free inside the phase diagram region of its applicability, can be employed for the acquisition of extensive quasiexact thermodynamic and structural results as well as for the first reconstruction of collective excitation spectra.

FIG. 2 .
FIG.2.PIMC results for the dependence of the average sign S on the system size N at the electronic Fermi temperature, Θ = 1.Blue diamonds (green crosses): UEG at rs = 2 (rs = 3.23); red circles (black stars): hydrogen at rs = 2 (rs = 3.23).Dashed black: exponential fit to the data at rs = 3.23.The UEG data for rs = 2 are partially taken from Ref.[114].

FIG. 13 .
FIG. 13.Ab initio PIMC results for the electronic density in the y-z-layer.Panels a) and b) have been computed for a layer with two ions in it, whereas panels c) and d) are taken from the interstitial region.The left and right columns show exact, direct fermionic PIMC calculations and corresponding ξ-extrapolated results based on input data from the FSP free region of ξ ∈ [0, 1], respectively.The dashed and dotted horizontal lines show scan lines that are investigated in more detail in Fig. 14.

FIG. 14 .
FIG.14.Scan lines of the electronic density for the dashed (top) and dotted (bottom) horizontal lines in Fig.13.The left and right columns correspond to the layer with and without ions in it.

ACKNOWLEDGMENTS
This work was partially supported by the Center for Advanced Systems Understanding (CASUS), financed by Germany's Federal Ministry of Education and Research (BMBF) and the Saxon state government out of the State budget approved by the Saxon State Parliament.This work has received funding from the European Research Council (ERC) under the European Union's Horizon 2022 research and innovation programme (Grant agreement No. 101076233, "PREX-TREME").Views and opinions expressed are however those of the authors only and do not necessarily reflect those of the European Union or the European Research Council Executive Agency.Neither the European Union nor the granting authority can be held responsible for them.Computations were performed on a Bull Cluster at the Center for Information Services and High-Performance Computing (ZIH) at Technische Universität Dresden, at the Norddeutscher Verbund für Hoch-und Höchstleistungsrechnen (HLRN) under grant mvp00024, and on the HoreKa supercomputer funded by the Ministry of Science, Research and the Arts Baden-Württemberg and by the Federal Ministry of Education and Research.