Nuclear spin coupling crossover in dense molecular hydrogen

One of the most striking properties of molecular hydrogen is the coupling between molecular rotational properties and nuclear spin orientations, giving rise to the spin isomers ortho- and para-hydrogen. At high pressure, as intermolecular interactions increase significantly, the free rotation of H2 molecules is increasingly hindered, and consequently a modification of the coupling between molecular rotational properties and the nuclear spin system can be anticipated. To date, high-pressure experimental methods have not been able to observe nuclear spin states at pressures approaching 100 GPa (Meier, Annu. Rep. NMR Spectrosc. 94:1–74, 2017; Meier, Prog. Nucl. Magn. Reson. Spectrosc. 106–107:26–36, 2018) and consequently the effect of high pressure on the nuclear spin statistics could not be directly measured. Here, we present in-situ high-pressure nuclear magnetic resonance data on molecular hydrogen in its hexagonal phase I up to 123 GPa at room temperature. While our measurements confirm the presence of ortho-hydrogen at low pressures, above 70 GPa, we observe a crossover in the nuclear spin statistics from a spin-1 quadrupolar to a spin-1/2 dipolar system, evidencing the loss of spin isomer distinction. These observations represent a unique case of a nuclear spin crossover phenomenon in quantum solids.

C hanges in electronic spin-statistics under changing thermodynamic conditions are an established physical crossover phenomenon [1][2][3] . It has direct applications for spintronics 4 and enables the understanding of the stabilisation of magnetospheres of rocky (Earth-like) planets 5 as well as gas-and ice-giants 6 . The degrees of freedom of the nuclei spins, however, are widely approximated as fixed within the analysis of experiments, due to large shielding by core electrons and the extremely short atomic distances necessary to induce such fundamental changes.
Hydrogen, on the other hand, exhibits no core electrons and when bound contributes its electron to the molecular bond. Furthermore, due to the low mass of the hydrogen nuclei, quantum nuclear effects are considered to be significantly more pronounced compared to all other elements. The combination of both effects results in a number of fascinating physical phenomena in molecular H 2 7-9 .
One property intrigued physicists in particular: the nature of the nuclear spin of the H 2 molecule and the resulting coexistence of the spin isomers ortho-(ortho-H 2 ) and para-hydrogen (para-H 2 ). Following Pauli's exclusion principle, in order for the total H 2 molecular wave function to be antisymmetric under exchange of atomic positions, demands for the rotational ground state J = 0, that the corresponding total nuclear wave function is antisymmetric (singlet state of I = 0, i.e. para-H 2 ). Analogously, for the J = 1 rotational state, the total nuclear wave function is required to be symmetric (triplet state of I = 1, i.e. ortho-H 2 ). Therefore, the spin allotropic isomerism of the H 2 molecule originates in the coupling of both rotational state and nuclear spin. It has been argued 10,11 that at high pressure (P) ortho-and para-hydrogen spin isomers remain stable up to the dissociative Wigner-Huntington transition at P > 400 GPa 12,13 . This can only be assumed for weak or moderate intermolecular interactions, i.e. when nearest neighbour distances (r n ) are much shorter (≈0.7Å at ambient conditions) than next-nearest neighbour distances (r nn ≈3Å at ambient conditions), allowing for sufficient intramolecular wave function overlap (left side of Fig. 1a).
Under high enough densities, however, intermolecular interactions increase significantly as r nn decreases by~70% within 100 GPa 14,15 . At these pressures, r nn approaches r n and collective nuclear quantum fluctuations increase rapidly 16 .
For decades, theoretical 8,17,18 and experimental 9 studies indicated that under such extreme pressures odd values of J become unstable, rapidly decaying into even states, which leads to potential indistinguishability of the hydrogen spin isomers on experimental timescales.
The only experimental technique in high-pressure research to directly study the H 2 nuclear spin states is nuclear magnetic resonance (NMR) spectroscopy, detecting the linear response of the nuclear spin system upon radio-frequency excitation in a magnetic field B 0 . The spin-singlet state of para-hydrogen is NMR Non-Pairing regime ( I = ) 1   silent, whereas application of B 0 lifts the three-fold degeneracy of the ortho states and allows for an excitation of nuclear spin transitions (Fig. 1b). Nuclear spin-pairing in ortho-hydrogen leads, furthermore, to a finite electric quadrupole moment, eQ, interacting with the local charge distribution defined by the structural arrangement of hydrogen molecules. Thus, the quadrupolar coupling can be considered the dominant spin interaction, resulting in characteristic NMR line shapes 19 (Fig. 1c).
Here, we present 1 H-NMR data of dense molecular hydrogen up to 123 GPa at room temperature and find a distinct crossover in the nuclear spin statistics of molecular hydrogen indicating a loss of ortho-para spin isomer distinction. Details on experimental conditions, spectral simulations, as well as data analysis, are provided in the Methods section.

