Laboratory advances have created new avenues of investigation into the warm dense matter (WDM) condensed-phase-state regime. Dynamic compression methods^1–5 and high-power laser ablation^6–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 are generally of 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 fs x-ray pulses have been reported.⁹ WDM produced by dynamic compression is usually in short-lived, spatially inhomogeneous states, and therefore detailed knowledge of such nonequilibrium (transient) WDM states is important.^10–13

An interesting and important category is the two-temperature case: electrons at temperature T_e and ions at T_i, with T_e ≠ T_i. Over 40 years ago, Anisimov et al.¹⁴ proposed a two-temperature model for the interpretation of laser-pulse-excited electron emission from a metal surface. Hohlfeld et al.¹⁵ proposed three regimes to describe the response of a system to subpicosecond 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.

A particularly pertinent description can be found in Ref. 18. 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–22 at “…different electronic and ionic temperatures. Such simulations correspond to separate equilibrium states of ions and electrons, but to an out-of-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-of-equilibrium state of decoupled temperatures.” Other examples of the use of models with differing ion and electron temperatures to interpret the equation of state, optical, and transport properties of WDM include those in Refs. 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 hydrogen in the universe makes precise knowledge of its behavior in WDM state conditions crucial for modeling the interiors of giant planets.²⁸ Its nuclear properties make hydrogen crucial for inertial confinement fusion experiments.^29,30 A free-electron x-ray laser pump–probe measurement of dense cryogenic hydrogen is an important example.^31–33 In those experiments, a pair of 300 fs, 92 eV pulses were separated by a time delay 0 < Δt ≤ 5 ps. Depending on the pump intensity (19 TW/cm² or 27 TW/cm²), 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.^31,32

In AIMD, Born–Oppenheimer forces from DFT drive classical point ions. The ionic configuration space is then sampled with molecular dynamics. Applications to finite-T systems should use free-energy DFT calculations,³⁴ but it is not unusual to see ground-state approximate functionals used with finite-T densities. Such was the case in Refs. 31 and 32. Although 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 hydrogen.

The most important of these three challenges concerns nuclear quantum effects (NQEs) on structure and thermodynamic properties.³⁵ NQEs are particularly significant for low-mass elements such as hydrogen and its isotopes.³⁶ 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 hydrogen or deuterium, with path integral molecular dynamics (PIMD) and PIMD-centroid approaches being frequent choices.^37–46 As reported in the literature,^37–39 NQEs have significant impacts on phase transitions between solid molecular hydrogen at megabar pressures, as well as the melting behavior of dense hydrogen.^41,42 Furthermore, NQEs play an important role in the accurate description of the electronic structure of dense hydrogen. 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.^36,43–48 Those differences cannot be described sufficiently even with the quasiharmonic approximation.³⁷

Second, and more generically, the cost of standard Kohn–Sham (KS) DFT⁴⁹ 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 is therefore the orbital-free 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 first-principles MD” scheme.⁵⁰ This 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) noninteracting free-energy functional and the ground-state local spin density approximation for the exchange-correlation (XC) free energy.⁵¹ By contrast, some state-of-the-art orbital-free noninteracting free-energy functionals have been developed recently and have been applied successfully to WDM.^52–54 These substantially improve the achievable accuracy of OFDFT calculations for WDM and 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 XC functionals evaluated with finite-T densities. This “ground-state approximation” is now known to introduce nontrivial errors in both the equation of state and the principal Hugoniot.^55–58 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. This 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 much-studied complications of the molecular H₂ to atomic H liquid–liquid transition, we consider ion temperatures and densities that, from the best available evidence,⁴¹ should, for the most part, correspond to the atomic liquid. Details are given below.

In this paper, we first briefly describe the PIMD methods for protons and finite-temperature OFDFT methods for electrons in Sec. II. Then, in Sec. III, we describe the computational details for two-temperature warm dense hydrogen. The results and a discussion are presented in Sec. IV in which the calibrating comparisons with other methods, NQEs on radial distribution functions, and equation of state are discussed in detail. Finally, we give a brief conclusion in Sec. V.

