# Dual quantum locking: Dynamic coupling of hydrogen and water sublattices in hydrogen filled ice

Loan Renaud, Tomasz Poreba, Simone Di Cataldo, Alasdair Nicholls, Léon Andriambariijaona, Maria Rescigno, Richard Gaal, Michele Casula, A. Marco Saitta, Livia Eleonora Bove

► To cite this version:

Loan Renaud, Tomasz Poreba, Simone Di Cataldo, Alasdair Nicholls, Léon Andriambariijaona, et al.. Dual quantum locking: Dynamic coupling of hydrogen and water sublattices in hydrogen filled ice. 2025. hal-05405284

HAL Id: hal-05405284

https://hal.sorbonne-universite.fr/hal-05405284v1

Preprint submitted on 8 Dec 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Dual quantum locking: Dynamic
coupling of hydrogen and water
sublattices in hydrogen filled ice

Loan Renaud<sup>a,1,2</sup>

Tomasz Poreba<sup>b,1</sup>

Simone Di Cataldo<sup>c</sup>

Alasdair Nicholls<sup>b</sup>

Léon Andriambariarijaona<sup>d</sup>

Maria Rescigno<sup>b,c,e</sup>

Richard Gaal<sup>b</sup>

Michele Casula<sup>a</sup>

A. Marco Saitta<sup>a</sup>

Livia Eleonora Bove<sup>a,b,c,2</sup>

October 2025

<footnote><sup>a</sup>Institut de Minéralogie, de Physique des Matériaux et de Cosmochimie (IMPMC), Sorbonne Université, CNRS UMR 7590, MNHN, 4, place Jussieu, Paris, France</footnote>

<footnote><sup>b</sup>Laboratory of Quantum Magnetism, Institute of Physics, École Polytechnique Fédéerale de Lausanne, CH-1015 Lausanne, Switzerland</footnote>

<footnote><sup>c</sup>Dipartimento di Fisica, Sapienza Università di Roma, Piazzale Aldo Moro 5, 00185 Roma, Italy</footnote>

<footnote><sup>d</sup>Laboratoire pour l'Utilisation des Lasers Intenses (LULI), École Polytechnique, 91120 Palaiseau, France</footnote>

<footnote><sup>e</sup>Institut Laue-Langevin, 71 Avenue des Martyrs, Cedex 9, Grenoble, France</footnote>

<footnote n="1"><sup>1</sup>L.R. and T.P. contributed equally to this work.</footnote>

<footnote n="2"><sup>2</sup>Corresponding authors. E-mail: loan.renaud@sorbonne-universite.fr, livia.bove@sorbonne-universite.fr</footnote>
Abstract

Hydrogen hydrates (HH) are a unique class of materials composed of hydrogen molecules confined within crystalline water frameworks. Among their multiple phases, the filled ice structures, particularly the cubic C2 phase, exhibit exceptionally strong host-guest interactions due to ultra-short H₂-H₂O distances and a 1:1 stoichiometry leading to two interpenetrated identical diamond-like sublattices, one comprised of water molecules, the other of hydrogen molecules. At high pressures, nuclear quantum effects involving both hydrogen molecules and the water lattice become dominant, giving rise to a dual-lattice quantum system. In this work, we explore the sequence of pressure- and temperature-driven phase transitions in HH, focusing on the interplay between molecular rotation, orientational ordering, lattice symmetry breaking and hydrogen bond symmetrization. Using a combination of computational modeling based on classical and path-integral molecular dynamics, quantum embedding, and high pressure experiments, including Raman spectroscopy and synchrotron X-ray diffraction at low temperatures and high pressures, we identify signatures of quantum-induced ordering and structural transformations in the C2 phase. Our findings reveal that orientational ordering in HH occurs at much lower pressures than in solid hydrogen, by inducing structural changes in the water network and enhancing the coupling of water and hydrogen dynamics. This work provides new insights into the quantum behavior of hydrogen under extreme mechanochemical confinement and establishes hydrogen-filled ices as a promising platform for the design of hydrogen-rich quantum materials.
The rotational dynamics of guest gas molecules confined in porous materials are profoundly influenced by their local environment [1]. In large, nearly spherical cages, guest species often behave as nearly free rotors [2–4]. In contrast, low-symmetry or anisotropic environments impose steric constraints that hinder rotation, leading to librational motions or even complete orientational freezing [5–7].

Such effects have been extensively studied in clathrate hydrates, where the host–guest potential not only quantizes the center-of-mass motion but can also induce partial orientational ordering at low temperatures and high pressures [5, 8–10]. Zeolite frameworks, with their inherent anisotropy and stronger host–guest interactions, exhibit similar behavior: small guest molecules often display restricted rotational freedom and develop orientational preferences within the channels or cages [11–15]. Endohedral fullerenes offer a particularly clean example of quantum confinement: molecules such as $H_2$ and $CH_4$ confined inside $C_{60}$ exhibit discrete quantum rotational states that are exquisitely sensitive to cage symmetry and weak van der Waals interactions [6, 16, 17]. These systems provide compelling evidence of how confinement and symmetry shape quantum rotational behavior across a wide range of host architectures.

In gas-filled ice structures under high pressure [18], the extreme confinement imposed by the ice hosting lattice leads to the suppression of molecular rotations even at high temperatures, enabling the onset of orientational ordering of guest molecules [19–21]. In certain cases, this ordering is accompanied by molecular distortion [20], and can even trigger symmetry-lowering transitions in the host lattice itself [22, 23]. These phenomena underscore the unique interplay between quantum confinement and structural flexibility in filled-ice systems—an interplay rarely observed in classical porous materials.

Among these, hydrogen hydrates (HH)—where hydrogen molecules are encapsulated within an extended water lattice—stand out as a particularly rich system for exploring host–guest interactions under extreme thermodynamic conditions. HH comprises five well-characterized phases. The first, i.e., the clathrate structure sII, contains hydrogen molecules trapped in large polyhedral cages within a hydrogen-bonded water network [3, 24]. The other four—C0, C1, C2, and C3—are filled-ice structures in which hydrogen molecules occupy interstitial voids within compact ice-like frameworks. Specifically, C0 adopts the open-channel ice XVII topology [25], C1 corresponds to the layered ice II structure [26], while C2 and C3 are based on the cubic ice (Ic) framework [27].

While sII is only stable at cryogenic temperatures, both C1 and C2 persist at room temperature under high pressure. C1 is stable between 0.9 and 2.7 GPa (and metastable up to 5.2 GPa), whereas C2 stabilizes above 2.7 GPa and transforms into C3 at temperatures above 670 K and pressures above 30 GPa. C3 remains stable up to at least 90 GPa and can be recovered metastably down to 7 GPa [23].

Among these phases, C2 is particularly intriguing. It adopts a topology analogous
to that of ice VII, consisting of two interpenetrating but non-interconnected diamond-
like sublattices. However, unlike ice VII—where both sublattices are composed of water
molecules—in C2, one sublattice is entirely formed by hydrogen molecules. This breaks
the structural equivalence, lowers the overall symmetry, and results in a unique dual-
framework architecture. Remarkably, C2 exhibits the shortest known H<sub>2</sub>-H<sub>2</sub>O distances
among hydrogen hydrates, and its 1:1 stoichiometry gives rise to exceptionally strong
host–guest interactions.

Under such extreme densities confinement effects predominate, affecting not only
the hydrogen molecules within their sublattice but also the protons constituting the ice
framework. As pressure increases, the distinction between host and guest progressively
diminishes, and the system is more appropriately described as two equivalent, interpene-
trating quantum lattices of distinct species: H<sub>2</sub> and H<sub>2</sub>O.

The rotational dynamics of hydrogen in hydrogen-filled ices are strongly pressure- and
temperature-dependent. Increasing pressure and lowering temperature progressively
restrict guest rotation, inducing stepwise orientational ordering [28], analogous to effects
seen in solid methane [29, 30] and solid hydrogen [31–33], where intermolecular inter-
actions produce roton band splittings. Similar but stronger effects occur in hydrogen
clathrates [10].

In the C2 phase, first-principles calculations predicted a pressure-induced cubic-to-
tetragonal distortion at low temperature via hydrogen alignment along the c-axis [34],
confirmed in D<sub>2</sub>O-D<sub>2</sub> by XRD peak splitting and roton band splitting above 20 GPa at 10 K
[35]. Additional high-pressure phases (HH-HP1 at 40 GPa and HH-HP2 at 60 GPa) involve
further orientational ordering and possible hydrogen-bond symmetrization [35–37].

Our recent work refined this phase diagram: we observed tetragonal distortion of C2
above 30 GPa at room temperature, amorphization above 60 GPa, and transformation to
C3 upon laser heating [23, 27]. We also determined that hydrogen-bond symmetrization
occurs at 26 GPa—much lower than in ice VII—due to shorter O–O distances in the cubic
ice Ic framework [38].

In the present work, we go beyond these macroscopic observations to elucidate their
microscopic origin: we demonstrate that hydrogen-bond symmetrization in the water
lattice triggers the nematic alignment of the H<sub>2</sub> sublattice, and that this collective orienta-
tional ordering, in turn, drives the tetragonal distortion of the host framework. By disen-
tangling this sequence of coupled events, we establish the microscopic mechanism linking
proton symmetrization, guest orientational order, and host lattice deformation in C2. We
probe these phase transitions via combined *ab initio* simulations, low-temperature Ra-
man spectroscopy, and synchrotron XRD [39–43], highlighting how enhanced host–guest
coupling under pressure suppresses rotational freedom and stabilizes orientational order
at much lower pressures than in pure hydrogen.
<figure><img src="image_2.png" /><figcaption>FIGURE 1. Pressure-Temperature phase diagram of the C2 hydrogen hydrate as derived by our quantum embedded calculations and experimental data. Some experimental pattern are indicated in Figure 1: XRD data are shown as circles, colored red or black depending on whether they belong to the temperature- or pressure-ramp series; Raman measurements are indicated by squares. The summary of experimental thermodynamic path is reported in Table 1. Symbols are filled in the cubic phase and open in the tetragonal phase, while half-filled symbols indicate the transition between them. The background color represents S, the orientation factor (defined in Methods) of H₂ computed in the quantum embedded framework. The white triangles with their associated error bars are derived from our molecular dynamics simulations, utilizing the orientation factor as detailed in the Methods section. The dashed lines serve as a guide to the eye to locate the cubic to tetragonal transition. The three insets display the optimized geometries obtained via DFT (0 K) for the quantum plastic phase (upper left), herringbone phase (lower left), and nematic phase (right) of H₂ orientational states.</figcaption></figure>
<table><caption>TABLE 1. Summary of experimental thermodynamic paths for X-Ray Diffraction (XRD) and Raman.</caption><thead><tr><th></th><th>XRD</th><th>Raman</th></tr></thead><tbody><tr><td><strong>Isobar</strong></td><td>2.45(3) GPa<br>300 → 30 K</td><td>3.5 GPa<br>300 → 100 K<br>14.7 GPa<br>300 → 5 K<br>35 GPa<br>300 → 20 K</td></tr><tr><td><strong>Isotherm</strong></td><td>300 K<br>2 → 46 GPa</td><td>300 K<br>10 → 40 GPa</td></tr></tbody></table>

