Toward first principles-based simulations of dense hydrogen Open Access

Hydrogen was the first element that formed after the Big Bang and it has remained the most abundant species in the Universe. Its properties, therefore, determine the structure and evolution of astrophysical objects including stars, many planets, and interstellar gas clouds. These objects exist in a huge range of temperatures—spanning from a few Kelvin to billions of Kelvin—and pressures that extend from zero to trillions of atmospheres in neutron stars. Astrophysical observations have revealed extensive though indirect information about the properties of hydrogen. For example, observations of the oscillations of the sun's surface are employed to constrain its interior properties (“helioseismology”) (e.g., Refs. 1–3). More recently, a major source of knowledge has come from laboratory experiments that reach ever higher pressures. Examples include static compression using diamond anvil cells (DAC) (e.g., Ref. 4), dynamic compression using shock waves that are generated by impactors,^5,6 explosives (e.g., Refs. 7 and 8), or high-intensity lasers.⁹ Extracting data from these measurements and translating them into static, dynamic, or optical properties of hydrogen can be challenging. The accuracy of the diagnostics is usually severely limited due to the short observation time and relaxation processes, unknown intrinsic properties of the equipment, and so forth. These gaps in knowledge are often bridged with simple models for the equation of state (EOS) or hydrodynamics, the validity of which may be questionable. All of this renders the accuracy of many measurements unclear. Moreover, high-pressure experiments are challenging and costly, and available only at a limited number of laboratories. Therefore, there is a very high demand for theoretical analysis and computer experiments that would allow for an improved interpretation of experimental data but also independent and reliable predictions.

Aside from experiments that focus on basic science aspects of the properties of hydrogen at high pressure, there is a rapidly increasing interest in technological applications of hydrogen. This includes inertial confinement fusion (ICF)^10–12 and the use of hydrogen in green energy applications. In particular, the ignition campaign at the NIF has recently reported major breakthroughs. Remarkably, the fusion gain exceeded unity (not accounting for wasted energy).^13–16 Technological advances, in particular ICF, rely heavily on theoretical modeling that allows one to understand the relevant physical processes, make reliable predictions, and enable the optimization of experiments. Moreover, the direct measurement even of basic parameters such as the temperature or density is often not possible; instead, they have to be inferred indirectly from other observations, which, in turn, requires theoretical results for different observables.¹⁷ 

The theory and modeling of dense hydrogen has a long history. The hydrogen atom and molecule were the focus of the early developments of quantum mechanics, and hydrogen was the first many-body system for which the consequences of the Fermi statistics of the electrons and the Pauli principle were explored by Fowler.¹⁸ Pressure ionization of atoms was already predicted by Hund.¹⁹ At about the same time, Wigner and Huntington predicted that atomic hydrogen under high pressure would be metallic.²⁰ Since then the questions about the liquid and solid phases of hydrogen, as well as of transitions between insulating, conducting and, possibly, a superconducting phase predicted by Ashcroft,²¹ have been a driver of both experimental and theoretical progress. There have been many reports of the experimental observation of hydrogen metallization at room temperature (see Refs. 22–24), but there is no consensus that metallization has been reproducibly observed.

On the other hand, it is well accepted that the reverberating shock experiments by Weir et al.⁶ achieved the metallization of hydrogen at high temperature and pressure. The conditions are broadly consistent with those in the interior of Jupiter. There is indeed strong evidence that a thick layer of metallic hydrogen exists in the planet's interior because Jupiter has a strong magnetic field, even though the precise location of the dynamo active layer is still being debated.²⁵ 

The equation of state of hydrogen (and helium) is of crucial importance^26,27 when measurements of Jupiter's and Saturn's gravity fields by orbiting spacecraft like Juno and Cassini are interpreted.^28,29 Interior models typically start from the isentropic pressure–density relationship for hydrogen (see Fig. 1) before helium and heavier elements are introduced. The gravity measurements have shown that Jupiter's interior cannot be perfectly homogeneous and almost all models have introduced a step in composition into the planet's interior, which requires a physical justification. Earlier models thus invoked a first-order plasma phase transition (PPT) in hydrogen³⁰ until first-principles computer simulations could not confirm its existence and instead predicted a liquid–liquid phase transition (LLPT).^31–33 However, the latter is believed to occur at too low temperatures to matter for Jupiter's interior. Modern Jupiter models now invoke the phase separation of hydrogen–helium mixtures^34–36 or a dilute core^37,38 to justify an inhomogeneous interior structure. Still, a major challenge in understanding Jupiter's interior is remaining. Measurements of the planet's atmosphere have revealed a heavy element abundance of $∼$ $4%$⁠. If one assumes these measurements to represent the bulk of the planet and if one invokes common EOSs for hydrogen and helium,^39–41 then gravity measurements of the Juno spacecraft cannot be reproduced.^42–44 Various hypotheses have been proposed to reconcile this discrepancy, which include lowering the density of hydrogen,⁴⁵ increasing the planet's interior temperature,⁴⁶ or reducing the deep abundance of the heavy elements.⁴⁷ In addition to the equation of state, further thermodynamic (e.g., compressibility, Grüneisen parameter) and transport properties (e.g., electrical and thermal conductivity, viscosity) are required to model the thermal evolution (cooling) and dynamo processes (magnetic field generation) in giant planets like Jupiter and Saturn; for recent work (see Refs. 48 and 49).

The relevance to planetary science makes it important to characterize how hydrogen transitions from a molecular, insulating state to an atomic and conducting state. As noted above, an early conjecture was the plasma phase transition (PPT), predicted by Schramm⁵⁰ and later by Norman and Starostin.⁵¹ Even though no longer favored today, this issue has led to considerable development in theory and simulations. The first models were based on the thermodynamics of chemically reacting gases applied to many-particle systems of electrons, protons, atoms, and molecules (“chemical models,” e.g., Ref. 51). They were systematically extended to take into account interactions between all constituents (“nonideality effects”) by Ebeling and co-workers,^52,53 Mihalas et al.,^54–56 Chabrier and coworkers,^39,57,58 and the Rostock school around Kremp et al.^59,60 and Röpke et al.^61–63 This has led to important advances in the many-body theory of Coulomb systems in the framework of kinetic theory and the theory of fluctuations,⁶⁴ statistical physics,^65–67 nonequilibrium Green functions,^68,69 density operator theory,⁷⁰ and linear response theory.^71,72

