## New Research In

### Physical Sciences

### Social Sciences

#### Featured Portals

#### Articles by Topic

### Biological Sciences

#### Featured Portals

#### Articles by Topic

- Agricultural Sciences
- Anthropology
- Applied Biological Sciences
- Biochemistry
- Biophysics and Computational Biology
- Cell Biology
- Developmental Biology
- Ecology
- Environmental Sciences
- Evolution
- Genetics
- Immunology and Inflammation
- Medical Sciences
- Microbiology
- Neuroscience
- Pharmacology
- Physiology
- Plant Biology
- Population Biology
- Psychological and Cognitive Sciences
- Sustainability Science
- Systems Biology

# Evidence for a first-order liquid-liquid transition in high-pressure hydrogen from ab initio simulations

Contributed by David Ceperley, May 26, 2010 (sent for review February 11, 2010)

### Related Articles

- The insulator-metal transition in hydrogen- Jul 20, 2010

## Abstract

Using quantum simulation techniques based on either density functional theory or quantum Monte Carlo, we find clear evidence of a first-order transition in liquid hydrogen, between a low conductivity molecular state and a high conductivity atomic state. Using the temperature dependence of the discontinuity in the electronic conductivity, we estimate the critical point of the transition at temperatures near 2,000 K and pressures near 120 GPa. Furthermore, we have determined the melting curve of molecular hydrogen up to pressures of 200 GPa, finding a reentrant melting line. The melting line crosses the metalization line at 700 K and 220 GPa using density functional energetics and at 550 K and 290 GPa using quantum Monte Carlo energetics.

Since the pioneering work of Wigner and Huntington (1), on the metallization of solid molecular hydrogen by pressure, there has been a great effort to understand the molecular dissociation process in high-pressure hydrogen from both experiment and theory. In the solid, at low temperatures, metallization has been expected to occur in conjunction with a transition to a solid atomic state, although a transition to exotic phases such as quantum fluids (2) or metallic molecular phases may also be possible (3–5).

For dense hydrogen in the liquid phase, metallization (probably accompanied by molecular dissociation) can occur either through a continuous process (a crossover) or through a sharp, first-order transition, often called the plasma phase transition. Numerous experiments have been performed in the liquid phase using dynamic compression techniques to measure both the principal Hugoniot in hydrogen and deuterium [using for example: gas guns (6), laser-driven compression (7–9), magnetically driven flyers (10, 11), and converging explosives (12)] and to measure off-Hugoniot properties at lower temperatures [electrical conductivity measurements using shock reverberation (13) and multiple shocks (14), compressibility measurements using explosive-driven generators (15), to mention a few]. The conductivity measurements by Nellis and coworkers (13) produced the first evidence of minimum metallic conductivity in fluid hydrogen at a pressure of 140 GPa and temperatures on the order of 3,000 K. Until recently, there were no experimental indications of a first-order liquid-liquid transition (LLT) in hydrogen. Even though results could not rule out the existence of the transition and most studies had been performed at fairly high temperatures, there was no sign of a sharp, first-order behavior. The only direct experimental evidence of a LLT is from the work of Fortov et al. (15) where reverberating shocks produced with high explosives were used to ramp compress hydrogen, presumably reaching temperatures in the range of 3–8 × 10^{3} K. Using highly resolved flash X-ray diagnostics, they were able to measure the compressibility of the liquid and found a 20% increase in density in the regime where the conductivity increases sharply by approximately 5 orders of magnitude. It is important to note that in these experiments the temperature is estimated using a “confined atom” model whose accuracy is difficult to assess. In general, accurate temperature measurements under these extreme conditions are notoriously difficult to achieve. Moreover, the existing experimental data is sparse and a rapid, yet continuous, change in conductivity with increasing pressure cannot be completely ruled out from the data.

The majority of the theoretical techniques that have been used to study dense liquid hydrogen fall into two groups: methods based on effective models (the chemical picture) and simulation methods starting from protons and electrons (first-principles approach). In general, the effective models exhibit a clear first-order LLT that persists for temperatures well above 10,000 K (16–19). However, the presence of a first-order phase transition in the effective models is now recognized to be spurious (20) and the most recent implementations now specifically include ad hoc terms to smooth out the transition region to ensure positive compressibility (21).

