Bound-State Breaking and the Importance of Thermal Exchange–Correlation Effects in Warm Dense Hydrogen

Hydrogen at extreme temperatures and pressures is of key relevance for cutting-edge technological applications, with inertial confinement fusion research being a prime example. In addition, it is ubiquitous throughout our universe and naturally occurs in a variety of astrophysical objects. In the present work, we present exact ab initio path integral Monte Carlo (PIMC) results for the electronic density of warm dense hydrogen along a line of constant degeneracy across a broad range of densities. Using the well-known concept of reduced density gradients, we develop a new framework to identify the breaking of bound states due to pressure ionization in bulk hydrogen. Moreover, we use our PIMC results as a reference to rigorously assess the accuracy of a variety of exchange–correlation (XC) functionals in density functional theory calculations for different density regions. Here, a key finding is the importance of thermal XC effects for the accurate description of density gradients in high-energy-density systems. Our exact PIMC test set is freely available online and can be used to guide the development of new methodologies for the simulation of warm dense matter and beyond.


Introduction
The properties of hydrogen under extreme temperatures and pressures are of prime relevance for numerous applications in both fundamental and applied sciences.It is ubiquitous throughout our universe and constitutes the predominant material in stars and giant planets [1,2].In addition, it is of key importance for technological applications, most notably inertial confinement fusion (ICF) [3,4].Despite its apparent simplicity, hydrogen offers a plethora of interesting physical effects at high-energy density (HED) conditions, including the infamous insulator-to-metal phase transition at high pressure [5,6] and a potential roton-type feature in the spectrum of density fluctuations in hydrogen jets [7,8].
A particularly important property of HED hydrogen is given by the subsequent breaking of molecular and atomic configurations with increasing densities.This bound state breaking plays an important role in astrophysical models, such as heat transport in Jupiter atmospheres through H 2 dissociation and recombination [9,10].Moreover, the breaking of hydrogen bound states at high pressures has been proposed as a possible mechanism of high-temperature superconductivity [11,12].Finally, we mention the importance of pressure ionization for the properties of a deuterium fuel capsule on its compression path in ICF experiments [3,13].
As a consequence, the properties of HED hydrogen are routinely probed in experiments at large research facilities such as the National Ignition Facility (NIF) [14,15], the Omega Laser Facility [16], and the Linac coherent light source (LCLS) [17].Yet, the extreme conditions make the rigorous diagnostics of such experiments challenging.Therefore, our understanding of the physical and chemical processes at HED parameters heavily relies on simulations.In practice, the combination of a thermal Kohn-Sham density functional theory (KS-DFT) [18][19][20] description of the quantum degenerate electrons with a semi-classical molecular dynamics propagation of the heavier ions constitutes the most widely used simulation method to support and explain experimental works in HED science.It is well known that the accuracy of a KS-DFT calculation decisively depends on the utilized exchangecorrelation (XC) functional, and the performance of different functionals has been analyzed extensively at ambient conditions [21].Further, it is common practice to assess the quality of different XC functionals by benchmarking against test sets based either on experiments or on highly accurate theoretical results [22].
Unfortunately, the situation is substantially more difficult at HED conditions.Firstly, the construction of the first thermal XC functionals [23][24][25] that are based on finite-temperature quantum Monte Carlo (QMC) simulations [25][26][27][28][29] is a relatively recent achievement, and the development of more advanced functionals that occupy higher rungs on Jacob's ladder of functionals [30] remains in its infancy [31][32][33].Moreover, the rigorous benchmarking of existing functionals has been hindered so far by the near complete lack of a reliable test set that would need to be based on exact simulation data for relevant properties of a real HED system.
In this work, we aim to fundamentally change this unsatisfactory situation by introducing a canonical test set for the simulation of warm dense hydrogen based on exact ab initio path integral Monte Carlo (PIMC) calculations [34,35].This is particularly motivated by ICF applications [36,37], where ignition-a net energy gain with respect to the energy that has been used to compress the fuel capsule-has been demonstrated in a recent mile stone experiment [4].Here, the hydrogen capsule was shock compressed and heated keeping electrons partially degenerate with T ≃ T F [3] (T F being the usual Fermi temperature [38]) before the final heating stage to the nuclear fusion regime takes place.At NIF, the compression of the hydrogen capsule drives it from ρ ≃ 10 −1 g/cc to ignition conditions with ρ > 10 2 g/cc.During this process, the mean-inter particle distance decreases from r s ≈ 4 to r s ≪ 1 (in atomic units).
Using our new ab initio PIMC results, we study the breaking of H and H 2 bound states upon increasing the density from r s = 4 to r s = 1 along the electronic Fermi temperature T = T F ∼ r −2 s (T = 3.13 − 50.1 eV).This has been achieved by combining information about the electron density with the reduced density gradient (RDG) to unambiguously identify the signature of a hydrogen bound state in bulk hydrogen.Moreover, we provide a generalisation of the RDG that allows us to study the interstitial electronic structure, which is of major interest for the exploration of new materials, e.g, with so-called interstitial quasiatoms [39,40].
Finally, we use our PIMC based test set to provide the first rigorous assessment of a variety of widely used XC functionals for KS-DFT calculations of warm dense hydrogen.In particular, we study the ability of functionals on the level of the local density approximation (LDA), generalized gradient approximation (GGA), and meta-GGA to capture the manifestation-and eventual breaking-of electron-proton bound states across a wide range of densities.Overall, we find that the inclusion of thermal XC effects even on the level of the LDA leads to an improved description in all cases, which has important implications for the future development of improved XC functionals for warm-dense matter (WDM) applications.Our test set is freely available online [41] and can be used to unambiguously benchmark both existing and novel tools for WDM theory.

