Two-temperature warm dense hydrogen as a test of quantum protons driven by orbital-free density functional theory electronic forces

by orbital-free density functional theory electronic forces Dongdong Kang, a) Kai Luo, b) Keith Runge, c) and S. B. Trickey d) Department of Physics, National University of Defense Technology, Changsha, Hunan 410073, P. R. China Earth and Planets Laboratory, Carnegie Institution of Washington, 5251 Broad Branch Rd. NW, Washington, DC 20015, USA Department of Materials Science and Engineering, University of Arizona, Tucson, Arizona 85721, USA Quantum Theory Project, Department of Physics and Department of Chemistry, University of Florida, Gainesville, Florida 32611, USA


I. INTRODUCTION
Laboratory advances have created new avenues of investigation into the warm dense matter (WDM) condensed phase state regime. Dynamic compression methods [1][2][3][4][5] and highpower laser ablation [6][7][8] allow the exploration of WDM over wide temperature-pressure ranges. Typical experimental creation of WDM deposits energy predominantly into either electrons or ions, on time scales that generally are on the order of nanoseconds or shorter. In fact, for laser excitation, pulse durations can be in the femtosecond time domain. Concurrent heating and measurement of the electron subsystem temperature during 25 femtosecond x-ray pulses has been reported 9 . WDM produced by dynamic compression usually is in shortlived, spatially inhomogeneous states, thus detailed knowledge of such non-equilibrium (transient) WDM states is important [10][11][12][13] .
An interesting and important category is the twotemperature case: electrons at temperature T e , ions at T i . Over forty years ago, Anisimov et al. 14 proposed a two-temperature model for the interpretation of laser-pulse-excited electron emission from a metal surface. Hohlfeld et al. 15 proposed three regimes to describe the response of a system to suba) corresponding author:ddkang@nudt.edu.cn b) kluo@carnegiescience.edu c) krunge@email.arizona.edu d) trickey@qtp.ufl.edu picosecond laser pulses. The first regime is characterized by ballistic transport of the excited electrons, which lasts until electron-electron scattering leads to a T e that, in the second regime, is not equal to T i . In the third regime, ion-electron coupling results in an eventual thermal equilibrium state of the whole system, described by a single temperature, T e =T i . Meanwhile, Murillo et al. 16,17 predicted that heating the electrons creates an instantaneous response that leads to spontaneous ionic heating before thermal equilibration in the plasma regime.
Ref. 18 gave a particularly pertinent description. Those authors considered ". . . ions at a temperature T i much less than the measured electron temperature T e of 10 eV" treated by orbital-free density functional theory (OFDFT) driving ab initio molecular dynamics (AIMD) [19][20][21][22] at ". . . different electronic and ionic temperatures. Such simulations correspond to separate equilibrium states of ions and electrons, but to an outof-equilibrium state of the whole system. Since the electron system is treated in the simulations within density functional theory at T e , the ion-electron system is frozen in the out-ofequilibrium state of decoupled temperatures." Other examples of use of models with differing ion and electron temperatures to interpret the equation of state, optical, and transport properties of WDM include 23-27. Here we use the same two-temperature approach to treat hydrogen but with far more advanced density functional approximations (see below). The majority abundance of H in the universe makes precise knowledge of its behavior in WDM state conditions crucial for modeling giant planetary interi-ors 28 . Its nuclear properties make H crucial for inertial confinement fusion experiments 29 . A free electron x-ray laser pump-probe measurement of dense cryogenic H is an important example [30][31][32] . In those experiments, a pair of 300 fs, 92 eV pulses were separated by a time delay 0 < ∆t ≤ 5 ps. Depending upon the pump intensity (19 or 27 TW/cm 2 ), from roughly ∆t = 1 ps onward to 5 ps, the electron temperature T e (from hydrodynamic analysis) was about 1.75 eV while the ion temperature T i was about 0.35 -0.4 eV. The Spitzer equilibration time for the two temperatures to stabilize was calculated to be about 0.9 ps, a value consistent with observation. Several of the conclusions were confirmed by AIMD simulations driven by free-energy DFT 30,31 .
In AIMD, Born-Oppenheimer forces from DFT drive classical point ions. The ionic configuration space then is sampled with molecular dynamics. Applications to finite-T systems should use free-energy DFT calculations 33 but it is not unusual to see ground state approximate functionals used with finite-T densities. Such was the case in Refs. 30 and 31. Though by now the approach has been applied extensively, there are three substantial challenges to such free-energy-DFT-based AIMD simulations that are highly relevant to understanding two-temperature H.
The most important of those three challenges is nuclear quantum effects (NQEs) upon structure and thermodynamic properties 34 . NQEs are particularly significant for low-mass elements such as H and its isotopes 35 . An especially significant issue is the variation of NQEs as a system of low-mass ions traverses a wide range of state conditions. Several studies have been conducted using path integral techniques for the description of warm dense H or D, with path integral molecular dynamics (PIMD) and PIMD-centroid approaches being frequent choices [36][37][38][39][40][41][42][43][44][45] . As reported in the literature [36][37][38] , NQEs have significant impacts on phase transitions between solid molecular hydrogen at megabar pressures, as well as the melting behavior of dense H 40,41 . Furthermore, NQEs play an important role in the accurate description of the electronic structure of dense H. With the inclusion of the quantum nature of ions, the metallization, dissociation, and transport properties of hydrogen at high pressures exhibit strikingly different behaviors from what is found in conventional calculations with classical ions 35,[42][43][44][45][46][47] . Those differences cannot be described sufficiently even with the quasi-harmonic approximation 36 .
Second, and more generically, the cost of standard Kohn-Sham (KS) DFT 48 calculations scales with the cube of the number of occupied KS orbitals. As the temperature increases, the computational cost of conventional KS-DFT inexorably becomes prohibitive. A promising tool for material modeling at high temperatures therefore is the orbitalfree form of KS theory (OFDFT) because of its linear scaling of computational cost with system size. An approximately linear scaling-method that also does not depend on the system state to achieve sparsity is the "extended firstprinciples MD" scheme 49 . It replaces high-lying, comparatively low-occupation KS states with plane waves, thereby sacrificing orthogonality but retaining atomic shell structure. It has been used to attempt to delineate the validity of OFDFT but only in the context of the outmoded and fundamentally flawed Thomas-Fermi (TF) non-interacting free energy functional and the ground state local spin density approximation for the exchange-correlation free energy 50 . In contrast, some state-of-the-art OF non-interacting free-energy functionals have been developed recently and applied successfully to WDM [51][52][53] . Those substantially improve the achievable accuracy of OFDFT calculations for WDM, hence go well beyond the TF approximation that was used, for example, in Ref. 18.
Third, as remarked above, the majority of finite-T AIMD calculations have used ground-state exchange-correlation (XC) functionals evaluated with finite-T densities. This "ground-state approximation" now is known to introduce nontrivial errors in both the equation of state and principal Hugoniot [54][55][56][57] . XC thermal effects need to be included in the WDM regime by use of genuine XC free energy functionals.
The advance of OFDFT approximations for both noninteracting and XC free energies provides the opportunity to explore whether NQEs can be treated adequately by path integral simulations driven by OFDFT calculations. If successful, this would open new possibilities for large-scale, detailed simulations of important WDM systems. Here we investigate precisely that potential payoff. We study the structure and thermodynamic properties of two-temperature warm dense hydrogen using PIMD simulations for protons driven by electronic forces obtained from finite-temperature OFDFT. The effects of the OFDFT functional approximation are assessed by comparison with conventional KS-DFT calculations. With that comparison in hand, we then focus on the physics of the Hydrogen problem by systematic assessment of NQEs. That is done by comparisons between orbital-free Born-Oppenheimer molecular dynamics (which for brevity we denote as OFMD) simulations and orbital-free-driven path integral PIMD (PI-OFMD) simulations. To avoid the well-known and muchstudied complications of the molecular H 2 to atomic H liquidliquid transition, we consider ion temperatures and densities that from the best available evidence 40 should, for the most part, correspond to the atomic liquid. Details are below.
In this paper, we first briefly describe the PIMD methods for protons and finite-temperature OFDFT methods for electrons in section II. Then, in section III we describe the computational details for two-temperature warm dense hydrogen. The results and discussions are presented in section IV, in which the calibrating comparisons with other methods, NQEs on radial distribution functions (RDFs) and equation of state are discussed in detail. Finally, we give a brief conclusion.

