Theory of high pressure hydrogen, made simple

Phase I of hydrogen has several peculiarities. Despite having a close-packed crystal structure, it is less dense than either the low temperature Phase II or the liquid phase. At high pressure, it transforms into either phase III or IV, depending on the temperature. Moreover, spectroscopy suggests that the quantum rotor behaviour disappears with pressurisation, without any apparent phase transition. Here we present a simple thermodynamic model for this behaviour based on packing atoms and molecules and discuss the thermodynamics of the phase boundaries. We also report first principles molecular dynamics calculations for a more detailed look at the same phase transitions.


Introduction
In 1935 Wigner and Huntington proposed a metallic modification of hydrogen as an atomic phase at high pressure. They calculated that metallisation would occur at about tenfold compression, which is remarkably close to modern estimates. They also quoted a transition pressure of 25GPa, which has proved over an order of magnitude too low, and has been widely ridiculed. This should stand as testament that calculations of volumes are more reliable than pressures.
The experimental search for crystalline metallic hydrogen continues, and with numerous nonmetallic phases reported, it has become clear that a different physical picture from Wigner and Huntingdon's "atomisation begets metallisation" is required.
The structure of solid hydrogen up to 250GPa is now generally agreed upon. Phase I is These phases can all be reproduced by considering classical protons interacting via forces calculated in the density functional theory. The calculated pressure and temperature phase boundaries depend on the details of the calculation which typically involve trade offs between system size, exchange correlation treatment, pseudopotential fidelity, quantisation of nuclear motion and k-point sampling.
It is widely assumed that quantum protons and advanced treatment of exchange-correlation for van der Waals bonding are essential to describe bonding in solid hydrogen [8]. So it is particularly curious that the "quantum rotor" phase I is well explained by MD using classical rotors [1]. Given that the melting point maximum was predicted by classical MD simulation with remarkable accuracy, one has to believe that if those quantum effects are essential, they are exactly compensated by some incorrect piece of classical physics.
In this paper we examine what is the minimal physics required for an understanding of the hydrogen phase diagram. We do not aspire to obtain accurate pressures or temperatures, rather the questions to address are: (i) Is there a simple explanation for the mixed atomic-molecular phase IV? (ii) Does Phase I comprise free rotors at all pressures? (iii) Why is the melt denser than the close packed solid phase I?

A simple packing model
Despite the wide range of candidate structures reported from ab initio structure search techniques [4,9,10], three motifs seem to be sufficient to capture the physics of the phase diagram: rotating molecules, non-rotating molecules, and layers comprising short-lived molecules where the time-averaged indistinguishable atomic positions form hexagonal layers. Structure search reveals numerous possible small distortion which MD or path integral methods show are eliminated by elevated temperatures and/or quantum behaviour of the protons [5,6,11]. Our simple model considers three objects: • S: Spherical molecules, which correspond to a J = 0 quantum rotor ground state, or a time-averaged classical free rotor. S objects have high volume and high internal entropy. The zero of energy is taken as the covalent bond strength. • R: Rod-like molecules, corresponding to the non-rotating type, with length similar to the diameter of S. R molecules have smaller volume than S, and zero internal entropy or energy. • A: Spherical "atoms", about half that of S. They have no internal entropy, and high internal energy due to the broken bond.
In our model, phase I comprises S-type objects, phase II R-like, phase III R-type (better packed, but elongated), phase IV is a binary mixture of A and S. The liquid contains S, R and A objects according to their Boltzmann weights. A packing fraction and an internal energy is assigned to each crystal structure. The total free energy of each phase is then the sum of the energies of the crystal and its component objects. We have considered other crystalline arrangements within the model, and the only competitive structure is close-packing involving a single H atom (labeled "Atomic" in figure 1).
We can model the hydrogen molecule as two spherical atoms with centers separated by about 0.74Å. In the absence of external pressure, the density is 23 cm 3 /mol, (about 39Å 3 per molecule) indicating that the molecules are far apart compared to their bond length. At low pressure, hydrogen is extremely compressible, which implies that in this regime close-packing of a rigid spheres is not a good approximation. The close-packed density for such S-objects, i.e. where the intermolecular distance becomes close to double the bohr radius, is reached around 200GPa [7,12,13] where phase transitions start to happen. Code is available online [14], where parameters can be adjusted for fit or sensitivity analysis.
It may seem curious that there is no explicit treatment of compressibility within a phase. This is because to calculate phase boundaries using the geometric model presented here requires only that the compressibility of competing phases is similar at the phase boundary (c/f figure 2), and the additional terms serve only to complicate the model.
The hcp structure has a packing fraction of 0.7405, which is maintained under a affine deformation to nonideal c/a by transforming spheres into oblate spheroids. Naively, one might expect a transformation of all molecules directly from the S to A at a pressure which favours the smaller objects. In fact, the model predicts that a mixed S-A structure would have lower free energy for a range of pressure provided it adopts the AB 2 , MgB 2 structure, one of only two structure which has a denser packing than hcp for binary hard spheres. The time-averaged structure of Phase IV is precisely this close-packed AB 2 structure, space group P6/mmm [5,6]. Therefore the stability of this apparently complex structure can be attributed to an exceptionally simple packing model. Thus our model for Phase IV regards it as close packing of binary hard spheres: S and A objects. This picture of Phase IV as an entropically stabilised phase is consistent with classical and quantum molecular dynamics the "G" site atoms temporarily pair up to form weak molecules, but over time sample the fully symmetric phase space. If the proton motions are treated quantum mechanically, either by path integral sampling or through consideration of tunneling probabilities, the many-proton wavefunction attains P 6/mmm symmetry directly and the weak molecules manifest as a correlation in the proton wavefunction.
The transition between phases I and II occurs at much lower densities, indicating that it is driven by the interplay between quadrupole interactions and entropy, as well as packing.
A phase diagram corresponding to this model is shown in figure 1, and the thermodynamic explanation for the phase transitions in this picture is: • I → II. Energy is gained by the quadrupole-quadrupole interaction. Volume decreases due to the rods being smaller than spheres, despite inefficient packing. Entropy favours free rotor phase I. • I → IV. Volume decreases due to denser packing of two sphere types. Entropy favours free rotor phase I. • I → liquid. Conventional hard sphere melting transition at low pressure, with an increased population of type R and A objects at high pressure causing the liquid to become denser than the solid. • III→IV. Entropy gain from rotation and volume decrease due to packing atoms.