Bound state breaking in warm dense hydrogen
For isolated molecular dimers, the dissociation process is theoretically studied by performing simulations at different values of the distance between the atoms constituting a molecule [42]; see the Supplemental Material [43] for additional information.In contrast, the compression induced breaking of bound states is caused by the increasingly close average distance of the particles in bulk hydrogen.To gain insight on how this effect works for H 2 in the HED regime, we consider a disordered configuration of protons with a single H 2 molecule in the center of the simulation cell.The corresponding PIMC results for the electronic density are shown in Fig. 1, with the left panel corresponding to the lowest considered value of the density, r s = 4.In addition to its relevance for ICF compression experiments, such dilute hydrogen can also be realized experimentally in hydrogen jets [8].The comparably strong electron-electron coupling makes it a potentially challenging test bed for different methods, and might give rise to exotic, hitherto unobserved phenomena such as the emergence of a roton-type feature in the dynamic structure factor [7]. Due to the large interatomic distance, the H 2 molecule can easily be identified with the bare eye from the electron density at r s = 4.
To observe the expected change in the H 2 bond due to compression, we shrink the entire system by re-scaling the position vectors of all atoms relative to the central H 2 molecule.Specifically, a re-scaling of the atomic position vectors by a factor of one half and one quarter leads to the decrease of the mean-inter particle distance to r s = 2 and r s = 1, respectively.Note that the distance between the two reference atoms that make up the H 2 molecule at the center of the simulation cell always remains fixed at d = 0.74 Å(white bars).
The highest considered density with r s = 1 is shown in the right panel of Fig. 1.In this case, d is comparable to the average interatomic distance, and it is clear that the molecular bond has been broken.Indeed, all electron-proton bound states are broken in this regime as a direct consequence of the large Fermi temperature of T F ≈ 50 eV [44].
A particularly interesting picture can be seen for r s = 2 (the center panel of Fig. 1), where we clearly observe the localisation of the electrons both around the central H 2 molecule and the surrounding hydrogen atoms and, at the same time, a stronger spreading of the electron density into the inter-atomic space compared to the dilute case of r s = 4.Such a situation is very typical for WDM [20], where one cannot clearly distinguish bound and free electrons [45].
Comparing the electron localisation in H 2 at r s = 2 to that at r s = 4, we observe that the electronic cloud in H 2 is somewhat more extended at the higher density.Evidently, the electronic density by itself does not provide sufficient information about the occurrence of either a molecular bond between two adjacent ions or even the formation of a proper electron-proton bound state in contrast to a screening cloud of a free electron around a nucleus.
To overcome this diagnostic bottleneck, we use the dimensionless reduced density gradient (RDG) measure [46], where n(⃗ r)) is the density, and q F (⃗ r) = (3π 2 n(⃗ r)) 1/3 it the local Fermi wavenumber.The RDG has been shown to be an effective tool to identify bonding at ambient conditions [47].Secondly, the RDG has an appealing physical meaning as the local representation of the electronic momentum [48] within quantum kinetic energy.In particular, s[n] is the ratio of the local electron momentum to twice the local Fermi momentum 2q F [n].This ratio also naturally appears in the density gradient expansion of XC functionals in KS-DFT [49][50][51] and of kinetic energy functionals in orbital-free DFT [52][53][54][55].In fact, s[n] is a key ingredient to the construction of XC functionals beyond the local density approximation.Finally, we note that the reduced density gradient is related to the ionization potential of atoms [56], which is an input quantity in multi-scale simulation models of ICF [37].In Fig. 2, we show the RDG s[n] as a function of the electronic density at r s = 4 (left), r s = 2 (center), and r s = 1 (right).In addition, we provide the corresponding distribution of s[n] within the 3D simulation cell in the subplots.Clearly, the dependence of the RDG on the density is non-monotonic.At large density values n ≫ n 0 (with n 0 = 3/(4πr 3 s ))-i.e., in close proximity to the protons-s[n] decreases with the increase in the density since the Fermi momentum scales as q F [n] ∼ n 1/3 while the density gradient amplitude |∇n| is limited by Kato's cusp condition [56].In the opposite limit of low densities, the density gradient |∇n| also decreases since the Coulomb field is effectively screened in the inter-atomic space.The RDG s[n] attains its largest values at a distance to a proton for which the local electron quantum momentum exceeds the Fermi momentum.To identify this region, we spatially resolve the top 40% values of s[n] in the subplots.Specifically, these values correspond to the data points above the horizontal dashed-line in the respective s[n] plots on the left.
At r s = 4, we observe that a sharp peak in s[n] at small n identifies the electronic shell structure of a bound state both of a hydrogen atom and, in particular, in the H 2 molecule.Indeed, this shell structure is characteristic for an isolated H 2 molecule [43], and we have confirmed this observation by performing KS-DFT simulations with the Fermi-Löwdin orbital selfinteraction correction [57][58][59][60] for a single molecule at different nuclei separations; see the Supplemental Material [43] for additional information.These calculations further substantiate our identification of the sharp peak in s[n] for low n with a signature of bound states in both H and H 2 .
Upon increasing the density to r s = 2, we observe a significant flattening of the s[n] distribution.Moreover, the maximum of s[n] is reduced by more than a factor of two, and the sharp signal of bound states disappears entirely.This constitutes a clear observation of the breaking of bound states in warm dense hydrogen, both for the disordered atoms and the single molecular configuration around the center.Further, it can be seen clearly that the distribution of s[n] in real space does not have a clear shell structure.Finally, we remark that any qualitative difference between H and H 2 disappears at r s = 2, which, too, is a direct consequence of the disappearance of the electronic bound states.In the right panel of Fig. 2, we show corresponding results for the highest considered density, r s = 1.Here the main trend is given by a further flattening of the distribution of s[n], and a further overall decrease of its magnitude.Looking at the depiction of s[n] in real space in the subplot, we find an even more clear suppression of any electronic shell structures compared to the previous case of r s = 2.
We note that the analysis presented in Fig. 2 has been carried out for a set of comparably small synthetic snapshots.A corresponding RDG analysis of real bulk hydrogen based on DFT-MD simulations of N = 500 hydrogen atoms is shown in Fig. 6 below and leads to the same trends.