The vast majority of first-principles simulations of extended systems are performed using density functional theory (DFT). Although realizations of DFT contain many approximations, current implementations provide a reasonable compromise between accuracy and computational requirements. At high pressures, where experimental results are scarce, more accurate electronic structure methods are needed to assess the quality of DFT predictions. In this work, we use the accurate quantum Monte Carlo (QMC) method to test the predictions of DFT in the dissociation and metallization regime of the phase diagram of hydrogen at high pressures. QMC has been shown to give very accurate results in extended systems and currently has the capability to reach near chemical accuracy in many systems (22).

Although there are numerous publications of DFT simulation studies of hydrogen, there is currently no consensus about its predictions regarding the existence of a first-order transition in the liquid. Two of the earlier studies using Car–Parrinello molecular dynamics (MD), by Scandolo (23) and Bonev et al. (24), showed evidence of a first-order transition, for (*P* = 125 GPa, *T* = 1,500 K) in the former and (*P* = 52 GPa, *T* = 2,800–3,000 K) in the latter. The use of the Car–Parrinello method in these simulations put in question the validity of their results because hydrogen metallizes at these conditions, which can potentially led to problems with this approach. As a result, the majority of the subsequent first-principles simulations of hydrogen have relied on the use of Born–Oppenheimer molecular dynamics (BOMD) to describe the onset of metallization. For example, the BOMD simulations reported by Vorberger et al. (25) and Holst et al. (26) exhibited a gradual crossover with no evidence of a sharp transition, and a recent BOMD study by Tamblyn and Bonev (27) indicated that a first-order transition is likely but still an open question. In addition a previous exploration by coupled electron-ion Monte Carlo method (CEIMC) (35) also suggested the occurrence of a continuous molecular dissociation. These studies have led to the general belief that the results reported in the earlier simulations were influenced by the use of the Car–Parrinello method and that the transition region as described by DFT is indeed smooth (36).

In this work we consider again the onset of molecular dissociation in liquid hydrogen at high pressure. We employ first principle simulations based on DFT within BOMD and QMC within the CEIMC approach (37, 38) with improved trial wave function (39, 40). To pinpoint a possible transition we trace isotherms using a much finer grid in density than in previous studies. In addition, we calculate the melting line of the molecular solid using free energy calculations in both liquid and solid phases with DFT. Fig. 1 shows the resulting phase diagram. A sharp, first-order phase transition, is observed at both levels of theory. At high temperatures, the transition ends at a critical point around *T* ≈ 2,000 K above which metallization is a continuous process with pressure. The values of the transition pressure depend on the simulation method, with the more accurate predictions coming from QMC simulations. At low temperature the DFT transition line meets the melting curve in a liquid-liquid-solid multiphase coexistence point at approximately 700 K whereas the extrapolation of the QMC transition line and the melting line at higher pressures should meet at 290 GPa and 550 K. The results presented here resolve the discrepancies among previous DFT calculations of hydrogen, reconcile the theoretical understanding of the metallization in hydrogen and the results of the dynamical compression in the liquid phase (13, 15) and the static compression* in the solid phase.

## Results and Discussion

### Liquid-Liquid Phase Transition.

Fig. 2 shows a comparison of the 1,000 K isotherm in hydrogen as a function of density for the different simulation methods used here. A plateau in the pressure over a narrow range of densities (on the order of 0.02 g/cm^{3}, which corresponds to a volume change of approximately 2%) is clearly seen in both CEIMC and BOMD simulation results, indicating the existence of a first-order transition. In this region (*dP*/*dV*)_{T} = 0, which implies a diverging isothermal compressibility, *κ*_{T} = -*V*^{-1}(*dV*/*dP*)_{T}. We observe similar behavior at *T* = 1,500 K for CEIMC and at *T* = 1,500 K, 800 K, and 600 K for BOMD. The transition occurs systematically at higher pressures in CEIMC as compared to BOMD. This is reasonable given that most DFT exchange-correlation functionals have band closure at densities that are too low and tend to favor delocalized electronic states over localized ones, presumably due to self-interaction errors (42). The estimated latent heat (enthalpy change per atom) associated with the transition at *T* = 1,000 K is 440 (90) K and 570 (160) K for the DFT and QMC transition lines respectively. The corresponding values obtained from the Clausius–Clapeyron equation Δ*H* = *T*Δ*V*(*dP*_{LLT}/*dT*) are 380 (90) K and 410 (90) K. In this case, the derivatives of the transition pressure with temperature are estimated numerically. There is good agreement between the two independent estimations, which suggest that our results are thermodynamically consistent. The overall location of the transition region found here is in agreement with a recent BOMD study reported in ref. 27. However, in this case a (*dP*/*dV*)_{T} = 0 region was not found, most likely due to inadequate sampling across the transition region.