Chemical models rely on effective interaction potentials between the considered “species.” Their density and temperature dependence are usually treated within simple approximations such as static screening and perturbation theory. With the emergence of modern computers, many-body simulations such as quantum Monte Carlo⁷³ (QMC) and density functional theory⁷⁴ (DFT) became feasible. These methods are based on the fundamental properties of electrons and nuclei and work in the “physical picture.” They avoid any artificial classification into different chemical species. QMC simulations of hydrogen and helium, based on Feynman's imaginary-time path integral representation of quantum mechanics (PIMC) were pioneered by Fosdick and Jordan^75,76 and Filinov et al.^77,78 However, when applied to fermions, these simulations, even though being free of systematic errors, are severely hampered by the fermion sign problem^79–81 (FSP) confining them to regions of low electron degeneracy. To address the FSP, Ceperley introduced the restricted PIMC (RPIMC) method,⁸² which is completely sign-problem free, as long as diagonal density matrix elements are evaluated, which is sufficient to compute the internal energy, pressure, pair correlation functions, and structure factor. However, some negative contributions need to be included when the momentum distribution is computed.⁸³ While being formally exact, RPIMC requires information about the nodal structure of the thermal density matrix, which has to be approximated in practice. This has allowed Ceperley and coworkers to perform simulations for temperatures where hydrogen is dominated by atoms and molecules.^84–88 

The second approach within the “physical picture” is Kohn–Sham (KS) DFT coupled to molecular dynamics for the ions (DFT-MD or Born–Oppenheimer MD, BOMD, also known as ab initio MD or AIMD) which have dramatically extended the scope of problems that can be studied while yielding a higher accuracy compared to simpler models that require the introduction of assumptions specific to hydrogen. While DFT-MD and RPIMC methods rely on some approximations, these do not depend on a particular material. Early DFT-MD simulations of hot, dense hydrogen were performed by Lenosky,⁸⁶ Scandolo,³¹ Desjarlais,⁸⁹ and Bonev et al.⁹⁰ 

The “first principles” character (formally exact quantum mechanical approach within the physical picture) of RPIMC and DFT simulations⁹¹ rests on two fundamental approximations: the choice of the nodal surface of the N-particle density matrix and the exchange–correlation (XC) functional, respectively. The exact versions of these quantities are, in general, not known and the accuracy of a particular choice is difficult to assess from within the methods, and it will depend on the phase of hydrogen. Of course, the ultimate accuracy arbitrator is given by experimental measurements; however, in the case of highly compressed matter, such comparisons are difficult and afflicted by large uncertainties and error bars. An instructive example provides the laser-driven shock wave measurements of liquid deuterium of Da Silva et al.^9,92 They used the Nova laser and reported a very high compression ratio, $ρ/ρ0$⁠, of up to 5.9 [cf. gray circles in Fig. 2(a)] significantly larger than predicted by the SESAME model.⁹³ This discrepancy sparked an intense discussion and led to a series of additional experiments and many theoretical investigations. Already the first two articles that reported results from first-principles simulations with RPIMC⁹⁴ and DFT-MD⁸⁶ consistently predicted a much smaller compression ratio of $∼$ 4.3 at 200 GPa pressure (right panel of Fig. 2). It was argued that the uncertainty in the simulation results was smaller than an energy increase in $∼3 eV$/atom that would be required to bring the simulation results in agreement with the Nova measurements. That conclusion was supported by many following articles that reported simulation results. Using chemical models, on the other hand, it is possible to find agreement with the Nova measurements. Examples are the linear mixing model of Ross⁹⁵ and the equation of state of Saumon et al.,³⁹ as shown in Fig. 2. Note that both models predated the experiments and were not fit to the results.

This discrepancy and the resulting uncertainty in the EOS led to a number of new measurements, such as the Z-pinch experiments by Knudson et al.¹¹² and shock experiments at the Omega laser^97,98 [see Fig. 2(a)]. The difference in maximum compression (4- to 4.5-fold¹¹² vs 5- to 6-fold^9,92,100) is striking. Interestingly, the theoretical results may also be sorted into two groups: high compression, close to a maximum of 6 has been predicted by chemical models, such as Ebeling's Padé formulas,⁵³ in addition to the two models mentioned above. On the other hand, there are the first-principle RPIMC and DFT-MD simulations, and also the SESAME tables,⁹³ that predicted significantly lower compressibilities. It is one aim of the present article to resolve these discrepancies. Testing the accuracy of the different simulations for dense partially ionized hydrogen allows us to confirm and quantify the reliability of RPIMC and DFT-MD.

Still, more work will be needed to resolve the remaining disagreements between various experimental and theoretical predictions for the Hugoniot state in Fig. 2. For example, measurements at the Z-machine indicate a lower compression ratio at 50 GPa than is predicted by coupled electron–ion Monte Carlo (CEIMC) simulations which, on the other hand, are in good agreement with the laser measurements by Sano et al.¹⁰⁰ Compared to the Nova results, the more recent laser experiments at the $Ω$ facility have yielded much lower compression ratios, which brought them into better agreement with predictions from the Z-facility and from first principles simulations.

Finally, we note that a comparatively small discrepancy has remained between first principles simulations and the recent laser shock experiments by Fernandez-Panella et al.¹⁰¹ on cryogenic liquid deuterium at very high pressures up to $550 GPa$ [see Fig. 2(c)]. Rygg et al.¹¹³ suggested that collective modes (plasmons) would be missing from existing first principles calculations and might be the reason for the discrepancy. They performed a Debye-type calculation of the specific heat of electron plasma waves¹¹³ and simply added their contributions to the RPIMC results of Hu et al.¹¹⁴ This brought the theoretical results within the one-sigma error bars of the measurements above 300 GPa but could not explain the deviations at lower pressures. The central question is whether such plasmon waves are excited under such conditions and whether they would introduce a sufficiently large correction to the computed PIMC energies. Note that PIMC simulations take into account the entire Hamiltonian, so that plasmons are automatically accounted for, as long as the simulation cell is sufficiently large. By comparing results featuring 32 and 64 atoms, Militzer et al.⁹⁴ determined that the compression ratio shifts by less than 0.01 at $∼580 GPa$ and renders it unlikely that finite size effects or plasmons can explain the deviation of $∼$0.1 between the measurements by Fernandez-Panella et al. and the RPIMC predictions. In addition, in this paper, we present novel independent tests of the RPIMC data of Ref. 94. A comparison with the first-principle fermionic PIMC (FPIMC) results of Filinov and Bonitz¹¹⁵ shows that the deviations in energy and pressure between RPIMC and FPIMC, for the parameters of interest, are less than $1%$ (cf. Sec. IV A 3), which rules out the plasmon hypothesis. Note that the FPIMC data have undergone a finite size extrapolation, an issue that will be discussed in Sec. III A.