Generalised dimensionless RDG
The dimensionless RDG s[n] plays an important role in the construction of XC functionals beyond LDA [46,50].To quantify the quality of thermal KS-DFT for the description of electronic density gradients, we analyse the deviation in s[n] between KS-DFT and PIMC data.In this context, we note that-as we have seen earlier-s[n] attains a maximum around the outer layer of an atom or molecule, and rapidly decays both in the interatomic region and around the protons.Hence, the usual definition of s[n] puts the focus on certain density regions, which potentially limits its utility as a benchmark for the quality of different XC functionals.To avoid this undesirable feature, we introduce a density re-weighting of the form α , thereby generalising the dimensionless RDG defined above.For example, by setting α = 1/3, we compensate for the fast decay of s[n] with density at n/n 0 ≫ 1 due to the q F ∼ n −1/3 dependence in Eq. ( 1).Exploring the behavior of s[n] (n/n 0 ) α , we find that the variation of α leads to a shift in the maximum of s[n] (n/n 0 ) α into the inter-atomic region for α < 0 and closer to protons for α > 0 compared to the usual case of α = 0.In practice, the modified RDG thus allows us to perform a scanning of the quality of density gradients across different regions.For benchmarking purposes, it is convenient to condensate the information about the generalised RDG into a single scalar number by defining the integrated measure with N being the total number of particles in the integration volume.Eq. ( 2) thus plays a central role for our assessment of a variety of XC functionals in thermal KS-DFT simulations below.In Fig. 3  to larger densities, and the RDG distribution is broadened.To get a more intuitive picture of the effects of the density re-weighting, we show the corresponding distribution of the generalised RDG in panels (c) and (d).Indeed, setting α = −1/3 reveals a hitherto hidden but remarkably rich structure of density gradients between the protons, whereas α = 1/3 puts the focus on the close vicinity of the latter.