Results
Two NMR-DACs equipped with diamond anvils of 250 and 100 µm culets were loaded with molecular H 2 . At low pressure (below 60 GPa), intense 1 H resonances of roughly 500 kHz width were detected. With increasing P, the resonance signals broadened significantly approaching 750 kHz at 68 GPa (Fig. 2). Above 68 GPa, we observed a resonance narrowing accompanied by the emergence of two Pake doublets 20 with increasing splitting upon compression.
For the quadrupolar nature expected for ortho-H 2 (I = 1), the electric field V(r), defined by the local charge distribution based on the crystal structure of phase I, should influence the shape of the observed resonance lines. Calculated line shapes for an I = 1 spin system are shown in Fig. 2 at pressures of up to 68 GPa. The order of magnitude of the quadrupolar interaction energy was considered small relative to the nuclear Zeeman energy 19 and consequently treated as a first-order perturbation (see "Methods" section for computational details). Up to 68 GPa, the measured 1 H-NMR spectra are well described by calculated line shapes broadened by first-order quadrupole interaction. The line shape is mainly controlled by two parameters: (i) the quadrupole coupling constant C q describing the coupling between eQ and V(r) as well as (ii) the electric field gradient asymmetry parameter η accounting for the geometry of V(r). Figure 3a (top panel) shows estimated values of C q , which increase from 28.1(6) kHz at 20 GPa to 61.9(7) kHz at 58 GPa. This increase is likely originated in the high compressibility and rapidly reducing next-nearest neighbour distances between molecular H 2 units, enhancing quadrupolar coupling. The asymmetry parameter η (Fig. 3a, bottom panel) was found to be almost constant within experimental errors varying between 0.44 (6) at 20 GPa and 0.49(9) at 58 GPa. Based on the hexagonal crystal structure of phase I 14,21 , η can be expected to be close to 0.5, which is in excellent agreement with values derived from the analysis of our NMR measurements.
Above P ≈ 68 Pa, however, we observed a sudden decrease in both C q and η coinciding with a resonance peak splitting. As no structural rearrangement of H 2 molecules is reported by diffraction methods 21 or Raman spectroscopy 22 at room temperature in this P-range, effects based on a modification of the H 2 nuclear spin system should be considered.
This pressure indicates a turning point in the behaviour of the 1 H spin system, as the observed peak splitting devolves from having dominantly quadrupolar characteristics towards a system controlled by nuclear dipole-dipole coupling, resulting in pronounced I = 1/2 line shapes 20 with a frequency difference between spectral density function singularities directly correlated to the distances between hydrogen atoms. Considering that in this case both, the interaction with the nearest and next-nearest neighbours will result in a dipolar NMR pattern, respectively, a superposition of signals as shown in Fig. 1c can be expected.
Computationally obtained values of nearest and next-nearest neighbour distances 15 are r n ≈ 0.731Å and r nn ≈ 1.342 Å and at 120 GPa. Calculating the distances from the parameter-set obtained through analysis of the NMR spectra for such a mixed scenario resulted in r n ≈ 0.727 Å and r nn ≈ 1.27(8)Å at 123 GPa, in excellent agreement with the computational estimates by Labet et al. 15 . Figure 2 shows the comparison between experimental spectra and calculated I = 1/2 line shapes between 71 and 123 GPa. The respective values derived for both nearest (top panel) and next-nearest (middle panel) neighbour distances can be found in Fig. 3b. Additionally, the comparison between the equation of state derived from ab initio computations 23 and diffraction data 14,21 , along with the unit cell volumes (blue dots) derived from r n and r nn extracted from the analysis of the NMR spectra are shown in the bottom panel of Fig. 3b.
Homonuclear Lee-Goldburg decoupling sequences 24 have been used to suppress quadrupolar and dipolar line broadening in order to resolve isotropic chemical shifts, δ iso . Figure 3c shows the evolution of δ iso : initially decreasing from 8.6 to 5.9 ppm between 20 and 59 GPa, δ iso has an inflection point at~60 GPa and raises under further compression to 22.7 ppm at 123 GPa. Comparison with Raman data 25 suggests that the minimum in δ iso coincides Spin Pairing regime I=1 with the well-known turn-over in the Raman shift of the H 2 vibron caused by a weakening of intramolecular and increased intermolecular interactions 26 .