As an alternative to laboratory measurements, new insights into the quality of theoretical models might be obtained with novel simulation methods that involve fewer or no uncontrolled approximations, such as exact diagonalization for model systems. Another tool for benchmarks that was already mentioned is fermionic PIMC simulations (featuring the FSP) under conditions for which they are feasible. Indeed, progress has been made recently with first-principles PIMC simulations for jellium (the uniform electron gas model,^116,117 UEG) by Bonitz and coworkers who were able to avoid the FSP by a combination of two complementary approaches—configuration PIMC (PIMC simulations in occupation number representation, CPIMC)^118–121 and permutation blocking PIMC (an advanced coordinate space approach, PB-PIMC).^122–124 Simulations turned out to be feasible for all densities and temperatures exceeding half of the Fermi energy.^125–128 Both were employed to benchmark RPIMC simulations^121,129 and popular chemical models,¹³⁰ including the Padé formulas of Ebeling,⁵³ the STLS parameterization of Ichimaru et al.,^131–133 and the models of Vashishta and Singwi¹³⁴ and Perrot and Dharma-wardana.¹³⁵ Another result of these simulations is the highly accurate parameterization of the exchange–correlation free energy—the functionals of Groth et al. (GDSMFB)^117,128 and of Karasiev et al. (KSDT)¹³⁶—which constitute explicitly thermal exchange–correlation functionals for DFT simulations on the level of the local density approximation.^137,138 This demonstrates the high potential of parameter-free computer experiments for theory and model development in the field of warm dense matter.

Similar benefits can be expected from recent first-principles fermionic PIMC simulations—further improvements of the PB-PIMC approach—for partially ionized dense hydrogen by Filinov and Bonitz^115,139,140 that cover temperatures above 15 000 K and densities corresponding to values of the Wigner–Seitz (Brueckner) coupling parameter $rs≳3$ (see Sec. II A). Reference 115 presented tables of benchmark data and comparisons with several alternative simulations. Another step forward is fermionic PIMC simulations for the static density response of hydrogen for fixed ionic configurations by Böhme et al.,^141,142 which extend previous first principles results for jellium.^143–147 Finally, we mention very recent works by Dornheim et al.^148,149 who have presented the first dependable PIMC results for a number of structural, imaginary-time, and density-response properties for hydrogen at the electronic Fermi temperature considering the cases of $rs=1$⁠, $rs=2$⁠, and $rs=3.23$⁠.

When considering hydrogen at temperatures well below the electronic degeneracy temperature, $TF=EF/kB$⁠, fermionic PIMC (FPIMC) becomes highly inefficient because of the FSP. RPIMC simulations have been applied to temperatures as low as $0.1 TF$⁠, but nodal restriction makes sampling new paths and moving the nuclei inefficient below such temperatures (see below). On the other hand, electronic temperature effects become less relevant here, and one can resort to a finite-temperature treatment of the nuclei, combined with a ground state description of the electrons. This is the realm of coupled electron ion Monte Carlo (CEIMC). Usually, the electronic problem is solved by DFT methods (BOMD), but for hydrogen, even more controlled electronic ground state QMC methods have been applied: the CEIMC method uses electronic energies from ground state QMC to sample nuclear degrees of freedom using Monte Carlo (MC),^150,151 while a similar method using QMC-derived forces in a Langevin dynamics has also been proposed.¹⁵² Despite being heavy in terms of computer resources, those methods are more controlled than DFT since their accuracy can be judged based on the variational principle.

The first goal of the present paper is to summarize the current knowledge of the phase diagram of hydrogen, in particular at low temperatures. Here, many new results have been obtained recently but many details still remain open. We update the review of Ref. 156. The analysis is also extended to higher temperatures of partially ionized hydrogen over a wide range of densities.

The second goal of this paper is to give an overview of the diverse arsenal of simulation methods that can be applied to dense hydrogen and range from first-principles methods, such as QMC, over Green functions and DFT, to simpler models, such as chemical models and hydrodynamic equations. We discuss which quantities are accessible by the different methods and compare them with respect to their accuracy and applicability range. To this end, we perform comparisons for the equation of state, pair distributions, and degree of ionization between recent first-principles QMC results and DFT simulations, semiclassical MD with quantum potentials, average atom codes, and chemical models.

The third goal of this article is to present an overview of recent developments in simulations of dense hydrogen. These include the influence of finite-temperature XC-functionals in DFT simulations for hydrogen^{107,136,138,157} and the development of wave packet molecular dynamics simulations. Today, ab initio approaches can be used to analyze x-ray Thomson scattering (XRTS) signals from warm dense matter experiments.¹⁵⁸ It is, for instance, possible to infer the experimental plasma parameters by comparing DFT-MD data for the dynamic structure factor to the XRTS scattering signal.¹⁵⁹ A novel first principles approach to accurately analyze XRTS experiments was developed by Dornheim et al.^{17,148,160–164} Thus, it has become possible to compare not just integrated quantities like the EOS but also dynamic (frequency dependent) quantities from first principles simulations to experimental measurements. PIMC results exist for jellium for the dynamic structure factor,^165,166 the dynamic local field correction,¹⁶⁷ the plasmon spectrum,¹⁶⁸ the dielectric function, and the conductivity.¹⁶⁷ These can be extended in a semi-quantitative way to hydrogen and predicted a similar roton-type feature in hydrogen as in strongly correlated jellium that might be observable.^166,169

