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

## Abstract

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 - direct to the atomic phase is found at GPa. The vibrational contribution lowers the transition pressure by 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-anvil-cell (DAC) techniques (7); (8); (9); (10); (11); (3); (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 graphene-like 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 – 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 high-pressure 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 - phase is insulating up to GPa in the approximation (28), although proton zero-point and finite-temperature effects are expected to lower the metallization pressure (29). The molecular -4 phase and the atomic phase of symmetry (the structure of Cs-IV) are both metallic (30); (31); (32). DFT with the PBE functional predicts that -12 is stable up to 385 GPa, -4 is stable in the range 385–490 GPa, and is stable from 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,

(1) |

where is a -dimensional vector of the positions of the electrons, is the position of the ’th spin-up electron, is the position of the ’th spin-down electron, is a Jastrow factor, and and 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 exchange-correlation 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 protons. We used twist-averaged 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 enthalpy-pressure relations are shown in the upper plot of Fig. 1. At the static lattice level we find a transition from -12 to -4 at 431 GPa and a transition from the molecular -4 to atomic structure at about 465 GPa. DMC calculations predict that the - 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 -, while the energies were lowered by about 17 meV per proton, approximately independently of density. The energy reduction in the molecular - phase is slightly larger than in the atomic phase, but the energy reductions of the -4 and 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 - and structures had stable modes at the supercell sizes considered. For - 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 importance of anharmonicity in hydrogen, especially in the high-density regime (29); (2). 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 eV/proton. All calculations were performed with the PBE functional and supercells containing and atoms. Supercell-size convergence tests indicate that the anharmonic ZP energy correction is accurate to within meV/proton for all three phases. Further details of the VSCF calculations are given in the Supplemental Material (39). We performed calculations for the atomic phase at pressures of , , and GPa, obtaining anharmonic corrections of , , and meV/proton, respectively. Similar calculations for the molecular - phase at pressures of and GPa give anharmonic corrections of and meV/proton, respectively. Calculations at GPa for the - structure lead to an anharmonic correction of 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 structure, we plot in Fig. 2 the Born-Oppenheimer energy surface and the corresponding anharmonic wave function density at GPa for a -point optical phonon of the structure, and a comparison with the harmonic quantities. The 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 potential are separated by about Å in real space. Adjacent minima correspond to equivalent structures connected by this in-plane proton motion. This is the mode with the largest anharmonicity in the 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 phase at GPa is Å, which is much smaller than the nearest neighbour distance of Å. 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 , but for quantum melting a value of is considered more accurate (52); (53); (54). From our data (39), at 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 GPa (24), which suggest a melting temperature of 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 phase increases with pressure from to GPa over the pressure range of the inset in the lower panel of Fig. 1, while the ZP pressure of the molecular -4 and -12 phases increases from to GPa and to GPa, respectively. As shown in the dynamical phase diagram of Fig. 1, the transition from the -12 to the atomic phase occurs at GPa, and there is no stability region for the - 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 . 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 - to atomic . Therefore, the molecular-to-atomic phase transition is fairly isotope-independent.

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 -4, -12 and 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 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 -12 molecular phase to the -4 molecular phase at GPa, and a transition to the atomic phase at GPa. At the dynamical level the molecular -12 phase transforms directly to the atomic phase at 374 GPa. The limited precision of our calculations prevents us from stating categorically that the - 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 - and - phases and a decrease in that of the 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 Cambridge High Performance Computing Service. We acknowledge support from the Thomas Young Centre under grant TYC-101

### References

- E. Wigner and H. B. Huntington, J. Chem. Phys. 3, 764 (1935).
- J. M. McMahon, M. A. Morales, C. Pierleoni, and D. M. Ceperley, Rev. Mod. Phys. 84, 1607 (2012).
- A. F. Goncharov, R. T. Howie, and E. Gregoryanz, Low Temperature Physics 39, 402 (2013).
- N. W. Ashcroft, Phys. Rev. Lett. 21, 1748 (1968).
- N. W. Ashcroft, Phys. Rev. Lett. 92, 187002 (2004).
- E. Babaev, A. Sudbo, and N. W. Ashcroft, Nature (London) 431, 666 (2004).
- M. I. Eremets and I. A. Troyan, Nature Mater. 10, 927 (2011).
- R. T. Howie, C. L. Guillaume, T. Scheler, A. F. Goncharov, and E. Gregoryanz, Phys. Rev. Lett. 108, 125501 (2012).
- R. T. Howie, T. Scheler, C. L. Guillaume, and E. Gregoryanz, Phys. Rev. B 86, 214104 (2012).
- C.-S. Zha, Z. Liu, and R. J. Hemley, Phys. Rev. Lett. 108, 146402 (2012).
- P. Loubeyre, F. Occelli, and P. Dumas, Phys. Rev. B 87, 134101 (2013).
- Chang-sheng Zha, Zhenxian Liu, Muhtar Ahart, Reinhard Boehler, and Russell J. Hemley, Phys. Rev. Lett. 110, 217402 (2013).
- C. J. Pickard, M. Martinez-Canales, and R. J. Needs, Phys. Rev. B 85, 214114 (2012), C. J. Pickard, M. Martinez-Canales, and R. J. Needs, Phys. Rev. B 86, 059902(E) (2012).
- Ross T. Howie, Eugene Gregoryanz, and Alexander F. Goncharov, J. Appl. Phys. 114, 073505 (2013).
- D. M. Ceperley and B. J. Alder, Phys. Rev. Lett. 45, 566 (1980).
- W. M. C. Foulkes, L. Mitas, R. J. Needs, and G. Rajagopal, Rev. Mod. Phys. 73, 33 (2001).
- C. Attaccalite, and S. Sorella, Phys. Rev. Lett. 100, 114501 (2008).
- V. Natoli, Richard M. Martin, and D. M. Ceperley, Phys. Rev. Lett. 70, 1952 (1993).
- S. Deemyad and I. F. Silvera, Phys. Rev. Lett. 100, 155701 (2008).
- M. I. Eremets and I. A. Trojan, JETP Lett. 89, 174-179 (2009).
- N. Subramanian, A. F. Goncharov, V. V. Struzhkin, M. Somayazulu, and R. J. Hemley, Proc. Natl. Acad. Sci. U.S.A. 108, 6014 (2011).
- S. Bonev, E. Schwegler, T. Ogitsu, and G. Galli, Nature 431, 669 (2004).
- I. Tamblyn and S. A. Bonev, Phys. Rev. Lett. 104, 065702 (2010).
- J. Chen, X. Li, Q. Zhang, M. I. J. Probert, C. J. Pickard, R. J. Needs, A. Michaelides, and E. Wang, Nature Communications 4, 2064 (2013).
- Bartomeu Monserrat, N. D. Drummond, and R. J. Needs, Phys. Rev. B 87, 144302 (2013).
- J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
- C. J. Pickard and R. J. Needs, J. Phys.: Condens. Matter 23, 053201 (2011).
- Sam Azadi, W. M. C. Foulkes, and T. D. Kuhne, New Journal of Physics 15, 113005 (2013).
- Miguel A. Morales, Jeffrey M. McMahon, Carlo Pierleoni, and David M. Ceperley, Phys. Rev. B 87, 184107 (2013), Phys. Rev. Lett. 110 065702 (2013).
- K. A. Johnson and N. W. Ashcroft, Nature 403, 632 (2000).
- C. J. Pickard and R. J. Needs, Nature Physics 3, 473 (2007).
- J. M. McMahon and D. M. Ceperley, Phys. Rev. Lett. 106, 165302 (2011).
- Sam Azadi and W. M. C. Foulkes, Phys. Rev. B 88, 014115 (2013).
- R. E. Cohen, Ivan I. Naumov, and Russell J. Hemley, 110, 13757 (2013).
- R. J. Needs, M. D. Towler, N. D. Drummond, and P. López Ríos, J. Phys.: Condens. Matter 22, 023201 (2010).
- P. Giannozzi et al., J. Phys.: Condens. Matt. 21, 395502 (2009).
- The Cmca-12 phase is metallic in the LDA but insulating up to 373 GPa in calculations. To check the accuracy of our metallic LDA-based trial functions for Cmca-12, we also used trial functions from hybrid PBE0 calculations with a much larger band gap. The difference between the two DMC energies was less than 2 meV per proton at all pressures considered.
- S. Azadi, C. Cavazzoni, and S. Sorella, Phys. Rev. B 82, 125112 (2010).
- See Supplemental Material.
- D. Alfè and M. J. Gillan, Phys. Rev. B 70, 161101 (2004).
- C. J. Umrigar, K. G. Wilson, and J. W. Wilkins, Phys. Rev. Lett. 60, 1719 (1988).
- N. D. Drummond and R. J. Needs, Phys. Rev. B 72, 085124 (2005).
- C. Lin, F. H. Zong, and D. M. Ceperley, Phys. Rev. E 64, 016702 (2001).
- S. Chiesa, D. M. Ceperley, R. M. Martin, and M. Holzmann, Phys. Rev. Lett. 97, 076404 (2006).
- N. D. Drummond, R. J. Needs, A. Sorouri, and W. M. C. Foulkes, Phys. Rev. B 78, 125106 (2008).
- P. López Ríos, A. Ma, N. D. Drummond, M. D. Towler, and R. J. Needs, Phys. Rev. E 74, 066701 (2006).
- B. Monserrat, N. D. Drummond, C. J. Pickard, and R. J. Needs, Phys. Rev. Lett. 112, 055504 (2014).
- Joon O. Jung and R. Benny Gerber, J. Chem. Phys. 105, 10332 (1996).
- Joel M. Bowman, J. Chem. Phys. 68, 608 (1978).
- Stewart J. Clark et al., Z. Kristallogr. 220, 567 (2005).
- F. Lindemann, Z. Phys. 11, 609 (1910).
- P. A. Whitlock, D. M. Ceperley, G. V. Chester, and M. H. Kalos, Phys. Rev. B 19, 5598 (1979).
- S. T. Chui, Phys. Rev. B 41, 796 (1990).
- C. A. Burns and E. D. Isaacs, Phys. Rev. B 55, 5767 (1997).
- C. Lee, W. Yang, and R. G. Parr, Phys. Rev. B 37, 785 (1988).
- L. Dubrovinsky, N. Dubrovinskaia, V. B. Prakapenka, and A. M. Abakumov, Nat. Commun. 3, 1163 (2012).