Discussion
The presented data analysis leads to the following interpretation of the observed effects: At P < 60 GPa, 1 H-NMR data is characteristic for an I = 1 quadrupolar spin system as expected for ortho-H 2 . In this regime, individual nuclear spin angular momenta couple with their nearest neighbours (separated on average by r n ), leading to a significant wave function overlap within the molecular units and a stabilisation of the spin isomers. The excellent agreement between NMR derived values for the electric field gradient asymmetry parameter η and values inferred from X-ray diffraction data 14,21 strengthens this assessment. Following the theoretical study of Strzhemechny et al. 27 , this compression driven enhancement of the quadrupolar coupling constant C q in this P-regime may be interpreted as experimental evidence for the mechanism of orthopara conversion through electric quadrupole interaction. At P > 60 GPa, quadrupolar coupling rapidly diminishes despite the absence of a structural rearrangement of the molecular H 2 units. Starting from about 70 GPa, spectral features characteristic of homonuclear dipole-dipole coupling between nearest and nextnearest neighbours become apparent. Provided the good agreement between NMR data with DFT 15 and experimentally 14,21 derived intermolecular and interatomic distances, this shift in behaviour implies intramolecular coupling of nuclear spins to become increasingly perturbed. The inflection point in the isotropic chemical shift δ iso strengthens this hypothesis as the increasing nuclear de-shielding above 60 GPa indicates a shift of electron density away from individual molecular centres towards intermolecular regions.    15 . Blue circles are extracted values of r n and r nn from the NMR spectra in the nonpairing I = 1/2 regime. Middle panel: nearest neighbour distances r n (blue circles) and DFT calculations 15 (red circles). The discontinuity at P > 125 GPa in the DFT calculations indicates a transition from the hcp based to a monoclinic structure. Bottom panel: comparison between the equations of state derived from ab initio computations 23 (yellow line) and diffraction data 14,21 along with the unit cell volumes (blue dots) derived from r nn and r n extracted from the NMR experiments. c Extracted isotropic chemical shift values δ iso after homonuclear Lee-Goldberg decoupling. Error bars are within the symbol size. The orange dotted line shows the room temperature Raman shift of the H 2 vibron 25 . d Selected Raman spectra of the H 2 vibron at increasing pressure. Spectra up to 64 GPa are taken from ref. 25 . Error bars in a and b were taken from spectral simulation and comparison with experimental data. Error bars in c are taken from Voigtian spectral line shape fitting to the experimental data. Error bars of the NMR data after Lee-Goldburg decoupling (blue points in c) were within the symbol size.
14 N-NMR on molecular nitrogen at P = 3 GPa (see "Methods" section) supports this argument, as the nitrogen spin system shows clear characteristics of a nuclear spin-triplet state anticipated within the non-pairing regime contrary to the quintuplet state stabilised in the molecular spin-pairing regime.
In this work, in-situ high-pressure nuclear magnetic resonance was used to investigate the nature of the nuclear spin statistics of molecular hydrogen up to 123 GPa in Phase I at room temperature. It was found that even at moderately high pressures (<100 GPa) intramolecular nuclear spin coupling broke down and the hydrogen spin system adopted an average dipolar I = 1/2 value. Crossovers of the nuclear spin statistics of a quantum solid such as hydrogen have so far not been observed and given the large compressibility of hydrogen in conjunction with strong nuclear quantum effects, this crossover phenomenon might only be experimentally observable in molecular H 2 . Nuclear spin statistics of similar diatomic molecules (e.g. N 2 ) are likely to be best described as non-pairing nuclear spins due to enhanced atomic masses as well as reduced compressibilities due to the presence of core electrons.
This nuclear spin-crossover may have far-reaching consequences for understanding different phenomena such as the stabilisation of magnetospheres of gas and ice giant planets containing large quantities of molecular H 2 .