This subsection 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,⁵⁹ and 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_BT_i (k_B is the Boltzmann constant) takes the well-known discretized path integral form^60,61ZP=mP2πβℏ23NP/2∫∏s=1P∏i=1Ndri(s)×exp−β∑s=1P∑i=1N12mωP2(ri(s)−ri(s+1))2+1PV({ri}(s)),(1)where P is the number of discretized points, ωP=P/βℏ, and V({ri}(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 interionic 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⁶² 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.^63,64 The approximation becomes a true isomorphism in the limit P→∞, i.e., Z=limP→∞ZP.

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 byHeff=∑s=1P∑i=1Npi(s)22m′+Veff,(2)where m′ is the fictitious mass that determines the efficiency of sampling of the polymer configurations, pi(s) is the fictitious momentum conjugate to the position ri(s), and the effective potential isVeff=∑s=1P∑i=1N12mωP2ri(s)−ri(s+1)2+1P∑s=1PV{ri}(s).(3)The equations of motion then follow:m′r̈i(s)=−mωP22ri(s)−ri(s+1)−ri(s−1)−1P∂V{ri}(s)∂ri(s).(4)

The foregoing formulation comprises the primitive PIMD algorithm. Because the harmonic confinement coupling grows as P, this algorithm is known to have practical limitations with respect to configuration space sampling.⁶⁵ Craig and Manolopoulos⁶⁶ 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 have any impact on the generation of thermodynamic equilibrium average quantities.⁶⁷ Thus, we adopt thermostatted RPMD (TRPMD)⁶⁸ with normal mode transformation⁶⁹ for the bead coordinates to perform 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⁶¹ϵkin=32NPkBTi−∑s=1P∑i=1N12mωP2ri(s)−ri(s+1)2,(5)the potential-energy (electronic free-energy) estimatorϵpot=1P∑s=1PV{ri}(s),(6)and the entropy contribution from the nuclear configurations. Here, we neglect the difference in the ionic configuration entropy between the classical and quantum treatment of ions.

The thermodynamic estimator for the pressure p isp=NPkBTiV−13V∑s=1P∑i=1N[mωP2(ri(s)−ri(s+1))2+1Pri(s)⋅∂V({ri}(s))∂ri(s)],(7)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 some 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.

In ab initio PIMD, the interpolymeric interaction V({ri}(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.³⁴ In this, 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⁵²Ω[n]=F[n]+∫drv(r)−μn(r),(8)where v(r) is the external potential on the electrons corresponding to electron density n and μ is the chemical potential. The free-energy functional F[n] is composed of the noninteracting free energy F_s[n], the classical Coulomb repulsion energy (i.e., Hartree energy) F_H[n], and the XC free energy F_xc[n],F[n]=Fs[n]+FH[n]+Fxc[n].(9)The Hartree energy has the same form as in zero-temperature DFT,FH[n]=12∫n(r,Te)n(r′,Te)|r−r′|drdr′.(10)The noninteracting free energy isFs[n]=Ts[n]−TeSs[n],(11)where T_s[n] and S_s[n] are the noninteracting 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 differences in the kinetic free energy and entropy between the fully interacting system and the noninteracting reference system, to wit,Fxc[n]=(T[n]−Ts[n])−Te(S[n]−Ss[n])+(Uee−FH[n]),(12)where 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 noninteracting 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 noninteracting free-energy functionals T_s[n] and T_eS_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.

In contrast to 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 these, minimization of the grand canonical potential in Eq. (8) with respect to the electron density n(r, T_e) gives the Euler–Lagrange equationδTs[n]δn−TeδSs[n]δn+δFH[n]δn+δFxc[n]δn=μ−v(r).(13)The computational cost of solving this equation scales linearly with system size and is essentially temperature-independent. 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 nonempirical orbital-free generalized gradient approximations for the noninteracting free-energy functional has been explored,⁵² and two distinct noninteracting free-energy functionals, VT84F⁵³ and LKTF,^54,70 have been developed. More or less concurrently, several explicitly temperature-dependent XC functionals have also been proposed at both the local density approximation (LDA) level of refinement^71,72 and the generalized gradient approximation (GGA) level.⁷³ They have been shown to have significant thermal effects on equations of state, principal Hugoniots, and optical properties of WDM.^56,57 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.

To investigate the impact of NQEs on 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 based on the ionic thermal de Broglie wavelength λ=h/(2πmkBTi)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, and thus it provides a measure of the degree of ionic quantum nature. The densities, ion temperatures, and α values used are displayed in Fig. 1. For the density of 1 g/cm³, we chose T_i = 300 K, 1000 K, and 5000 K, corresponding to α = 0.681, 0.373, and 0.167, respectively. At 2.5 g/cm³, we used T_i = 553 K, 1845 K, and 9204 K, respectively, and at 5 g/cm³, 879 K, 2929 K, and 14 610 K, respectively.

The phase diagram of the real physical system is, of course, at local thermodynamical equilibrium, T_i=T_e. This 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 nonequilibrium (two-temperature) system gives pressures for which the equilibrium physical system is predicted to be in a molecular liquid state^41,43 or a molecular solid.⁷⁴ This distinction from the results of the equilibrium calculations does not persist at the higher temperatures.

The PIMD simulations used the i-PI code⁷⁵ combined with a locally modified version of the PROFESS package^76,77 for OFDFT calculations. In terms of Eq. (4), the discretized Trotter beads run on the potential surface generated by the harmonic interaction between the neighboring beads and the interpolymeric interaction calculated by OFDFT. The accurately parametrized finite-temperature local density XC free-energy functional KSDT⁷¹ was used. (Use of the corrKSDT LDA functional would not alter the results for the state conditions studied.⁷⁸) The nonempirical constraint-based generalized-gradient-approximation (GGA) noninteracting free-energy functional LKTF⁵⁴ was used. This is the finite-temperature extension of the noninteracting kinetic energy functional LKT.⁷⁰ 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,⁵³ an earlier nonempirical constraint-based noninteracting free-energy approximation.

We used a local pseudopotential (LPP) for all of the OFDFT calculations. Details regarding this can be found in Ref. 52. For comparison, 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 free-energy and pressure convergence.

Periodic supercells including 128–250 atoms were employed, depending on the bulk density, and the body-centered 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. 79 for details. After equilibration of 5000 MD steps, 15 000 configurations with time steps of 0.2 fs–0.3 fs were generated for taking the thermodynamical averages.

For comparison with conventional KSMD, we performed simulations using the Quantum-ESPRESSO package.⁸⁰ To achieve consistency of comparison with the orbital-free calculations, the KSDT XC free-energy functional was used. The projector augmented wave procedure was employed with a 120 Ry plane-wave energy cutoff. A sufficient number of energy bands were considered in order to make the highest occupied band energy higher than the chemical potential by at least 20k_BT_e. Γ-point Brillouin-zone sampling was used.

Insight into the accuracy of the noninteracting free-energy functionals is provided by two preliminary studies. First, we carried out KSMD 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^54,70 and VT84F⁵³ GGA noninteracting functionals and finite-temperature Thomas–Fermi functional (TTF)⁵² are shown in Table I. It can be seen that LKTF is somewhat better in terms of 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 is also important to note 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 hydrogen.⁵³

Next, we investigated the behavior of the three noninteracting functionals (LKTF, VT84F, and TTF) for the two-temperature situation, again with OFMD calculations compared with KSMD 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³, with 300 kK ≤ T_i ≤ 20 kK. Comparisons of the electronic pressure are shown in Fig. 2. The gross failure of TF is now obvious, as is the improvement of LKTF compared with 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 decreases slowly and evenly with increasing ion temperature, from about 19% at T_i = 300 K to about 11% at T_i = 20 kK.

On the basis 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 hydrogen total free energy and pressure with corresponding results from coupled electron–ion Monte Carlo (CEIMC) calculations⁸¹ 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 bohrs, 1.10 bohrs, and 1.25 bohrs. (Note that since the system is hydrogen, 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,⁸¹ 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 nonzero 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.

To assist in understanding NQEs on the structure of two-temperature hydrogen, we calculated the radial distribution functions (RDFs) from PI-OFMD and compared these with their classical counterpart OFMD results. In PIMD, the RDF can be calculated viag(r)=〈ni(r)4πni,0r2Δr〉PI,(14)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, and 〈〉_PI is the PI ensemble average.

Figure 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 the 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. This enhances the anharmonic vibrational effects, thereby lowering the melting temperature.^42,82 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 will inevitably affect the x-ray scattering properties of light-atom matter such as hydrogen under high pressures and at low temperatures.⁴⁰

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, it was shown in Ref. 52 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 has also been confirmed by comparison of KSMD using LPP results with results from KSMD using PAW; see Refs. 52 and 53. The occurrence of the atomic solid is thus diagnosed as arising from the orbital-free approximate noninteracting free-energy functional. This 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³, T_i = 300 K, and for ρ = 5 g/cm³, T_i = 879 K, i.e., α = 0.681. We also carried out KSMD calculations with standard Quantum-Espresso PAW data sets for the same conditions. Figure 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) are obviously liquid in character, while the OFMD RDFs show the characteristic solid peak sequence.

Despite these detailed differences, the RDFs provide an unmistakable semiquantitative criterion for the onset of meaningful 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 this plot is that α ≈ 0.30 is a lower bound for meaningful NQEs. Independent calculations by some of us⁴⁷ on the insulator–metal transition in hydrogen 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 semiquantitative lower bound for meaningful NQEs is universal for systems irrespective of different ionic compositions, not just for hydrogen. Obviously, this requires further investigation. The accuracy limitations of the LKTF functional used here do not alter these deductions about a roughly universal α.

We can examine the differing electron temperature effects on 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, we present in Figs. 6 and 7 comparisons of the RDFs only for α = 0.373 and 0.681, respectively. For α = 0.373 at the higher density, 5 g/cm³, 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³, 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 that of 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.

To obtain 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, this is defined asRig=1P∑s=1Pri(s)−ric,(15)where ric is the centroid of the ith polymer. Figure 8 shows changes in the radius-of-gyration distribution of the two-temperature 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 a larger quantum “size” and stronger relative quantum fluctuation of particles. This behavior is also the reason why the RDFs from quantum simulations are quite different from the classical RDFs.

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 interproton 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 neighboring 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 on the specific ion temperature and density conditions. The only exception is the case of small dimensionless size, α = 0.167. Even for this, variation of T_e has almost no impact on ionic pressures. These behaviors are a consequence of the definition of Piclassical 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 Piclassical calculated with those two temperatures, indicate the thermostat quality. In our calculations, the difference between Piclassical 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 their classical treatment. We note that the ionic configuration is significantly affected by NQEs, as shown in Fig. 3. This corresponds to a different external potential for the electrons to which the electron density must adapt. Inclusion of NQEs gives the protons a 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, and hence they should be considered when highly accurate 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 nonequilibrium process through extrapolations from those somewhat limited PIMD simulations.

As remarked at the outset, in this study, we used the TRPMD⁶⁸ method to sample the ionic configuration space of the isomorphic polymer system rather than the primitive algorithm.⁶¹ For clarity about methodological effects, in Table IV, we show the difference between the electronic pressures obtained from the primitive algorithm and from TRPMD for the case ρ = 1 g/cm³. 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 with TRPMD results. This is consistent with the known sampling adequacy of ionic configurations for the primitive algorithm.^67,83 The TRPMD technique can provide a solution to this problem with rapid, proper canonical sampling. This is simply a reminder of the effects of the choice of algorithm for the study of NQEs on equation-of-state or other thermodynamic properties of WDM.

In summary, we have performed PI-OFMD simulations for two-temperature warm dense hydrogen and have 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 currently available (LKTF) combined with the LDA finite-temperature XC functional improves upon the ground-state approximation by the inclusion of an explicit T dependence. The combination gives electronic pressures that reasonably reproduce the trends in the KSDFT calculations but are offset. However, the offset is substantially smaller than with other noninteracting free-energy functionals, especially TF, although further improvement in F_s[n] is still needed.

We have 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 that needs to be taken into account in the development of a better F_s[n] approximation.

NQEs calculated with this methodology are substantially larger than those from the CEIMC methodology. We have suggested the very different (reversed) two-temperature relationships in the two methods as the most obvious plausible source of the discrepancy. 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 ionic quantum effects alters the ionic RDFs perceptibly when the parameter α (defined as the ratio of the ionic thermal de Broglie wavelength to the mean interionic 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 that needs further investigation.

Importantly, NQEs raise both the ionic and electronic pressures. This is because quantum protons have a larger effective size than classical point ions. When highly accurate equation-of-state data for warm dense hydrogen are required, NQEs should not be ignored, especially under conditions of relatively low ionic temperature and high density. The extent to which this is true of other light elements needs to be ascertained on a case-by-case basis.

In addition, we have found that the use of the primitive algorithm of PIMD leads to an overestimation of electronic pressures of two-temperature warm dense hydrogen. Any attempt to use this primitive algorithm as a computational economy is therefore ill-advised.