It is interesting to note that as the transition region is approached, the simulations exhibit quite strong size effects. To obtain converged results in the BOMD simulations (sampled at the Γ-point), we had to use a system with 432 atoms. Our liquid simulations with 432 atoms produced radial distribution functions that were in excellent agreement with larger system sizes, but showed significant differences with 128 and 256 atom simulations. A similar observation of strong size effects was made in ref. 27, which indicated that the orientational order of the liquid is particularly sensitive as well. For a much smaller system size of 64 atoms, we observed first-order behavior during the metallization transition only when a 3 × 3 × 3 grid of *k*-points was used; simulations at the Γ-point showed no evidence of significant discontinuities. In the case of the CEIMC simulations we use an integration over the boundary conditions of the simulation cell with 64 (4 × 4 × 4) *k*-points. Simulations with 54 atoms produced radial distribution functions that were indistinguishable from simulations with 128 atoms and 64 *k*-points, for both the metallic and insulating sides of the transition.

The simulation results presented so far have assumed that the protons are classical particles. To test the validity of this assumption and to study its influence on the nature of the transition region, we performed path integral molecular dynamics (PIMD) simulations with DFT at *T* = 1,000 K. The path integrals were discretized with a time step of (8,000 K)^{-1}, which corresponds to 8 beads at *T* = 1,000 K, enough to obtain converged results. A pressure isotherm from the PIMD simulations is shown in Fig. 2. The inclusion of zero point motion in hydrogen produces a decrease in the transition pressure by approximately 40 GPa, but does not change the nature of the transition; a clear plateau, indicating a first-order transition still exists in the pressure isotherm.