A classical-proton AIMD model for the free and arrested rotor phases
We have carried out extensive ab initio MD simulations on the phase boundaries around Phase I using the CASTEP code with PBE exchange and classical protons [14]. We carry out NPT simulations at 100GPa, starting at 0K in P 63m. NPT enables us to determine the c/a ratio of hexagonal phases, and the volume changes on transition. The rapid heating rates and finite system size mean that the transition temperatures will be an overestimate of their true  In Phase I, it is higher due to rotations. Above 1100K, a linear increase in MSD indicates melting. (d) Raman vibron frequency obtained via projection [15], dropping with increasing T, and further at the melting transition. This is consistent with experimental data [3].
On heating at 100GPa, (figure 2) we find that the symmetry-broken phase II transforms into the molecular rotor phase I at 100K, ultimately melting at 1100K. The close-packed phase has the lowest density of the three. The phase transitions are evident to movie visualisation. II → I most clearly shown in the disappearance of AAF (figure 3a), while melting becomes evident in the mean squared displacement (MSD figure 3c).
Direct measurement of quadrupole interaction energy from the electronic structure is swamped by fluctuations in the MD, but we can measure it via orientational correlations. The quantity EQQ (figure 4) is the interaction energy of two linear quadrupoles pointing in directions, k, k', separated by a vector r. The molecular quadrupole moment for hydrogen is around Q = 0.26DÅ [16], and the ensemble average of EQQ is the quadrupole-quadrupole energy. A series of simulations of pressurisation at 600K (Fig 5) showed phase I persisting until about 200GPa at which point it transformed into a disordered structure. This is not the expected melting transition, as evidenced by the enhanced AAF, indicating the arrest of the rotors, and the disappearance of the strong vibron. The structure may be similar to the chain-like structure reported elsewhere in this volume. Angular autocorrelation function (AAF), radial distribution function (RDF), and mean square displacement (MSD) graphs confirm the broken symmetry -free rotor -liquid sequence with temperature and the weakening of the vibron with increased pressure, followed by extreme broadening and further softening, plus suppression of molecular rotation in Phase III.

k-point convergence
Previous work has reported both P ca21 and P 63m as the most stable structure for phase II in this model. Preliminary calculations with a single k-point found a transition to P 63m, however at full k-point convergence P ca21 is more stable. This conclusion was tested in static relaxations of MD snapshots, and is the likely resolution to the historical discrepancy. The description of phase I appears to be converged with Γ-point sampling, but k-point sampling issues also appear in the melting. Sampling with a single k-point we observed melting at 1100K; this calculation is essentially a repeat of previous work [1]. However, using a denser k-point mesh systematically raised the melting point by several hundred degrees. It is unusual to require dense k-point sampling in molecular crystal, but careful convergence is needed to properly resolve the quadrupole moments in phase II, which are crucial to the ordering. Meanwhile, in the liquid phase, free diffusion allows the atoms to arrange themselves in a way most favourable for minimising the Γ-point electron energy at the expense of unsampled contributions from other parts of the Brilloiun Zone.

Lattice parameters, local and long ranged order
Our simulation spontaneously transformed from P ca21 to hcp for phase I. We find the c/a ratio to be essentially independent of temperature, but to decrease with increased pressure from ideal to about 1.57 at 150GPa. This, and the PV equation of state, are in good agreement with the experiment [17].
Phase I has no long range order, but the negative EQQ (Fig. 4) indicates short ranged order and increasing coupling between the rotors with pressure of a few meV, increasing approximately as V − 5 3 with pressure. As the correlation increases, there will be a gradual breakdown in J as a good quantum number and consequent disappearance of the distinctive roton bands. Although there will be a distinctive signature in spectroscopy, this correlation between rotors does not become long-ranged and so does not constitute a thermodynamic phase transition.
At low temperature, EQQ is sufficient to cause long-range ordering in the broken-symmetry phase II: The binding due to EQQ is up to 30meV, enough to stabilise P ca21 below 100K.

Conclusions
We have presented a model which distills the essential physics determining the phase stability in hydrogen. Phases II and III are energetically favored, through quadrupole interaction and electron delocalisation respectively. Phases I and IV are entropically favoured, primarily due to the high entropy of rotating molecules, which can be understood classically or quantum mechanically. The J=0 quantum rotor is an intrinsically large object: the higher compressibility of the liquid compared to phase I can be understood by the increased population of "smaller" molecules in higher energy states.
The success of the classical MD models can be attributed to this near equivalence of quantum uncertainty described by the wavefunction, and classical uncertainty described by the entropy. Phase II is unstable with Γ-point sampling (red), but stabilised by slower heating and more kpoint sampling (green). The blue dot shows EQQ for P ca21. The free rotors in phase I become increasingly correlated as pressure increases, due to the strengthening quadrupole moment and steric effect. The many-body wavefunction for this behaviour is not presented here, but clearly the excited states observable in spectroscopy should not have the characteristic overtones of a free rotor. This change in behaviour does not indicate a phase transition.