Given the high interest in ICF, we also reevaluate the open issues for the accurate theoretical modeling of the compression path of the deuterium–tritium fuel. The density temperature trajectory of the fuel and ablator materials traverses the warm dense matter into the high dense matter regimes as the ICF capsule implodes and then undergoes thermonuclear ignition, burns its fuel, and ultimately decompresses (cf. Fig. 3). Pressures range from around a Megabar, at a density of 10²⁴ cm⁻³ and temperatures of 1 eV, to tens of Gigabars of pressure, during ignition and burn, where densities are on the order of 10²⁶ cm⁻³ and temperatures are many kiloelectron volt. It is expected that the novel results obtained from first principles simulations presented in this article will also be valuable for ICF modeling.

The structure of this article is as follows: In Sec. II, we summarize the current knowledge as well as open questions on the phase diagram of hydrogen, including a discussion on metallization, the plasma phase transition, and the liquid–liquid phase transition(s). In Sec. III, we present an overview of important simulation approaches for dense partially ionized hydrogen including the regimes of low and high temperatures. In Sec. IV, we present our numerical results and comparisons of methods, including thermodynamic properties of hydrogen, density response functions, and transport properties. We also make suggestions for the future developments of the different methods and conclude with a summary and outlook in Sec. V where we also outline strategies for accurately combining different simulation approaches.

In this section, we present an overview on the most important properties and phases of hydrogen at high pressures. This includes the hypothetical plasma phase transition (PPT) (Sec. II B 1), the insulator–metal (liquid–liquid) transition (Sec. II B 2), and the proton crystal phase (Sec. II B 3). First, we introduce the basic parameters that are being used to describe highly compressed hydrogen. Since we discuss a large variety of properties and quantities as well as simulation methods and approximations, we provide a list of common acronyms in  Appendix in Table IX.

Despite its chemical simplicity, warm dense hydrogen exists in a broad variety of phases, ranging from solid to liquid, gas, and plasma (for an overview, see Figs. 3 and 4). For systematics, it is useful to introduce dimensionless parameters. In the following, we will use CGS units, i.e., set $4πε0=1$⁠.

The properties of the electrons and protons in thermal equilibrium are characterized by a combination of standard jellium parameters and Coulomb one-component plasma (OCP) parameters:¹⁷⁰ 

• The Wigner–Seitz (Brueckner) coupling parameter of electrons

  $rs=daB,$

  (1)
  is the ratio of the mean interparticle distance (or Wigner–Seitz radius), d, given by

  $4π3d3=1n,$

  (2)

  to the Bohr radius [Eq. (18)], where n denotes the number density.

• The degeneracy parameter

  $Θ=kBTEF=TTF,$

  (3)
  is the ratio of thermal energy to the Fermi energy, $EF=(ℏqF)2/2m$⁠, where the Fermi wave number is $qF=(3π2n)1/3$⁠. The parameter $Θ$ is closely related to the degeneracy parameter $χa$⁠, given in terms of the thermal de Broglie wavelength $Λ$ of electrons and protons as

  $χa=naΛa3,$

  (4)

  $Λa =h(2πmakBT)1/2 , a=e,p.$

  (5)

  This expression for $Λ$ holds in thermal equilibrium for ideal particles, its modification due to interaction effects is discussed in Sec. III G 2.

• The generalized coupling parameter

  $Γ=e2dkBT2σ+1nΛ3I3/2(βμ),$

  (6)

  which takes into account the quantum kinetic energy via the Fermi-integral $I3/2$ and the spin degeneracy factor $2σ+1$⁠.¹⁷¹ It is shown as the red line in Fig. 3. For non-degenerate particles in the gas (plasma) phase, i.e., for $χa≪1$⁠, this simplifies to the well-known classical Coulomb coupling parameter, $Γ=e2/(dkBT)$⁠. For high degeneracy, $Γ$ becomes constant when the zero-temperature quantum coupling parameter, $rs$⁠, is constant.

• Screening and dynamic properties of the free electrons are characterized by the Debye length, $rD$ and the plasma frequency, given as

  $rD=(kBT4πne2)1/2,$

  (7)

  $ωp =(4πne2m)1/2,$

  (8)

  $rTF =ℏ2m1/2(π3n)1/6.$

  (9)

At low temperatures, in the condensed phase, the atoms or molecules of hydrogen form an ordered lattice. Then pressure ionization and dissociation is a many-particle effect. Compression of the crystal leads to a change of the band structure and, eventually, to a closure of the energy gap between the highest occupied and the lowest unoccupied band: an insulator to metal transition (IMT). This effect was first studied by Wigner and Huntington,²⁰ predicting that molecular hydrogen (an insulator) would transition to metallic atomic hydrogen because of band filling. Mott^174,175 further considered the effects of electron correlation and computed critical values for the electron density beyond which metallization occurs. These values strongly depend on the material and the crystal symmetry, and the corresponding values for the coupling parameter are significantly higher than the hydrogen plasma result (27), for a recent overview (see Ref. 176). The IMT in hydrogen will be discussed in detail in Sec. II B 2.

Hydrogen has a very rich phase diagram, many details of which are under active investigation. Figure 4 shows a simplified overview. At high temperatures and low pressures (top left corner) hydrogen behaves as a classical gas of charged particles—a two-component electron–ion plasma (TCP). Upon cooling below the binding energy, $T≤Tion$⁠, hydrogen atoms and molecules form, the fraction of which increases as the temperature decreases.

At higher pressures, cooling potentially gives rise to a second fluid–fluid phase transition in partially ionized hydrogen. There have been extensive discussions in the literature about whether there are two independent phase transitions. Landau and Zeldovich,¹⁷⁷ Schramm,⁵⁰ and Norman and Starostin⁵¹ predicted that ionization and recombination of the plasma proceeds via a “plasma phase transition” (PPT) (see Fig. 4). We discuss theoretical and experimental predictions in Sec. II B 1 and conclude that there is no evidence for the PPT. A second liquid–liquid phase transition (LLPT) at megabar pressures and temperatures around $1000 K$ likely occurs between molecular hydrogen and metallic fluid hydrogen.⁶ 

Lowering the temperature further, hydrogen freezes. At these lower temperatures, nuclear quantum and isotopic effects become important. We discuss the location of the melting line in Sec. II B 2. There, we also consider phenomena predicted to occur upon compression, in particular the insulator-to-metal transition in solid hydrogen.²⁰ Whether this transition coincides with the transition from molecular hydrogen to atomic hydrogen is still an open question. Ashcroft et al.^21,178 have suggested that the solid atomic hydrogen phase will be a room-temperature superconductor because of the large phonon energies and electron–proton coupling. Finally, in Sec. II B 3, we concentrate on the compression of hydrogen gas up to TPa and PPa at relatively high temperatures (see Fig. 4). These pressures give rise to the formation of an atomic (proton) crystal that is bounded by a classical proton fluid at lower pressures, and a quantum proton liquid at higher pressures.