Fig. 3 shows the proton-proton distribution function and electronic conductivity at the two ends of the transition for the DFT isotherm at *T* = 1,000 K. The transition has strong effects on electronic properties, in particular, those related to electron localization. The molecular fraction changes slowly across the transition; in fact, partial dissociation is seen at densities well below the transition. The molecular state in these conditions is very different from that found at lower densities where there is no overlap between molecules. Close to the transition, the molecular state is transient and very weakly bound; molecules form and break on small time scales of the order of the collision frequency (25, 27). Whereas the description of molecular and atomic fractions based on somewhat arbitrary bonding state criteria has been used to describe the transition region, and has led to conclude in favor of a continuous dissociation transition (35), we believe that it is more indicative to use the discontinuity of the electrical conductivity: It coincides precisely with the onset of the (*dP*/*dV*)_{T} = 0 plateau.

Fig. 4 shows the (time averaged) DC conductivities as a function of pressure for CEIMC, BOMD, and PIMD simulations. The calculated conductivities are in rough agreement with experimental observations, as also shown by previous simulations (26, 43). At 2,000 K and above, the conductivity increases smoothly as a function of pressure, but for *T* ≤ 1,500 K there is a sharp and discontinuous increase in the conductivity at the transition. The size of the discontinuity increases with decreasing temperature in both BOMD and CEIMC simulations. By extrapolating to the pressure where the size of the discontinuity in conductivity goes to zero, we can estimate the critical point of the transition. Above this temperature, the conductivity and the dissociation as a function of density is continuous. At lower temperatures, a first-order transition is encountered as we go across the dissociation regime with a volume discontinuity and a sharp metallization of the liquid. The transition corresponds to a sharp change from semiconducting to metallic behavior caused by the sudden collapse of the molecular state. Even though strong short range correlations persist after the transition is crossed, the resulting atomic liquid is metallic. Our findings are in agreement with a recent qualitative analysis of the hydrogen phase diagram by Tamblyn et al. (44), where the authors propose that the semiconductor-to-metal and dissociation transitions may merge in a critical point. Notice that as mentioned before, the location of the transition is different when computed with CEIMC and BOMD. We expect the CEIMC results to be more accurate because they avoid the approximate exchange-correlation functional. Nonetheless, both methods produce the same qualitative picture of the transition.

### Melting Line.

As indicated in Fig. 1, at low temperatures the predicted LLT is expected to meet the melting line of the molecular solid. In this work, we determine the melting line of the system by comparing the free energy of the liquid phase and the free energy of the molecular solid phase (phase I, hcp rotationally disordered). For computational reasons, we limit our study to the DFT level of theory. It is well known that the semilocal exchange-correlation functionals in DFT underestimate the band gap in most semiconductors and favor delocalized states. This explains why the LLT is predicted at lower pressure with DFT than with QMC. For pressures below the metallization, the band gap is finite and the ground state properties should be accurately reproduced. As soon as the DFT-band gap closes, the nature of the ground state in this theory changes significantly and the predictions from DFT become inaccurate. This discrepancy will continue until the true band gap of the system closes. At higher densities, DFT will again produce reliable results. This is consistent with our finding here and in previous work (40) and limits the range of pressure for which DFT can be used to predict the melting line (Fig. 1).

Using thermodynamic integration with BOMD (45), we performed free energy calculations in the solid and liquid phases to determine the melting line at high pressures. We neglected quantum effects on the nuclei for these calculations. The melting line of hydrogen should be well represented by our calculations with classical protons for pressures below 200 GPa, because the system remains insulating. Unfortunately, we were unable to study the influence of nuclear quantum effects on the melting line due to the large computational demands of the PIMD method. We used 432 atoms for the calculations in the liquid and 360 atoms in the solid, assuming an hcp molecular crystal with no orientational order. Variable cell MD simulations at selected pressures were used to equilibrate the hcp cell as a function of density in the solid. These were subsequently used in simulations at constant volume. Using coupling constant integration with BOMD, we calculated the free energy difference between the system described by DFT and a reference system described by a suitable effective potential. In the solid, we used a reference Einstein crystal (46) to calculate the free energy of the reference solid. In addition, for both phases, we performed BOMD simulations in a dense grid of points for pressures between 15 and 200 GPa, and temperatures between 500 and 1000 K. The resulting set of pressures and energies, together with the results from the coupling constant integration, were used to fit a functional of the Gibbs free energy in each phase. The melting line shown in Fig. 1 was obtained from the fitted functions and is in good agreement with previous DFT simulations that were calculated with a combination of Car–Parrinello MD and the two-phase coexistence approach (32).

The melting line intersects the DFT-LLT line at 700 K and 220 GPa and should result in a liquid-liquid-solid multiphase coexistence point^{†}. Unfortunately, the use of DFT to extend the melting curve to higher pressures is likely to be highly inaccurate as it would involve an insulator-metal transition. However, the maximum in the melting line and the subsequent change in slope is a robust prediction of DFT because at this point the system is far from the LLT and the thermodynamic properties of the liquid and the solid are in excellent agreement with QMC calculations in this regime. Given this, we have chosen to extrapolate the DFT melting line to higher pressure using Kechin’s equation, *T*_{m} = 14.025 K(1 + *P*/*a*)^{b} exp(-*P*/*c*) (*a* = 0.1129 GPa, *b* = 0.7155, *c* = 149 GPa) (32, 33), and extrapolate the QMC-LLT line to higher pressures assuming that the shift between the QMC- and DFT-LLT lines are systematic; they meet at a pressure of 290 GPa and a temperature of 550 K. The uncertainty in this prediction is admittedly large because the two lines have similar slopes and such an extrapolation can be unreliable. However, this gives a QMC-based estimate of the location of the multiphase coexistence point^{‡}. This prediction is also compatible with a reported liquid state at *T* = 364 K and *P* = 335 GPa, obtained via an alternative QMC-based ab initio method (47).

In summary, we have used both DFT- and QMC-based simulations to examine the onset of molecular dissociation in hydrogen under high-pressure conditions. Our simulation results clearly indicate that for temperatures below 2,000 K, there is a range of densities where (*dP*/*dρ*)_{T} = 0 and the electrical conductivity of the fluid increases in a discontinuous fashion. These observations are consistent with the existence of a first-order phase transition in liquid hydrogen that ends in a critical point at approximately 2,000 K. By determining the melting line of hydrogen with highly accurate DFT-based free energy calculations we predict that there is a liquid-liquid-solid multiphase coexistence point in the hydrogen phase diagram at a pressure of approximately 290 GPa and a temperature of 550 K that corresponds to the intersection of the liquid-liquid, metal-insulator, phase transition with the melting curve.

## Methods