Benchmarking XC functionals and the role of thermal effects
The integrated RDG measure S(α), cf.Eq. ( 2) above, constitutes a convenient scalar quantity to analyse the capability of different XC functionals to capture density gradients, and the rich physics they entail, in real WDM systems.We analysed the quality of eight XC functionals commonly used in material science, including the LDA as it has been parametrized by Perdew and Zunger [61], the GGA level functionals PBE [62], PBEsol [63] and the Armiento-Mattsson functional (AM05) [64], and the meta-GGA approximations TPSS [65], revTPSS [66], and SCAN [67].It is important to note that these functionals have exclusively been constructed for applications at ambient conditions, where the electrons are in their respective ground state.However, using such a functional for applications in the WDM regime can potentially lead to substantial inaccuracies, as it has been reported by independent groups [31,33,[68][69][70].
To rigorously assess the importance of thermal XC effects that are, by definition, not included into the aforementioned ground-state functionals, we also study the finite-temperature LDA (T-LDA) XC functional by Groth et al. [24].For completeness, we note that more sophisticated thermal XCfunctionals on the level of GGA and meta-GGA have been developed very recently [32,33,71].
Applying such novel functionals to the present canonical test set of warm dense hydrogen constitutes an important project for future work, which can easily be done once they are made openly available e.g.via the Libxc library of exchangecorrelation and kinetic energy functionals for DFT [72].
In the left column of Fig. 4, we compare our exact PIMC reference results (black) for S(α) to analogous KS-DFT results based on ground-state LDA (blue), PBE (green), and thermal LDA (red); the top, center, and bottom rows correspond to r s = 4, r s = 2, and r s = 1.Note that all curves are normalised to the result for S(0) of PIMC at a given r s .
Let us postpone the discussion of S(α) for now and proceed to the right column of Fig. 4 that shows the deviation of KS-DFT simulations to the PIMC results in the electron density projected along the y-axis.Evidently, all considered XC functionals exhibit a similar level of accuracy with deviations of a few percent at r s = 4 and less than one percent at r s = 1.On average, thermal LDA performs slightly better than the ground state LDA and PBE for the two highest densities, but similar to PBE at r s = 4. Therefore, we cannot clearly resolve the role of thermal XC effects from this analysis.
In contrast, our results for S(α)/S(0) that are shown in the left column of Fig. 4 reveal a substantial, systematic improvement due to the thermal LDA functional over a wide range of α-values for all considered densities.This clearly demonstrates the role of thermal effects in capturing the correct RDG.
From a physical perspective, the inverted dome shape of S(α) is a specific feature of bulk systems.In the Supplementary Material [43], we show that S(α) decays exponentially for large α in the case of isolated atoms and molecules due to the decay of the electron density at large distances to the nuclei.For the case of warm dense hydrogen that we consider in the present work, we observe that the value of α at which S(α) attains its minimum changes its sign with the breaking of bound states with increasing density.The possibility that this finding might constitute a universal sign of bound state breaking in monatomic materials constitutes an interesting route for future studies.
To get a broader and more complete picture of the performance of a variety of XC functionals on various rungs of Jacob's ladder of functional approximations, we consider the three representative cases of α = −1/3, α = 0, and α = 1/3 in Fig. 5.More specifically, we analyse the difference in S(α) between the exact PIMC benchmark results and the different KS-DFT data sets via which we visualise as a bar plot.This leads us to the following conclusions: (i) thermal LDA performs better than groundstate LDA and all other considered GGA level functionals (PBE, PBEsol, AM05); (ii) ground state PBE improves the quality of the RDG description compared to groundstate LDA; (iii) comparing the results for meta-GGA functionals with each other, we see that SCAN performs significantly better than TPSS and revTPSS; (iv) thermal LDA and SCAN provide results for the RDG of a similar quality.Despite the fact that SCAN is built on top of the ground-state LDA, the observed similarity between the quality in the RDG to that of thermal LDA might be a result of the information about the orbital kinetic energy densities in SCAN that automatically take into account the thermal spreading of the occupation numbers.
From a physical perspective, the particular importance of thermal XC effects for the accurate description of density gradients can be understood heuristically as follows.A substantial contribution to the XC free energy is given by the entropy term, which is completely absent in the zero-temperature limit.Specifically, this term contributes to the kinetic energy of the system, which is related to the density gradient.Therefore, it makes sense that thermal XC effects affect the density gradient measure more substantially than the actual single-electron density itself.
In general, we observe that the quality of the RDG from KS-DFT simulations is worse at r s = 4 compared to r s = 2 and r s = 1.This is expected since the degree of electron localization around the protons is much stronger at low densities.From Fig. 5, we see that the quality of the RDG from thermal KS-DFT with the thermal LDA and ground-state SCAN functionals is extremely good at r s = 2.This is relevant information for WDM research since WDM is often generated using metals [73] in experiments.Interestingly, the quality of the analysed thermal KS-DFT results for the RDG is somewhat worse at r s = 1 compared to r s = 2 despite the significantly less pronounced electronic structure in this regime.This is a direct and somewhat artificial consequence of the small variations in the amplitude of the RDG.While being a subtle feature that is hard to capture from a theoretical perspective, it is of near negligible importance for any physical properties of the system.This is similar to the role of the local field correction (or, equivalently, the static XC-kernel) at high densities and temperatures [74], which is not accurately captured by common dielectric theories; at the same time, physical observables such as the linear density response function are described with high accuracy even on the mean-field level in this regime.

Utility of the RDG measure for DFT-MD simulations of bulk hydrogen
To unambiguously demonstrate the general utility of the RDG for the detection of bound state breaking in warm dense hydrogen, we have carried out additional extensive simulations with N = 500 hydrogen atoms.The snapshots are generated by performing KS-DFT-MD simulations for r s = 1, r s = 2, and r s = 4 at the electronic Fermi temperature, T = T F .In contrast to the previous example, we use snapshots that are self-consistently computed using KS-DFT-MD independently for each set of parameters.The results for the dependence of the RDG on the density are presented in Fig. 6 for r s = 4 (left), r s = 2 (center) and r s = 1 (right).Clearly, we find the same qualitative trends as for the smaller synthetic snapshots analyzed in Fig. 2 above.Specifically, we observe a sharp peak of the RDG at r s = 4, which is a signal of the presence of bound states in the system.When the density parameter is decreased to r s = 2 and r s = 1, the sharply peaked structure disappears as the bound states are being broken.This is accompanied by a significant decrease of the maximum value of the RDG distribution with respect to the density.The presented results for 500 particles thus further demonstrate the utility of the RDG for the detection of bound state breaking in warm dense hydrogen.1)] with respect to density for warm dense hydrogen computed using KS-DFT data for 500 atoms using the T-LDA [24] XC functional.Snapshots for rs = 4, rs = 2, and rs = 1 were generated by self-consistent KS-DFT-MD simulations at the corresponding densities and temperatures.