We start by considering partially ionized hydrogen at temperatures $10 000 K≲T≲Tion$ [Eq. (21)] and densities below the Mott density, $rs≳rsMott$ [Eq. (27)]. The degree of ionization $α$ increases monotonically with temperature (thermal ionization) and density (pressure ionization) which is accompanied by a strong increase in electrical conductivity. This parameter range attracted particular attention in the 1970s to 1990s, due to the prediction of the so-called plasma phase transition (PPT)—a hypothetical first-order phase transition that was thought to be analogous to the van der Waals gas–liquid transition. The idea was^50,51 that the combination of overall attractive Coulomb forces in the plasma and short-range electron repulsion (due to quantum effects), below a critical temperature may lead to a fluid with two phases of different degree of ionization. Additional arguments for the existence of the PPT were similar effects in electrolytes, alkali metals, and electron–hole plasmas in semiconductors (e.g., Ref. 68). The critical point of the PPT was closely linked to the ionization potential of the hydrogen atom, and most predictions settled around $Tcr≈0.1 Tion∼$ 15 800 K; for a representative collection of results (see Table I). The majority of predictions of the PPT originated from chemical models (CM, the top part of Table I) (for details, cf. Sec. III B).

The PPT was also investigated with first principles fermionic PIMC simulations.^77,185,186 However, these simulations are severely hampered by the sign problem (cf. Sec. III E 2) and could not achieve converged results.¹⁸⁷ On the other hand, restricted PIMC simulations with free-particle nodes reported a PPT⁸⁵ but were not confirmed by more accurate later simulations with variational nodes⁸⁷ (see also Refs. 188 and 189, cf. bottom part of Table I).

Finally, on the experimental side (bottom of Table I) hydrogen metallization and a density jump at high pressure were reported in shock compression experiments using explosives.⁷ On the other hand, Weir et al.⁶ observed metallization of hydrogen and deuterium in shock experiments, around $p=140 GPa$ at temperatures where hydrogen is a fluid, but did not observe a phase transition. Fortov et al. also used explosives-driven shock compression of deuterium⁸ and claimed proof of the PPT, based on observation of a “density jump”' in the range of p = 127–150 GPa. However, the evidence is not convincing,¹⁹⁰ and the studied temperature range of 2550–4000 K (as deduced from models) is confined to the fluid phase. Thus, the observations could rather be an indication of the liquid–liquid transition (see Sec. II B 2).

The topic of the PPT has been frequently discussed in reviews and textbooks (see, e.g., Refs. 156, 191, and 192). Also, Refs. 193 and 194 contain a recent overview of the history of the PPT, from the point of view of their Russian authors. Table I summarizes the extensive literature on the PPT. Note that chemical models become questionable if the ionization fraction changes significantly with density. This is the case, in particular, if the binding energies vanish near the Mott density, where atoms (or molecules) transform to short-lived transient states, as typical for WDM. The PPT may be an artifact that is “built into” these models.^195,196 For more details, see Sec. III B. State-of-the-art first-principles simulations in the physical picture show no indications of a PPT. DFT-MD simulations of Lorenzen et al.³³ and coupled electron–ion (CEIMC) simulations of Morales et al.^32,197 have produced evidence of a phase transition at much lower temperatures in the fluid phase (LLPT), which differs qualitatively from the PPT (for more details, see Table II and Sec. II B 2).

The phase diagram of hydrogen under extreme conditions is still uncertain, mainly because of experimental difficulties in producing and controlling the extreme physical conditions required¹⁹⁹ and the limited information experiments obtained about the system. At low temperatures in the solid phase, the most important missing information is the crystal structure as a function of temperature and pressure. In Fig. 5, we report the present low-temperature phase diagram emerging from experimental information. Up to six different crystal structures have been detected but the transition lines are traced based on discontinuities and changes, respectively, of the vibrational frequencies of the infrared (IR) and Raman spectra, while measurements of Bragg peaks are mostly missing. A notable exception is along the 300 K isotherm where special x-ray spectroscopy has been developed and applied to confirm that the structure of phase I is m-HCP up to 250 GPa.^200,201 Candidate structures have been inferred by the ab initio Random Structure Search (AIRSS) method²⁰² based on DFT calculations with phonon corrections.

A relevant and still partially unanswered question concerns the mechanism by which solid hydrogen metallizes upon increasing pressure.¹⁹⁹ 

Important advances toward the metallization of solid hydrogen at low temperature were announced by different groups in recent years. In 2017, the group at Harvard reported the observation of metallic hydrogen in diamond anvil cell (DAC) experiments at 495 GPa below 80 K.²⁰³ This conclusion was based on the sudden appearance of a reflective sample that was interpreted as evidence for the metallization of hydrogen. This interpretation has been criticized by others,¹⁹⁹ and so far the results have not been reproduced. Almost simultaneously the group in Mainz reported evidence of the formation of a semi-metallic, still molecular, phase at around 350 GPa and below 100 K. They employed both optical probes and direct electric measurements in a DAC.²⁰⁴ In 2020, Loubeyre et al.²⁴ reported results using a toroidal DAC with synchrotron radiation. They measured the IR absorption profile over a wide range of pressure and detected complete absorption at 427 GPa and 80 K, which was interpreted as a sudden closure of the direct energy gap, a strong indicator for the state of a “good metal.” Although sample visual inspection and reversibility of the transition upon pressure release suggest that this metal is still molecular, strong experimental evidence, such as the persistence of the vibron signal, is missing. Thus, it is conceivable that the observed collapse of the infrared gap signals a metallization through molecular dissociation.

