Quantum simulation of low-temperature metallic liquid hydrogen

Experiments and computer simulations have shown that the melt-ing temperature of solid hydrogen drops with pressure above about 65 GPa, suggesting that a liquid state might exist at low temperatures. It has also been suggested that this low temperature liquid state might be non-molecular and metallic, although evidence for such behaviour is lacking. Here, we report results for hydrogen at high pressures using ab initio path-integral molecular dynamics methods, which include a description of the quantum motion of the protons at finite temperatures. We have determined the melting temperature as a function of pressure by direct simulation of the coexistence of the solid and liquid phases and have found an atomic solid phase from 500 to 800 GPa which melts at<200 K. Beyond this and up to pressures of 1200 GPa a metallic atomic liquid is stable at temperatures as low as 50 K. The quantum motion of the protons is critical to the low melting temperature in this system as ab initio simulations with classical nuclei lead to a considerably higher melting temperature of \sim300 K across the entire pressure range considered.

1 Ever since Wigner and Huntington's prediction in 1935 that solid molecular hydrogen would dissociate and form an atomic and metallic phase at high pressures [1], the phase diagram of hydrogen has been the focus of intense experimental and theoretical studies.
Advanced experimental techniques, notably diamond anvil cell approaches, mean that it is now possible to explore hydrogen at pressures up to about 360 GPa [2][3][4][5], and new types of diamond anvil cell may be able to access much higher pressures [6]. These experiments, along with numerous theoretical studies, have revealed a remarkably rich and interesting phase diagram comprising regions of stability for a molecular solid, a molecular liquid and an atomic liquid, and within the solid region four distinct phases have been detected [4,5,7,8]. Furthermore, metallic hydrogen has been observed in high temperature shock wave experiments [9][10][11][12] and it is accepted to be a major component of gas giant planets such as Jupiter and Saturn [13]. Despite the tremendous and rapid progress, important gaps in our understanding of the phase diagram of high pressure hydrogen remain, with arguably the least well understood issue being the solid to liquid melting transition at very high pressures. Indeed the melting curve is only established experimentally and theoretically up to around 200 GPa [14][15][16][17][18]. However, from 65 GPa up to about 200 GPa the slope of the melting line is negative (i.e., the melting point drops with increasing pressure), which suggests that at yet higher pressures a low temperature liquid state of hydrogen might exist or, as suggested by Ashcroft, perhaps even a metallic liquid state at zero K [19].
Further interest in hydrogen at pressures well above 200 GPa stems from other remarkable suggestions such as superfluidity [20] and superconductivity at room temperature [21][22][23], all of which imply that hydrogen at extreme pressures could be one of the most interesting and exotic materials in all of condensed matter.
In the current study we use computer simulation techniques to probe the low temperature phase diagram of hydrogen in the ultra-high 500-1,200 GPa regime to try and find this potential low temperature liquid state of hydrogen. We describe the proton motion using ab initio path-integral molecular dynamics (PIMD) methods [24][25][26][27][28][29][30], which are based on forces computed "on the fly" as the dynamics of the system evolves, and can account for bond formation and breaking in a seamless manner. To compute the melting line we simulate solid and liquid phases in coexistence [18,[31][32][33]. The coexistence approach minimises hysteresis effects arising from super-heating or super-cooling during the phase transition, and although it has been used with model potentials and ab initio approaches before, the current study is, to the best of our knowledge, the first to combine it with the ab initio PIMD method. With this combination of approaches we have found a low-temperature metallic atomic liquid phase at pressures of 900 GPa and above, down to the lowest temperature we can simulate reliably of 50 K. The existence of this low temperature metallic atomic liquid is associated with a negative slope of the melting line between atomic liquid and solid phases at pressures between 500 and 800 GPa. This low temperature metallic atomic liquid is strongly quantum in nature since treating the nuclei as classical particles significantly raises the melting line of the atomic solid to ∼300 K over the whole pressure range. It also destroys the negative slope of the melting line, and consequently it does not predict a low temperature liquid phase. Our study provides a finite temperature phase diagram of hydrogen under ultrahigh pressures and gives direct numerical support for Ashcroft's low temperature metallic liquid.
Extensive computational searches for low enthalpy solid structures of hydrogen have been performed using density functional theory (DFT) methods [34][35][36][37][38] and recently a systematic analysis of the evolution of the low-enthalpy phases has been provided by Labet et al. [39][40][41]. These searches found a metallic phase of I4 1 /amd space group symmetry to be stable from about 500 GPa to 1,200 GPa [35], when quasi-harmonic proton zero-point motion was included. We use this phase, therefore, as the starting point for our finite temperature exploration of the phase diagram and melting line. With the coexistence method we have performed a series of two phase solid-liquid simulations at different temperatures (from 50 K to 300 K) which are then used to bracket the melting temperature from above and below.
We begin by considering the 500-800 GPa pressure regime and show an example of the data obtained from the coexistence simulations at 700 GPa in Fig. 1. At this pressure we find that for T ≥125 K the system transforms in to a liquid state, while for T ≤100 K, it ends up as solid. The phases on either side of the melting line (at this and other pressures) were characterized by their radial distribution function (RDF) and as can be seen in Fig. 1 (d) upon moving from 100 to 125 K the system clearly transforms from solid to liquid. The phases were also characterised by the variations in mean square displacement (MSD) of the nuclei of the particles over time. Since PIMD rigorously provides only thermally averaged information we used the adiabatic centroid molecular dynamics (MD) approach within the path-integral scheme to obtain real-time quantum dynamical information [42]. Again, as shown in Fig. 1 (e), the distinction between the solid phase at 100 K and liquid phase at 125 K is clear.

3
The same coexistence procedure was used to locate the melting point at 500 and 800 GPa, leading to the melting line shown in Fig. 2. The up (down) triangles indicate the highest (lowest) temperatures at which the solid (liquid) phases are stable, bracketing the melting temperatures within a 25 K window. From this we see that the melting temperature is 150-175 K at 500 GPa and that it drops rapidly with increasing pressure yielding a melting temperature of only 75-100 K at 800 GPa. Thus the melting curve has a substantial negative slope (dP /dT < 0) in this pressure range. Across this entire pressure range the molten liquid state is atomic and the solid phase which grows is the original atomic I4 1 /amd phase that was used as the starting structure. Given that molecular phases have been observed at pressures lower than 360 GPa in both experimental and theoretical studies [3,5,43], we suggest that a molecular-to-atomic solid-solid phase transition should occur between 360 and the lowest pressure of 500 GPa considered here.
The negative slope of the melting line up to 800 GPa suggests that at even higher pressures a lower temperature liquid phase might exist. Motivated by this, we carried out simulations at 900 and 1,200 GPa. However, in this pressure range we have to consider nuclear exchange effects, which are neglected in the PIMD simulations performed here, but which could potentially become significant. Indeed, analysis of our simulations reveals that at these pressures the dispersion of the beads in the path integral ring polymer becomes comparable to the smallest inter-atomic separations when the temperature is below ∼40 K (see Fig. S6). This is the so-called quantum degeneracy temperature, below which the exchange of nuclei will be important. With this in mind we performed all simulations in this very high pressure regime at T ≥ 50 K. Interestingly we find that at 50 K at both 900 and 1,200 GPa the systems are already in the liquid state, revealing that the melting temperate at these pressures is < 50 K. Whether the liquid phase is the 0 K ground state of hydrogen at these pressures is not something we can establish at this stage. However, the large negative slope of the melting line at lower pressures and the observation of a liquid phase at temperatures as low as 50 K present strong support for Ashcroft's low temperature liquid metallic state of hydrogen [19], and it implies that any room temperature superconductor in this regime would have to be a liquid.
It is instructive to compare the results of the ab initio PIMD simulations with those obtained from the ab initio MD approach in which the nuclei are approximated by classical point-like particles. To this end we performed a second complete set of coexistence 4 simulations with ab initio MD across the entire 500-1,200 GPa range. The ab initio MD melting line is shown by the red data in the inset of Fig. 2, where it can be seen that the melting temperatures obtained from the MD simulations are much higher than those from the fully quantum PIMD simulations. The ab initio MD melting temperature is well above 200 K across the pressure range 500-1,200 GPa, and the slope of the melting line is small. These remarkably large differences between the ab initio MD and PIMD melting lines clearly demonstrate that quantum nuclear effects play a crucial role in the stability of the low temperature liquid.
In conclusion, our results are consistent with a melting curve with a negative slope between an atomic solid and a low temperature (50 K or below) metallic atomic phase. A comparison of the MD and PIMD results shows that the quantum nature of the nuclei is responsible for the large negative slope of the melting line and the consequent existence of the low temperature metallic liquid. A natural conclusion from this is that the phase diagram for hydrogen in the pressure regime considered should exhibit rather large isotope effects. It is not clear from our studies whether the liquid phase survives down to 0 K, and this needs to be studied using other methods. The occurrence of the I4 1 /amd solid above about 500 GPa indicates that a molecular-to-atomic solid-solid phase transition must occur at lower pressures, which could be accessible to experiments in the near future.

Methods
Density-functional theory (DFT) calculations with the PBE exchange correlation functional were performed with the VASP code [46][47][48][49][50]. Projector augmented wave potentials were used along with a 500 eV plane wave cutoff energy [51]. A Monkhorst-Pack k-point mesh of spacing 2π × 0.05Å −1 was used along with a unit cell that typically contained 200 atoms. Convergence tests with larger unit cells and higher plane wave cutoff energies are reported in the SI.
MD and PIMD simulations were of the "Born-Oppenheimer" type within the NVT ensemble, using a Nosé-Hoover chain thermostat [52]. A 0.5 fs time step was used and all simulations were run for at least 20,000 steps (10 ps). The PIMD simulations made use of a recent implementation of this method in VASP. 32 beads were used to sample the imaginary-time path-integral at each temperature, as this was found to be the minimum number of beads needed to obtain reliable estimates of the melting temperature (SI). When calculating real-time atomic diffusion in the liquid phase using the adiabatic centroid MD method ( Fig. 1 (e)), a smaller time step of 0.05 fs was used after thermalization and 20,000 steps were used to characterize the atomic diffusion.
In the coexistence simulations, each run was initiated from a structure comprising a solid