Discussion
The formation and, conversely, breaking of bound states in hydrogen and heavier elements has an important role in HED physics; it is directly relevant both for experiments, e.g. in the context of ICF, and for our understanding of astrophysical objects such as Jupiter.Yet, the extreme conditions in HED experiments constitute a significant obstacle for the experimental study of this process.In this regard, ab initio simulations are indispensable to understand the physics and chemistry at these conditions.In practice, the most widely used first-principles method for such studies is KS-DFT.Having originally being developed for applications at ambient conditions, the extension of KS-DFT to high temperatures remains significantly less developed compared to the ground-state case.In particular, the accuracy of a KS-DFT simulation decisively depends on the utilized XC functional.Being of an a-priori unknown quality, new XC functionals require the rigorous benchmarking against test sets that are based on either experimental observations or higher level simulation methods.In the present work, we present the first suitable test set for a real warm-dense matter system based on exact PIMC calculations.It is freely available online [41] and can be used to benchmark both existing and newly developed tools for HED theory.
Among the presented findings, we first highlight the utility of the dimensionless RDG s[n] in identifying the pressure induced ionisation in the medium.Second, the generalized RDG s[n] (n/n 0 ) α that has been introduced in this work can be used as a tool for the study of the interstitial electronic structure.Third, our analysis clearly highlights the importance of explicitly thermal contributions to the XC functional for the description of the density gradients.It is most likely a direct consequence of the entropic contribution to the kinetic part of the total XC free energy that is completely missing from groundstate functionals.In this context, we reiterate the relevance of information about the distribution of the RDG for the further development of thermal XC functionals, for example for ICF applications.
We are convinced that our work opens up a number of new avenues for impactful future research.First and foremost, our findings will inform the development of novel functionals that are specifically designed for applications in the HED regime.In this regard, a particularly promising route is given by a new, fully nonlocal class of thermal XC functionals based on a combination of the adiabatic connection formula with the fluctuation-dissipation theorem [75].This line of research may eventually lift thermal KS-DFT calculations onto the same level of accuracy as KS-DFT calculations at ambient conditions, where its success regarding the description of real materials arguably remains unrivaled.
In addition to its demonstrated utility for the assessment of XC functionals for KS-DFT, the presented test set constitutes an unassailable benchmark for any theoretical method that is available in the warm-dense matter regime.This is particularly interesting to understand the accuracy of the fixed-node approximation in PIMC simulations.On the one hand, this restricted PIMC method [76] allows one to completely circumvent the exponential computational bottleneck due to the fermion sign problem [77], which has allowed Militzer et al. [78,79] to present extensive numerical results for a gamut of WDM systems including second-row elements and composite materials [80].On the other hand, this paradigm constitutes a de-facto uncontrolled approximation in practice, and rigorous benchmarks have been hitherto limited to the comparably simple uniform electron gas model [28,29,81]; here systematic errors of the order of ∼ 10% have been reported by Schoof et al. [29].Our new test set thus presents the first opportunity to assess the accuracy of the fixed-node approximation in PIMC for a real warm-dense matter system.
Finally, we mention that the present study of the RDG can be extended beyond hydrogen, which has a particular relevance for HED science.More specifically, the complex interplay of numerous physical effects at these conditions leads to interesting effects such as partial ionization, in particular of heavier elements.Checking the capability of the RDG to capture the nontrivial superposition of different charge states thus constitutes a natural follow-up project of this work.Indeed, this line of thought might, for example, help to shine light onto the delocalization of atomic orbitals in Be at extreme temperature and density that has very recently been observed by Döppner et al. [82] at the NIF.