In Fig. 5, we sketch the current low-temperature phase diagram up to the IMT region. Below 200 K in phase III, we report the transition to the semimetallic state at 350 GPa²⁰⁴ as a blue shaded area and the later transition to the “metallic” state at 425 GPa²⁴ as the red shaded area. At 495 GPa, we also report the results of the Harvard group²⁰³ as the dark red-filled area. No experimental information is available for metallization at temperatures higher than ∼100 K. Candidate structures for the low temperature phase III, emerging from AIRSS^202,216 and a subsequent QMC analysis,²¹⁷ are C2/c-24 and Cmca-12. Gap closure of those structures has been detected with QMC methods²¹⁸ finding that the C2/c-24 structure undergoes an indirect gap closure around 360 GPa, in quantitative agreement with the experimental transition to the semimetallic phase. Moreover, it was shown that the direct gap of the C2/c-24 phase follows the experimental absorption edge behavior with pressures from Ref. 24. We emphasize that this agreement can only be obtained by taking into account the zero point motion of the protons.²¹⁸ Moreover, modeling of the correlated electrons requires going beyond common DFT approximations and using QMC methods.^218–220 A recent study,^221,222 based on the Stochastic Self-Consistent Harmonic Approximation (SSCHA) and QMC-corrected energies suggests that the observed collapse of the direct gap at 425 GPa²⁴ is related to a structural transition from the C2/c-24 to the Cmca-12 metallic structure, in both hydrogen and deuterium, as verified by experiments.²²³ 

Above 200 K, phase III transforms into a new phase IV when the molecular vibron line splits into two characteristic frequencies suggesting that different molecules may experience two different environments.²²⁴ For this phase, theory predicted a layered structure in which regular molecular layers are intercalated with planar layers where molecules form almost regular hexagons, precursors of a transition to a hexagonal monatomic layer at higher pressure.^216,225 Not much is known about metallization in this higher temperature phase (see, however, Refs. 225 and 226).

The black lines in Fig. 5 are experimental phase boundaries. However, experimental information about the crystalline structures of those phases are missing, and the experimental determination of the melting line is not based on the evidence of the vanishing of long-range order (Bragg peaks) above the melting line. We also report the melting line from two different theories. The orange line is from BOMD simulations with the DFT-PBE approximation.^32,90 The green line is a recent prediction from a Machine Learning force field trained on QMC energies and forces (ML-QMC).²⁰⁶ Within DFT, the location of the melting line and LLPT line strongly depends on the approximation for the exchange–correlation functional. The DFT-PBE melting line seems to agree better with the experiments, although it was observed that, considering quantum nuclei within this theory, leads to unphysical results.^197,227 At pressures between 50 and 200 GPa, ML-QMC predicts a different solid structure that is much more stable than in PBE-PIMD calculations. This prediction has yet to be confirmed by experiments.

Above the melting line, fluid hydrogen can be either molecular or atomic. First-principle simulations, both by Born–Oppenheimer molecular dynamics (BOMD) and by coupled electron–ion Monte Carlo (CEIMC),^{31–33,197,214,228,229} predict the existence of a weakly first order LLPT between the molecular insulating fluid and the mostly monatomic metallic fluid, below a critical temperature estimated to be $TLLPTcr∼$ 2000 K. This is signaled by a small discontinuity in the specific volume at a given pressure and temperature. The precise value of $TLLPTcr$ is not known but, above $TLLPTcr$⁠, molecular dissociation and metallization increase progressively with pressure. In Fig. 6, we report CEIMC predictions (cf. blue and red shaded lines) for the LLPT line of hydrogen and deuterium between 600 and 1500 K.¹⁹⁷ Other lines from BOMD are in qualitative agreement, but their precise location depends on the specific XC approximation adopted.^197,230 Unequivocal experimental evidence of the occurrence of a first-order transition is presently missing.

Experiments on fluid hydrogen are performed either by static compression, with DAC,²⁰⁷ or by dynamic compression, with shock wave techniques^6,8,212,213 (cf. Fig. 1). During the pressure increase, the sample first becomes opaque when the electronic gap (cf. Sec. IV C 3) becomes comparable to the energy of the probe laser (absorption coefficient $μ≥1μm−1$⁠), and later, at higher pressure, it turns reflective, signaling the occurrence of the metallic state (reflectivity $R=0.3$⁠).

In Fig. 6, we show static compression data at $R=0.3$ (DAC-r) from Refs. 207 and 231 (short-pulse DAC, sp-DAC) and Ref. 211 (long pulse DAC, lp-DAC) on hydrogen and deuterium, together with shock wave data on deuterium from the Sandia group²¹³ and from the Livermore group.²¹² The inception of absorption (⁠ $μ≥1μm−1$⁠) in those experiments is represented by open symbols while closed symbols represent the occurrence of reflective samples (R = 0.3). Also, data from sp-DAC show their temperature plateau (DAC-p). Together with experimental data, Fig. 6 also reports theoretical predictions (obtained within the Kubo–Greenwood approach, cf. Sec. IV E 2) for the inception of an opaque sample and inception of a metallic sample, from Ref. 214. While we see an essential agreement between different experiments and with the theory, for the absorption threshold, reflectivity measurements are not in agreement between the two different shock wave experiments. Note that in those experiments, temperature is not directly measured but it is inferred from a model EOS, and predictions are sensitive to the adopted model. Theoretical predictions for the reflectivity threshold are closer to the NIF-r data and are also in agreement with sp-DAC-r measurements.

In Ref. 214, it was shown that $R=0.3$ corresponds roughly to a conductivity value of $σ=$ 2000 S/cm, although the correspondence depends on the refractive index of the considered interface (see Ref. 214 for details). This value of the conductivity of metallic hydrogen was observed as a saturation value in early shock wave experiments.⁶ In Fig. 6, we also report theoretical predictions of the occurrence of metallic hydrogen (⁠ $σ=2000 S/cm$⁠) from PIMD with a modern XC approximation (SCAN+rVV10).²¹⁵ While below the critical point, those predictions are in agreement with other estimates, at higher temperatures this behavior deviates toward the absorption threshold.

Before closing this section, we should note that the new ML-QMC melting line has a potentially large impact on interpreting and understanding this part of the phase diagram. For instance, the thermodynamic path followed during the compression in the shock wave experiments might be influenced by the occurrence of crystallization and subsequent re-melting of the sample. We observe that, except for the $μ=1 μm−1$ data, all other data show a “kink” as the temperature is lowered including that of the sp-DAC-p data. The relative location of the new melting line and the LLPT line from CEIMC suggests that this kink might be the signature of a reentrant melting transition—a speculation that needs confirmation. Also, we note that a higher melting line restricts the LLPT region to a much smaller portion of the phase diagram.