## Results and Discussion

Figure 1 presents the pressure-temperature phase diagram constructed from our combined computational and experimental approaches, with the thermodynamic paths probed experimentally summarized in Table 1. Within the C2 phase, H₂ exhibits three distinct regimes, dictated by temperature and the degree of coupling with the surrounding water lattice.

At high temperatures and pressures below 30 GPa (above the dashed line) the H₂ sublattice in H₂ hydrate realizes a quantum plastic crystal. In this state, the H₂ molecules are ordered on a diamond-type lattice, defining well-localized translational equilibrium sites, yet their orientations remain dynamically disordered [44]. Each H₂ behaves as a nearly free quantum rotor, i.e., the rotational degrees of freedom are not frozen but remain active down to relatively low temperatures due to quantum fluctuations. This combination of translational order and orientational quantum mobility is the defining feature of a quantum plastic phase. In contrast, the water sublattice is proton-disordered but not plastic, since the orientations of H₂O molecules are constrained by the ice rules and do not undergo free rotational dynamics. The overall symmetry of the crystal is cubic, consistent with previous DFT predictions and experimental observations [23, 27, 34].

Upon crossing the dashed boundary, either through temperature reduction or compression, H₂ molecules undergo orientational ordering. At high pressure this results in a phase where H₂ are aligned along the c axis of the crystal. This ordering transition can be quantified by the orientation factor *S* (See Methods) computed based on our quantum embedded calculations as described in the following, and visualized as a color gradient in Fig. 1. The factor *S*, which ranges from 0 (complete orientational disorder) to 1 (perfect alignment), measures the collective alignment of the hydrogen molecules along the c axis of the crystal, capturing both the orientational order imposed by the lattice geometry and
the intrinsic thermal quantum fluctuations encoded in the density matrix.

Within the orientationally ordered phase—occurring below 30 GPa and 100 K–H₂ molecules adopt a herringbone-like arrangement, as shown in the inset of Fig. 1. This configuration favors alignment along the crystalline c axis. The tilt angle of the molecules relative to the c axis varies continuously with pressure, reflecting the energetically favorable response to lattice compression at low temperatures. As pressure increases, the tilt angle gradually decreases, approaching full alignment along the c axis, i.e the nematic phase. This progressive reorientation is also captured by the pressure dependence of the orientation factor *S* in Fig. 1.

The low-temperature, high-pressure phases (< 50 GPa) of H₂ and D₂ are unique molecular systems in which free rotations are governed by both quantum effects and temperature [45]. Molecular hydrogen occurs in two nuclear-spin isomers, both partially occupied at ambient conditions. When the proton spins couple to the singlet state with total nuclear-spin quantum number *J* = 0, defined by *J* = *J*₁ + *J*₂, where *J*₁ = *J*₂ = ½ are the proton resolved nuclear spins, the resulting para-hydrogen possesses an antisymmetric spin wavefunction and, to satisfy overall fermionic antisymmetry, its spatial part is restricted to even rotational quantum numbers *l*. Coupling the same proton spins to the triplet state (*J* = 1) produces ortho-hydrogen, whose symmetric spin wavefunction permits only odd *l* values.