Methods
Diamond anvil cell preparation. Two diamond anvil cells, equipped with pairs of diamond anvils with a culets size of 250 and 100 µm, were prepared. Rhenium gaskets were pre-indented to 25 and 10 µm, respectively, and 80 and 40 µm diameter holes were laser drilled in the centre of the indentation to form the sample cavities, resulting in sample volumes of about 125 and 13 pl, respectively.
The diamond anvils were coated with a 1-µm-thick layer of copper using physical vapour deposition 28 . Double 29 (in the case of the 250 µm diamonds) and triple 30 (for the 100 µm diamonds) stage Lenz-lens radio-frequency resonators were produced by using focused ion beam milling. To ensure electrical insulation and avoid hydrogen diffusion into the rhenium, the gaskets were coated by 500nm-thick layers of Al 2 O 3 . Radio-frequency excitation coils were made from 100 µm thick, Teflon insulated, copper wire and arranged such that a Helmholtz coil pair is formed.
Hydrogen loading was conducted at the ESRF and pressure was increased at cryogenic temperatures to avoid rapid hydrogen diffusion into the diamond anvils. The pressure was calibrated by means of the diamond edge Raman scale 31,32 . Comparison of the vibron frequencies of the H 2 samples at elevated pressures shows a slight systematic offset of less than 5 GPa at the highest pressures where Raman data was collected 33 .
NMR experiments. All NMR experiments were conducted using a solid-state NMR spectrometer from Tecmag Inc. (Redstone) equipped with a 100 W pulse amplifier. To polarise the nuclear spin system, we used a sweepable electromagnet with an average magnetic field of 1 T and sufficiently high homogeneity. Intense 1 H-NMR signals were recorded at frequencies of 45.26 MHz, corresponding to an external magnetic field strength of about 1063 mT. Using nutation experiments, we found optimal excitation pulses between 1 and 1.2 µs for both cells, in reasonable agreement with earlier experiments [28][29][30]34 .
Free induction decays were excited using a single pulse of 833 kHz to 1 MHz bandwidth. The spectrometer was blanked off for 1 µs to avoid damage to the preamplifier. Supplementary Figs. S2 and S3 show all 1 H-NMR spectra recorded by this method. 25,000 scans were accumulated for each spectrum (Fig. 2).
In order to resolve isotropic chemical shifts, δ iso , a Lee-Goldburg pulse for homonuclear decoupling was used 24 . The resulting narrowed NMR spectra had an FWHM linewidth of about 3 ppm, thus the resolution accuracy of δ iso was found to be in the order of 0.1 ppm (Supplementary Fig. S4). Resonance frequencies were referenced towards an aqueous solution of tetramethylsilane in a similar DAC at ambient pressure conditions.
Computation of NMR line shapes and asymmetry parameters of the electric field gradient. Calculation of the NMR line shapes was carried out following the analytical method outlined by Bloembergen and Rowland 35 , Pake 20 , and Hughes and Harris 36 : Using the standard expressions for the resonance frequency distribution ω for both first-order quadrupole interaction as well as homonuclear dipole-dipole interaction: with where the Euler angles α and β describe the orientation of the crystallites with respect to the external magnetic field. γ n is the gyromagnetic ratio of the hydrogen nuclei, m the nuclear spin quantum number (m = 1, 0, −1) and r i the average distance between interacting hydrogen nuclei, r n or r nn , respectively. η describes the asymmetry of the electric field gradient tensor (V ij ) in the principal axis system as: The line shape function, P(ω), for quadrupolar spin interactions, is given by: where β(ω, α, m) denotes the inverse function of Eq. (1) with respect to β, and µ accounts for the multiplicity of spectral functions. For the dipolar interaction P(ω), is given by where β(ω i ,α) denotes the inverse function of Eq. (2) with respect to β and µ accounts for the multiplicity of spectral functions. Cut-off frequencies of the resulting spectral line functions were chosen according to the Heaberlein convention for NMR shift tensors 37 . Spectral line broadening was accounted for by convolution of the total line shape function with a Voigtian line of defined Lorentzian and Gaussian widths. In order to fit the experimental data, the respective line shape function P(ω) is optimised by varying C q and η for quadrupolar coupling and r n and r nn for dipolar coupling. The corresponding Python scripts are available from the authors upon request. Table 1 summarises all fit parameters.
In order to calculate the asymmetry parameter η of the electric field gradient tensor in the spin-pairing regime, we used the second derivative of the electric potential, V(r), defined by the molecular centre of gravity positions from diffraction measurements 14 : Using Eq. (6) under consideration of the ordering of the components of V ij in the principal axis system allows the computation of η from crystallographic data. 14 N-NMR of molecular nitrogen at 3 GPa. Molecular nitrogen was measured using natural isotopic composition, where the majority of molecules can be expected to be pairs of 14 N nuclei. As 14 N nuclei have a nuclear spin of I = 1, one can expect a spin-pairing scenario similar to molecular D 2 : the para-N 2 states consist of a quintuplet subsystem with I = 2, whereas the ortho-N 2 states are a triplet subsystem.
The electric field gradient asymmetry parameter η was estimated according to diffraction data 14 to be around 0.23. Recorded 14 N-NMR spectra (Supplementary Fig. S2; right panel) do not show pronounced shoulder, expected for an I = 2 quadrupolar powder pattern in absence of m −2→−1 and m 1→2 transitions. In fact, the spin system is well described by an I = 1 spin system using the estimated value for η (Supplementary Fig. S2; left panel).
According to structural data 38 , r n can be estimated to be around 1.2 Å at this pressure; four times longer than the thermal de-Broglie wavelength of a single 14 N atom. Therefore, the wave function overlap should be negligible in molecular nitrogen at these pressures and nuclear spins considered unpaired.

Data availability
The data supporting the findings of this study are publicly available from the corresponding author upon request.

Code availability
The custom Python code used for simulating the NMR spectra is available from the corresponding author upon request.
Received: 8 August 2020; Accepted: 6 November 2020; C q is the quadrupole coupling constant, η the asymmetry parameter of the electric field gradient tensor in the principle axis system, r n and r nn are the nearest and second nearest neighbour distances, respectively. The isotropic chemical shift, δ iso , was derived after homonuclear Lee-Goldburg decoupling.