As noted above, the crystallization of heavy charged particles is also relevant for other systems. Prominent examples are white dwarf stars where crystallization of carbon ions (⁠ $C6+$⁠) and oxygen ions (⁠ $O8+$⁠) is expected to occur in the core. Due to the increased ion mass, the inequality (28) leads to a much broader density range and higher maximum temperature. Equations (30)–(32) yield, for carbon (oxygen), $n(1)=2×1026cm−3(6.6×1026cm−3),n(2)=3.7×1033cm−3(2.7×1034cm−3)$⁠, and $T*=109K(4.2×109K)$⁠. As a final note, liquid or crystal formation is also relevant for certain semiconductors. In fact, hole crystallization was predicted, among others, by Abrikosov^239,240 who predicted a critical hole to electron mass ratio of $Mcr≈100$ for CuCl. Fermionic PIMC simulations for electron–hole plasmas yielded^172,241,242 $Mcr≈83$⁠, in three dimensions (3D), and $Mcr≈60$⁠, in 2D [cf. Eq. (33)]. Such mass ratios are feasible in intermediate valence semiconductors, such as Tm[Se, Te].²⁴³ 

We conclude the analysis of the phase diagram of dense hydrogen by outlining open problems and indicating possible solutions. As we have seen before, the phase diagram of hydrogen is surprisingly complex and identification of phases and phase transitions or crossovers very sensitively depends on the analysis tools. As we have seen in the example of the shock Hugoniot (cf. Sec. I and Fig. 2), various experiments and simulation methods may lead to strongly differing predictions which requires a careful analysis.

• The high-pressure phases of condensed hydrogen and the LLPT are particularly complex, with many competing possible orders. The simulation results presented in Sec. II B 2 and summarized in the phase diagrams of Figs. 4–6 include extensive predictions that still await experimental tests. For completeness, we mention that a different picture has been put forward by Fedorov et al. (see Ref. 244 and references therein). Based on restricted open shell Kohn–Sham simulations they predict the existence of excitons in the low-temperature fluid. The dissociation of excitons would represent an efficient energy transfer from nuclear to electronic degrees of freedom beyond the BO approximation which, however, should be clearly visible in the optical data. They hypothesize the existence of a plasma phase and PPT, at low temperatures.

• Accurate experiments and first principle simulations are needed for definitive findings. Reliable simulation results are of great value for benchmarking other less fundamental methods. Examples are fermionic PIMC simulations for the jellium model,^117,121,127 as well as for hydrogen¹¹⁵ as discussed in Sec. IV A.

• The DFT and QMC simulations are currently limited to treating a few hundred atoms. While pressure and internal energy seem to converge fairly quickly with increasing system size, this limitation may lead to biases in the vicinity of phase transitions. Extrapolation to the thermodynamic limit (cf. Sec. III A 4) is necessary to establish an accurate long-range order and critical exponents. Methods that can treat larger system sizes are helpful. A promising technique is to combine DFT or QMC results with machine learning techniques (e.g., Refs. 206, 245, and 246) but the accuracy and reliability to reproduce in general the underlying DFT or QMC simulations still needs to be established.

• The high-temperature and high-pressure range is of particular interest for ICF (cf. Fig. 3). This wide range of material conditions places a severe set of constraints on computational techniques used to compute equations of state and transport properties. The radiation-hydrodynamic codes (cf. Sec. III H 5) use physics models such as the equation of state, electrical and thermal conductivity, and plasma viscosity as input to the Navier–Stokes equations. It is important that the phenomena are well understood theoretically and modeled accurately so that the ICF design simulations are as accurate as possible.

A major challenge for experiments with matter at extreme densities, temperatures, and pressures, in general, and dense hydrogen, in particular, is given by accurate diagnostics. Indeed, very often, basic parameters such as density and temperature cannot be directly measured and have to be inferred from other observations and theories. The velocity interferometer for any reflector (VISAR) diagnostics is a standard tool to obtain pressure and density in a shocked state of the sample.^247,248 If the pressure standard is known to have good accuracy and the shock is planar enough, it will deliver reliable results. Streaked optical pyrometry (SOP) is a diagnostic used to extract the temperature of the sample.^248,249 However, it usually measures (near) surface temperatures and the results depend on the (model-dependent and often unknown) reflectivities of the warm dense matter state in a gray body model.

In practice, however, XRTS experiments with hydrogen are notoriously difficult. Hydrogen has a small scattering cross section, leading to a low photon count and, consequently, a noisy intensity signal (e.g., Ref. 255). This is particularly problematic for XRTS experiments with fusion plasmas, which have a very low repetition rate (e.g., one to three shots per day at the NIF²⁵⁶). XRTS experiments with hydrogen jets^257,258 at modern x-ray free electron laser (XFEL) facilities, such as LCLS in the USA or the European XFEL in Germany, on the other hand, allow one to average over thousands of shots. Here, a crucial problem is given by the uncertain degree of inhomogeneity and nonequilibrium of the sample.²⁵⁹ 

Let us conclude by outlining two additional, general challenges of XRTS experiments with warm dense matter. On the experimental side, the accurate interpretation of the measured intensity requires detailed knowledge of the source-and-instrument function $R(ωs)$⁠, which is a nontrivial task for both backlighter sources²⁶⁰ and XFEL facilities.²⁵³ From the theoretical perspective, understanding the measured signal usually requires model or simulation results for $See(q,ω)$⁠, which is a difficult task. In our opinion, the most flexible tool for this purpose is given by time-dependent DFT (TDDFT), which is discussed in Secs. III H 4 and IV F 1. An alternative strategy has recently been suggested by Dornheim et al.,^17,160,164 who have proposed to evaluate Eq. (34) in the imaginary-time domain (see Sec. IV D 2). Remarkably, this allows for model-free diagnostics of parameters such as the temperature^160,161 and the static structure factor $See(q)$⁠,¹⁶² and allows for direct comparisons with quasi-exact PIMC simulations^148,163 (see also Sec. IV D 3).

Given these recent developments, we are confident that future XRTS experiments with hydrogen will give valuable new insights into the equation of state components and other important properties.

