Dissociation of high-pressure solid molecular hydrogen: Quantum Monte Carlo and anharmonic vibrational study

A theoretical study is reported of the molecular-to-atomic transition in solid hydrogen at high pressure. We use the diffusion quantum Monte Carlo method to calculate the static lattice energies of the competing phases and a density-functional-theory-based vibrational self-consistent field method to calculate anharmonic vibrational properties. We find a small but significant contribution to the vibrational energy from anharmonicity. A transition from the molecular Cmca-12 direct to the atomic I4_1/amd phase is found at 374 GPa. The vibrational contribution lowers the transition pressure by 91 GPa. The dissociation pressure is not very sensitive to the isotopic composition. Our results suggest that quantum melting occurs at finite temperature.

In 1935, Wigner and Huntington [1] predicted that solid molecular hydrogen would dissociate at high pressure to form a metallic atomic solid. The properties of atomic hydrogen have fascinated high-pressure scientists and astrophysicists ever since [2,3]. Various exotic predictions have been made, such as the stability of atomic metallic hydrogen in a superfluid state or as a room-temperature superconductor [4][5][6], but neither an insulator-to-metal transition nor a molecular-to-atomic transition has yet been observed unambiguously at low temperatures.
The nature of hydrogen at high pressure is currently the subject of intense interest. Experimental studies of hydrogen and deuterium have been performed up to pressures above 300 GPa using static diamond-anvilcell (DAC) techniques [3,[7][8][9][10][11][12]. A new high-pressure phase IV of hydrogen and deuterium was recently observed, which is believed to consist of alternate layers of strongly bonded molecules and weakly bonded graphenelike sheets [8,13]. The precise pressures achieved in these experiments may have been overestimated and are still controversial [14]. Even more controversial is the suggestion that conductive dense hydrogen has been produced in room-temperature experiments [7]. However, the discovery of weak bonding in phase IV suggests that static DAC experiments could probe the conditions under which full molecular dissociation of hydrogen and deuterium take place. The results of our work corroborate this suggestion.
We have studied hydrogen in the pressure range 300-650 GPa, within which the transition from molecular to atomic structures is thought likely to occur. The most important contribution to the structural energy is the static lattice energy. The energy differences between competing phases in hydrogen are small and a very accurate description of the electronic energy is required to resolve them. We have therefore calculated static lattice energies using the diffusion quantum Monte Carlo (DMC) method, which is the most accurate method known for evaluating the energies of large assemblies of interacting quantum particles [15][16][17][18].
Experimental measurements [19][20][21] and classical molecular dynamics (MD) simulations using density functional theory (DFT) methods [22,23] suggest that the melting temperature of hydrogen increases with pressure and reaches a maximum value of roughly 1000 K at a pressure in the region of 100 GPa, whereafter it declines with pressure. Path-integral molecular dynamics (PIMD) simulations have suggested that the inclusion of the zero-point (ZP) energy of the protons reduces the melting temperature to about 160 K at 500 GPa and 100 K at 800 GPa [24]. This suggests that the inter-atomic bonding becomes very weak at these pressures, and anharmonic effects could become important.
We have performed vibrational self-consistent field (VSCF) calculations within DFT to calculate the anharmonic vibrational ZP energies [25]. We used the Perdew-Burke-Ernzerhof (PBE) generalized gradient approximation density functional, which is well suited for very highpressure studies, as the charge density is more uniform than at low densities, and it obeys the uniform limit and gives a good account of the linear response of the electron gas to an external potential [26].
Static lattice DFT calculations using ab initio random structure searching (AIRSS) [27] indicate that there are three energetically competitive structures in the range of interest. The molecular Cmca-12 phase is insulating up to 373 GPa in the GW approximation [28], although proton zero-point and finite-temperature effects are expected to lower the metallization pressure [29]. The molecular Cmca-4 phase and the atomic phase of I4 1 /amd symmetry (the structure of Cs-IV) are both metallic [30][31][32]. DFT with the PBE functional predicts that Cmca-12 is stable up to 385 GPa, Cmca-4 is stable in the range 385-490 GPa, and I4 1 /amd is stable from 490 GPa up to pressures beyond 1 TPa.
DFT studies of high-pressure phases of hydrogen have been performed using several approximate density functionals [13,24,29,33,34] and a significant dependence of the results on the functional has been noted. The enthalpy differences between phases are so small that changes of only a few meV per proton can make a noticeable difference to the phase diagram. It is therefore important to use an accurate approach for calculating the energies of the competing phases. We have chosen to use the DMC method [15,16] to calculate the static lattice energies. This method solves the many-electron Schrödinger equation, is in principle exact for the ground states of the hydrogen atom and molecule, and rigorously excludes self-interaction errors. It is likely that DMC provides a considerably more accurate description of the energetics of hydrogen than the currently available exchange-correlation density functionals.
We used the casino code [35] to perform fixed-node DMC simulations with a trial wave function of the Slater-Jastrow (SJ) form, where R is a 3N -dimensional vector of the positions of the N electrons, r ↑ i is the position of the i'th spinup electron, r ↓ j is the position of the j'th spin-down electron, exp[J(R)] is a Jastrow factor, and det[ψ n (r ↑ i )] and det[ψ n (r ↓ j )] are Slater determinants of spin-up and spin-down one-electron orbitals. These orbitals were obtained from DFT calculations performed with the plane-wave-based quantum espresso code [36], employing a norm-conserving pseudopotential constructed within DFT using the local density approximation (LDA) exchange-correlation functional. The choice of exchangecorrelation functional used to generate the orbitals has almost no effect on the DMC energies of solid hydrogen phases [28,37]. Earlier work also suggests that using a pseudopotential has only a small impact on results for high-pressure solid hydrogen [32]. We chose a very large basis-set energy cut-off of 300 Ry to approach the complete basis set limit [38], as detailed in the Supplemental Material [39]. The plane-wave orbitals were transformed into a localized "blip" polynomial basis [40]. Our Jastrow factor consists of polynomial one-body electron-nucleus and two-body electron-electron terms, the parameters of which were optimized by minimizing the variance of the local energy at the variational Monte Carlo (VMC) level [41,42]. The QMC calculations were performed with simulation cells containing N = 128 protons. We used twistaveraged boundary conditions with 24 randomly chosen twists to reduce the single-particle finite-size effects [43]. We have corrected our results for the effects of using finite simulation cells, employing the approach described in Refs. [44] and [45]. The residual finite-size effects are estimated to lead to errors in the enthalpy differences between phases of less than 5 meV per proton. The finite-size corrections are detailed in the Supplemental Material [39]. The statistical errors in our data are smaller than the reported accuracy [39].
The enthalpy was evaluated by fitting a polynomial to the finite-size-corrected DMC energy as a function of volume and differentiating the fit. The resulting enthalpypressure relations are shown in the upper plot of Fig. 1. At the static lattice level we find a transition from Cmca-12 to Cmca-4 at 431 GPa and a transition from the molecular Cmca-4 to atomic I4 1 /amd structure at about 465 GPa. DMC calculations predict that the Cmca-4 phase is significantly less stable than in PBE DFT.
To explore the accuracy of the DMC results further we performed calculations with trial wave functions incorporating an inhomogeneous backflow (BF) transformation [46], which modifies the nodal surface of the wave function and can introduce additional correlation effects. The BF wave functions give energies lower than the SJ wave functions by about 19-13 meV per proton with increasing density for Cmca-4, while the I4 1 /amd energies were lowered by about 17 meV per proton, approximately independently of density. The energy reduction in the molecular Cmca-4 phase is slightly larger than in the I4 1 /amd atomic phase, but the energy reductions of the Cmca-4 and I4 1 /amd phases are almost the same in the region of the phase transition. We conclude that the introduction of BF correlations does not significantly alter our results. Further details of the effects of BF are described in the Supplemental Material [39].
The above results are based on static lattice calculations in which the vibrational motion of the protons has been neglected. We have also performed calculations of the ZP enthalpy arising from the proton motion using (a) the quasiharmonic approximation and (b) a VSCF approach that enables the calculation of anharmonic vibrational energies [25]. The quasiharmonic phonon calculations were performed with the PBE functional using both the supercell finite displacement method and density-functional perturbation theory (DFPT) as implemented in quantum espresso [36].
Previous calculations of quasiharmonic proton ZP energies in solid hydrogen have encountered significant numbers of unstable phonon modes [31,32] at high pressures. We found that the Cmca-12 and I4 1 /amd structures had stable modes at the supercell sizes considered. For Cmca-4 we found a small unstable region around the Γ point which was further reduced, but not entirely eliminated, by using cells with up to 256 protons. This small unstable region does not affect our estimates of the ZP energy, as shown in the Supplemental Material [39]. As illustrated in the inset of the lower panel of Fig. 1, the proton ZP enthalpy of all three phases increases with pressure.
Systems of light atoms with weak bonding often exhibit large vibrational amplitudes, which are likely to give rise to anharmonic vibrations. There is evidence for the im- portance of anharmonicity in hydrogen, especially in the high-density regime [2,29]. Utilizing our recently developed variational VSCF scheme [25,47], we have calculated the anharmonicity of the proton ZP motion of both the molecular and atomic phases. We use the principal axes approximation to map the Born-Oppenheimer energy surface along independent but anharmonic vibrational modes [25,48], and solve the resulting equations within a VSCF scheme [25,49]. We also calculate the contribution from phonon-phonon two body coupling in the most anharmonic modes to estimate the effects of these terms on the anharmonic vibrational energy [39]. The Born-Oppenheimer energy surface is mapped within plane-wave DFT using the castep code [50]. By comparing the energies of the highest-and lowest-energy modes with those of the static lattice, we estimate that our choice of computational parameters leads to energy differences between frozen phonon configurations that are converged to within 10 −4 eV/proton. All calculations were performed with the PBE functional and supercells containing 96 and 108 atoms. Supercell-size convergence tests indicate that the anharmonic ZP energy correction is accurate to within 1 meV/proton for all three phases. . We performed calculations for the atomic I4 1 /amd phase at pressures of P = 400, 500, and 600 GPa, obtaining anharmonic corrections of −7.2, −8.1, and −7.3 meV/proton, respectively. Similar calculations for the molecular Cmca-4 phase at pressures of P = 400 and 500 GPa give anharmonic corrections of +8.7 and +8.3 meV/proton, respectively. Calculations at P = 400 GPa for the Cmca-12 structure lead to an anharmonic correction of +4.0 meV/proton. The anharmonic corrections to the proton ZP energy lower the energy of the atomic phase and raise the energy of the molecular phases.
As an example of the vibrational properties of the I4 1 /amd structure, we plot in Fig. 2 the Born-Oppenheimer energy surface and the corresponding anharmonic wave function density at P = 500 GPa for a Γ-point optical phonon of the I4 1 /amd structure, and a comparison with the harmonic quantities. The I4 1 /amd structure can be viewed as a sequence of four stacked planes with square lattices as shown in the lower panel of Fig. 2. The mode corresponds to an in-plane motion of the protons, where alternate stacked planes oscillate in opposite directions. The minima of the anharmonic po-tential are separated by about 0.427Å in real space. Adjacent minima correspond to equivalent I4 1 /amd structures connected by this in-plane proton motion. This is the mode with the largest anharmonicity in the I4 1 /amd structure.
We note that the fermionic nature of the protons has not been taken into account at either the harmonic or anharmonic levels. In order to estimate the effects of this approximation, we consider the real space amplitude of the motion of the protons about their equilibrium positions. The root-mean-square atomic amplitude in the I4 1 /amd phase at P = 500 GPa is u 2 = 0.099Å, which is much smaller than the nearest neighbour distance of a = 0.99Å. This indicates a small overlap between the proton wave functions and justifies the neglect of their quantum statistics. The Lindemann criterion [51] for the melting of a solid is usually taken to be u 2 /a > ∼ 0.1, but for quantum melting a value of u 2 /a > ∼ 0.25 is considered more accurate [52][53][54]. From our data [39], u 2 /a = 0.1 at T = 0 K, which is not in the regime of a zero-temperature quantum liquid [6]. The inclusion of quantum statistics would increase the kinetic energy of the liquid and make it even less favorable at zero temperature. The Lindemann criterion applied to the corresponding atomic phase for the heavier deuterium would lead to a higher melting temperature. The Lindemann criterion cannot definitely answer the question of whether the ground-state atomic phase is a metallic solid or a quantum liquid, but it suggests that the melting temperature is higher than zero. These conclusions should be compared with recent PIMD results at P = 500 GPa [24], which suggest a melting temperature of 160 K.
We estimate the dynamical enthalpy as the sum of the static DMC enthalpy and the ZP enthalpy calculated using the quasiharmonic approximation corrected by the VSCF scheme to account for the effects of anharmonicity. The contribution to the total pressure from the ZP motion for the atomic I4 1 /amd phase increases with pressure from 25 to 45 GPa over the pressure range of the inset in the lower panel of Fig. 1, while the ZP pressure of the molecular Cmca-4 and Cmca-12 phases increases from 19 to 24 GPa and 10 to 13.5 GPa, respectively. As shown in the dynamical phase diagram of Fig. 1, the transition from the Cmca-12 to the atomic I4 1 /amd phase occurs at 374 GPa, and there is no stability region for the Cmca-4 phase.
At the static level, the phase diagrams of hydrogen and deuterium are identical. At the dynamic level and using the quasiharmonic approximation, the deuterium ZPE is calculated by dividing the hydrogen ZPE by √ 2. As illustrated in the Supplemental Material [39], the dynamical phase diagram of deuterium shows that the molecular-to-atomic phase transition happens at a pressure of 390 GPa, also through a structural transformation from molecular Cmca-12 to atomic I4 1 /amd. Therefore, the molecular-to-atomic phase transition is fairly isotopeindependent.
We find that our DMC results are essentially independent of the exchange-correlation (XC) functional used to calculate the orbitals for the trial wave function. However, the value of the proton ZP energy depends on the choice of XC functional. To investigate the effect of this we have recalculated the harmonic ZP enthalpy for the Cmca-4, Cmca-12 and I4 1 /amd structures using the BLYP XC functional [55], as detailed in the Supplemental Material [39], which gives significantly different results from the PBE functional at the static lattice level. The phase diagram including the effects of the ZP harmonic enthalpy calculated with the BLYP functional leads to a reduction of 28 GPa in the transition pressure for hydrogen dissociation, compared to the PBE-based phase diagram. The differences in dissociation pressure due to the flavour of XC functional used for the treatment of atomic vibrations do not affect the qualitative results presented in this work, and only have a limited quantitative effect.
In conclusion, we have studied the dissociation of solid molecular hydrogen at the static lattice and dynamical lattice levels. At the static lattice level our calculations give a transition from the Cmca-12 molecular phase to the Cmca-4 molecular phase at P = 431 GPa, and a transition to the I4 1 /amd atomic phase at 465 GPa. At the dynamical level the molecular Cmca-12 phase transforms directly to the atomic I4 1 /amd phase at 374 GPa. The limited precision of our calculations prevents us from stating categorically that the Cmca-4 phase does not exist, but the pressure range over which it might exist is very narrow. The atomization pressure is close to being within range of DAC experiments [56]. Therefore the low temperature molecular-to-atomic phase transition of high pressure hydrogen might be observable experimentally. By comparing the dynamical phase diagrams of hydrogen and deuterium, we predict that the molecular-to-atomic phase transition is almost isotope-independent. The proton ZP vibrational energies increase with pressure and the anharmonic contribution leads to an increase in the vibrational energy of the molecular Cmca-4 and Cmca-12 phases and a decrease in that of the I4 1 /amd atomic phase. Our results suggest that quantum melting of hydrogen would occur at finite temperature. Since metallic hydrogen is thought to be present in large amounts in the interiors of Jupiter, Saturn, and some extra-solar planets, planetary models should consider incorporating our prediction of the existence of an atomic metallic state at lower pressures than previously assumed.
We acknowledge the financial support of the UK Engineering and Physical Sciences Research Council under grants EP/I030190/1 and EP/I030360/1, the use of the HECToR computing facilities (the UK national supercomputing service), the Imperial College London High Performance Computing Centre, and the Cam-bridge High Performance Computing Service. We acknowledge support from the Thomas Young Centre under grant TYC-101