PIMC simulations
We have carried out ab initio PIMC calculations to get exact results for the 3D single-electron density profile using the set-up described in Refs.[34,35].We use P = 500 imaginary-time propagators within the pair approximation [83], and the convergence with P has been carefully checked.Note that we do not impose any nodal restrictions [76] on the thermal density matrix.Therefore, our calculations are computationally involved (we use O 10 5 CPUh for a single calculation) due to the fermion sign problem [77], but exact within the given Monte Carlo error bars.The PIMC data are freely available online [41], and can be used for a variety of applications.

Thermal KS-DFT simulations
The results from thermal KS-DFT simulations were computed using the GPAW [84,85] electronic structure code.For r s = 4, we used a standard PAW setup provided by GPAW to accelerate the convergence of our results [86].For r s = 4, we set the k-points sampling to 10 × 10 × 10, the energy cutoff to 440 eV and the main simulation cell size to 15.54 a 0 (with a 0 being the first Bohr radius).For r s = 2 and r s = 1, we have performed all-electron calculations with the field of each ion being represented by a bare Coulomb interaction.This is done to avoid potential problems related to the finite effective cutoff range at small distances in PAW setups.For r s = 2 and r s = 1, we used the energy cutoff 2700 eV and k-points sampling 10 × 10 × 10.We set the number of bands for N = 14 particles to be N b = 280 at all considered densities.The main simulation cell size at r s = 2 (r s = 1) is 7.77 a 0 (3.885 a 0 ).For occupation numbers smearing we have T = 3.132 eV at r s = 4, T = 12.528 eV at r s = 2, and T = 50.112eV at r s = 1.As a convergence criterion for the selfconsistent field cycle we required the change in energy in the last three iterations to be less than 0.5 meV per electron, the change in density to be less than 10 −4 per electron, and for the eigenstates the integrated value of the square of the residuals of the KS equations to be less than 4 × 10 −8 eV 2 per electron.These simulation parameters have been rigorously tested with respect to convergence of the total energy.
For the RDG calculations of the system with N = 500 particles at the Γ point and using T-LDA [24], we set the number of bands to N b = 6000 at r s = 1 and T = 50.112eV, to N b = 3700 at r s = 2 and T = 12.528 eV, and to N b = 2700 at r s = 4 and T = 3.132 eV.We performed KS-DFT-MD simulations using the GPU implementation of the Vienna ab initio simulation package (VASP) [87][88][89].To accelerate the equilibration part of the MD simulations, we used the hybrid approach reported in Ref. [90].Accordingly, the KS-DFT-MD simulations where initialised using ionic coordinates from snapshots generated by performing orbitalfree DFT-MD using the DFTpy package [91].Equilibrated snapshots were obtained from KS-DFT-MD simulations with 480 MD steps at r s = 1, 700 MD steps at r s = 2, and 3200 MD steps at r s = 4 (with each MD time step set to 0.01 fs).The KS-DFT-MD simulations were performed using a Nose-Hoover thermostat with the coupling parameter of the system to the heat bath Q (effective mass) set to Q = 0.5.