A. Path integral molecular dynamics for protons
This sub-section provides a brief description of the PIMD treatment for protons. Note that the temperatures under consideration in this work are higher than the quantum degeneracy temperature 58 , hence proton exchange effects are neglected in the following formulation. The canonical partition function of N identical ions of mass m at inverse temperature β = 1/k B T i (k B is the Boltzmann constant) takes the wellknown discretized path-integral form 59,60 where P is the number of discretized points, ω P = √ P/β , and V({r i } (s) ) is the potential energy of the system at the imaginary time slice s for system configuration {r i }. For the two-temperature system, V is conventionally noted as the electronic free energy although the inter-ionic potential is included.
In Eq. 1, the canonical partition function of a quantum system has been expressed as a configurational integral of Boltzmann-weighted continuous paths. Those closed paths are discretized through the Trotter approximation 61 into P "beads" circularly connected via harmonic springs, Thus, as is well-known, an N-particle quantum system can be made approximately isomorphic to a classical system consisting of N ring polymers, in which each quantum particle is mapped onto a closed, flexible polymer of P beads 62,63 . The approximation becomes a true isomorphism in the limit P → ∞, i.e., Z = lim P→∞ Z P .
To sample the configuration space of the classical polymer system using MD simulations, and thereby obtain the thermodynamic properties of the isomorphic quantum system, an effective Hamiltonian is defined by where m ′ is the fictitious mass which determines the efficiency of sampling of the polymer configurations, p (s) i is the fictitious momentum conjugate to the position r (s) i , and the effective potential is The equations of motion then follow: The foregoing formulation comprises the primitive PIMD algorithm. Because the harmonic confinement coupling grows as √ P, that algorithm is known to have practical limitations with respect to configuration space sampling 64 . Craig and Manolopoulos 65 extended the algorithm to ring-polymer molecular dynamics (RPMD) for simulating the fictitious quantum dynamics approximately. In RPMD, the fictitious bead mass m ′ is chosen to be the physical mass m and the harmonic bead-coupling and potential energy terms are taken to be P times larger than their counterparts in Eq. 4. It can be shown that these choices do not impact the generation of thermodynamic equilibrium average quantities 66 . Thus we adopt the thermostatted RPMD (TRPMD) 67 with the normal mode transformation 68 for the bead coordinates to do the PIMD simulations. At the end we give a brief comparison of results from the primitive PIMD algorithm.
The ensemble averages for thermodynamic quantities can be derived from the path integral representation of the partition function in Eq. 1. In particular, the thermodynamic estimator for the total free energy is contributed by the quantum kinetic energy estimator 60 the potential energy (electronic free energy) estimator and entropy contribution from the nuclei configurations. Here, we neglect the difference of the ionic configuration entropy between the classical and quantum treatment of ions. The thermodynamic estimator for the pressure p is where V is the system volume. The first and second terms are the contributions of the ionic quantum motions to the total pressure, while the third term is the electronic contribution (in parts of the literature the electron pressure is known as the excess pressure). Averaging these estimators along the MD trajectory yields the corresponding observable thermodynamic quantities.

