Structural and electronic properties of solid molecular hydrogen from many-electron theories

We study the structural and electronic properties of phase III of solid hydrogen using accurate many-electron theories and compare to state-of-the-art experimental findings. The atomic structures of phase III modelled by C2/c-24 crystals are fully optimized on the level of second-order perturbation theory, demonstrating that previously employed structures optimized on the level of approximate density functionals exhibit errors in the H$_2$ bond lengths that cause significant discrepancies in the computed quasi particle band gaps and vibrational frequencies compared to experiment. Using the newly optimized atomic structures, we study the band gap closure and change in vibrational frequencies as a function of pressure. Our findings are in good agreement with recent experimental observations and may prove useful in resolving long-standing discrepancies between experimental estimates of metallization pressures possibly caused by disagreeing pressure calibrations.

The seminal work of Wigner and Huntington -that first predicted a metallization of hydrogen [1] in 1935 at a pressure of about 25 GPa -has sparked a continuous interest in the pressure-temperature phase diagram of hydrogen. However, state-of-the-art experiments [2][3][4] have not been able to conclusively detect metallic behaviour with the exception of some recent experimental studies [5][6][7] that are still under debate [8,9]. Until today, one of the most reliable experimental estimates for the metallization pressure range is approximately 425 GPa-450 GPa [7]. The lower value was obtained by the discontinuous pressure evolution in the infrared absorption, assuming a structural phase transition to the atomic structure, whereas the higher value was obtained by extrapolation of the band gap, assuming hydrogen remains in phase III. Determining the metallization pressure accurately is extremely challenging. This is partly reflected by the disagreement of the measured H 2 vibron frequency peaks as a function of the pressure, which is crucial for pressure calibration in many experiments [9]. In addition to the electronic structure, questions concerning the atomic structure are also difficult to address. Using Xray scattering to determine the crystal structure experi-mentally is hampered by the low scattering cross section of hydrogen. Depending on pressure and temperature, hydrogen has been predicted to condense in different orientationally ordered molecular crystals [10][11][12][13][14][15][16] or (liquid) metallic [1,2,5,[17][18][19][20][21][22] phases.
Accurate theoretical predictions of the equilibrium phase boundaries and other properties of high pressure hydrogen require an appropriate treatment of quantum nuclear and many-electron correlation effects [23][24][25][26][27][28], which can only be achieved using state-of-the-art ab initio methods. Hitherto, most ab initio studies of solid hydrogen are based either on density functional theory (DFT) [15,18,29,30] or quantum Monte Carlo calculations [24][25][26][31][32][33][34]. DFT employing approximate exchange and correlation (XC) energy functionals can be applied to compute infrared and Raman spectra as well as equilibrium phase boundaries, facilitating a direct comparison between theory and experiment [10-12, 16, 35-38]. However, different parameterisations of the XC functional in DFT yield inconsistent predictions [26,32,39]. Diffusion Monte Carlo (DMC) produces more reliable pressure temperature phase diagrams [24-26, 33, 34]. Furthermore DMC can also be used to compute quasi arXiv:2011.09529v1 [cond-mat.mtrl-sci] 18 Nov 2020 particle gaps including nuclear quantum effects [40]. Recently we have shown that coupled cluster singles and doubles (CCSD) theory predicts static lattice enthalpies of solid hydrogen phases with high accuracy and computational efficiency [41]. CCSD results for the most stable model phases including phase II and III are in good agreement with those obtained using diffusion Monte Carlo. However, these studies are based on structures optimised using approximate XC functionals, causing uncontrollable errors when comparing computed transition pressures or band gaps to experiment. Here, we employ accurate many-electron theories to predict the atomic structure of crystalline molecular hydrogen phases and related properties, enabling a more rigorous study of band gaps and vibrational frequencies.
Methods. We optimise the atomic structure of model phase III using nuclear gradients calculated on the level of second-order Møller-Plesset (MP2) perturbation theory and a plane wave basis set [42]. We note there are some earlier implementions of MP2 forces in periodic solids using Gaussian basis set [43][44][45]. All periodic calculations have been performed using the Vienna ab initio simulation package (VASP) [46] in the framework of the projector augmented wave method [47], interfaced to our coupled cluster code [48] that employs an automated tensor contraction framework (CTF) [49]. We use Hartree-Fock orbitals in all post Hartree-Fock methods [50]. Computational details are discussed in Ref. [42]. Although MP2 theory can be considered a low-order approximation to CCSD theory, it predicts lattice constants for a wide range of solids with higher accuracy than DFT-PBE when compared to experiment [51]. Due to the manyelectron nature of the employed Ansatz, CCSD theory is exact for two-electron systems. The coupling between electron pairs is, however, approximated by truncating the many-body perturbation expansion in a computationally efficient manner and performing a resummation to infinite order of certain contributions only.
Phase III is modelled by C2/c-24 crystals [15] initially predicted by ab initio simulations and random structure searches [15,24]. The structure is labelled by its symmetry followed by the number of atoms in the primitive cell. C2/c-24 consists of layered hydrogen molecules whose bonds lie within the plane of the layer, forming a distorted hexagonal shape. We note that previous DMC studies employed structures that have been optimized using a range of approximate density functionals, indicating that an appropriate choice is crucial [26]. In this work we employ supercells containing up to 96 atoms for the relaxation of the atomic positions. The convergence with respect to computational parameters such as number of virtual orbitals, k-meshes for the Hartree-Fock energy contribution and energy cut offs for the employed plane wave basis set have been checked carefully and are summarized in the supplemental materials [42].
Beyond the static lattice model, the T -dependent band gap renormalization of the single particle excitation energy due to electron-phonon interactions (EPIs) was also studied, using a dynamical extension of the static EPIs theory originally proposed by Heine, Allen, and Cardona (HAC) [52,53]. The quasiparticle approximation (QPA) was used to correct the DFT-PBE eigenvalues based on the EPIs self-energies. These calculations are performed using QUANTUM ESPRESSO (QE) [54] and YAMBO [55,56]. The excitonic effects were obtained by solving the Bethe-Salpeter equation (BSE), as implemented in VASP. The EPIs and the excitonic effects are calculated using DFT-PBE optimized primitive cell structure.
Results. We fully relax the internal degrees of freedom of DFT-PBE structures by minimizing the atomic forces computed on the level of MP2 theory, while keeping the lattice vectors fixed and maintaining the space group symmetry. The MP2 structures are published in the supplemental information alongside additional results, demonstrating that further effects resulting from the relaxation of the lattice vectors can be disregarded [42]. For the purpose of the following discussion we will focus on the shortest hydrogen bond length in these structures, which represents the most striking difference between MP2 and DFT-PBE results. At a pressure of 250 GPa, the shortest hydrogen molecule bond length in the DFT-PBE structures for phase III is 0.75 Å , whereas MP2 theory predicts 0.72 Å. Similar findings apply to the structures at other pressures. In passing we note that the shortest hydrogen molecule bond length obtained using the vdW-DF functional [57] is 0.72 Å, which is fortuitously close to our MP2 findings and agrees with findings reported in Ref. [26]. However, it is important to assess the reliability of these newly optimised structures further by comparing to CCSD results. Fig. 1 illustrates that the total MP2 energy per atom of phase III at a volume of 1.57 Å 3 /atom (corresponding to a DFT-PBE pressure of 250 GPa) is lowered by about 5 meV/atom during the structural relaxation. The initial 11 steps of the relaxation were carried out using a 72 atom supercell only, whereas all further optimization steps have been performed using a 96 atom supercell, indicating that finite size effects become negligible. The shortest bond length is only changed by about 0.01 Å betweeen the 11th and the final step. After 14 steps the remaining forces on the atoms are smaller than 0.05 eV/Å. Fig. 1 also depicts that the CCSD energy is lowered in total by 11 meV/atom during the full MP2 relaxation trajectory, which is similar to the change in MP2 theory. The latter observation is important because it demonstrates that MP2 and CCSD equilibrium structures are expected to deviate only slightly. This justifies the main assumption of the present work which states that MP2 structures for phase III are very accurate. To further substantiate this claim, we note that MP2 theory predicts lattice constants for a wide range of solids with significantly higher accu- racy than DFT-PBE when compared to experiment [51].
As a first demonstration for the far-reaching consequences of the structural changes, we discuss its impact on the quasi particle band gap of model phase III (C2/c-24). Fig. 2 depicts the electronic band structure for phase III at a pressure of 250 GPa employing the atomic structures optimised using DFT-PBE and MP2 theory. The Kohn-Sham band structures are computed using the PBE functional, exhibiting an indirect band gap with the valence band maximum at X and the conduction band minimum at L. The direct gap is located at Γ. The direct and indirect PBE band gaps for the MP2 structure are 2.97 eV and 1.9 eV, respectively. However, due to the reduced hydrogen bond length, the direct and indirect band gaps are about 1 eV larger in the MP2 structure compared to the DFT-PBE structure. We note that this increase in the band gap persists for the more accurate quasi particle band gaps computed on the level of the G 0 W 0 approximation. We stress that due to the strong dependence of the electronic gap on the pressure, an underestimation of the band gap by 1 eV results in a decrease in the predicted metallization pressure by more than 50 GPa. We note that the previously employed vdW-DF structures in Refs. [26,40] yield band gaps that agree with our findings obtained using the MP2 structures to within about 0.1 eV. The direct and indirected PBE band gaps computed using the vdW-DF structures are 2.88 eV and 1.74 eV, respectively.
We now turn to the comparison between computed G 0 W 0 band gaps and experimental findings. As shown in Ref. [40], the inclusion of zero point vibrational effects to the quasi particle gaps is crucial. At 0 K, this is termed as zero-point renormalization (ZPR). At finite T s, T -dependent band gap renormalization also exists, originating from the Fan and Debye-Waller terms as described in the dynamical HAC theory. More details can be found in Ref. [56,59]. Unfortunately a seamless inclusion of the electron-phonon coupling contributions to the band gap on the level of MP2 theory would be computationally too expensive at the moment. Therefore we estimate these renormalizations using DFT-PBE phonons and include them in the G 0 W 0 quasi particle band gaps [42]. Our calculations yield a ZPR of the direct and indirect gap of about -1 eV, which is by coincidence a similar magnitude as the band gap increase caused by structural relaxation but significantly smaller in magnitude than the -2 eV ZPR reported previously [40]. The difference of the T -dependent indirect band gap renormalizations in Ref. [40] between 200 K and 300 K is about 0.2 eV, which is an order of magnitude larger than our estimate of 0.02 eV [42]. The difference in the experimental indirect band gaps between 100 K [4] and 300 K [16] is about 0.02 eV, which agrees much better with our result. The computed G 0 W 0 gaps with EPIs are depicted in Fig. 3 for a range of pressures alongside experimental findings [4,7,16,40]. We note that the direct G 0 W 0 gaps includes ≈-0.12 eV exciton binding energy [42] in order to enable a direct comparison to the optical measurements from Ref. [4,7]. Furthermore we plot the G 0 W 0 gaps with respect to the CCSD pressures computed from the enthalpy versus volume curves, enabling an accurate and direct comparison to experimental findings. Compared  this contribution is the remaining leading order error in our study. However, the experimental metallization pressure of about 450 GPa lies within our theoretical uncertainties.
For a deeper understanding of the comparison with experiments, we also assess the reliability of the experimental pressure calibration. This is done by analyzing the dependence of the H 2 vibron peak frequency as a function of the pressure. As pointed out in Refs. [7,9,60] and depicted in Fig. 4, the currently available experimental estimates for the H 2 vibron peak frequency vary significantly at high pressures, questioning the reliability of experimentally determined pressures. Possible reasons for the experimental uncertainties are summarised in Ref. [9]. However, theoretical estimates of the vibron peak frequency with respect to pressure also vary significantly with respect to the employed XC parametrization on the level of DFT [7,26]. We have estimated the vibrational frequency for the MP2 structures by computing the MP2 and CCSD energies as a function of H 2 bond lengths around the equilibrium. Molecular orientations, locations of the centers of mass and volumes are fixed while changing the bond lengths in accordance with Ref. [26]. The change of the harmonic frequency with respect to the pressure can be used as a reliable calibration for pressures depicted in Fig. 4. We find that both the MP2 and CCSD frequencies have a similar slope as the H 2 vibron frequency peak measured by Loubeyre et. al. in Ref. [7,61]. From this we conclude that the experimental band gaps depicted Fig. 3  The approximate vdW-DF (taken from the supplemental material of Ref. [26]), MP2 and CCSD harmonic frequencies are shown by the pink, blue and green solid lines, respectively. The brown solid line shows the experimentally measured relation between the H2 vibron frequencies and pressures from Silvera (2018) [60] and Zha (2012) [35]. Another two experimental data lines are from Loubeyre (2020) [7] (red) and Loubeyre (2017) [61] (purple).
that are in good agreement with our most accurate estimates. In passing we note that despite the good agreement of vdW-DF structures with our MP2 structures, vdW-DF vibrational frequencies are in better agreement with experimental results of Ref. [35,60]. However, we argue that this agreement is most likely fortuitous, because both MP2 and the more accurate CCSD vibrational frequencies exhibit a very similar and steeper slope with respect to pressure. From the above findings, we conclude that the vibrational frequencies of high pressure hydrogen phases are very sensitive to the structural parameters and the corresponding electronic structure method. This has potential implications for estimates of the zero point motion energy contribution to the lattice enthalpies of accurate ab initio calculations of transition pressures [24]. Having established the good agreement between our pressure estimates and those reported in Ref. [7], we can also comment on the observed evidence of a phase transition at 425 GPa. As predicted by both DMC and CCSD calculations [24,26,41] at low temperature, phase III (C2/c-24) transforms into Cmca-12 at this pressure. However, these calculations have been performed using DFT optimized structures. We have investigated the lowering of static lattice enthalpies resulting from MP2 lattice relaxations for both structures at a selected volume corresponding to the DFT pressure of 450 GPa, finding that changes to the previously calculated transition pressures are negligible. This is surprising given the relatively large changes to the H 2 bond lengths.
Conclusions. Our work demonstrates the strengths and weaknesses of widely-used approximate DFT methods for simulating high pressure phases of hydrogen by comparing to more accurate results obtained using manyelectron methods including coupled cluster theory. Although approximate density functional theory is a computationally efficient tool for performing random structure searches [15], further structural optimization is required to achieve good agreement of band gaps and vibrational frequencies with experimental findings in solid hydrogen. Here, we demonstrate that periodic manyelectron perturbation theory calculations using plane wave basis sets have become increasingly efficient in previous years [41,58], making such optimisations feasible for systems with an increasing number of atoms. Our findings show that compared to MP2 theory, DFT-PBE structures exhibit too large hydrogen bond lengths causing too small band gaps. Although vdW-DF calculations predict structures that are closer to MP2 theory, vibrational frequencies that agree with experiment for a wide range of pressures can only be obtained on the level of CCSD. Furthermore we have demonstrated that the remaining leading order error of ab initio band gaps in solid hydrogen crystals is likely to originate from approximations used to estimate the EPI contributions. Nevertheless, it is worth pointing out that T -dependent fundamental band gap renormalization based on DFT-PBE structure is in better agreement with the experimental data. Combining accurate benchmark results with hybrid or non-local XC functionals using adjustable parameters could be useful for materials modeling in this case. Alternatively, machine-learning from MP2 forces or even more accurate ab initio data could be used to produce accurate potential energy surfaces and corresponding vibrational entropy contributions. Future work will focus on a seamless integration of electron-phonon interaction on the level of many-electron theories to further improve the accuracy of such ab initio simulations.