Data Availability
The data supporting the findings of this study are available on the Rossendorf Data Repository (RODARE) [41].

Figure 1
Figure 1 Ab initio PIMC results for the electron density in warm dense hydrogen for a snapshot of 14 protons at rs = 4 (left), rs = 2 (center), and rs = 1 (right) for T = T F .In the central region of the simulation cell, two protons are positioned with a distance of d = 0.74 Å to each other (white bars); this molecular configuration is not changed for different rs.The positions of the surrounding protons are re-scaled while keeping the angular orientation towards to the center fixed.Note that the density is normalised by the mean density n 0 , and we cut out values above n = 0.2n 0 for better visibility.

Figure 2
Figure 2 Distribution of the RDG [Eq.(1)] with respect to density for the warm dense hydrogen system shown in Fig. 1.The subplots show the corresponding distribution of the RDG in real space.More specifically, we show the RDG values above the horizontal dashed line in the s[n] plot.This allows us to identify regions where density gradients are particularly important.
, we show s[n] (n/n 0 ) α for α = 1/3 and α = −1/3 computed from our exact PIMC results for hydrogen at r s = 4. Comparing these data with the results for α = 0 shown in the left subplot of Fig.2above, we find that setting α to negative values leads to a shift of the maximum of s[n] (n/n 0 ) α to smaller densities, i.e., to the inter-atomic region.Moreover, the distribution s[n] becomes substantially more narrow.In contrast, setting α = 1/3 shifts the maximum

Figure 3
Figure 3 Generalized dimensionless RDG s[n] (n/n 0 ) α at α = −1/3 and α = 1/3 for rs = 4 and T = T F .We show that the variation of α allows emphasizing the rich behaviour of the density gradient at different distances from protons.Panels (a) and (b) show the distribution of s[n] (n/n 0 ) α as a function of the density at α = −1/3 and α = 1/3, respectively.Panels (c) and (d) show corresponding results for s[n] (n/n 0 ) α in real space, including data points above the dashed horizontal lines in panels (a) and (b).

Figure 6
Figure 6 Distribution of the RDG [Eq.(1)] with respect to density for warm dense hydrogen computed using KS-DFT data for 500 atoms using the T-LDA[24] XC functional.Snapshots for rs = 4, rs = 2, and rs = 1 were generated by self-consistent KS-DFT-MD simulations at the corresponding densities and temperatures.