Onset of metallic transition in molecular liquid hydrogen

Liquid-liquid phase transition of hydrogen is at the center of hydrogen phase diagram as a promising route towards emergent properties such as the Wigner-Huntington metallization, superconductivity, and superfluidity. Here we report a study on the liquid-liquid phase transition of hydrogen using the state-of-the-art diffusion quantum Monte Carlo and density functional theory calculations. Our results suggest that the metallization process happens at lower pressures and temperatures compared to the structural phase transition of molecular to atomic hydrogen. The consequence is that metallized molecular hydrogen is stable at a wide range of pressures and temperatures. Our study breaks the conventional assumption that metallization coinciding with dissociation of hydrogen molecule, and the molecular metallic hydrogen liquid phase is likely to become the frontier of studying hydrogen phase transitions.

Hydrogen, the simplest and lightest element in nature, is famous for a complex phase diagram and a rich variety of phase transitions. Among the many phases of hydrogen, the liquid-liquid phase transition (LLPT) is of particular interest ever since metallic hydrogen was predicted in 1935, and LLPT is considered as a path to metallic hydrogen [1,2]. Experimentally, the LLPT of hydrogen has been achieved under dynamic and static compression in the fluid phase, and the evidence of metallization was mostly provided by the change in optical reflectivity [3][4][5][6][7][8][9]. Nevertheless, challenges and controversies still exist in experiments especially at low temperatures, including the uncertainty in determining the extreme pressure and temperature conditions, the difficulties in structure characterizations, and the lack of accurate transport measurements of conductivity [1,10,11].
Theoretically, LLPT of hydrogen is often described as a transition from the insulating molecular phase to the metallic atomic phase, supported by computational simulations [12][13][14][15][16][17][18][19][20], indicating or assuming that the molecular-atomic transition (MAT) and the insulatingmetallic transition (IMT) occur simultaneously. In recent years, new attention has been paid to the supercritical behavior and the location of liquid-liquid critical point in the computational community, involving studies using density functional theory (DFT), machine learning potentials and quantum Monte Carlo (QMC) [21][22][23][24][25][26]. Based on these studies, it is understood that below the critical point, if the temperature is above the melting line, LLPT is one phase transition where H 2 molecule dissociates and metallizes simultaneously; but above the critical point, the dissociation and metallization happen in a smooth and continuous manner. In literature, both the molecular-atomic transition and the insulatingmetallic transition have been analysed in detail using the state-of-the-art computational methods (see e.g. Refs. 22 and 27). However, to the best of our knowledge, there is no theoretical guarantee or experimental evidence that the metallization of molecular hydrogen should exactly coincide with the dissociation, regardless if LLPT is first order or not. On the contrary, experimental studies performed in the solid regime of dense hydrogen has suggested that metallization can happen without a molecular-atomic transition, motivating a re-examination of the nature of LLPT [28].
In this paper, we study the molecular-atomic and insulating-metallic phase transition of liquid hydrogen, to answer the basic question of whether they coincide with each other or one occurs before another. To provide microscopic insights, we propose the molecular fraction of liquid hydrogen as the key order parameter to compare the two kinds of phase transitions on  the same foot. We use the state-of-the-art diffusion quantum Monte Carlo (DMC) method to benchmark DFT functionals and find that vdw-DF can accurately describe the structural phase transition of dense liquid hydrogen. We also study the electronic phase transition of liquid hydrogen by calculating the fundamental gap as a function of molecular fraction using DMC and DFT methods. We find that the insulating-metallic transition begins when around 20% hydrogen molecules dissociate into hydrogen atoms, which requires much lower dissociation than the half molecular-atomic phase transition. By mapping such results onto the pressure-temperature phase diagram, we estimate that the pressure of IMT is tens of GPa lower than structural MAT.
The interactions between the valence electrons were treated with Hamann-Schlüter-Chiang-Vanderbilt pseudopotentials [36,37]. Each simulation was performed in the canonical (NVT ) ensemble, and the temperature was controlled by a Nosé-Hoover thermostat [38]. DFT total energy calculations with the PBE, the vdW-DF, the BLYP-D2, the SCAN, the HSE06 [39] functionals were performed with both the QE package and the VASP code [40]. DMC calculations were performed using the CASINO package [41], with the size-consistent DMC algorithm ZSGMA [42]. The recently developed energy consistent correlated electron pseudopotentials (eCEPPs) [43] were used. Further computational details can be found in the Supplemental Material [44].
To investigate the molecular-atomic phase transition we use an order parameter, namely the molecular fraction. As suggested in previous studies, a molecular fraction of 0.5 is a good criterion to characterize the structural molecular-atomic phase transition [17,23,45].
Specifically, we follow the definition of Cheng et al. [23], where the molecular fraction x was defined as the fraction of atoms with one neighbour within a smooth cutoff function starts from 0.8Å and decays to 0 at 1.1Å (see Supplementary Material II for more details and further validation using an agnostic classification scheme [44]). There are other definitions of molecular fraction [17,45] and other order parameters, e.g. density/volume, heat capacity, compressibility, and diffusion coefficient, that can be used to describe the molecular-atomic transition on the temperature pressure phase diagram [21,23,24]. A detailed discussion on the thermodynamic properties involved in the molecular-atomic LLPT is beyond the scope of this study. In this study we focus on the molecular fraction, which provides the clearest microscopic insight to the molecular-atomic transition. Fig. 1 (a) and (b) show snapshot of two structures with a high molecular fraction of (x=0.96) and a low molecular fraction (x=0.48), respectively.
In Fig. 1 (c) we plot the molecular fraction as a function of temperature by performing systematic AIMD simulations using different DFT exchange correlation functionals at the density of r s = 1.50 Bohr with a system size of 128 H atoms. Tests on e.g. size effects and density dependence are further presented and discussed in Fig. S2 and S3 [44]. The simulations confirm that the molecular fraction can be used to monitor the molecular-atomic LLPT as a function of temperature. The structural transition is further illustrated by the pair correlation functions in Fig. 1 [44]). The relative binding energy is Among the various functionals considered, we find that PBE predict the lowest phase transition temperature and vdW-DF2 has the highest phase transition temperature. The other three functionals, namely SCAN, vdW-DF, and BLYP-D2, show behaviors that are in between the cases of PBE and vdW-DF2. Overall, these computational observations of the qualitative difference among different DFT functionals agree with the comparisons reported in literature [17,46].
The different results from AIMD simulations suggest that it is necessary to further evaluate the accuracy of DFT functionals for liquid hydrogen. Therefore, we calculate the binding energies of a set of structures via DMC and use them as the benchmark to evaluate the performance of DFT functionals, a strategy that has been adopted successfully for solid hydrogen [46][47][48]. In Fig. 2, we plot the relative binding energies of liquid hydrogen structures as a function of the molecular fraction with respect to the DMC calculations ). The binding energy is defined as E liquid is the total energy computed using DFT or DMC, and E H is the energy of a single H atom. As one can see, vdw-DF and BLYP+D2 describe the interactions in liquid hydrogen most accurately among all the functionals. From the atomic side to the molecular side, the relative binding energy differs by less than 10 meV per atom, and the minor slope of the curve also benefits the prediction of phase boundary between the molecular phase and the atomic phase. We also find that PBE and SCAN tend to overestimate the stability of liquid hydrogen, and a large positive slope means an severe overestimate of the relative stability of the atomic phase. In contrast, vdw-DF2 underestimates the stability of liquid hydrogen and the underestimation of the stability of the atomic phase is stronger. Based on the benchmark, we conclude that vdw-DF and BLYP+D2 describe the structural phase transition of liquid hydrogen most accurately, whereas PBE and SCAN tend to facilitate the molecular-atomic transition and vdw-DF2 hinders the hydrogen molecules from dissociating into atoms. The benchmark conclusions are in line with the AIMD results: the vdw-DF and BLYP+D2 results coincide with each other within a small variation; PBE and SCAN predict lower phase transition temperatures while vdw-DF2 has a higher phase transition temperature. Note that the DFT methods that we used are different in both the nature of the underlying exchange correlation functionals and the van der Waals correction scheme.
The importance of van der Waals interactions in dense hydrogen have been mentioned in some studies, however further analyses show that the differences in describing the molecular and atomic liquid hydrogen are mainly dominated by the underlying exchange correlation functional while the van der Waals interactions do not play a significant role (Fig. S4 [44]).
We now consider the electronic insulating-metallic phase transition of liquid hydrogen.
The fundamental gap of liquid hydrogen as a function of molecular fraction calculated by DMC is shown in Fig. 3, computed according to ∆ = E(N + 1) + E(N − 1) − 2E(N). We find that the fundamental gap ∆ closes at molecular fraction x = 0.8 (red dashed line), which is far from the half molecular-atomic phase transition boundary x = 0.5 (green dashed line).
A fraction of x = 0.8 means liquid hydrogen becomes conductive when about 20% hydrogen molecules dissociate into atomic hydrogen. Note that the structures used in computing the fundamental gap are collected from simulations covering a wide range of temperatures and pressures [44]. Therefore, despite that hydrogen metallizes easier under higher pressures, the molecular fraction can be a robust order parameter to characterize the metallization of hydrogen.
We can also compute the fundamental gap using DFT methods. As show in Fig. 3 HSE06 and vdW-DF underestimates the gap by up to 30%, whereas PBE may underestimate by more than 50%. However, an interesting fact is that the molecular fraction where the fundamental gap closes is independent of the calculation methods. All DFT functionals predict consistently a threshold of molecular fraction x = 0.8 for gap closing, therefore we can directly calculate the fundamental gap using more economic methods to test the size effects in fundamental gap calculations. So we select more configurations from vdw-DF AIMD simulations with 512 hydrogen atoms at different densities (r s = 1.35 Bohr,r s = 1.50 Bohr and r s = 1.65 Bohr respectively) and calculate their fundamental gaps using PBE, as shown in Fig. S5 [44]. We find that the insulating-metallic phase transition boundary is always around molecular fraction x = 0.8.
A clear prediction from the above results is that dense liquid hydrogen can be metallized before the molecular-atomic phase transition, which implies that experiments can find metallic hydrogen at lower temperatures/pressures than atomic hydrogen. To provide a more quantitative prediction, we draw the insulating-metallic (x = 0.8) and molecular-  [4,49] are presented by the dark red points.
atomic (x = 0.5) phase transition boundaries in the pressure-temperature phase diagram, using AIMD simulations with the benchmarked vdW-DF functional (See Supplementary Material I and Fig. S6 for more details [44]). As shown by the pink shaded area in Fig. 4, at a constant temperature the insulating-metallic transition pressure is shifted downward  GPa from the molecular-atomic transition pressure. At a constant pressure, the insulatingmetallic transition temperature is lower than the molecular-atomic transition temperature 8 by a range of 300 to 1000 K. The two phase boundaries mean that liquid hydrogen can be divided into three phases, namely the molecular insulating hydrogen (MIH), the molecular metallic hydrogen (MMH), and the atomic metallic hydrogen (AMH), as labelled in Fig. 3 and Fig. 4. Molecular semi-metallic phase has been observed in the low temperature solid regime of dense hydrogen [28], and partially dissociated molecular phases of solid hydrogen have also been discussed experimentally and theoretically [50,51]. Here our prediction suggests that there is a wide range of molecular metallic hydrogen in the liquid regime, where we can further explore the intriguing electronic transitions of hydrogen. We note that a different definition of molecular fraction was used in Ref. 17 to determine the MAT boundary, and our estimates agree quite well with theirs, suggesting that the conclusions reached are insensitive to the definition of the molecular fraction.
In experiments of liquid hydrogen, most of the measures are monitoring the metallization instead of the structural phase transition. Therefore, the separation of the two phase transition boundaries allows us to make further quantitative comparison with experimental data, some of which are also plotted in Fig. 4. Apart from the original results of Knudson et al., which was later re-interpreted by Celliers et al. [9], our new prediction of the insulating-metallic transition boundary are in better agreements with experimental data [9,49]. Vibrational measurements using e.g. Raman and Infrared spectroscopy have been performed to identify structural transitions of dense solid hydrogen [52], and further extension of such techniques are desired to determine the molecular-atomic LLPT. On the contrary, theoretical studies of LLPT have mostly focused on structural phase transitions, and there is space for further developments of theoretical methods, e.g. employing other electronic structure calculations. Nuclear quantum effects can also lead to a shift of both boundaries by a small amount [16,20], but it remains to be investigated whether nuclear quantum effects on MAT and IMT are effectively the same, and how sensitive are nuclear quantum effects to the choice of DFT functional.
To conclude, we have performed DMC calculations and AIMD simulations of liquid hydrogen at a wide range of pressures and temperatures. Our main finding is that the LLPT of hydrogen is in fact two separate transitions that do not coincide with each other. Specifically, the insulting-metallic transition occurs before the molecular-atomic phase transition as the temperature and pressure increases, leading to a regime of molecular metallic hydrogen with a width of 30-60 GPa and 300-1000 K. Our results provide an important addition to the current understanding of the phase diagram of hydrogen. Last but not the least, the onset of metallic transition in molecular liquid hydrogen is encouraging for the experimental chasing of the holy grail of metallic hydrogen and the seeking of high temperature superconductor in dense hydrides.