B. Finite-temperature orbital-free density functional theory for electrons
In ab initio PIMD, the inter-polymeric interaction V({r i } (s) ) at the imaginary-time slice s has two parts, namely the ion-ion electrostatic interaction and the electron-ion interaction energy. The latter contribution can be calculated via free-energy DFT 33 . In it, the electronic free energy is obtained by minimizing the grand canonical potential with respect to the electron density n(r, T e ). The grand canonical potential has the form 51 where v(r) is the external potential on the electrons corresponding to electron density n and µ is the chemical potential.
The Hartree energy has the same form as in zero-temperature DFT The non-interacting free energy is where T s [n] and S s [n] are the non-interacting kinetic free energy and entropy, respectively. The XC free energy is the excess electron-electron interaction energy with respect to the Hartree energy plus the difference for the kinetic free energy and entropy between the fully interacting system and the noninteracting reference system, to wit: Here U ee [n] is the full electron-electron Coulomb interaction free energy.
In conventional KS-DFT, a sophisticated scheme exploits the one-electron orbitals of the non-interacting system to construct the electron density of the real system and thereby, the total free energy. The advantage of conventional KS-DFT is that the non-interacting free energy functionals T s [n] and TS s [n] can be constructed exactly from the one-electron orbitals and their Fermi-Dirac occupations, thereby giving an explicit Euler equation once a suitable approximate F xc is provided.
Distinct from conventional KS-DFT, in OFDFT the noninteracting functionals T s [n] and S s [n] are formulated directly in terms of the electron density rather than the KS orbitals. Doing so requires approximations. Given them, minimization of the grand canonical potential in Eq. 8 with respect to the electron density n(r, T) gives the Euler-Lagrange equation The computational cost for solving this equation scales linearly with system size and is essentially temperatureindependent. The challenge is that both accurate noninteracting and XC free energy functionals are essential ingredients for successful application of OFDFT to WDM. Recently, a framework for generating constraint-based, non-empirical orbital-free generalized gradient approximations for the non-interacting free energy functional was explored 51 , and two distinct non-interacting free energy functionals, VT84F 52 and LKTF 53,69 , were developed. More or less concurrently, several explicitly temperature-dependent XC functionals also have been proposed at both the local density approximation (LDA) level of refinement 70,71 and the generalized gradient approximation (GGA) level 72 . They have been shown to have significant thermal effects on equations of state, principal Hugoniots, and optical properties of WDM 55,56 . The opportunity opened by these functionals that is explored here is to broaden both the scope and the accuracy of ab initio PIMD calculations of NQEs.