There exist a few methods that allow us to solve the many-body problem with the Hamiltonian (35) exactly (i.e., without approximations or statistical noise from sampling). These include exact diagonalization (configuration interaction, CI) and density matrix renormalization group (DMRG) approaches which are, however, limited to small particle numbers N and/or small number of basis functions and have not been used for dense hydrogen. They are valuable for benchmarking of other methods in model cases that might be relevant to dense hydrogen.

Before discussing simulation approaches that are being applied to dense hydrogen, we briefly summarize the main many-particle properties that are of interest, both for comparison with experiments and for basic physical understanding.

The main quantities of interest can be classified into the following groups for which we also indicate the sections where they are discussed in this paper.

• The first relevant set of observables is given by the thermodynamic properties of dense hydrogen. This includes the equation of state (EOS, pressure) (Sec. IV A 3), total energy, free energy etc. More detailed information on the structural properties and interaction effects is obtained from the species-resolved static structure factors $Sab(q)$ and pair distribution functions $gab(r)$ which will be studied in Sec. IV A 5. Complementary microscopic information is contained in both the electronic and ionic momentum distributions $ne(q)$ and $ni(q)$ (Sec. IV C), which are important, e.g., for the characterization of the effective charge state or, in the case of the ions, for the estimation of ion impact reaction or nuclear fusion rates (e.g., Refs. 261 and 262).

• The second group of important observables describe electronic dynamic properties. These include the dynamic structure factor, $See(q,ω)$⁠, which is the key property in XRTS experiments [cf. Eq. (34)], and the imaginary-time correlation function (ITCF), $Fee(q,τ)$⁠, which is connected to the former via a two-sided Laplace transform [Eq. (121)]. Interestingly, the ITCF has emerged as an important observable for XRTS diagnostics in its own right,^{160–163,263} as it gives one straightforward and model-free access, e.g., to the temperature.

• The third group of observables are electronic transport and optical properties, such as the conductivity $σ(ω)$⁠, the heat conductivity $λ(ω)$ (cf. Sec. IV E 3), and the opacity $κ(ω)$ (Sec. IV E 4). The latter is a key property in optical transmission experiments, e.g., at the NIF.²⁶⁴ 

• A related fourth group of properties includes spectral information, such as the density of states (DOS), as well as the single-particle spectral function, $A(q,ω)$⁠, and the Matsubara or nonequilibrium Green functions (Sec. III H 3). Note that the spectral function can be directly measured in photoemission experiments. Closely related properties are the energy gap that is important to characterize the low-temperature phases of hydrogen (cf. Sec. IV C 3) and the ionization potential depression (Sec. IV B 1), which characterizes the ionization of high density hydrogen in the gas phase.

• The fifth set of relevant observables is given by the ion dynamic properties, that have traditionally been estimated on the basis of molecular dynamics simulation. In particular, we consider the ionic structure factor, $Sii(q,ω)$⁠, which is important for an analysis of the collective modes (ion acoustic modes) and sound speed of the heavy particles (Sec. IV D 6). Other important observables are the diffusivity, $Dion$⁠, viscosity, $ηion$⁠, and the ionic thermal conductivity $λion$⁠.

• The sixth group of observables are related to important non-equilibrium properties which characterize the equilibration of hydrogen after an external excitation. We list three examples: the relaxation time, $trel$⁠, during which the electron momentum distribution, $ne(q,t)$⁠, thermalizes to a Fermi function (in general, to a correlated equilibrium distribution), the equilibration time between electrons and protons, $teq$⁠, and the stopping power, $Sx$ (Sec. IV F 3) (and related energy loss properties).

The above list does not pretend to be complete but concentrates on quantities that have been in the focus recently. Other observables that are potentially interesting for future investigations, but which are not considered here, include various nonlinear response functions,^265–269 effective forces,²⁷⁰ and potentials.^271,272

Now, obvious questions arise: which simulation approaches are available to compute the above quantities? And, in case a computation is possible, how accurate and reliable are the results? How can the accuracy be estimated or verified? Finally, if several methods are available—which one is preferable? In the following, we attempt answers to these questions, with a focus on applications to dense hydrogen. A comparison of accuracy and limitations of various methods is presented in Tables III–V. A summary of the relation between methods and observables is presented in Table VI.

For dense hydrogen in the ground state, the most accurate results, so far, have been obtained using diffusion or variational Monte Carlo (DMC, VMC). Ground state properties are not the focus of the present paper, for an overview (see Ref. 156). For finite temperature properties, one has to switch to a mixed ensemble description based on the N-particle density operator $ρ̂$⁠. Here, the most accurate approach is provided by fermionic path integral Monte Carlo (FPIMC) which is free of systematic errors, and convergence to the exact thermodynamic limit can be achieved by proper extrapolation with respect to the particle number N and the number of high-temperature factors P (see Sec. III E 1). FPIMC is the only method which is capable of producing unbiased accurate predictions without requiring benchmarks by other simulations. The limiting factor of FPIMC is the fermion sign problem (FSP, cf. Sec. III E 2) manifesting itself by an exponentially small signal to noise ratio, which currently limits simulations to temperatures $T≳0.5TF$⁠. There are various concepts to avoid this limitation by formulating PIMC in second quantization (CPIMC) (see, e.g., Sec. III E), but the accessible parameter range remains restricted by the FSP.

To extend the simulations of hydrogen to the parameters of interest (or to more complex systems), as well as to increase the computational efficiency, a large arsenal of approximate methods has been developed, both for equilibrium and nonequilibrium situations. An overview is presented in Tables IV–VI, and a more detailed discussion is given in Secs. III A 4 and III A 5. The purpose of these tables is to provide a compact comparison of important methods which includes their applicability range and the quantities that they are capable of providing. The methods are listed on the order of their accuracy that is generally expected, based on their construction and the involved basic approximation schemes. Note that the actual accuracy of different methods, when applied to dense hydrogen, may differ from these expectations. While the ultimate tests are, of course, experiments, they are extremely challenging and complex in the case of dense hydrogen and currently afflicted by large uncertainties or disagreements between different experimental groups and techniques. Therefore, we base our tests on computer experiments for which the accuracy has been convincingly established and which serve as benchmarks. These are fermionic PIMC (FPIMC) results for the uniform electron gas (e.g., Refs. 117 and 129) and hydrogen.¹¹⁵ Presently, only a few such benchmarks are available; therefore, in the present paper, we produce additional comparisons which are all listed in Tables IV and V.