High pressure hydrogen and the Potts model on a triangular lattice

We present Monte Carlo studies and analysis of the frustrated antiferromagnetic Potts model of a triangular lattice. This Potts model shows a remarkably rich range of structures, and striking similarities to the high pressure phases of hydrogen which are typified by hexagonal close packed layered structures [1]. There are four known H2 molecular phases, all of which are isostructural to within the resolution of x-ray diffraction. Experimentally, the phase lines have been mapped by spectroscopy, which cannot reveal the structure. Study by density functional theory (DFT) has suggested a large number of candidate structures, based on the hexagonal-close packing of H2 molecules. The Potts model exhibits structures similar to DFT candidate hydrogen phases I, II and III: the range of different Potts model structures suggests that the hydrogen system in the ‘phase II’ region, may exhibit more than a single phase. It also suggests reorientational excitations which may be detectable in spectroscopy.


Introduction
In this paper we investigate an intriguing similarity between the Potts model and the behaviour of solid hydrogen.
Solid molecular hydrogen exists at low temperature and high pressures. X-ray diffraction suggests that the structure comprises molecules arranged on an hcp lattice, and four different phases have been identified based on spectroscopy [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15]. Density functional theory (DFT) calculations based on structure search have revealed several proposed phases, a common factor of which is that the molecular centres are close to an hcp structure, while the molecular axes lie in the basal plane [16][17][18][19][20]. DFT calculations for phases I-III have proposed many possible structures structures which are based on different stackings of such layers. Both calculation and x-ray diffraction show that the c/a ratio is below the ideal close-packed value of 8/3, further reducing with pressure.
The generally accepted picture is of a phase I, comprising rotating molecules, phase II with molecular ordering determined by quadrupole interactions, and phase III with molecule ordering determined by efficient packing. Raman and DFT studies suggest that the higher temperature and pressure phase IV, has a hexagonal structure with alternating layers [14,[20][21][22]. It appears to have two types of H 2 molecule with different bondlengths and strengths. A recent x-ray diffraction [1] identified three diffraction peaks which were assigned to an hcp lattice for phase IV, although the x-ray data is also compatible with various alternating-layer structures observed in ab initio molecular dynamics [23].
Density functional theory has an impressive record for determining candidate crystal structures either by structure search [24,25] at zero temperature, or molecular dynamics at higher temperatures [21,23,[26][27][28][29]. However, free energy calculations accurate enough to determine phase boundaries remain challenging. Simpler models, guided by the DFT results, can provide insight into the phase diagram, and more intuitive explanations for the structures. For example, the panoply of DFT candidates for phase II is because the orientation of quadrupoles on an hcp lattice is highly frustrated (unlike fcc) [30]. Similarly, in phase III, the 2D arrangement of quadrupoles on a triangular lattice has a single, simple minimum, but there is no efficient stacking of such layers in 3D [22].

Methodology
From studying the candidate DFT structure for phases II and III we observe that the in-plane molecules: • Are located on sites of a triangular lattice.
• Are orientated to point away from neighbouring sites, in one of three equivalent in-plane directions, or directly out of the plane.
Furthermore, the large number of near-degenerate different possible stackings of these layers suggests that the dominant contributions to the ordering hydrogen lies in the 2D in-plane arrangements.
This suggests that insights into high pressure hydrogen can be gained by analogy with the four-state Potts model, in which one of the states is non-equivalent. Under pressure, there is an increasing tendency for molecules to lie in-plane, evidenced by molecular dynamics and by the reduction of the c/a ratio with pressure [1,9]. Thus we can assume that the non-equivalent state is unfavourable at high pressure. By contrast, in the quantum limit, the four-state Potts model remains relevant, but the para state is energetically favoured. To emphasize the analogy with hydrogen, the different states are represented by molecular orientation. This highlights how elongated molecules pack in this arrangement. The arrangement is not centrosymmetic, emphasized by the touching ellipses. The hcp structure involves extended to three dimensions by stacking layers above the sites labelled B. This layered arrangement is characteristic of the low-energy phases found for phase III of hydrogen [16,19], which differ mainly in the stacking sequence.
The phase III appears as a high pressure phase: unlike phase II, the molecular orientations do not minimise the quadrupolequadrupole energy, rather they are arranged for the most effective packing, with the molecular pointing away from other in-plane molecular centres, with some out of plane tilting. By symmetry, this means a given molecule on a triangular lattice has three possible orientations. We model these as the three states of a Potts model s i .
We further observe that in all the energy-minimising DFT phases, adjacent molecules have different orientations. This Typical results from the MC run for 4-state Potts model system sizes 96 2 at various temperatures with near neighbour antiferromagnetic interactions. All data comes from calculations with 10 9 attempted flips. The noise in the B = 0 data comes from the fact that the ground state is multiply degenerate, including states in the range 0 < n 4 < 1/3, but these states are not simply connected by single spin flips.
suggests mapping the problem onto the simplest three state Potts model, which has enthalpy.
with J (assumed positive) being an antiferromagnetic coupling constant and δ being the Kronecker delta, which equals one for s i = s j and zero otherwise. i runs over all sites and j 1 over near-neighbour pairs. The fourth orientation is not symmetrically equivalent to the others, and a practical issue with incorporating a fourth state is that the ground state becomes superdegenerate unless second neighbour interactions are included. We postulate an extension of equation (1) to account for these effects: The first term is the same as the three state model, the second term extends it to second neighbours and the third term distinguishes the out-of-plane fourth state which contributes and the excess enthalpy B(V). Our interpretation is that this state increases the enthalpy via the PV term, by increasing the width of the layer. We anticipate that it will become less favoured in denser phases in the hydrogen application. Specifically, at low pressure the enthalpy (H = U + PV) is dominated by the quadrupole-quadrupole interaction (internal energy), for which we expect B(V) < 0, whereas at high pressure the packing (PV term) is important, so ∂B(V) ∂P > 0. In the limit where B(V) → ∞ and J 2 → 0 we recover the model from equation (1).

Three-state Potts model
We implemented a three state Potts model on a triangular lattice, and ran a single-flip metropolis Monte Carlo algorithm [45,46]. The ground state of the model is a hexagonal pattern with three sites per unit cell. This is illustrated in the lowest panel of figure 5 (n 4 = 0).
The MC simulations (figure 1) show that the system has a well defined order-disorder transition which is easily identified by the divergence in the specific heat. An interesting feature of this calculation comes from the study of finite size effects: the divergence is far more pronounced where the system size is a multiple of three primitive cells. Only in this case is a perfect arrangement of the s i compatible with the boundary conditions [48,49]. For other system sizes, an antiphase boundary or array of point defects must be created somewhere in the system, and because it could be anywhere, the 'ground state' is multiply degenerate. These defects can move, but their existence is topologically protected by the boundary conditions. The requirement for system sizes tuned to the expected ground state has also been noted in molecular dynamics of hydrogen [23]. Figure 2 shows how this Potts ground state relates to a molecular system. The DFT predictions for phase III are all based on layers with orderings of this type. The arrangement breaks centrosymmetry which, in molecular physics, would mean that the molecular vibration becomes infrared active-the characteristic signature of the hydrogen phase III [50,51]. In the figure this is emphasized by drawing touching  [30], which is also the ground state of the four-state Potts model, and the motif from which several candidates for phase II hydrogen are built. The spin equivalent of this structure is shown in figure 5. ellipses, whose centres are thereby necessarily displaced from the triangular lattice sites.
Considering the sites where another layer could be placed above or below the central molecular in a close packed stacking, we see the molecule has six adjacent interstitial sites one has three molecules pointing to it, three have one molecule pointing in, and two have none. We refer to these and 3, 1 and 0-interstices respectively. Note that despite the inequivalent interstices, the molecules themselves are all in equivalent environments, with one end pointing into the a 1-interstice and the other to a 3-interstice.
If one considers hcp as ABA stacking, the shown layer is A, while the B layer has both 0-and 3-interstices, meaning that it is impossible to stack this kind of layer without breaking symmetry. This gives a strong hint that there is no optimal 3D stacking. The fcc-like 'C' layer sites (1-interstices, not labelled) can be seen to be obscured and distorted by the molecular relaxation implied by the ellipses-this helps to explain why hcp-like stacking is favoured over fcc-like.

Four-state Potts model
See figure 3.

Four-state Potts model with near neighbours
The general form of the hydrogen-inspired four-state model considered here has three parameters (equation (2)).
The classic four-state model [52,53] is recovered if we set B = 0 and J 2 = 0. This has a massively degenerate ground state of little physical significance [36]. However, a small positive value of B breaks the degeneracy, and the zero-entropy ground state has n 4 = 0 for B > 0, exactly the same as for the 3-state model (figure 2). A small negative value of B 0 reduces the degeneracy with the n 4 = 1/3 becoming ordered, similar to candidate structures for hydrogen phase II (figure 4) but there are still many degenerate options for the other three states (see figure 5).
The consequence of this as we show in figure 3 is that for B 0 there are two phase transitions, the first corresponding to the 'impurity' n 4 causing a loss of long range order in the 3state-like ground state, and the second to a fully paramagnetic state with equal amounts of all states. An example of a zeroenergy point-defect impurity is shown in figure 7.
The B 0 case has only one transition: since the ground state only has long range order for the triangular arrangement of the n 4 's, there is no disorder transition in the other sublattice. The transition to a paramagnetic state arises when the n 4 sublattice disorders.
This demonstrates that the 4-state Potts model can be regarded as having a multicritical transition at T = 0, as three distinct phase boundaries terminate there.

Four-state Potts model with second nearest neighbours
Another simple case is to set J 2 = J 1 , B = 0, penalising second neighbours [35,54]. This is more physically interesting, since it has a zero-entropy ground state (non-degenerate aside from a trivial permutation of the states). This state (figure 4 and n 4 = 1/4 in figure 5) is equivalent to a layer of phase II hydrogen. The model also exhibits a phase transition to a disordered phase analogous to hydrogen phase I, the same transition occurring for any value of 0 < J 2 < J 1 . Varying the value of J 2 > 0 has no qualitative effect-the ground state is the same, and a single order-disorder transition occurs. For our finitesized system the order of the transition cannot be determined. In hydrogen, the symmetry change is expected to couple to a lattice distortion, making the change first-order. This is not the Table 1. Perfect crystal and conventional mean-field free properties at T = 0, related to the structures shown in figure 6. Enthalpy shown for molecules in state 4, and states 1-3 separately, P is pressure, B is defined in equation (2)  Phase diagram for J 1 = 1 = 3 5/2 J 2 . Data was collected for Monte Carlo simulations with for each (T, P) combination. Data was then aggregated into isotherms and isobars and the highest peaks in the heat capacity ( H 2 − H 2 ) or susceptibility (< n 2 4 > −< n 4 > 2 ) as a function of temperature and pressure were found. The (P, T) conditions of these peaks for each isotherm and isobar are shown here as black and red circles respectively. The isobars were then re-examined for other peaks, which are shown in green circles. These points trace out the phase boundaries of the structures which are labelled according to their ground state, and shown in figure 5. case in the Potts model and in fact, such a symmetry change has yet to be resolved in hydrogen.

Four state Potts model with pressure
The final model to consider is our best approximation to hydrogen. It has weak second-neighbour interactions (enough to lift degeneracy) and considers the effect of varying B. This can be related to a PT phase diagram because B, an energy penalty for molecules perpendicular to the layer, is then a proxy for pressure which reduces the c/a ratio.
The full three-dimensional parameter space of J 2 , J 1 , B is large, so we concentrate on the case relevant to interacting molecules. For H 2 the leading term with orientational dependence is the electrostatic quadrupole-quadrupole interaction, which drops off with distance as R −5 . The second neighbour distance on a triangular lattice is √ 3 so we choose J 2 /J 1 = 3 −5/2 . Figure 6 shows the phase diagram for this model as calculated with Monte Carlo, with boundaries detected from fluctuations in energy and/or n 4 . Only single spin-flips were allowed in the Monte Carlo.
Compared with the next-neighbour model, adding this additional term to the Hamiltonian has the surprising effect of reintroducing degenerate ground states. Figure 6 shows no fewer than seven separate low-T phases, depending on the value of B. For large positive B, the ground state is the same as the three-state model. B = 0 is identical to the four-state model, and this ground state is stable for small values of B. There are then a series of ground state structures with increasingly large n 4 , as shown in figure 5. Of these, only 1/4 and 2/3 are non-degenerate.
The n 4 = 2/3 state can be readily identified as equivalent to an expanded three-state model, with the interactions between RGB states being determined by J 2 . The difference in energy scale between B and J 2 mean that there are two temperaturedriven transitions, firstly to a disordering on the RGB sublattice, then to fully paramagnetic disorder. In the simulations, the Black circles denote state 4, RGB denote equivalent states 1-3, unit cell outlined in black, with a single localised exchange defect, ringed. The defect is localised and has zero excess energy which means the groundstate has non-zero entropy.
sublattice disorder is only picked up in the energy fluctuations, since n 4 remains constant on its sublattice.

Conclusions
We have examined the three and four state antiferromagnetic Potts models on a triangular lattice. Motivated by the orientation of diatomic molecules on such a lattice in hydrogen, we consider three equivalent states corresponding to inplane orientation, and a fourth inequivalent state representing out-of-plane orientation.
This system has a range of different ground states, often with large unit cells. We find phase transformations between these as a function of temperature and the relative enthalpy of the 4th spin state. Care must be taken to ensure that the boundary conditions are compatible with the ground state, otherwise spurious topological defects are introduced which broaden and obscure the phase transition in Monte Carlo simulation.
We find no fewer than 11 distinct states, as shown in figure 6. These include seven ordered ground states: ferromagnetic and the six shown in figure 5, three states with second-neighbour sublattice disorder and the paramagnetic state.
The analogy with the phase diagram for hydrogen comes from assigning phase 1 and the paramagnetic phase, which have the same symmetry, to the free rotor/quantum rotor phase I of hydrogen. The transformation from phase n 4 = 1 to paramagnet resembles the quantum-classical rotor, where state 4 maps to the symmetric J = 0 quantum rotor state, while states 1-3 map to the three J = 1 states.
The phase n 4 = 0, containing only three types of molecules in plane, maps to hydrogen phase III.
The remaining phases all lie in the region of PT space designated phase II in hydrogen. We do not claim that there is an exact match between the simple Potts model and the hinderedquantum-rotor phase [55][56][57][58][59], but we regard the complexity of the Potts diagram in this regions as an indication for why there is so much confusion about the structure of hydrogen phase II [18,25,30].
There is an interesting feature of the relationship between the excited states of molecular systems and the Potts model ordered ground states with single site flips or two-siteexchange localised defects such as figure 7. The excited states involve singly reoriented sites, and these local reorientations become thermally excited approaching the phase line. In a molecular material, excited states can be probed by spectroscopy, and there is an ongoing discrepancy in hydrogen phase II between the theoretical and measured Raman and IR signal [16,59]. DFT calculations of Raman spectra typically assume the excitations are quantum oscillators or rotors [5,12,59,60]. The Potts model analogy suggests that spectroscopic signals may also come from reorientation events. The idea has only very recently been investigated qualitatively [60].