At low pressures, Raman spectra of solid hydrogen and clathrates [4, 10, 27, 46, 47] display free-rotor peaks with energies *E*<sub>*l*</sub> = *ħ*²/(2*I*)*l*(*l* + 1) where *I* is the moment of inertia of the rotor and *l* the angular momentum, with selection rule *Δ*(*l* = ±2. In the free molecule, rotational levels are (2*l* + 1)-fold degenerate, and the main gas-phase Raman bands are *S*₀(0) (*l* = 0 → 2) (354 cm⁻¹) and *S*₀(1) (*l* = 1 → 3) (587 cm⁻¹) [46, 48]. In the clathrate phase and in low-pressure filled ice C0 and C1 phases, the rotor bands remain close to the free-rotor values and their degeneracy is preserved. In the high-pressure C2 phase, however, the anisotropic crystal field lifts this degeneracy and shifts the main ortho-para transition band by 4 meV—about 30% of its value [49]. With increasing pressure, these bands broaden, making individual rotational levels increasingly difficult to resolve.

Figure 2a,b shows the evolution of the Raman intensity distribution among the *S*<sub>*x*</sub>(0) transitions in a gas-loaded hydrogen hydrate sample (see Methods) as a function of temperature and pressure. XRD confirmed phase-pure C2 after gas loading, whereas cryo-loaded samples contained 83 mol% H<sub>2</sub>O due to decomposition of the sII substrate (6:1 H<sub>2</sub>O:H<sub>2</sub> stoichiometry).

At low temperatures, the *S*₀(2) and *S*₀(3) bands weaken and vanish below 150 K due to thermal depopulation of higher rotational levels; in pure hydrogen, these disappear at 120 K (3 GPa) [50]. In the compressed hydrate, the Raman evolution is more complex: a sudden broadening and asymmetry around 30 GPa suggest a change in the rotational
<figure><img src="image_3.png" /><figcaption>FIGURE 2. Raman spectra as a function of temperature at P= 14.7(1) GPa (a) and as a function of pressure at T=300 K (b). Represnetative X-ray diffraction data as a function of temperature at P=2.45(3) GPa (left panel,c-e) and as a function of pressure at T=300 K (right panel, f-h). The extracted evolution of the lattice parameters (in the tetragonal setting) and the degree of asymmetrization upon cooling (i,j) and compression (k,l), are plotted in the bottom panels, respectively.</figcaption></figure>

state of encapsulated H₂ [46, 51], likely driven by modifications to the crystalline field
from the surrounding water framework during its transition to a less compressible phase
[23]. The instability of C2 above 20 GPa toward decomposition into ice VII and molecular hydrogen [37] may also contribute to inhomogeneous Raman signals from free H₂. The intensity oscillations reported near 20 GPa [37] were not observed here; spectra were consistently collected from the same spot after each pressure step, with a 30-minute equilibration, ensuring that the observed increases in the S₀(1)/S₀(0) intensity ratio upon cooling (14.7 GPa) and compression (300 K) reflect enhanced ortho-para conversion. Both high pressure and low temperature are known to facilitate this process in solid hydrogen [46, 51].

XRD measurements reveal that either cooling below 150 K at 2.45 GPa or compressing above 27 GPa at RT induces peak splitting (Fig. 2c-h), signaling a cubic ($Fd\bar{3}m$) to tetragonal (I4₁/amd) transition. Although the splitting pattern is similar in both cases, it is markedly sharper and an order of magnitude larger under compression. The pressure-driven transition begins abruptly at around 30 GPa, coinciding with the onset of Raman anomalies—evidence of coupling between H₂ rotational dynamics and framework distortion. In both scenarios (cooling and pressurization), the c axis stiffens while the orthogonal directions contract (Fig. 2i,k). The distortion, quantified as $c/\sqrt{2}a$, increases smoothly, reaching 0.5% at 30 K (2.45 GPa) and 4.6% at 48 GPa (300 K) (Fig. 2j,l).

To probe the molecular-scale mechanism linking proton alignment to the lattice distortions seen in XRD, we performed simulations of a crystal-embedded H₂ molecule by solving the Schrödinger equation with an effective potential corresponding to the crystal field generated by the water framework and surrounding H₂ molecules. Figure 3(a-b) shows the evolution of H₂ density isosurfaces with temperature and pressure, along with the angular component $P_{r_{eq}}(\theta, \phi)$ of the external potential $V_{ext}(r, \theta, \phi) = R(r) + P_r(\theta, \phi)$ at the equilibrium H₂ bond length $r = r_{eq}$. $\theta$ and $\phi$ are respectively the polar and azimuthal angles of the molecular orientation in a cubic lattice reference. As described in Methods and Supporting Information, densities at a given $T$ are obtained from the para and ortho thermal density matrices mixed together with a fixed 1:3 para-ortho ratio.

At low pressure ($P = 3$ GPa, panel 3(a)), the external angular potential exhibits minima that force the molecules into a herringbone-like arrangement. As $T$ increases within the C2 stability range, para- and ortho-wavefunctions mix and the density approaches spherical symmetry. With increasing pressure, intermolecular interactions intensify and potential barriers rise, while the minima shift toward $\theta = 0$ [$\pi$], ultimately favoring full nematicity. At room temperature (panel 3(b)), the density evolves from nearly spherical shells at 10 GPa—characteristic of the quantum-plastic phase—to well-defined lobes aligned with the c-axis at 50 GPa.

The eigenvalue evolution shown in Fig.3(c) demonstrates how the anisotropic crystal field in the C2 phase lifts the free-rotor degeneracy and induces a pressure-dependent giant splitting, in agreement with earlier findings [49]. This quantum mechanical picture
<figure><img src="image_4.png" /><figcaption>FIGURE 3. Isosurfaces of the density ρ(⃗r) and angular potential P<sub>req</sub>(θ, φ) evaluated at the equilibrium H<sub>2</sub> bond length, shown as a function of temperature at P = 3 GPa (a), and as a function of pressure at T = 300 K (b). The density is projected onto the three planes (xy, yz, and xz) and normalized with respect to the maximum density at 0 K, 3 GPa for panel (a), and at 300 K, 50 GPa for panel (b). Pressure dependence of the first nine eigenvalues at T = 0 K. Ortho and para levels are shown as circles and squares, respectively (c). Boltzmann weights vs. temperature at P = 50 GPa for the first five eigenvalues (d). Isosurfaces of the wavefunctions squared modulus for the ground state (para, l = 0) and the first three excited states (ortho, l = 1) of the H<sub>2</sub> molecule at P = 50 GPa (e).</figcaption></figure>

explains the emergence of nematic ordering at high pressures and its persistence to
elevated temperatures. At low pressure, within the C2 stability range, the small energy
gap between the fundamental para state (l = 0) and the three ortho states (l = 1, λ₁, λ₂, λ₃)
enables thermal mixing above ~100 K, yielding nearly spherical densities. Increasing
pressure drives a sharp rise in the splitting between λ₁ and λ₂, bringing the λ₁ level close
<figure><img src="image_5.png" /></figure>

FIGURE 4. (a) Comparison of the experimental Raman points and all simulated S₀(0) and S₀(1) rotons. The blue and orange lines are theoretical whilst the red and green point indicate the S₀(0) and S₀(1) rotons respectively. The rotons were modelled with three components each, as indicated by the square, crosses and diamonds (b) Pressure dependence of the S₀(1) frequency ratio ωH₂/ωD₂ derived from the experimental Raman data. The rigid rotor ratio is indicated by the dotted line at ωH₂/ωD₂ = 2, whilst the quantum harmonic oscillator ratio is indicated by the dashed line at ωH₂/ωD₂ = √2. HH crosses into the harmonic oscillator regime at around 27 GPa, the same pressure a structural change is previously reported [23] (c) Experimental Raman frequencies derived from the fitting procedure described in Methods (squares) and computed Raman frequencies obtained via DFPT (circles) for the hydrogen vibron upon compression of the hydrogen hydrate. The change in slope, highlighted by the shadow area, corresponds to the pressure where the hydrogen molecules order.

to *l* = 0. This restricts significant mixing to *l* = 0 and *l* = 1, λ₁, both oriented along the *c* axis, producing a nematic density distribution. The fixed 1:3 para-ortho ratio further enhances the statistical weight of λ₁, stabilizing orientational order up to ~500 K (Fig. 3d-e).

From these computed eigenvalues, we can determine the Stokes rotational transitions S₀(0) and S₀(1), plotted in Fig. 4(a) as orange and blue curves, respectively. The comparison with the experimental Raman data in the same panel shows good agreement, particularly in the pressure evolution of the main features, thus validating the simulation framework. The pronounced change in slope of both S₀(0) and S₀(1) near 30 GPa coincides with the cubic-tetragonal distortion observed in XRD and with the calculated onset of dominant λ₁ - *l* = 0 mixing (Fig. 3c-d), supporting the interpretation that the structural transition is driven by a pressure-induced orientational ordering of the guest molecules.

Panel 4(b) further supports this picture by showing the pressure dependence of the ratio between the S₀(1) Raman frequencies measured for an H₂-hydrate (ωH₂) and for a D₂-hydrate (ωD₂). At low pressure, the ratio remains close to the free-rotor limit of 2, reflecting the 1/I scaling of rotational energies with the molecular moment of inertia. Upon compression, the increasingly strong crystalline field exerted by the rigid water lattice progressively hinders the molecular rotations, effectively locking the H₂/D₂ orienta-
<figure><img src="image_6.png" /><figcaption>FIGURE 5. (a) Crystal structure of C2 hydrogen hydrate at low pressures ($P \lesssim 30$ GPa), showing the two interpenetrating cubic sublattices of $H_2O$ and $H_2$. The water-proton network is disordered, while the $H_2$ molecules are orientationally disordered at high temperature and adopt a herringbone-like arrangement at low temperature. (b) Pressure dependence of the first (solid lines) and second (dotted lines) peak positions in the radial distribution functions $g(r)$ for the $H_w$~O and O~O pairs, from PIMD simulations at $T = 300$ K. The background color gradient represents the orientation factor $S_{PIMD}$ extracted from those simulations. (c) Crystal structure of C2 hydrogen hydrate in the symmetrized (water-proton) state at high pressures ($P \gtrsim 30$ GPa). At the highest pressures, $H_2$ molecules become nematic along the c-axis, coinciding with the tetragonal distortion direction.</figcaption></figure>

tion relative to the host cages. As rotational freedom is quenched, the excitations acquire librational character and the ratio approaches the harmonic mass-scaling limit $\sqrt{2}$, characteristic of an oscillator with frequency $\propto m^{-1/2}$. A qualitatively similar evolution has been reported for pure solid hydrogen and deuterium at much higher pressures (above ~140 GPa) [46, 48], and a direct comparison between the dispersion curves for the hydrate and for the elemental solids is provided in the Supplementary Information (Fig. A6). In analogy with the dense elemental phases, if the hindering potential in the hydrate is sufficiently anharmonic — i.e. weaker than purely harmonic at large angular displacements — the ratio can transiently overshoot below $\sqrt{2}$ before converging, reflecting an enhanced softening of the librational mode. Notably, the roton-libron crossover in the hydrate occurs at pressures an order of magnitude lower than in pure hydrogen/deuterium, highlighting the much stronger orientational constraints imposed by the host lattice. The observed trend thus constitutes direct spectroscopic evidence for the progressive locking of the guest-rotor orientation due to the strengthening of $H_2/D_2-H_2O$ coupling under pressure.

This progressive hindering of the guest rotations revealed in panel 4(b) has a direct counterpart in the behavior of the intramolecular H-H stretching mode, since changes in rotational freedom modify the coupling between rotational and vibrational degrees of freedom and thus the vibron pressure dependence. Consequently, from the analysis of the H-H stretching mode in Fig. 4(c) comparing Density Functional Perturbation Theory
(DFPT) calculations with experimental Raman frequencies it turns out that both simulation and experiment display a clear change in slope near 30 GPa. This is consistent with the signatures seen in the rotational modes and further linking the vibron behavior to the onset of pressure-induced orientational ordering predicted by our quantum embedded model.

The structure of the C2 hydrogen hydrate at low pressure is shown in Fig.5(a). The two interpenetrating hydrogen and water cubic sublattices are fully orientationally disordered: the protons in the water framework are randomly distributed between equivalent hydrogen-bond configurations, and the H₂ molecules adopt random orientations within their cages. Panel 5(b) reports the evolution with pressure of the positions of the two oxygen atoms closest to a given water proton Hw—the intramolecular oxygen covalently bound to Hw and the oxygen of the neighboring molecule to which Hw is hydrogen-bonded—shown in blue, as extracted from our Path-Integral Molecular Dynamics (PIMD) trajectories (see Methods). With increasing pressure, these two nearest oxygen atoms approach to each other and their corresponding distances to Hw converge at ~26 GPa [38], indicating proton symmetrization and the formation of symmetric O–H···O bonds. The first peak position of the oxygen-oxygen (O–O) radial distribution function is shown in black, while the background color encodes the orientation factor S<sub>PIMD</sub> of the H₂ molecules, computed from path-integral molecular dynamics as an ensemble average over beads and time steps, as defined in Methods. This quantity implicitly incorporates both the geometrical arrangement of the molecules and their thermal-quantum fluctuations relative to the c axis. Before symmetrization, compression primarily shortens the hydrogen bonds, with little ordering of the H₂ sublattice. Proton symmetrization abruptly stiffens the O–H···O network, producing a sharp increase in the bulk modulus B [23, 27] and making the ice sublattice much less compressible. Beyond this point, further isotropic compression is energetically disfavored; instead, it becomes more favorable to reorient the H₂ molecules with their bond axes parallel to the c direction [panel 5(c)], reducing the volume along a and b owing to the increased compressibility along these axes. This collective alignment of the H₂ sublattice drives the tetragonal elongation of the H₂O framework and produces a second, more gradual rise in B [23, 27]. In other words, the lattice first complies with increased pressure through hydrogen-bond rigidification in the water network, and then through anisotropic reorientation of the hydrogen network, rather than by further reducing the O–O separation. Figure 5 thus reveals a two-stage process in which the quantum-mechanical behavior of the water and hydrogen sublattices becomes tightly interlocked, their structural and dynamical degrees of freedom evolving in unison and merging into a single quantum lattice.

Our results demonstrate that C2 hydrogen hydrate undergoes a pressure-driven, two-stage transformation in which the hindering of the rotational dynamics of the hydro-
gen sublattice becomes intertwined with the distortion of the ice sublattice. The first stage—proton symmetrization in the host lattice—abruptly stiffens the hydrogen-bond network, while the second—anisotropic alignment of the H₂ molecules—relieves further compression and induces a tetragonal distortion. This collective alignment is directly manifested in the structural transition of the lattice, marking the point where the guest reorientation and host deformation proceed in lockstep. The coupling is evident across structural, elastic, rotational, and vibrational observables, linking the roton-libron crossover of the guest molecules to the proton ordering of the host. In this regime, the motions of the host lattice and the guest molecules become strongly interconnected, forming collective quantum modes rather than distinct subsystems. This provides an exquisite example of coupled nuclear quantum dynamics in a crystalline solid, with potential implications for other hydrogen-rich materials under high pressure.

## Methods

### Experimental Details

Samples were prepared in two different ways, depending on the specific experiment, cryo-loaded or gas-loaded in the diamond anvil cell (DAC). In the first case, polycrystalline hydrogen hydrate samples with the sII clathrate structure were prepared by exposing H₂O ice (Ih) to gas H₂ at 0.28 GPa for 20–30 min, following the method described in ref. [3]. The ice was made of spheres with typical diameters of a few tens of micrometers. After preparation, the samples were recovered at ambient pressure and stored at liquid nitrogen temperature. Then, a small amount of sample was loaded at liquid nitrogen temperature into a diamond anvil cell partially immersed in a liquid nitrogen bath and in an atmosphere saturated with nitrogen vapor to avoid water or oxygen condensation on the culet. The loaded samples were then compressed to pressures of a few GPa before being warmed to room temperature. This is a similar loading procedure to what we had previously done for methane hydrate [20, 22]. This procedure inevitably produces a mixture of C₂ hydrogen hydrate and excess pure ice, due to the difference in molar ratios between the clathrate sII phase and the C₂ phase. Clathrate sII samples prepared following our method have a molar ratio H₂O:H₂ between 4 and 5, as we have verified by neutron diffraction in the past [3], and the C₂ phase is characterized by a molar ratio of 1. Two samples were prepared directly in the DAC by loading room-temperature liquid H₂O and high-pressure H₂ using a gas loading setup. Briefly, liquid water was loaded together with an air bubble, consisting of approximately 70% of the sample chamber volume. The air bubble was replaced with 99.99+% H₂ at 1400 bar, and the DAC was subsequently pneumatically closed and pressurized. These samples were measured by both Raman spectroscopy and by XRD. Some of the cryoloaded samples contained a small amount of
pure N₂ and/or pure O₂, both of which can easily be identified from their respective vibron Raman peaks. They were trapped in the sample chamber during cryoloading and would have acted as a pressure-transmitting medium. In the gas-loaded samples, H₂ itself was acting as the pressure transmitting medium. Culet diameters ranging between 200–500 µm for the IIas quality diamond anvils were used. Re-foil gaskets were used to contain the sample whilst ruby and/or gold spheres were used as pressure calibrants in Raman and X-ray experiments respectively.

Angle-dispersive X-ray powder diffraction patterns were acquired at ID15b at ESRF (Grenoble, France) using a monochromatic X-ray beam (λ=0.41 Å) and an Eiger 2 9M CdTe flat-panel detector, with a typical acquisition time of a ten seconds.[43] The beam spot size was 6×6 µm, and the DAC was rocked 5° during the acquisition. The 2D diffraction patterns were treated to mask the Bragg peaks of the diamond anvils, then integrated into one-dimensional patterns, and a smooth polynomial background was subtracted. Le Bail refinements were performed using Topas Academic version 6 [52]. In the XRD measurements, the pressure was determined using the equation of state of gold from ref. [53]. Raman spectra at ambient temperature were acquired using a commercial Horiba Jobin-Yvon LabRam HR800 Raman spectrometer in a backscattering geometry, equipped with a COBOLT SambaTM green 532 nm laser, with a nominal power of 1000 mW. Low-temperature Raman spectra were collected using an Argon laser (120 mW output power) tuned to 514.32 nm and a Peltier-cooled HR460 spectrograph with 1500 gr/mm grating. The acquisition time was typically about a minute. Reference spectra were recorded at each pressure from the gasket close to the sample chamber, and particular care was taken to check the effect of the background subtraction on the fit results. In the Raman measurements, pressure was determined by the shift of the R1 ruby fluorescence line [54], from the edge of the stressed diamond signal (above 40 GPa) [55], or from the vibron frequency of pure solid H₂ [56]. A helium cryostat was used for low-temperature measurements.

## Computational Details

*Ab initio* molecular dynamics (AIMD) simulations were carried out on a 2 × 2 × 2 supercell of the C2 phase (space group Pna2₁) over a grid of pressures from 3 to 50 GPa and temperatures from 25 to 300 K. Calculations employed the Perdew-Burke-Ernzerhof (PBE) exchange-correlation functional [57] and were performed with both VASP [58] and Quantum Espresso [59, 60] packages. The Path-Integral molecular dynamics approach (PIMD) [61–65] was used to efficiently sample nuclear-quantum effects under the desired thermodynamic conditions. Simulations were conducted for temperatures above 25 K, by using the Path-Integral Ornstein-Uhlenbeck Dynamics (PIOUD) algorithm implemented in Quantum Espresso [66]. After performing convergence tests on both the virial and primitive estimators of the quantum kinetic energy, the number of beads (N<sub>b</sub>) was set to 40 at
300 K and increased to 160 at 25 K for PIMD runs. For both AIMD and PIMD, a plane-wave cutoff of 100 Ry and a 2 × 2 × 2 k-point mesh were sufficient to achieve convergence in the DFT solution of the electronic problem. Trajectories were obtained with a time-step of 0.5 fs and a total duration of approximately 10 ps for each (P, T) condition.

Spin isomers such as ortho- and para- hydrogen are not taken into account in molecular dynamics (MD) simulations. Indeed, to access the true low-temperature physics of the C2 phase while retaining quantum effects and accounting for spin isomers, we solved the Schrödinger equation at 0 K for an H₂ molecule in an effective potential modeling the induced interactions between H₂ and the water lattice. As detailed in the Supplementary Information section A, going beyond previous studies [67–70], we evaluated the full, non-separable potential V<sub>ext</sub>(r, θ, φ) by performing self-consistent DFT calculations (Quantum Espresso, 300 eV cutoff, 5×4×5 k-mesh) for each configuration generated on a dense spherical grid of bond lengths and orientations [70]. These 0 K DFT energies directly capture the interaction of H₂ with the water lattice. We used several exchange-correlation functionals and their comparison is shown in Fig. A1. For consistency with the MD calculations, the main results are obtained using PBE.

The Schrödinger equation has been solved using a basis set of the form $\psi_n^{lm}(r, \theta, \phi) = \chi_n(r)Y_l^m(\theta, \phi)$, where $\chi_n(r)$ is the radial basis of the Morse potential and $Y_l^m(\theta, \phi)$ are the spherical harmonics. This approach can be viewed as a mean-field model, in which a single H₂ molecule is embedded in the field generated by the surrounding water lattice.

The spatial density at finite temperature T was obtained by treating the two nuclear-spin manifolds of H₂ independently. Assuming that ortho ($\mathcal{J}=1$) and para ($\mathcal{J}=0$) populations do not interconvert on the time scale of the simulation, the total density reads

$$ \rho(\mathbf{r}) = \sum_{i \in \text{ortho}} \frac{f_o e^{-\beta (E_i - E_0^o)}}{Z_o} |\psi_i(\mathbf{r})|^2 + \sum_{j \in \text{para}} \frac{f_p e^{-\beta (E_j - E_0^p)}}{Z_p} |\psi_j(\mathbf{r})|^2 $$

where

$$ Z_o = \sum_{i \in \text{ortho}} e^{-\beta (E_i - E_0^o)}, \quad Z_p = \sum_{j \in \text{para}} e^{-\beta (E_j - E_0^p)}, $$

Here $\beta = (k_B T)^{-1}$, $|\psi_n|^2$ and $E_n$ are the eigenfunctions and eigenvalues obtained from the Schrödinger equation evaluated in the relaxed hydrate lattice; $E_0^o$ and $E_0^p$ are the ortho and para ground-state energies, and $f_p : f_o = 1 : 3$ is the fixed statistical spin ratio. This approach neglects cage deformation effects, as the potential is evaluated using the relaxed lattice geometry. A key advantage of the embedded quantum model is that it naturally incorporates the correct ortho/para mixture of H₂. In contrast, PIMD simulations describe only the para state (the ground rotational state), while in reality the ortho species is more populated than the para one in a 3:1 ratio. Because the ortho states are more localized than the para state, this difference has a major impact on orientational properties, particularly
at the lowest temperatures. Accounting for the ortho/para mixture is therefore essential for a quantitatively accurate description of the orientational regimes and phase transitions. Another clear advantage of the quantum embedded model is the possibility to explore low-temperature regimes at constant computational effort, compared to PIMD whose effort increases dramatically as temperature is lowered. We therefore adopted the quantum embedded model as viable approach to study the quantum phase diagram of hydrogen filled ice as a function of pressure and temperature.

To characterize the three distinct orientational regimes for H₂ (quantum plastic, herringbone-like configurations, and nematicity [27, 34]) from the quantum embedded model, we introduce the orientation factor S as follows [71, 72].

In analogy with diffraction theory, where the scattering amplitude can be factorized into a lattice contribution and a local form factor, we define the orientation factor as the product of a collective geometric term and an intrinsic quantum factor:

$$ (1) \qquad S(P, T) = S_g(P) S_q(P, T). $$

The geometric contribution $S_g$ is evaluated on the 0 K equilibrium geometry obtained from DFT calculations at a given pressure P. It is defined from the orientation of each molecule i, through the angle $\theta_i$ between the crystallographic c axis (preferential axis of alignment) and the molecular axis. Its definition coincides with the nematic order parameter [73, 74]:

$$ (2) \qquad S_g = \frac{1}{N} \sum_{i=1}^{N} \left( \frac{3}{2} \cos^2 \theta_i - \frac{1}{2} \right). $$

This factor thus encodes the collective orientational arrangement of the atomic sites. By construction, $S_g = 1$ when all H₂ molecular axes are perfectly aligned along the crystallographic c axis, while $S_g = 0$ corresponds to a completely random, isotropic distribution of the molecular orientations in the simulation cell.

We then apply a *quantum dressing* to this geometric factor, obtained from the probability density of the quantum embedded model at pressure P and temperature T. For this purpose, we define the alignment tensor

$$ (3) \qquad \mathbf{Q}_q = \frac{3}{2} \int \rho(\mathbf{r}) \mathbf{u}_r \mathbf{u}_r \mathrm{d}^3\mathbf{r} - \frac{1}{2} \mathbf{I}, $$

where $\mathbf{u}_r$ is the unit vector along $\mathbf{r}$, and $\rho(\mathbf{r})$ is the diagonal part of the thermal density matrix. We define the quantum contribution $S_q$ as the largest eigenvalue of $\mathbf{Q}_q$. By construction, $S_q = 1$ when the quantum distribution of the H₂ molecular axes is perfectly aligned along a given direction while $S_q = 0$ corresponds to an isotropic quantum distribution. This quantum dressing ensures that the final orientation factor incorporates both the
orientation linked to the underlying geometry ($S_g$) and the intrinsic thermal quantum fluctuations encoded in the density matrix ($S_q$), in the same way as diffraction separates lattice geometry from the internal structure of the scatterer.

With this definition, the orientation factor $S(P, T)$ can be mapped across the phase diagram, as shown in Fig. 1.

$S(P, T)$ can also be computed from PIMD simulations at pressure $P$ and temperature $T$, whenever available, instead of the quantum embedded model. In this case, the orientation factor reads $S_{\text{PIMD}} = \int dX \rho(X) \frac{1}{NN_b} \sum_{i=1}^{NN_b} \left( \frac{3}{2} \cos^2 \theta_i(X) - \frac{1}{2} \right)$, where the integral runs over the path configurations $X$, and the sum runs over all $N$ H$_2$ molecules and the beads. Here, $\rho(X)$ denotes the normalized configurational probability density sampled during the path-integral molecular dynamics, so that the integral represents an ensemble average over beads and time steps. Notice how the latter equation contains implicitly both a geometrical contribution, given by the orientational arrangement of the H$_2$ molecules in the supercell, and thermal-quantum fluctuations, provided by the path-integral formalism. The orientation factor obtained from the PIMD simulations $S_{\text{PIMD}}$ at 300K is plotted in Fig. 5, by neglecting spin isomers effects and relying on the similar density distribution of the $J = 0$ and $J = 1$ states at room temperature.

To locate the transition between the quantum plastic and orientationally ordered phases, we employ a complementary global order parameter $M_{\text{tot}}$. It is defined as a structure-factor measure of orientational coherence, constructed from angular phases chosen such that symmetry-related orientations (phase/antiphase configurations, e.g. in the herringbone pattern) interfere constructively. Randomly oriented H$_2$ molecules then cancel out ($M_{\text{tot}} \approx 0$), while herringbone or nematic H$_2$ molecules add coherently ($M_{\text{tot}} \approx 1$). The precise mathematical definition, including the angular cumulative distribution function transforms and rank-m order parameters, is given in Section B of the SI. In molecular dynamics trajectories sampled over the $(P, T)$ grid, orientational transitions are identified by the vanishing of $M_{\text{tot}}$, which signals the loss of orientational order in the H$_2$ sublattice. For each pressure, the temperature at which $M_{\text{tot}}$ vanishes is referenced as the transition point.

Since $M_{\text{tot}}$ cannot be evaluated within the quantum embedded model, which describes individually embedded H$_2$, we instead compute the global order parameter from classical MD simulations. At the classical transition points, we evaluate the corresponding orientation factors, where this time thermal quantum fluctuations in Eq. 3 are replaced by a thermal average performed over the classical canonical ensemble in the simulation cell. Within the approximation of dominant electrostatic effects, the resulting ensemble-averaged value of $S$ at the classical transition points provides a threshold, $S^{\text{th}}$, for the critical value of $S$ computed with the quantum embedded model. Within the embedded framework, the transition temperature $T_P^*$ at pressure $P$ is then obtained from the
condition

$$S(P, T_p^*) = S^{th}$$

Plotting $T_p^*$ as a function of pressure gives the points of the white dashed transition line shown in Fig.1. The error bars on $T_p^*$ are obtained by propagating the uncertainty in $S^{th}$: we evaluate the inverse relation at the two bounds $S = S^{th} \pm \sigma_S^{th}$. The temperature uncertainty is then defined as $\Delta T = \frac{1}{2}|T_{high}^* - T_{low}^*|$, where $T_{high}^* (T_{low}^*)$ corresponds to the temperature extracted from the inverse map at the $S^{th} + \sigma_S^{th} (S^{th} - \sigma_S^{th})$ bound (see Fig. A4 of SI).

Finally, because $S_g$ is computed with respect to the c axis in Eq. 2, the orientation factor also provides direct information on the nematic orientation of the system, the largest values of S being taken in the nematic aligned phase.

The phonon spectrum and Raman cross section of the system were computed using density functional perturbation theory (DFPT) as implemented in the Quantum ESPRESSO package [75, 76], employing a unit cell and a 4×4×4 k-point grid. This is the only calculation performed entirely within a classical nuclear framework, as no implementation currently exists for computing Raman cross sections from PIMD trajectories.

To investigate the relationship between H₂ orientational order, H₂O-H₂ interactions, and proton symmetrization within the water sublattice, the radial distribution functions $g(r)$ as a function of pressure were computed from our MD trajectories for Hw-O, Hh-O, and O-O pairs. Here, Hw denotes protons belonging to the water sublattice, Hh refers to protons from H₂ molecules, and O represents oxygen atoms.

## Acknowledgments

We acknowledge the European Synchrotron Radiation Facility for provision of synchrotron radiation facilities at the ID15B and ID27 beam lines and assistance from G. Garbarino and S. Gallego Parra. We acknowledge U. Ranieri, A. Falenti, W. Kuhs and L. Ulivi for samples preparation. L.E.B. and M.R. acknowledge the financial support by the European Union - NextGenerationEU (PRIN N. 2022NRBLPT), the ANR-23-CE30-0034 EXOTIC-ICE. L.E.B., T.P. and A.N. acknowledge the financial support by the Swiss National Fund (FNS) under Grant No. 212889. S.D.C. acknowledges computational resources from CINECA, proj. IsC90-HTS-TECH and IsC99-ACME-C, and the Vienna Scientific Cluster, proj. 71754 "TEST". L.R. and A.M.S acknowledge GENCI for providing computational resources on the IDRIS Jean-Zay supercomputing clusters under project numbers 2024-A0160901387, 2025-A0180901387. L.R. and M.C. acknowledge GENCI for providing computational resources on the IDRIS Jean-Zay supercomputing clusters under project number A0170906493. M. C. thanks the European High Performance Computing Joint Undertaking (JU) for the partial support through the "EU-Japan Alliance in HPC" HANAMI project (Hpc AlliaNce for Applications and supercoMputing Innovation: the Europe - Japan collaboration).
## Author Contributions

L.E.B. designed the research. S.D.C. initiated the theoretical investigation. L.R. and S.D.C. performed the MD simulations and data analysis with guidance from M.C., A.M.S., and L.E.B.. T.P., L.A., A.N., R.G., and L.E.B. performed the XRD experiments. T.P., A.N., and M.R. performed the Raman experiments. T.P. and L.A. performed the XRD data analysis, A.N., R.G. and M.R. performed the Raman data analysis, with guidance from L.E.B.. L.R., T.P., and A.N. prepared the figures, with inputs from L.E.B., M.C., and A.M.S.. L.E.B, L.R., and T.P. wrote the manuscript with input from M.C., and A.M.S. All authors discussed the results and the manuscript.

## References

[1] Saman Alavi, Livia E. Bove, Niall J. English, Donglin Jiang, Rocio Semino, and Amadeu K. Sum. Porous solids for energy applications. *Journal of Chemical Physics*, 161(20):200401, 2024.

[2] Cynthia J. Jameson and David Stueber. Nuclear magnetic resonance line shapes of xe in the cages of clathrate hydrates. *Journal of Chemical Physics*, 120(21):10200–10211, 2004.

[3] Umbertoluca Ranieri, Michael Marek Koza, Werner F. Kuhs, Richard Gaal, Stefan Klotz, Andrzej Falenty, Dirk Wallacher, Jacques Ollivier, Philippe Gillet, and Livia E. Bove. Quantum Dynamics of H2 and D2 Confined in Hydrate Structures as a Function of Pressure and Temperature. *The Journal of Physical Chemistry C*, 123(3):1888–1903, January 2019. Publisher: American Chemical Society.

[4] Umbertoluca Ranieri, Leonardo del Rosso, Livia Eleonora Bove, Milva Celli, Daniele Colognesi, Richard Gaal, Thomas C. Hansen, Michael Marek Koza, and Lorenzo Ulivi. Large-cage occupation and quantum dynamics of hydrogen molecules in sII clathrate hydrates. *The Journal of Chemical Physics*, 160(16):164706, April 2024. _eprint: https://pubs.aip.org/aip/jcp/article-pdf/doi/10.1063/5.0200867/19891268/164706_1_5.0200867.pdf.

[5] L. Ulivi, M. Celli, A. Giannasi, A. J. Ramirez-Cuesta, D. Bull, and M. Zoppi. Quantum rattling of molecular hydrogen in clathrate hydrate nanocavities. *Physical Review B*, 76:161401, 2007.

[6] Min Xu, Shuji Ye, Amanda Powers, Ron Lawler, Nicholas J. Turro, and Zlatko Bačić. Inelastic neutron scattering spectrum of h₂@c₆₀ and its temperature dependence decoded using rigorous quantum calculations and a new selection rule. *Journal of Chemical Physics*, 139(6):064309, 2013.

[7] Stefano Longo, Davide Lonigro, Giovanni Lerario, Alessandro Foti, Emilio Fiordilino, Giovanni Mazzola, and Rosario Pucci. Quantum states of h₂⁺ and h₂ in an icosahedral potential well. *European Physical Journal D*, 77(29):29, 2023.

[8] J. S. Tse, C. Gutt, W. Press, D. D. Klug, S. Janssen, I. P. Swainson, and V. P. Shpakov. Rotational dynamics in methane hydrates. *Applied Physics A*, 74:S1299–S1301, 2002.

[9] Hideki Hirai and Hidetaka Kadobayashi. Significance of the high-pressure properties and structural evolution of gas hydrates for inferring the interior of icy bodies. *Progress in Earth and Planetary Science*, 10(1):3, 2023.
[10] Timothy A. Strobel, E. Dendy Sloan, and Carolyn A. Koh. Raman spectroscopic studies of hydrogen clathrate hydrates. *The Journal of Chemical Physics*, 130, 1 2009.

[11] Surya Vellamarthodika and Sandeep Gautam. Role of orientational disorder in zsm-22 in the adsorption of so₂. *arXiv preprint*, 2022.

[12] H. Jobic, A. N. Fitch, and A. Tuel. Molecular motions of propane and propene in zsm-5: A quasielastic neutron scattering study. *Journal of Physical Chemistry B*, 103(4):660-664, 1999.

[13] F. Leroux, P. Boulet, and P. A. Jaffrès. Interaction of so₂ with zeolite materials studied by ftir spectroscopy and dft calculations. *Journal of Physical Chemistry C*, 112(43):16833-16839, 2008.

[14] Ishu Dhiman, Rounak Sarkar, Salman Alikhan, and Michael L. Klein. Anomalous orientational ordering and dynamics of so₂ in zsm-22 zeolite channels: A molecular dynamics study. *arXiv preprint*, 2023.

[15] U. Vietze, O. Kühn, and D. Michel. Molecular orientation and ordering of dyes in alpo₄-5 zeolite channels. *Physical Review Letters*, 81(21):4628-4631, 1998.

[16] Marina Carravetta, Odd G. Johannessen, and Malcolm H. Levitt. Beyond the motional narrowing limit: Nuclear spin relaxation of encapsulated h₂ in c₆₀. *Physical Review Letters*, 97(16):167401, 2006.

[17] A. J. Horsewill. Quantum dynamics of molecular rotation in the solid state. *Progress in Nuclear Magnetic Resonance Spectroscopy*, 35(4):359-389, 1999.

[18] L. E. Bove and U. Ranieri. Salt- and gas-filled ices under planetary conditions. *Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences*, 377(2146):20180262, 2019.

[19] Takuo Okuchi, Masashi Takigawa, Jinfu Shu, Ho kwang Mao, Russell J. Hemley, and Takehiko Yagi. Fast molecular transport in hydrogen hydrates by high-pressure diamond anvil cell nmr. *Physical Review B*, 75:144104, 4 2007.

[20] Sofiane Schaack, Umbertoluca Ranieri, Philippe Depondt, Richard Gaal, Werner F. Kuhs, Andrzej Falenty, Philippe Gillet, Fabio Finocchi, and Livia E. Bove. Orientational Ordering, Locking-in, and Distortion of CH₄ Molecules in Methane Hydrate III under High Pressure. *The Journal of Physical Chemistry C*, 122(20):11159-11166, May 2018. Publisher: American Chemical Society.

[21] Chengyang Yuan, Hongxiang Zong, Hongsheng Dong, Lei Yang, Yufei Gao, Zhen Fan, Lunxiang Zhang, Jiafei Zhao, Yongchen Song, and John S. Tse. Pressure-regulated rotational guests in nano-confined spaces suppress heat transport in methane hydrates. *Nature Communications*, 15(9477), 2024.

[22] Sofiane Schaack, Umbertoluca Ranieri, Philippe Depondt, Richard Gaal, Werner F. Kuhs, Philippe Gillet, Fabio Finocchi, and Livia E. Bove. Observation of methane filled hexagonal ice stable up to 150 gpa. *Proceedings of the National Academy of Sciences*, 116(33):16204-16209, 2019.

[23] L. Andriambariijaona, T. Poręba, S. Di Cataldo, R. Gaal, U. Ranieri, M. Santoro, T. C. Hansen, M. Mezouar, G. Tobie, and L. E. Bove. Hydrogen hydrates under extreme conditions: Insights into high-pressure phases and implications for planetary interiors. *Physical Review B*, 111:214109, 6 2025.

[24] Konstantin A. Lokshin, Yusheng Zhao, Duanwei He, Ho-Kwang Mao, Wendy L. Mao, Russell J. Hemley, Mikhail V. Lobanov, Martha Greenblatt, and Willem L. Vos. Structure and dynamics
of hydrogen molecules in the novel clathrate hydrate by high pressure neutron diffraction. *Physical Review Letters*, 93(12):125503, 2004.

[25] Leonardo del Rosso, Milva Celli, Daniele Colognesi, Francesco Grazzi, and Lorenzo Ulivi. Irreversible structural changes of recovered hydrogen hydrate transforming from c0 phase to ice xvii. *Chemical Physics*, 544:111092, 4 2021.

[26] Paulo H. B. Brant Carvalho, Amber Mace, Inna Martha Nangoi, Alexandre A. Leitão, Chris A. Tulk, Jamie J. Molaison, Ove Andersson, Alexander P. Lyubartsev, and Ulrich Häussermann. Exploring high-pressure transformations in low-z (h2, ne) hydrates at low temperatures. *Crystals*, 12(1), 2022.

[27] Umbertoluca Ranieri, Simone Di Cataldo, Maria Rescigno, Lorenzo Monacelli, Richard Gaal, Mario Santoro, Leon Andriambariijaona, Paraskevas Parisiades, Cristiano De Michele, and Livia Eleonora Bove. Observation of the most H2-dense filled ice under high pressure. *Proceedings of the National Academy of Sciences*, 120(52):e2312665120, December 2023. Publisher: *Proceedings of the National Academy of Sciences*.

[28] H. Hirai. Stabilizing of methane hydrate and transition to a new high-pressure structure at 40 gpa. *American Mineralogist*, 91:826-830, 5 2006.

[29] Roberto Bini and Gabriele Pratesi. High-pressure infrared study of solid methane: Phase diagram up to 30 gpa. *Physical Review B*, 55:14800-14809, 6 1997.

[30] Maxim Bykov, Elena Bykova, Chris J. Pickard, Miguel Martinez-Canales, Konstantin Glazyrin, Jesse S. Smith, and Alexander F. Goncharov. Structural and vibrational properties of methane up to 71 gpa. *Physical Review B*, 104:184105, 11 2021.

[31] Ho-kwang Mao and Russell J Hemley. Ultrahigh-pressure transitions in solid hydrogen. *Reviews of modern physics*, 66(2):671, 1994. Publisher: APS.

[32] Russell J. Hemley and Ho kwang Mao. Dense molecular hydrogen: Order, disorder, and localization. *Journal of Non-Crystalline Solids*, 205-207:282-289, 10 1996.

[33] Isaac F. Silvera and Ranga Dias. Phases of the hydrogen isotopes under pressure: metallic hydrogen. *Advances in Physics: X*, 6, 1 2021.

[34] Jingyun Zhang, Jer-Lai Kuo, and Toshiaki Iitaka. First principles molecular dynamics study of filled ice hydrogen hydrate. *The Journal of Chemical Physics*, 137, 8 2012.

[35] Hisako Hirai, Shingo Kagawa, Takehiko Tanaka, Takahiro Matsuoka, Takehiko Yagi, Yasuo Ohishi, Satoshi Nakano, Yoshitaka Yamamoto, and Tetsuo Irifune. Structural changes of filled ice ic hydrogen hydrate under low temperatures and high pressures from 5 to 50 gpa. *The Journal of Chemical Physics*, 137, 8 2012.

[36] Shin ichi Machida, Hisako Hirai, Taro Kawamura, Yoshitaka Yamamoto, and Takehiko Yagi. Structural changes of filled ice ic structure for hydrogen hydrate under high pressure. *The Journal of Chemical Physics*, 129, 12 2008.

[37] Shin ichi Machida, Hisako Hirai, Taro Kawamura, Yoshitaka Yamamoto, and Takehiko Yagi. Raman spectra for hydrogen hydrate under high pressure: Intermolecular interactions in filled ice ic structure. *Journal of Physics and Chemistry of Solids*, 71:1324-1328, 9 2010.

[38] Lorenzo Monacelli, Maria Rescigno, Alasdair Nicholls, Umbertoluca Ranieri, Simone Di Cataldo, and Livia Eleonora Bove. Accepted: Hydrogen bond symmetrization in high-pressure ice clathrates. *Physical Review B*, 2025.

[39] K. Aoki, H. Yamawaki, M. Sakashita, and H. Fujihisa. Infrared absorption study of the
hydrogen-bond symmetrization in ice to 110 gpa. *Phys. Rev. B*, 54:15673-15677, Dec 1996.

[40] A. F. Goncharov, V. V. Struzhkin, M. S. Somayazulu, R. J. Hemley, and H. K. Mao. Compression of ice to 210 gigapascals: Infrared evidence for a symmetric hydrogen-bonded phase. *Science*, 273(5272):218-220, 1996.

[41] Ph. Pruzan, E. Wolanin, M. Gauthier, J. C. Chervin, B. Canny, D. Häusermann, and M. Hanfland. Raman Scattering and X-ray Diffraction of Ice in the Megabar Range. Occurrence of a Symmetric Disordered Solid above 62 GPa. *The Journal of Physical Chemistry B*, 101(32):6230-6233, August 1997. Publisher: American Chemical Society.

[42] Kazuki Komatsu, Takanori Hattori, Stefan Klotz, Shinichi Machida, Keishiro Yamashita, Hayate Ito, Hiroki Kobayashi, Tetsuo Irifune, Toru Shinmei, Asami Sano-Furukawa, and Hiroyuki Kagi. Hydrogen bond symmetrisation in D2O ice observed by neutron diffraction. *Nature Communications*, 15(1):5100, June 2024.

[43] Tomasz Poręba, Davide Comboni, Mohamed Mezouar, Gaston Garbarino, and Michael Hanfland. Tracking structural phase transitions via single crystal x-ray diffraction at extreme conditions: advantages of extremely brilliant source. *Journal of Physics: Condensed Matter*, 35(5):054001, dec 2022.

[44] Maria Rescigno, Alberto Toffano, Umbertoluca Ranieri, Leon Andriambariijaona, Richard Gaal, Stefan Klotz, Michael Marek Koza, Jacques Ollivier, Fausto Martelli, John Russo, Francesco Sciortino, Jose Teixeira, and Livia Eleonora Bove. Observation of plastic ice vii by quasi-elastic neutron scattering. *Nature*, 640:662-667, 4 2025.

[45] J. VAN KRANENDONK and G. KARL. Theory of the rotational and vibrational excitations in solid parahydrogen, and frequency analysis of the infrared and raman spectra. *Rev. Mod. Phys.*, 40:531-555, Jul 1968.

[46] Peter I. C. Cooke, Ioan B. Magdău, Miriam Peña-Alvarez, Veronika Afonina, Philip Dalladay-Simpson, Xiao-Di Liu, Ross T. Howie, Eugene Gregoryanz, and Graeme J. Ackland. Raman signal from a hindered hydrogen rotor. *Physical Review B*, 102(6):064102, 2020.

[47] Timothy A. Strobel, Maddury Somayazulu, and Russell J. Hemley. Phase behavior of h<sub>2</sub>O at high pressures and low temperatures. *The Journal of Physical Chemistry C*, 115:4898-4903, 3 2011.

[48] Miriam Peña Alvarez, Veronika Afonina, Philip Dalladay-Simpson, Xiao-Di Liu, Ross T. Howie, Peter I. C. Cooke, Ioan B. Magdău, Graeme J. Ackland, and Eugene Gregoryanz. Quantitative rotational to librational transition in dense H₂ and D₂. *The Journal of Physical Chemistry Letters*, 11(16):6626-6631, 2020.

[49] Simone Di Cataldo, Maria Rescigno, Lorenzo Monacelli, Umbertoluca Ranieri, Richard Gaal, Stefan Klotz, Jacques Ollivier, Michael Marek Koza, Cristiano De Michele, and Livia Eleonora Bove. Giant splitting of the hydrogen rotational eigenenergies in the c2 filled ice. *Phys. Rev. Lett.*, 133:236101, Dec 2024.

[50] Miriam Peña-Alvarez, Philip Dalladay-Simpson, Xiao-Di Liu, Veronika Afonina, Hui-Chao Zhang, Ross T. Howie, and Eugene Gregoryanz. Intensity of raman modes as a temperature gauge in fluid hydrogen and deuterium. *Journal of Applied Physics*, 125(2):025901, 01 2019.

[51] Jon H. Eggert, Eran Karmon, Russell J. Hemley, Ho kwang Mao, and Alexander F. Goncharov. Pressure-enhanced ortho-para conversion in solid hydrogen up to 58 gpa. *Proceedings of the National Academy of Sciences*, 96(22):12269-12272, 1999.
[52] Alan A. Coelho. *TOPAS* and *TOPAS-Academic*: an optimization program integrating computer algebra and crystallographic objects written in C++. *Journal of Applied Crystallography*, 51(1):210-218, Feb 2018.

[53] Yingwei Fei, Angele Ricolleau, Mark Frank, Kenji Mibe, Guoyin Shen, and Vitali Prakapenka. Toward an internally consistent pressure scale. *Proceedings of the National Academy of Sciences*, 104:9182-9186, 5 2007.

[54] G. Shen and et al. Toward an international practical pressure scale: A proposal for an IPPS ruby gauge (IPPS-Ruby2020). *High Pres. Res.*, 40:299-314, 2020.

[55] Yuichi Akahama and Haruki Kawamura. Pressure calibration of diamond anvil Raman gauge to 310GPa. *Journal of Applied Physics*, 100(4):043516, August 2006. Publisher: AIP Publishing.

[56] Ross T. Howie, Christophe L. Guillaume, Thomas Scheler, Alexander F. Goncharov, and Eugene Gregoryanz. Mixed molecular and atomic phase of dense hydrogen. *Physical Review Letters*, 108(12):125501, March 2012. Publisher: American Physical Society (APS).

[57] John P. Perdew, Kieron Burke, and Matthias Ernzerhof. Generalized gradient approximation made simple. *Phys. Rev. Lett.*, 77:3865-3868, Oct 1996.

[58] G. Kresse and J. Furthmüller. Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set. *Phys. Rev. B*, 54:11169-11186, Oct 1996.

[59] Paolo Giannozzi, Stefano Baroni, Nicola Bonini, Matteo Calandra, Roberto Car, Carlo Cavazzoni, Davide Ceresoli, Guido L Chiarotti, Matteo Cococcioni, Ismaila Dabo, Andrea Dal Corso, Stefano de Gironcoli, Stefano Fabris, Guido Fratesi, Ralph Gebauer, Uwe Gerstmann, Christos Gougoussis, Anton Kokalj, Michele Lazzeri, Layla Martin-Samos, Nicola Marzari, Francesco Mauri, Riccardo Mazzarello, Stefano Paolini, Alfredo Pasquarello, Lorenzo Paulatto, Carlo Sbraccia, Sandro Scandolo, Gabriele Sclauzero, Ari P Seitsonen, Alexander Smogunov, Paolo Umari, and Renata M Wentzcovitch. Quantum espresso: a modular and open-source software project for quantum simulations of materials. *Journal of Physics: Condensed Matter*, 21(39):395502, sep 2009.

[60] Paolo Giannozzi, Oscar Baseggio, Pietro Bonfà, Davide Brunato, Roberto Car, Ivan Carnimeo, Carlo Cavazzoni, Stefano de Gironcoli, Pietro Delugas, Fabrizio Ferrari Ruffino, Andrea Ferretti, Nicola Marzari, Iurii Timrov, Andrea Urru, and Stefano Baroni. Quantum espresso toward the exascale. *The Journal of Chemical Physics*, 152(15):154105, 04 2020.

[61] RP Feynman and AR Hibbs Quantum Mechanics. Path integrals. *Lecture Notes Phys*, 106:64, 1979.

[62] David Chandler and Peter G. Wolynes. Exploiting the isomorphism between quantum theory and classical statistical mechanics of polyatomic fluids. *The Journal of Chemical Physics*, 74(7):4078-4095, 04 1981.

[63] Kenneth S. Schweizer, Richard M. Stratt, David Chandler, and Peter G. Wolynes. Convenient and accurate discretized path integral methods for equilibrium quantum mechanical calculations. *The Journal of Chemical Physics*, 75(3):1347-1364, 08 1981.

[64] D. M. Ceperley. Path integrals in the theory of condensed helium. *Rev. Mod. Phys.*, 67:279-355, Apr 1995.

[65] Thomas E. Markland and Michele Ceriotti. Nuclear quantum effects enter the mainstream. *Nature Reviews Chemistry*, 2(3):0109, February 2018.

[66] Félix Mouhat, Sandro Sorella, Rodolphe Vuilleumier, Antonino Marco Saitta, and Michele
Casula. Fully quantum description of the zundel ion: Combining variational quantum monte carlo with path integral langevin dynamics. *Journal of Chemical Theory and Computation*, 13(6):2400-2417, 2017. PMID: 28441484.

[67] Francesco Sebastianelli, Minzhong Xu, Zlatko Bačić, Ronald Lawler, and Nicholas J. Turro. Hydrogen molecules inside fullerene c70: Quantum dynamics, energetics, maximum occupancy, and comparison with c60. *Journal of the American Chemical Society*, 132(28):9826-9832, 2010. PMID: 20583809.

[68] Minzhong Xu, Francesco Sebastianelli, and Zlatko Bačić. Quantum dynamics of h2, d2, and hd in the small dodecahedral cage of clathrate hydrate: Evaluating h2-water nanocage interaction potentials by comparison of theory with inelastic neutron scattering experiments. *The Journal of Chemical Physics*, 128(24):244715, 06 2008.

[69] Minzhong Xu, Francesco Sebastianelli, Zlatko Bačić, Ronald Lawler, and Nicholas J. Turro. H2, hd, and d2 inside c60: Coupled translation-rotation eigenstates of the endohedral molecules from quantum five-dimensional calculations. *The Journal of Chemical Physics*, 129(6):064313, 08 2008.

[70] Simone Di Cataldo, Maria Rescigno, Lorenzo Monacelli, Umbertoluca Ranieri, Richard Gaal, Stefan Klotz, Jacques Ollivier, Michael Marek Koza, Cristiano De Michele, and Livia Eleonora Bove. Giant splitting of the hydrogen rotational eigenenergies in the c₂ filled ice. *Phys. Rev. Lett.*, 133:236101, Dec 2024.

[71] L Stryer. Fluorescence energy transfer as a spectroscopic ruler. *Annual Review of Biochemistry*, 47(Volume 47, 1978):819-846, 1978.

[72] Doros N. Theodorou and Ulrich W. Suter. Shape of unperturbed linear polymers: polypropylene. *Macromolecules*, 18(6):1206-1214, 1985.

[73] Simon C. McGrother, Dave C. Williamson, and George Jackson. A re-examination of the phase diagram of hard spherocylinders. *The Journal of Chemical Physics*, 104(17):6755-6771, 05 1996.

[74] P.G. de Gennes and J. Prost. *The Physics of Liquid Crystals*. International Series of Monographs on Physics. Clarendon Press, 1993.

[75] Stefano Baroni, Stefano de Gironcoli, Andrea Dal Corso, and Paolo Giannozzi. Phonons and related crystal properties from density-functional perturbation theory. *Rev. Mod. Phys.*, 73:515-562, Jul 2001.

[76] Michele Lazzeri and Francesco Mauri. First-principles calculation of vibrational raman spectra in large systems: Signature of small rings in crystalline sio₂. *Phys. Rev. Lett.*, 90:036401, Jan 2003.

[77] Raymond C. Clay, Jeremy Mcminis, Jeffrey M. McMahon, Carlo Pierleoni, David M. Ceperley, and Miguel A. Morales. Benchmarking exchange-correlation functionals for hydrogen at high pressures using quantum Monte Carlo. *Physical Review B*, 89(18):184106, May 2014. Publisher: American Physical Society (APS).
Computational supplementary

**Supplementary Information A. Quantum solution of the embedded H₂ molecule in the crystal field**

Building upon the method described in the Supplementary Material of Ref.[70], we evaluated the full, non-separable external potential $V_{\text{ext}}(r, \theta, \phi)$ and solved the corresponding Schrödinger equation for the quantum embedded H₂ molecule where r is the H₂ molecule elongation, θ is the polar and φ the azimuthal angles describing the molecule orientation (see Fig. A2).

Quantum embedding leads to a system of quantum rigid rotor in an effective potential representing an isolated H₂ in the C2 phase. The associated Schrödinger equation reads:

$$ \left(-\frac{\nabla^2}{2\mu} + V_{\text{ext}}\right) \psi(\vec{r}) = \epsilon\psi(\vec{r}), $$

where $\mu = \frac{m_H}{2}$ is the reduced mass, with $m_H$ the mass of a hydrogen atom. Since $V_{\text{ext}} \neq 0$, the energy levels deviate from those of a free rotor, given by $E_l = \frac{\hbar^2}{2I}l(l+1)$, where $I = \mu R^2$ is the moment of inertia and $l$ is the angular momentum quantum number.

The external potential takes the form $V_{\text{ext}} = V_{\text{DFT}}(\vec{r}) = V(r, \theta, \phi) = R(r) + P_r(\theta, \phi)$. The potential has been computed numerically from the DFT-relaxed structure of the C2 hydrogen hydrate, using several DFT functionals reported in Fig. A1. A script generated all configurations on a spherical coordinate grid with a fixed center of mass, varying the bond length $r$ and orientation angles $\theta$ and $\phi$. The grid used has dimensions $n_r \times n_\theta \times n_\phi = 64 \times 32 \times 32$.

As in Ref.[70], a self-consistent DFT calculation was performed in the primitive unit cell for each configuration. Since the simulations are carried out at 0K, the resulting total energies represent the potential experienced by the H₂ molecule around the classical equilibrium geometry. These calculations were performed with Quantum Espresso [59, 60], on a k grid of 5x4x5 and a plane-waves expansion cut at 300 eV.

The radial part of the potential $R(r)$ is fitted using a Morse function, while the angular component is interpolated linearly between the computed data points. After having determined numerically the effective potential $V_{\text{ext}}$, we solved the corresponding Schrödinger equation numerically using the Lanczos algorithm in shift-invert mode, following the approach of Ref.[70]. At variance with Ref.[70], in the present case we took into account the full 3D potential in a non-separable form, by developing the Hamiltonian matrix elements on a basis set of the type $\psi_n^{lm}(r, \theta, \phi) = \chi_n(r)Y_l^m(\theta, \phi)$. Here, $\chi_n(r)$ are the eigenfunctions of the Morse potential obtained by solving the one-dimensional Schrödinger equation for the fitted radial part, and $Y_l^m(\theta, \phi)$ are spherical harmonics, which form a complete basis
set for the angular part of a free quantum rotor.

Since the Hamiltonian can be decomposed as $\hat{H} = \hat{T} + \hat{V}$, with:

$$ \hat{T} = \frac{\hbar^2}{2\mu} \left( \frac{1}{r^2} \frac{\partial}{\partial r} \left( r^2 \frac{\partial}{\partial r} \right) + \frac{1}{r^2} \hat{L}^2 \right), \quad \text{and} \quad \hat{V} = R(r) + P_r(\theta, \phi), $$

the Hamiltonian matrix elements in this basis are given by the following terms:

$$ \text{a. Kinetic term: } \langle n'l'm' | \hat{T} | nlm \rangle = \frac{\hbar^2}{2\mu} \left[ \int_0^\infty \frac{d\chi_{n'}^*(r)}{dr} \frac{d\chi_n(r)}{dr} r^2 dr - \int_0^\infty \chi_{n'}^*(r) \chi_n(r) l(l+1) dr \right] \delta_{ll'} \delta_{mm'}. $$

$$ \text{b. Potential term: } \langle n'l'm' | \hat{V} | nlm \rangle = \int_0^\infty \chi_{n'}^*(r) R(r) \chi_n(r) r^2 dr \delta_{ll'} \delta_{mm'} + \int_0^\infty \chi_{n'}^*(r) \chi_n(r) r^2 dr. $$

$A_{l'm'lm}(r)$, where the angular coupling coefficient is defined as:

$$ A_{l'm'lm}(r) = \int d\Omega Y_{l'}^{m'*}(\theta, \phi) P_r(\theta, \phi) Y_l^m(\theta, \phi). $$

All matrix elements can be computed numerically based on the evaluated potential and the chosen basis. The angular integrals $A_{l'm'lm}(r)$ are evaluated on the r-grid and interpolated accordingly.

Our method is compared to that of Ref. [70] in Fig. A1. When the full, non-separable potential is taken into account, the energy level splitting is reduced compared to the separated approach (PBE separated vs. PBE), where the radial and angular parts are integrated independently. This difference is due to a more accurate treatment of the interaction potential in the present formalism.

We also investigated the impact of the exchange-correlation functional within our method. The use of vdW-DF leads to results that are in better agreement with experiment than those obtained with PBE for DFT energy estimation, consistent with previous findings [77]. Since vdW-DF tends to slightly overestimate the lattice parameters, we computed the eigenvalues using the experimental volume at 3 GPa. As expected, a slight increase in the splitting is observed compared to the standard vdW-DF calculation. This is consistent with the pressure dependence of the splitting, as shown in Fig. 3 (c).

Next, we evaluated the influence of proton disorder on the eigenvalues. To do this, we generated external potentials for a set of proton-disordered configurations that satisfy the ice rules. These potentials exhibit significant variations in magnitude. For each configuration, the corresponding eigenvalues were computed and the average is obtained using an arithmetic mean. As shown in Fig. A1, this leads to a substantial reduction in the energy level splitting, bringing it closer to the experimental values. This reduction is attributed to the decrease in anisotropy of the effective potential when proton disorder is included.

Incorporating all these effects—namely, the use of vdW-DF and averaging over multiple proton-disordered structures—significantly increases the computational cost. In this work, our aim is to provide insight into the orientational transitions and the coupling between the H₂ and H₂O sublattices, as well as to highlight the role of quantum effects. Furthermore, since both the *ab initio* molecular dynamics (AIMD) and path-integral molecular dynamics simulations were carried out using the PBE functional due to computational cost, all
quantum embedded calculations shown in the main paper have been performed with the same functional and in the proton-ordered configuration (i.e., the ice VIII sub-lattice) for consistency and direct comparison.

<figure><img src="image_7.png" /><figcaption>FIGURE A1. Comparison of the computed eigenvalues at P=3GPa using our method (with various exchange-correlation functionals: PBE, BLYP, and vdW-DF) is shown alongside results obtained with the approach presented in Ref.[70] (PBE with a separated potential). The labels vdW exp. lat. and vdW exp. lat. proton disordered refer to calculations using the vdW-DF functional on the C2 structure with experimental lattice parameters, with and without averaging over disordered proton configurations, respectively. Experimental energy levels from Ref.[70] are shown as black dotted lines.</figcaption></figure>

**Supplementary Information B. Order parameter**

**B.1. Symmetries**

The orientational order of the H₂ sub-lattice is monitored through a global order parameter. Figure A2 (c) displays a representative two-dimensional histogram of the polar angle θ and the azimuthal angle φ gathered along a *T* = 50 K, *P* = 20 GPa molecular-dynamics trajectory. Eight bright spots appear, revealing the eight symmetry orientations available to an H₂ molecule in the herringbone phase. As discussed in the main text, increasing pressure forces the sticks, representing the H₂ geometry, to align; in the fully ordered regime, i.e. nematic, the eight spots collapse into two. Since φ is not well defined for θ = 0 we see in practice two bands as shown in A2 (b).
<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /></figure>

FIGURE A2. (a) Definition of the spherical coordinates used to describe the orientation of the H₂ molecule; the *x*, *y*, and *z* axes are aligned with the crystallographic axes. Histogram of the orientation angles sampled by H₂ molecules during an ab initio molecular dynamics simulation at *T* = 100 K and *P* = 50 GPa (b), *T* = 50 K and *P* = 20 GPa (c) and *T* = 300 K and *P* = 3 GPa (d).

Although eight maxima are visible, only two *physically distinct* orientations remain once crystalline symmetries and the head-to-tail indistinguishability of the sticks are taken into account:

(i) **Polar symmetry.** Because the colatitude of an unoriented stick is defined modulo inversion, the two spots at $\theta_1$ and $\theta_2 = \pi - \theta_1$ correspond to the *same* configuration.

(ii) **Azimuthal symmetry.** In the longitudinal direction one has $\phi_2 = \phi_1 + \pi$; the pair $(\phi_1, \phi_2)$ therefore describes a single molecular orientation. The remaining spots at $\phi_3 = \pi - \phi_1$ and $\phi_4 = 2\pi - \phi_1$ are likewise equivalent to each other.

Accordingly, at any pressure the eight histogram maxima partition into two families of four spots, leaving only *two* genuinely different orientational states in the herringbone phase. The quantum plastic case, with no orientational order, is shown in A2 (d).
## B.2. Global order parameter

As explained in the Methods section, we define the order parameter to locate the transition between the quantum plastic and orientationally ordered phases. This parameter is designed to capture the orientational coherence of the H₂ sublattice. It is constructed as a structure factor, which vanishes when the H₂ molecules are randomly oriented in the quantum plastic phase, and interferes constructively when the molecules are in the herringbone or nematic phase.

The global order parameter, $M_{tot}$, can be separated into two components: $M_{\theta}$, which characterizes the order in $\theta$, and $M_{\phi}$, which characterizes the order in $\phi$:

$$M_{tot} = \sqrt{M_{\theta}^{2} + M_{\phi}^{2}}.$$

These components can be written individually as:

$$M_{\theta} = \langle |\psi_{\theta}(t)| \rangle_{t}, \quad M_{\phi} = \langle |\psi_{\phi}(t)| \rangle_{t}.$$

Here, $\psi_{\theta}(t)$ and $\psi_{\phi}(t)$ are the rank-$m$ order parameters with $m_{\theta} = 4$ and $m_{\phi} = 2$ at time $t$. They are defined as:

$$\psi_{\theta}(t) = \frac{1}{N} \sum_{\mu=1}^{N} e^{im_{\theta}\theta_{\mu}^{\text{corr}}(t)}, \quad \psi_{\phi}(t) = \frac{1}{N} \sum_{\mu=1}^{N} e^{im_{\phi}\phi_{\mu}^{\text{corr}}(t)},$$

The sum runs over the $N$ H₂ molecules in the system ($\mu = 1, \dots, N$). The angles $\theta_{\mu}^{\text{corr}}$ and $\phi_{\mu}^{\text{corr}}$ are the polar and azimuthal angles of the H₂ molecules, corrected to give the same phase in symmetrically equivalent configurations, as discussed in the previous section. The goal is to have an order parameter that vanishes when the molecules are disordered and equals 1 when the molecules are perfectly aligned, following the structure factor concept.

For each molecule, the bond vector $\mathbf{r}_{\mu}$ is expressed in spherical coordinates ($r_{\mu}$, $\theta_{\mu}$, $\phi_{\mu}$) in the crystal frame:

$$\theta_{\mu} = \arccos\left(\frac{y_{\mu}}{r_{\mu}}\right), \quad \phi_{\mu} = \arctan 2(z_{\mu}, x_{\mu}) \in [-\pi, \pi].$$

In order for the symmetric configurations mentioned above to yield the same phase in the herringbone-like situation, one possible estimator for the $\theta$ angle is $\sin\theta$. We flatten the distribution of $\sin\theta$ via the cumulative distribution function (CDF) transform:

$$\theta_{\mu}^{\text{corr}} = \frac{\pi}{2} (1 - \cos \theta_{\mu}) \in [0, \pi].$$
This transformation ensures that, for an isotropic rotor, the phase $exp(i\theta_{\mu}^{corr})$ is uniformly distributed on the circle, allowing the terms to cancel out statistically.

A similar approach is applied to $\phi$, since symmetric configurations have the same phase when using $|\sin \phi|$. Then, using the CDF to obtain a continuous uniform distribution, we have:

$$ \Phi_{\mu}^{\text{corr}} = 2 \arcsin (|\sin \Phi_{\mu}|) \in [0, 2\pi[ $$

**Baseline.** To account for the finite size of our simulation (both in time and length), we subtract a baseline: the same order parameter computed for a system of identical size and the same number of time steps, but with completely random H₂ positions.

We show in Fig. A3 the evolution of the global order parameter as a function of temperature for different pressure values. For example, at 3 GPa, we observe that $M_{tot}$ vanishes at 100 K.

Since $M_{tot}$ is a global parameter, it cannot be evaluated within the framework of the quantum embedded model, as it describes a single embedded H₂. As shown in the Methods section, in the quantum embedded model, the appropriate parameter is the orientation factor $S$. Thus, to predict the transition line between the disordered and ordered phases in the quantum embedded model, we need to find a suitable criterion, i.e., the critical value of $S$, that can indicate the transition from quantum plastic to ordered phases.

To do this, we compute the orientation factors $S$ at the transition points predicted by $M_{tot}$ from molecular dynamics simulations at different pressures. This provides a criterion to locate the phase transition in terms of $S$ in the quantum embedded framework.

In practice, to evaluate $S$ from molecular dynamics runs at the transition point (for example, at 100 K and 3 GPa, as shown in Fig. A3), we need an estimate of $S_q$ as defined in Eq. 3 in the main paper. In the quantum embedded model, $S_q$ is defined using the alignment tensor involving the thermal density matrix, which is not obtained in the classical molecular dynamics (MD) runs. Under the approximation of dominant charge distribution effects via electrostatic interactions, we can replace thermal quantum effects by classical thermal fluctuations, and estimate the critical value of $S_q$ for the quantum model as a thermal average performed over the classical canonical ensemble in the simulation cell at the classical transition point. Since the effective potential $V_{ext}$ of the quantum embedded model is obtained from the classical geometry at zero temperature, the geometric contribution $S_g$ to the orientation factor is the same for both quantum embedded and classical MD frameworks. Thus, finding the threshold for $S_q$ is sufficient. Moreover, we found that within the pressure and temperature resolution of our classical MD simulations, $S_q$ at the transition can be safely taken as pressure independent. By averaging the obtained estimations of $S_q$ at the transition point over pressure (see Tab. A1),
<figure><img src="image_12.png" /><figcaption>FIGURE A3. Global order parameter, $M_{tot}$, plotted as a function of temperature for AIMD trajectories.</figcaption></figure>

<table><caption>TABLE A1. Temperature of transition T estimated within our temperature resolution from Fig A3 and the associated $S_q$ at the transition point.</caption><thead><tr><th>P (GPa)</th><th>3</th><th>10</th><th>15</th><th>20</th></tr></thead><tbody><tr><td>T (K)</td><td>100</td><td>100</td><td>125</td><td>150</td></tr><tr><td>S<sub>q</sub></td><td>0.24 ± 0.08</td><td>0.25 ± 0.08</td><td>0.27 ± 0.07</td><td>0.27 ± 0.07</td></tr></tbody></table>

we determined the threshold $S_q^{th} = 0.26 \pm 0.04$ for the fluctuating part. Applying this to the datapoints plotted in Fig. A4, we can find the temperatures at which the transition occurs in the quantum embedded framework. We also displayed the upper and lower bounds $S_q^{th} \pm \sigma_q^{th}$ used to determine the errors on temperature in Fig. A4, as explained in Methods.

**Supplementary Information C. Free Rotor to Quantum Oscillator**

Figure A5 is a visualisation of how the rotons of $H_2 - H_2O$ and $D_2 - D_2O$ change as the pressure is increased to 40 GPa at 300 K. In panels A5(a) and (b), the $S_0(0)$, $S_0(1)$, $S_0(2)$ and $S_0(3)$ rotons are shown with increasing Raman shift frequency, from left to right, respectively. The $S_0(1)$ is the most intense. The intensity of all the rotons decrease as the pressure is increased, whilst the peaks significantly broaden.

The evolution of the rotons of $D_2 - D_2O$, shown in panels A5(c) and (d), is different since the Raman spectrometer that was used for the experiments could not well resolve the $S_0(0)$ roton, and they were exceedingly weak in comparison to the $H_2 - H_2O$ case as well, as indicated in the intensity scale of the different panels. Nonetheless, the $S_0(1)$
<figure><img src="image_13.png" /></figure>

FIGURE A4. Quantum contribution $S_q$ to the orientation factor evaluated from the quantum embedded model as a function of temperature for different pressures, together with the computed threshold $S_q^{th}$ used to determine the transition temperature. In dashed lines are displayed the upper and lower bounds $S_q^{th} \pm \sigma_q^{th}$ used to determine the error bars on the transition temperature.
roton was clearly seen and traced in pressure. This peak broadens significantly more than
the equivalent HH $S_0(1)$ roton.

Figure A6 shows the general broadening of the $S_0(0)$ and $S_0(1)$ rotons in $H_2 - H_2O$ and $D_2 - D_2O$ as pressure is increased to 40 GPa at 300 K. The rotons were modelled with three components, referring to the $S_0(x)_0$, $S_0(x)_1$ and $S_0(x)_2$ modes ($x = [0, 1]$). These were iteratively fitted, using the previous fitting to ensure that the peaks evolved naturally from pressure to pressure, and that there were no sudden changes to the fitting regime. As can be seen, our measured experimental data has been compared to the experimental data of pure $H_2$ and $D_2$ as reported by Pena-Alvarez et al. [48]. The rotons of HH shows a good agreement with the rotons of pure $H_2$ or $D_2$.

Note that the $S_0(1)$ roton progression in pressure is not shown for $D_2 - D_2O$. Again, this is due to the fact that the Raman spectrometer that was used did not have a good enough resolution below $180 \text{ cm}^{-1}$ to determine the weak roton from the large background from the elastic line.

Figure A7 compares the 'Centre of Mass' of the $S_0(1)$ rotons for $H_2 - H_2O$ and $D_2 - D_2O$ in pressure. The 'Centre of Mass' of the peak was determined by fitting a single peak to the data, instead of the three components of $S_0(1)_0$, $S_0(1)_1$ and $S_0(1)_2$. The grey band indicates the region where we observe the change in regime from a free roton to a quantum harmonic oscillator potential, and the enclathrated $H_2$ moves from a herringbone-like arrangement and the excitation of the molecules adopt a librational character.
<figure><img src="image_14.png" /></figure>

FIGURE A5. Raman shift for $H_2 - H_2O$ and $D_2 - D_2O$ as pressure is increased to 40 GPa at 300 K, visualised in a 3D-surface plot and a 2D-heat plot. Panels (a) and (b) show, with increasing frequency, the four $S_0(0)$, $S_0(1)$, $S_0(2)$ and $S_0(3)$ rotons of $H_2 - H_2O$ in pressure, whilst panels (c) and (d) only show the weak $S_0(1)$ roton of $D_2 - D_2O$.
<figure><img src="image_15.png" /></figure>

FIGURE A6. High-Pressure Raman measurements of the rotons of $H_2 - H_2O$ and $D_2 - D_2O$. The different roton components $S_0(x)_0$, $S_0(x)_1$ and $S_0(x)_2$, are shown by the blue, orange and magenta branches respectively. Open triangles represent measurements of pure $H_2$ or $D_2$ as reported by Pena-Alvarez et al. [48], filled circles represent measurements of $H_2 - H_2O$ or $D_2 - D_2O$ up to 40 GPa, and the thick grey line represents the position of the 'Centre of Mass' of the entire $S_0(x)$ peak. Panel (a) shows the $S_0(0)$ roton for $H_2 - H_2O$, panel (b) shows the $S_0(1)$ roton for $H_2 - H_2O$ and panel (c) shows the $S_0(1)$ roton for $D_2 - D_2O$.

<figure><img src="image_16.png" /></figure>

FIGURE A7. Comparison of the $S_0(1)$ roton 'Centre of Mass' of $H_2 - H_2O$ (blue) and $D_2 - D_2O$ (orange) in pressure. The grey shaded block indicates the change from a free rotor to a quantum oscillator regime, as detailed in the main text.