The CEIMC method was used to perform the QMC simulations. CEIMC is a first-principles simulation method that uses fixed-node reptation QMC to solve the electronic problem within the Born–Oppenheimer approximation. QMC methods for hydrogen do not rely on pseudopotential approximations and treat electron correlation explicitly by solving the many-electron problem directly. At present, it has the potential to be more accurate than DFT while being considerably less expensive than traditional quantum chemistry methods. We refer to previous publications for details about the CEIMC method (37–40). In this work, we used a more accurate wavefunction than in previous studies by using orbitals obtained from DFT with a backflow transformation. Both Jastrow and backflow functions contained, in addition to the analytic forms given in (48), short range components with free parameters optimized to minimize the energy of a set of configurations at any given density. This trial wave function provides a very accurate description of hydrogen, while maintaining high efficiency.

DFT-based first-principles simulations were performed using the Born–Oppenheimer molecular dynamics method, with the CPMD simulation package (49). The BOMD and PIMD simulations were performed with 432 atoms for the liquid and 360 atoms for the solid in the canonical ensemble, using a massive Nose–Hoover chain thermostat (50). The electronic structure was described with the Perdew–Burke–Ernzerhof exchange-correlation functional and a Troullier–Martins norm conserving nonlocal pseudopotential with a core radius of *r*_{c} = 0.5 a.u. to represent hydrogen. The equations of motion were integrated with a time step of 8 a.u. at the Γ-point with a plane-wave cutoff of 90 Ry in the DFT simulations. The path integrals in the PIMD simulations were discretized with a time step of (8,000 K)^{-1}, which corresponds to 8 imaginary time slices at *T* = 1,000 K and is sufficient to obtain converged results for this system. Electronic conductivities were calculated using DFT and the standard Kubo–Greenwood relation (51).

## Acknowledgments

This research was supported in part by the National Nuclear Security Administration under the Stewardship Science Academic Alliances program through US Department of Energy (DOE) Grant DE-FG52-06NA26170 and the Lawrence Livermore National Laboratory under Contract DE-AC5207NA27344. M.A.M. acknowledges support of a Stockpile Stewardship Graduate Fellowship; and C.P. thanks the Institute of Condensed Matter Theory at the University of Illinois at Urbana-Champaign for a short term visit, and acknowledges financial support from Ministero dell’Università e della Ricerca, Italy (Grant PRIN2007). Computer time was made available from the US DOE INCITE program, the National Center for Supercomputer Applications, Lawrence Livermore National Laboratory, and CASPUR (Italy) in the framework of Competitive HPC Grants 2009.

## Footnotes

^{1}To whom correspondence should be addressed. E-mail: ceperley{at}uiuc.edu.Author contributions: M.A.M., C.P., E.S., and D.M.C. designed research, performed research, analyzed data, and wrote the paper.

The authors declare no conflict of interest.

See Commentary on page 12743.

↵

^{*}In the solid phase the highest pressure reported so far in the diamond anvil cell experiments is of 320 GPa ; the sample still remains in the insulating state (41). This pressure value is compatible with our estimate of the pressure of metallization extrapolated from QMC data at melting.↵

^{†}This point could be either a triple point if the two liquid phases coexist with a single solid phase, or a quadruple point if the metal-insulator transition extends below freezing. Note that quadruple points are compatible with Gibb’s phase rule in a two component system. The characterization of this point requires a study of metal-insulator transition in the crystal phase, a work which is in progress.↵

^{‡}It should be mentioned that the rotationally invariant hcp-molecular phase (phase I) of hydrogen has been observed to transform into a symmetry broken phase (phase III) at*P*≲ 160 GPa (52–54). Despite a number of investigations (55, 56), the structure of phase III is still under debate. Moreover experiments have been limited to*T*≲ 250 K and it is not known whether phase III will extend in temperature up to the melting point.

Freely available online through the PNAS open access option.

## References

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- Collins GW,
- et al.

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- Ebeling W,
- Richert W

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- Scandolo S

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- Fortney JJ,
- Baraffe I,
- Militzer B

- ↵
- ↵
- Ferrario M,
- Ciccotti G,
- Binder K

- Pierleoni C,
- Ceperley DM

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- Morales MA,
- et al.

- ↵
- Frenkel D,
- Smit B

- ↵
- ↵
- ↵CPMD V3.13, Copyright IBM Corp. 1990–2008, Copyright MPI fuer Festkoerperforschung Stuttgart 1997–2001, http://www.cpmd.org/.
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵

## Citation Manager Formats

## Sign up for Article Alerts

## Article Classifications

- Physical Sciences
- Physics