III. COMPUTATIONAL DETAILS
To investigate NQE impacts upon the material properties of a sensitive two-temperature system, we have performed extensive PI-OFMD simulations of two-temperature warm dense hydrogen. The electron temperatures studied were T e =20 kK, 50 kK, and 100 kK. The densities and ion temperatures were specified in terms of a dimensionless scale parameter α := λ/2r s with the ionic thermal de Broglie wavelength λ = h/(2πmk B T i ) 1/2 and the ionic Wigner-Seitz radius r s = (3V/4πN) 1/3 , where h is the Planck constant, m is the ionic mass, T i is the ion temperature, and V is the system volume. The parameter α is the ratio of the effective quantum size of the ions to their mean separation. Thus it provides a measure of the degree of ionic quantum nature. The densities, ion temperatures, and α values used are displayed in Fig. 1. We chose T i = 300 K, 1000 K, and 5000 K for the density of 1 g/cm 3 , corresponding in order to α = 0.681, 0.373, and 0.167. At 2.5 g/cm 3 we used T i = 553 K, 1845 K, and 9204 K, and at 5 g/cm 3 , 879 K, 2929K, and 14610 K, both with the same order of α values.
The phase diagram of the real physical system is, of course, at local thermodynamical equilibrium, T i =T e . That phase diagram therefore provides only limited context for these α choices. At T i = 300 K and the lower densities chosen, the equation of state that we find for the non-equilibrium (twotemperature) system gives pressures for which the equilibrium physical system is predicted to be in a molecular liquid state 40,42 or molecular solid 73 . That distinction with the equilibrium calculations difference does not persist at the higher temperatures.
The PIMD simulations used the i-PI code 74 combined with a locally modified version of the PROFESS package 75,76 for OFDFT calculations. In terms of Eq. 4, the discretized Trotter beads run on the potential surface generated by the harmonic interaction between the neighbor beads and the interpolymeric interaction calculated by OFDFT. The accurately parametrized finite-temperature local density XC free energy functional KSDT 70 was used. (Use of the corrKSDT LDA functional would not alter the results for the state conditions studied 77 .) The non-empirical, constraint-based generalizedgradient-approximation (GGA) non-interacting free energy functional LKTF 53 was used. It is the finite-temperature extension of the non-interacting kinetic energy functional LKT 69 . LKT, in turn, is specifically designed to meet rigorous constraints when used with pseudopotentials, as is the case here. For comparative purposes, we also used VT84F 52 , an earlier non-empirical, constraint-based non-interacting free energy approximation.
We used a local pseudopotential (  ison conventional Kohn-Sham AIMD (KSMD) calculations were done with that same LPP unless otherwise noted. In the OFDFT calculations with PROFESS, the numerical grid for real-space integrations was set to 64 × 64 × 64 to ensure the free energy and pressure convergence.
Periodic supercells including 128-250 atoms were employed depending upon the bulk density, and the bodycentered cubic configuration was used as the starting ionic structure for all the MD simulations, both conventional KSMD and OFMD, in this study.
In PIMD, the Trotter bead number was P = 16. AIMD with classical ions but with the same algorithmic and code structure was achieved by setting P = 1. The stochastic path integral Langevin equation thermostat was used to improve the sampling efficiency. See Ref. 78 for details. After equilibration of 5000 MD steps, 15000 configurations with time steps of 0.2-0.3 fs were generated for taking the thermodynamical averages.
For comparison with conventional KSMD, we did simulations using the Quantum-ESPRESSO package 79 . To achieve consistency of comparison with the OF calculations, the KSDT XC free energy functional was used. The projector augmented wave procedure was employed with a 120 Ry plane-wave energy cutoff. Enough energy bands were considered in order to make the highest occupied band energy higher than the chemical potential by at least 20k B T . Γ point Brillouin zone sampling was used.

A. Comparison with other methods
Insight as to the accuracy of the non-interacting free energy functionals is provided by two preliminary studies. First, we did KS-AIMD and OFMD calculations (no NQEs) with the same LPP in local thermodynamic equilibrium at fixed, comparatively low temperature, T e = T i = 2 kK for three densities, r s = 1.05, 1.10, and 1.25. The results for the modern LKTF 53,69 and VT84F 52 GGA non-interacting functionals and finite-temperature Thomas-Fermi functional (TTF) 51 are shown in Table I. There one sees that LKTF is somewhat better in isothermal pressure error than VT84F or TTF in the sense that the range of LKTF fractional error is smallest of the three for the density range considered. It also is important that the TTF error is deceptive in that TTF does not give bound crystals in the static lattice limit and has other failures even for H 52 .
Next, we investigated the behavior of the three noninteracting functionals (LKTF, VT84F, TTF) for the twotemperature situation, again with OFMD calculations compared to KS-AIMD calculations (no NQEs), all using the same LPP. The two-temperature warm dense hydrogen had electron temperature T e = 20 kK and bulk density ρ = 1 g/cm 3 with 300 ≤ T i ≤ 20kK. Comparisons of the electronic pressure are shown in Fig. 2. The gross failure of TF now is obvious as is the improvement of LKTF compared to VT84F. However, as in the preceding case, LKTF again underestimates the electronic pressure relative to the KSMD values. As before, the underestimate is smooth and weakly varying. It de- On the grounds of these two comparisons, we make the quite plausible assumption that this smooth offset is inconsequential for reasonable estimates of NQEs. Therefore for the remainder of this study we use LKTF for the OFMD and PI-OFMD calculations.
Next we compare the PI-OFMD quantum nuclear corrections to the H total free energy and pressure with corresponding results from coupled electron-ion Monte Carlo (CEIMC) calculations 80 in the case of thermodynamic equilibrium states. Specifically we considered T i = T e = 2000 K and three densities corresponding to ionic Wigner-Seitz radii r s =1.05, 1.10, and 1.25 bohr. (Note that since the system is H, the ionic and electronic Wigner-Seitz radii are identical.) As shown in Table II, the quantum nuclear corrections from PI-OFMD are substantially larger than those from CEIMC for both the total free energy and pressure.
In the CEIMC calculation 80 , in order to assess NQEs on the thermodynamic properties, path integral Monte Carlo was performed for the protons on the potential-energy surface defined by the zero-temperature reptation quantum Monte Carlo method. In this work, the quantum protons are driven by the on-the-fly potential obtained from OFDFT calculations of electrons at non-zero temperatures. Different treatments of electronic structure lead to the much lower NQE estimates from CEIMC. It is important to note, however, that both sets of NQE corrections (to both the free energy and pressure) increase with density. See Table II. This accords with basic expectations about NQEs (specifically, a sort of excluded volume effect) and reinforces the conclusion that the two schemes are describing how the same nuclear effect is manifested in very different physical circumstances. See below for further discussion regarding that effect in the equation of state. TABLE II. Nuclear quantum effect (NQE) corrections to the total free energy and pressure of hydrogen from CEIMC and PI-OFMD calculations at the ion temperature T i = 2000 K. In PI-OFMD, T e = 2000 K, while in CEIMC the NQE corrections are obtained based on the zero-temperature potential-energy surface 80 . The total energy includes the ionic kinetic energy and the electronic free energy, F tot = ǫ kin + ǫ pot , as defined by Eqs. (5) and (6) respectively. Corrections are defined as ∆F = (F tot − F tot,classical )/N and ∆P = P − P classical , respectively. The ratio of the pressure corrections to the pressure obtained classically ∆P/P is presented. Statistical errors are reported in parentheses as the uncertainty on the last digit.

B. Radial distribution function
To assist in understanding NQEs on the structure of twotemperature hydrogen, we calculated the RDFs from PI-OFMD and compared those to their classical counterpart OFMD results. In PIMD, the RDF can be calculated via where n i (r) is the mean number of atoms in a shell of width ∆r at distance r, n i,0 is the mean atom density, PI is the PI ensemble average. Fig. 3 shows the RDF of two-temperature hydrogen at the ion temperature and density state points considered in this study when T e = 20 kK. For the three cases that correspond to α = 0.167, the RDFs from OFMD and PI-OFMD calculations completely overlap one another. When the parameter α is increased to 0.373, the first RDF peak is broadened by inclusion of NQEs. The degree of broadening is similar in the three cases.
When α is raised to 0.681, however, the RDFs from PI-OFMD calculations exhibit distinctly different behaviors from those from OFMD calculations. The structure delivered by the OFMD is the cubic solid from which the simulation started, with a well-defined rather narrow first peak, followed by a pair of structured peaks for all three densities. For all three densities, the PI-OFMD results lower the first peak and broaden it slightly while softening and smoothing the subsequent pair of peaks into a somewhat structured single peak. At the highest density, a third peak is pushed into the picture, again modestly broadened and lowered by NQEs relative to the OFMD result. The RDF behavior confirms that in quantum simulations the protons vibrate in a larger volume around their equilibrium positions than in the classical case. That enhances the anharmonic vibrational effects, thereby resulting in lowering the melting temperature 41,81 . Since the RDF is directly related (by Fourier transformation) to the structure factor, which is a key quantity in interpreting x-ray scattering measurements, NQEs on RDF inevitably will affect the x-ray scattering properties of light-atom matter such as hydrogen under high pressures and low temperatures 39 .
That said, there is a question. Why does the atomic solid appear in the OFMD calculations for α = 0.681? Someone unfamiliar with the LPP utilized here might speculate that it is implicated. However, Ref. 51 showed that it is realistic under a rather large range of state conditions. To the extent that one trusts a hard (small core) but ground-state PAW for extreme state conditions, the validity of the LPP also has been confirmed by comparison of KSMD with LPP results to results from KSMD with PAW; see Refs. 51 and 52. The occurrence of the atomic solid thus is diagnosed as arising from the orbital-free approximate non-interacting free energy functional. The diagnosis is confirmed by direct comparison of KSMD and OFMD calculations with the same LPP. We used T e = 20 kK for ρ = 1 g/cm 3 , T i = 300 K, and ρ = 5 g/cm 3 ,T i = 879 K, i.e. α = 0.681. We also did KSMD calculations with standard Quantum-Espresso PAW data sets for the same conditions. Fig. 4 shows clearly that the atomic solid is a consequence of the orbital-free approximate functional. The KSMD RDFs (which are identical for LPP or PAW on the scale of the Figure) obviously are liquid in character while the OFMD RDFs show the characteristic solid peak sequence.
Despite those detailed differences, the RDFs provide an unmistakable semi-quantitative criterion for the onset of mean- ingful NQEs. In Fig. 5 we display the ratio of the first peak height from OFMD (classical nuclei) to that for PI-OFMD as a function of the parameter α. What is qualitatively and at least roughly quantitatively indicated in that plot is that α ≈ 0.30 is a lower bound for meaningful NQEs. Independent calculations by some of us 46 on the insulator-metal transition in H and deuterium corroborate this interpretation. In those other calculations, PIMD driven by accurate conventional KS-DFT calculations and by OFDFT also exhibit notable NQEs above α ≈ 0.30. Moreover, because of the ionic mass dependence of α (smaller values for heavier ions means smaller NQEs), we suggest that α ≈ 0.30 as a semi-quantitative lower bound for meaningful NQEs is universal for systems irrespective of different ionic compositions, not just for Hydrogen. Obviously that requires further investigation. The accuracy limitations of the LKTF functional used here do not alter these deductions about a roughly universal α.
We may examine the differing electron temperature effects on the ionic structural properties by comparing the proton RDFs at the three electron temperatures T e = 20 kK, 50 kK, and 100 kK. Since the RDFs corresponding to the parameter α = 0.167 are indistinguishable, in Figs. 5 and 7 we present comparisons of the RDFs only for α = 0.373 and 0.681, respectively. For α = 0.373 at the higher density, 5 g/cm 3 , there are only slight RDF changes as T e is increased from 20 kK to 100 kK. The main effect is from NQEs, which, irrespective of T e , lower and broaden the first RDF peak modestly. At 1 g/cm 3 , the height of the first peak increases just enough to be perceptible as T e grows to 100 kK but still the stronger influence is NQEs. For α = 0.681, Fig. 7, the outcome is similar but more pronounced. Only slight changes in the RDF appear over the entire T e range, whereas the NQEs are strong enough to wash out the crystalline order that OFMD produces. In short, because even the smallest out-of-equilibrium temperature ratio is very large, T e /T i ≥ 6.8, the calculated RDFs are insensitive to growth in that out-of-equilibrium ratio.
As a further insight into NQEs, it is helpful to consider the spatial spread of the discretized path. A useful measure is the radius of gyration of the ith polymer. In PIMD that is defined   Fig. 6 for state points of (5 g/cm 3 , 879 K) and (1 g/cm 3 , 300 K), α = 0.681. as where r c i is the centroid of the ith polymer. Fig. 8 shows changes of the radius of gyration distribution of the twotemperature hydrogen at T e = 20 kK. Given a fixed density, the mean radius of gyration increases as the ion temperature is decreased, with great broadening at the lowest temperatures. This distinct temperature effect clearly shows the nature of quantum delocalization of protons, i.e., the lower temperature results in the larger quantum "size" and stronger relative quantum fluctuation of particles. That behavior is also the reason that the RDFs from quantum simulations are quite different from the classical RDFs.

C. Equation of state
The ionic, electronic, and total pressures calculated from OFMD and PI-OFMD simulations are summarized in Table III. The most obvious outcome is that the ionic pressure from quantum ions (PI-OFMD calculations) is significantly larger than from classical ions (OFMD calculations). The disparity ∆P i /P i grows as the ratio of quantum proton size to the mean inter-proton distance, i.e., the parameter α, is increased. This unsurprising behavior is a consequence of the inclusion in the PI-OFMD simulations of the quantum kinetic energy contribution to the ionic pressure (via the harmonic interactions between neighbor beads in the classical isomorphic polymer).
Observe that for a fixed α value, ∆P i /P i stays nearly constant with respect to T e , with only small variations. ∆P i /P i thus depends on the degree of quantum delocalization of protons rather than the specific ion temperature and density conditions. The only exception is the small dimensionless size, α = 0.167, case. Even for it, variation of T e has almost no impact on ionic pressures. These behaviors are a consequence of  the definition of P classical i as the ideal gas pressure for temperature equal to the MD-run average of the instantaneous thermostatted ionic temperature. For diverse technical reasons, that average temperature may differ slightly from the nominal (specified) ionic temperature. As a result, the ideal gas pressures calculated with the average and the nominal ionic temperatures may differ slightly. Close values of the average and nominal T i , or, equivalently, close values of P classical i calculated with those two temperatures, indicate the thermostat quality. In our calculations the difference between P classical i calculated from MD directly and the ideal gas pressure calculated with the nominal T i does not exceed 2%. We note that for α = 0.167 there are a few anomalous cases with negative pressure increments. Those result from statistical fluctuations related to the closeness of ionic pressures obtained from OFMD and PI-OFMD simulations.
More interestingly, the electronic pressures from quantum simulations are larger than those from classical treatments of the protons. This increase of P e by NQEs is below 5% at the temperature and density conditions considered here. This behavior should be attributed to the alteration of the electron density for quantum simulations of protons relative to the classical treatment of them. We note that the ionic configuration is significantly affected by NQEs as shown in Fig. 3. That corresponds to a different external potential for the electrons to which the electron density must adapt. Inclusion of NQEs gives the protons larger effective size than in the classical case, so the effective volume available for electrons is reduced and the electronic pressure goes up.
While the ionic pressure is dramatically raised by NQEs, the increment of total pressure with quantum protons is below 10% at the temperature and density conditions in this study. Nevertheless, the quantum nuclear corrections to the total pressure of two-temperature warm dense hydrogen are not negligible, hence should be considered when high-accuracy equation of state data are required, especially for low ion temperature and high density conditions. The quantum nuclear corrections to the free energy are shown in Fig. 9. First, note that, similar to the pressure corrections, the differing electron temperature has almost no impact on free-energy corrections. More importantly, the free energy corrections for protons increase approximately linearly with increasing density for all three cases, α = 0.167, 0.373, and 0.681. Therefore, the corrections to the free energy can be estimated in the entire non-equilibrium process through extrapolations from those somewhat limited PIMD simulations.
As remarked at the outset, in this study we used the TRPMD 67 method to sample the ionic configuration space of   the isomorphic polymer system rather than the primitive algorithm 60 . For clarity about methodological effects, in Table IV we show the difference of electronic pressures obtained from the primitive algorithm and TRPMD for the case ρ = 1 g/cm 3 . For the simulations with the number of beads P = 1, both algorithms yield almost the same results. But with the number of beads P = 16, we find that the electronic pressure is overestimated significantly by the primitive algorithm as compared to TRPMD results. This is consistent with the known sampling adequacy of ionic configurations for the primitive algorithm 66,82 . The TRPMD technique can provide a solution to this problem with rapid, proper canonical sampling. This simply is a reminder of the effects of algorithmic choice for study of NQEs on either equation of state or other thermodynamic properties of WDM.

V. CONCLUSION
In summary, we have performed PI-OFMD simulations for two-temperature warm dense hydrogen and systematically estimated nuclear quantum effects upon its structure and thermodynamic properties. The use of OFDFT with PIMD has provided a less computationally intensive option than conventional KSMD for studying two-temperature effects. The best OFDFT noninteracting free energy functional (LKTF) currently available combined with the LDA finite-temperature exchange-correlation functional used improve upon the ground state approximation by inclusion of explicit T-dependence. The combination gives electronic pressures which reasonably reproduce the trends in the KSDFT calculations but are offset. The offset is substantially smaller than with other non-interacting free energy functionals, especially TF. Further improvement in F s [n] is still needed.
We identified an unexpected limitation of LKTF, namely the appearance of spurious solid-like character in the RDF at the largest α = 0.681. Exploration of the cause is a consideration for development of a better F s [n] approximation.
NQEs calculated with this methodology are substantially larger than those from CEIMC methodology. We proposed the very different, (reversed) two-temperature relationship in the two methods as the most obvious plausible source of the difference. The use in CEIMC of clamped nuclei to generate the T e = 0 electronic potential surface would in principle amplify the difference. That proposed identification is a matter for further methodological exploration. What is learned unequivocally from the two very different sets of circumstances and methods is the ubiquitous importance of NQEs for the equation of state.
Inclusion of the ionic quantum effects alters the ionic RDFs perceptibly when the parameter α (defined by the ratio of the ionic thermal de Broglie wavelength to the mean inter-ionic distance) is as large or larger than about 0.30. This interpretation is confirmed by examination of first peak heights in the RDFs.
A significant physical finding is that variation of T e over a large range has only small to negligible effects on the ionic RDFs in the range of densities considered. This may be an indication of non-Born-Oppenheimer effects. Another possible reason is that the electronic structure does not change much at T e lower than the Fermi temperature. This is another issue about which investigation is needed.
Importantly, NQEs raise both the ionic and electronic pressures. This is because quantum protons have larger effective size than classical point ions. When high-accuracy equation of state data for warm dense hydrogen are required, NQEs should not be set aside especially at relatively low ionic temperature and high density conditions. The extent to which this is true of other light elements needs to be ascertained on a case-by-case basis.
In addition, we find that the primitive algorithm of PIMD induces the overestimation of electronic pressures of twotemperature warm dense hydrogen. Any attempt to use the primitive algorithm as a computational economy is ill-advised.