Large-cage occupation and quantum dynamics of hydrogen molecules in sII clathrate hydrates

Molecular hydrogen is only feebly soluble in ice at ambient pressure, but when pressure is increased, this gas can penetrate into the solid, giving rise to different phases of the mixed crystalline solid depending on pressure and temperature.¹ The solid H₂O–H₂ system has been largely studied in the past, and several phases of the mixture, stable at different temperatures and pressures, have been identified. Hydrogen clathrate hydrates have been first identified over 20 years ago^2,3 and, together with other hydrogen hydrate systems, are currently considered as candidate materials for efficient solid hydrogen storage.^4–9 Here, we are interested in the compound that is stable in the pressure range roughly below 400 MPa at temperatures lower than 0 °C. This is a non-stoichiometric clathrate hydrate, having the so-called cubic structure II (abbreviated as sII) common also to other water–gas solid compounds.¹⁰ In this structure, the H₂O molecules are connected by hydrogen bonds and form cavities of two different sizes, named small cages (SC) and large cages (LC), including either 20 (SC) or 28 (LC) H₂O molecules, arranged so that the oxygen atoms sit at the vertices of either a pentagonal dodecahedron (for SC) or a hexakaidecahedron (for LC), i.e., a solid having 12 pentagonal and four hexagonal faces.¹⁰ An atomistic view of this arrangement is shown in Fig. 1. Simple (i.e., with only one type of guest species) hydrogen clathrate hydrate forms at pressures above $≈100$ MPa, with the H₂ molecules confined both in the SC (one molecule per cage) and in the LC (up to four molecules per cage, depending on the thermodynamic synthesis conditions).^11–13 

The study of the dynamics of nanoconfined H₂ molecules, which is intrinsically quantal, is of great importance from both an applied and a fundamental point of view as it allows us to model the interaction of the guest hydrogen with the water cage. In the singly occupied SCs of sII clathrate hydrates, H₂ molecules perform an almost free rotation and a deeply anharmonic center-of-mass (CM) vibrational motion (also known as “rattling”). This motion has been experimentally investigated at low temperatures by spectroscopic means, making use of both inelastic neutron scattering (INS)^14–18 and Raman spectroscopy.^13,19–21 More recently, the CM motion of H₂ in the LCs has been probed by quasi-elastic neutron scattering²² and INS²³ experiments at low temperatures. Such a motion was found to give rise to a triplet at low energy transfer values (located between 1 and 3 meV) in the INS spectra.²³ In the present paper, we further investigate the same triplet at low energy transfer values by measuring sII clathrate hydrate samples containing various amounts of hydrogen guests, ranging from almost completely filled to mostly empty cages. First, detailed neutron diffraction (ND) measurements were carried out in order to characterize the structure of the samples, both pristine (i.e., as prepared) and undergoing an evacuation process. Second, INS data were collected at various hydrogen filling levels, providing additional information on the H₂ dynamics and showing some counterintuitive results, such as an increase in intensity of the mentioned rattling bands as the average H₂ level decreases. Finally, an interpretation of the band intensities as a function of momentum transfer making use of analytical models containing no fitting factors is reported, surpassing the previous²³ suggestion of a jump motion between four hypothetical equilibrium positions inside the LCs. In our view, this represents a substantial extension of the study contained in Ref. 23 and a clear interpretation of the spectroscopic bands at low energy transfer values in terms of the H₂ rattling motion in singly occupied LCs has now been achieved.

The rest of this paper will be organized as follows: Sec. II will deal with the sample preparation as well as with the performed ND and INS experiments, Sec. III will present the ND data analysis and results, and Sec. IV will be fully devoted to the INS data analysis procedure, the result presentation, and their scientific discussion. Finally, Sec. V will conclude the present study providing some perspectives for possible future theoretical work on this subject.

As illustrated in Table I, the samples investigated in the present study, i.e., sII clathrate hydrates containing H₂ and D₂, were prepared following two different production routes: one performed at Istituto di Fisica Applicata “N. Carrara” (IFAC) in Sesto Fiorentino, Italy, and the other performed at Helmholtz-Zentrum Berlin (HZB) in Berlin, Germany. In the former synthesis protocol, about 20 g of D₂O ice Ih powder, previously prepared by finely grinding a block of ice in an inert N₂ atmosphere, was pressurized in a Cu–Be autoclave with hydrogen gas up to a pressure of about 0.30 GPa at a temperature of 257 K for 6–12 h. After this operation, we recovered the sample at ambient pressure and a temperature of 77 K by quenching the autoclave in a liquid nitrogen bath and then releasing the pressure. The collected clathrate hydrate powder was inserted into cryogenic vials and subsequently stored in a liquid nitrogen Dewar suitable for shipment to Institut Laue-Langevin (ILL, Grenoble, France) for the neutron scattering measurements. The quality of the recovered clathrate hydrate samples was checked by inspecting the Raman spectra shown in Fig. 2. In the left panels of this figure, one can recognize the band originated from the lattice translational motion of the D₂O molecules in the spectral range 180–310 cm⁻¹ together with strong and well distinct bands originated from the rotational transitions S₀(j) (i.e., with ν = 0 and Δj = 2) of D₂ (three bands) and H₂ (two bands) molecules, with ν and j indicating the vibrational and rotational quantum numbers of the guest molecule, respectively. The width of these rotational bands is larger than in the free gas phase due to a more intense interaction of the confined guest molecules with the cages, both small and large. Moreover, from the Raman spectra of the vibrational mode of D₂ and H₂ (right panels of Fig. 2), one can observe the roto-vibrational transitions Q₁(j) (i.e., with ν = 1 and Δj = 0) of D₂ and H₂ molecules, verifying the presence of these guest molecules in both SCs and LCs, while no free gas phase is detected.^12,13 In the second protocol, sub-millimetric ice Ih spheres, obtained by spraying ultra-pure D₂O water into a liquid N₂ bath within a N₂-filled glovebox to prevent moisture condensation, were exposed to an H₂ pressure of 0.28 GPa at a temperature of 244 K for about 0.5 h. In this case, the hydrate samples were grown and collected in small aluminum vials and then stored in a Dewar suitable for shipment as in the case of the first protocol. It is important to mention that for samples prepared using both production routes, the obtained occupation level of the SCs is always close to 1. Meanwhile, LC occupation strongly depends on the synthesis conditions: samples prepared at HZB typically contain 1.5–1.8 H₂ or D₂ molecules on average in their LCs,²³ while for the samples prepared at IFAC, a mean occupation level of the LCs close to 3 is systematically obtained in both H₂ and D₂ cases. In particular, for the IFAC sample batches used in the present study, the estimated LC occupation level close to 3 was inferred via Raman spectroscopy measurements as described in Ref. 13 (see also Table I). The major difference in the LC occupation between the two mentioned sample sets is probably due to the reduced size of the clathrate crystallites one can obtain when the grinding method for the ice powder production is adopted (as done at HZB). The H₂ uptake under gas pressure is essentially comparable for the two production routes, but it is partially reversible when the sample is recovered at ambient pressure and low temperatures. The reduced size of the clathrate hydrate crystallites implies a faster (and, therefore, a larger) decrease in the LC occupation during such a recovery process (see also our results in Sec. III).

The first part of this study is devoted to determining the kinetics of the gas release from both types of clathrate hydrate samples. To this end, we performed an experiment on the high-intensity D20 neutron diffractometer at ILL. For all the ND-measured samples, a 30 mm-long and 8 mm-thick cylindrical aluminum vial containing the sample powder was inserted into an appropriate vanadium can, which was equipped with a capillary connected to the pumping system via a gas rig. After each loading, we verified the absence of the most intense reflections belonging to ice Ih and ice Ic in the collected diffractogram, in order to check that the starting ice powder was completely transformed into clathrate hydrate during the sample preparation procedure and no reverse transformation of the sample under measurement occurred. Sample no. 2 actually contained a very small amount (about 1%) of ice Ic. In order to probe the emptying of the deuterated clathrate hydrate cages, we applied dynamic vacuum to the sample can for several hours at the fixed temperature of 130 K, a procedure that guaranteed a significant out-diffusion of the guest molecules from the LC without exiting the thermodynamic metastability region of the sII clathrate hydrate structure.²⁴. Five-minute acquisitions with λ = 2.419 Å allowed us to obtain the best compromise between the count rate and d-space resolution. For calibration purpose, the diffraction pattern of an Na₂Ca₃Al₂F₁₄ powder sample, arranged in the same experimental configuration, was also measured at the end of the ND experiment.

The second part of the present work is devoted to the INS study of the CM vibrational dynamics of the guest molecules inside the clathrate hydrate cages as a function of their decreasing occupation level. Since we are interested in the study of the dynamics of a single guest molecule inside an LC, given the faster kinetics of the LC evacuation for H₂ (which requires 8 h as a maximum to get a mean N_LC ≤ 1.0) with respect to D₂, we chose to focus our INS experiments only on the H₂–D₂O clathrate hydrate system, namely on sample nos. 4 and 5 in Table I. The INS measurements were performed on the IN5 instrument²⁵ at ILL, a direct-geometry time-of-flight spectrometer, using λ = 3.50 Å. This wavelength value corresponds to an incoming neutron energy E_i = 6.68 meV, which allows us to attain a maximum energy transfer ℏω = 4.00 meV and a maximum wavevector transfer Q = 3.28 Å⁻¹. Making use of a chopper velocity v₀ = 12 000 rpm (revolution per minute), an elastic energy resolution Δℏω ≈ 240 μeV was achieved. This figure represents the full width at half maximum at ℏω = 0. As for the negative values of ℏω (neutron energy gain), only data with ℏω > −16.00 meV were considered. Both samples were contained in cylindrical aluminum vials, as in the D20 experiment, but in this case also the external can was made of aluminum. Sample no. 4 (belonging to the same preparation batch as the aforementioned ND sample no. 3) was measured as synthesized (i.e., with about three H₂ guests in an LC) and at different stages of the process of emptying the LCs from the guest H₂ molecules. The emptying procedure was carried out by heating the sample at T = 130 K under dynamic vacuum for ∼30 min. After each heating treatment, the sample was cooled down to base temperature [i.e., T = (1.5 ± 0.1) K] for a new INS measurement. This cycle was repeated four times on sample no. 4 in order to achieve a good single-occupation level of the LCs. Sample no. 5 coincided with the aforementioned sample no. 2, previously measured and partially emptied during the D20 experiment. It had an average of about 0.5 H₂ guests in each LC at the end of the D20 experiment and, subsequently, during the IN5 experiment. Raw experimental data were transformed from time-of-flight and scattering angle to ℏω and ℏQ using the appropriate instrument routines (i.e., the LAMP code). Then, the processed data were reduced to constant-Q spectra with the recorded intensity being a function of ℏω. They covered the experimentally accessible Q range, which stretches from 0.6 to 2.8 Å⁻¹, with a selected Q-spacing equal to 0.2 Å⁻¹. The samples always exhibited the crystalline structure of the type sII clathrate hydrates in both diffraction and spectroscopy experiments. Powder diffraction patterns were obtained directly on IN5 from the intensities of the elastic peak as a function of the scattering angle. In spite of the limited instrumental resolution, they were sufficient to identify the crystalline phase of the measured samples.

Comparing the diffraction patterns of the samples measured before and after the evacuation process (see Fig. 3), one can only see a difference in the relative intensities of the Bragg reflections, proving that no structural phase transition occurred, while there was a different occupation level of the crystallographic sites by the H₂ or D₂ guest molecules. As a function of time, for the D₂–D₂O sample (i.e., no. 1), we observe an increase in the intensity of the 111, 220, 311, 222, 400, and 331 reflections and a decrease in the intensity of the 511, 440, and 531 ones. This finding is in agreement with a previous report.²⁶ As for the H₂–D₂O samples (i.e., nos. 2 and 3), the situation is essentially reversed, with the reflections 111, 220, 311, and 222 losing intensity with time and the reflections 422, 511, 440, and 531 gaining it. In the case of the sample nos. 2 and 3, the removal of the H₂ guest molecules is also visually confirmed by the decrease in the background level [see panels (b) and (c) in Fig. 3], due to the much larger incoherent scattering cross section of the H₂ molecules with respect to D₂. Figure 3 also shows a small shift in the position of the peaks for sample nos. 1 and 3, particularly at high 2θ values, indicating a change (namely, a reduction) in the unit cell volume for these two samples along with time (at fixed temperature). The diffraction data were analyzed by means of the Rietveld refinement using the known sII clathrate crystal structure, which is common to many other gas hydrates for the water host lattice. The space group is Fd$3̄$m, and the cubic lattice constant a is about 17.2 Å. For the atomic positions of the guests, the structural model of Ref. 23 was applied. This model is a small adjustment to the model for the Ne-filled sII clathrate hydrate reported in Ref. 24 and has also been successfully employed in Ref. 27. It describes the orientationally disordered H₂ or D₂ molecule in the SCs with a 96g Wyckoff position (0, 0, ≈0.022), accounting for three orientations of a hydrogen molecule located in the SC center, which is at (0, 0, 0). The positionally and orientationally disordered H₂ or D₂ molecules in the LCs are also modeled with a 96g position (≈0.4, ≈0.43, ≈0.43), which yields a distribution of 12 H or D atoms around the LC center, which is at (3/8, 3/8, 3/8). The correlations among the parameters make it impossible to simultaneously refine atomic positions, thermal parameters, and crystallographic occupancies. So, in our refinements, the atomic positions and the thermal parameters for the guest molecules were restrained to the values of Ref. 23 and only the crystallographic occupancies were refined.

The obtained mean occupation level of the guest molecules in both cage types is shown in Fig. 4, where the following various results are worth noticing:

• The difference in the mean initial occupation number of the LC, N_LC, as measured by Raman spectroscopy at 77 K after the sample preparation (N_LC ≈ 3 as given in Table I) and by ND at 130 K at the beginning of the actual kinetic study (N_LC ≈ 2 as shown in panel (a) of Fig. 4 and reported in Table I) must be attributed to a very significant outgassing phenomenon happening at temperatures between 77 and 130 K, which appears to be particularly important in the case of hydrogenous guests. The sample transfer procedure could also have played a role in this effect.

• The kinetic curves in panel (a) of Fig. 4 show that H₂ molecules leave the LCs slightly faster than D₂ molecules, which could be a consequence of their lighter mass.^28,29 Whether a manifestation of residual quantum effects (analogous to the well-known “quantum sieving” in porous media^30,31), or, more simply, an expression of the classical Graham’s law for effusion,³² which justifies a difference in diffusion out of the cages proportional to the square root of the molecular mass ratio, this is still rather unclear. However, the latter scenario seems much more probable than the former considering the relatively high temperature of the evacuation procedure, i.e., T = 130 K.

• The H₂-containing clathrate hydrate sample produced at HZB exhibits a slightly faster decrease in the LC occupation if compared to that produced at IFAC [see panel (a) of Fig. 4]. This can be explained by the reduced size of the clathrate crystallites, resulting in a faster emptying of the samples. Such a peculiarity can be also held responsible for the lower starting filling rate for this sample.

• Interestingly, H₂ and D₂ molecules appear to leave the SCs at the same rate, which is also independent of the sample synthesis route: the three curves in panel (b) of Fig. 4 can be perfectly superposed by simply applying a vertical shift. In this context, the fact that some of our initial SC occupation values are slightly larger than 1.0 should be most likely interpreted as being due to the typical uncertainty level of our ND results for H₂.

• The occupation level of the SCs decreases by about 20% in 12 h. This rate happens to be much smaller than the evacuation rate of the LCs (namely, 42% and 63% in 12 h for D₂ and H₂ guests, respectively), and this is an expected behavior similar to that reported for the evacuation of the Ne guest atoms from the cages of the sII Ne clathrate hydrate.²⁴ This fact might be explained by the hindered diffusion of a guest molecule through the (small) pentagonal faces of an SC.³³ Conversely, hydrogen molecules can easily travel between large cages passing through the hexagonal faces.

• We can also notice a smaller rate of out-diffusion from the LCs for hydrogen molecules with respect to Ne atoms: in the latter case, the LCs resulted completely emptied after just 4 h at the same temperature of 130 K.²⁴ This discrepancy is probably due to the difference in the strength of the guest–host interactions in the case of Ne-clathrate as compared to H₂- and D₂-clathrate hydrates.

In Fig. 5, we report another result of the Rietveld refinement of the measured diffraction patterns: the lattice constant a at 130 K as a function of the total (weighted) hydrogen occupation in the small and large cages. For sample no. 2, which exhibited SCs that were always almost completely filled, a was rather constant throughout the evacuation process. However, in sample nos. 1 and 3, the lattice constant is significantly larger (by approximately a hundredth of an Å) for hydrogen fillings above 25 H₂ or D₂ molecules per unit cell, meaning that, not entirely surprisingly, the encaged hydrogen molecules expand the water host lattice to some extent. To the best of our knowledge, no such observation has been reported so far for the hydrogen clathrate hydrate or for clathrate hydrates of any other guest as small as H₂. However, Lokshin et al. in Ref. 11 described an anomalous decrease in the lattice constant of the sII hydrogen clathrate hydrate during heating from 90 to 130 K, in conjunction with a reduction of the D₂ occupation level of the LCs, roughly from four to three. This observation is in qualitative agreement with our results and has been attributed to the interaction between guest molecules and water framework.¹¹ Data on the link between cage fillings and lattice constants are, in general, scarce in the literature on gas clathrate hydrates, and the role of the guest–host attractive-vs-repulsive interactions is far from being understood. Some authors^34,35 previously suggested that there is a connection between cage filling and lattice constant of sI clathrate hydrates in order to reconcile their experimental results with the literature data. Specifically, this concerned an increase in the lattice constant of the Xe clathrate hydrate (at fixed temperature) for a high Xe occupation of the small and large cages, as well as in the case of the CO₂ clathrate hydrate for a high occupation of the SCs. An expansion of the lattice was also reported for the phase-C₀ hydrogen hydrate with increasing H₂ filling³⁶ and similarly for the isostructural Ne- and O₂-filled hydrates.³⁷ Finally, the empty sII clathrate remarkably exhibits a lattice constant larger than those in the N₂ and Ne sII clathrate hydrates^24,34 (e.g., with the Ne-filled one having 0.86 and 1.18 guest atoms in the small and large cages, respectively²⁴). This is, not unexpectedly, an effect exhibiting its largest magnitude at the lowest investigated temperatures.^24,34 By comparing our lattice constant values for hydrogen clathrate hydrate with that of the empty sII clathrate (also known as ice XVI), namely 17.142 Å at 130 K,³⁴ one finds that the values of a for the former system are larger than those for the latter at the same temperature. However, Ref. 34 also noticed that comparing absolute lattice constant values from independent investigations might lead to misleading conclusions because of possible systematic experimental errors.

The INS spectra are shown in Fig. 6 as a function of the emptying of the cages for sample no. 4, limited to the low-frequency range that is of interest in the present work, for some selected Q values. We did not repeat the same emptying procedure for sample no. 5 (also shown in Fig. 6). In this section, we will focus our attention on the low frequency zone (namely, from 0.5 to 3.5 meV) of the least-loaded spectra (e.g., those plotted in magenta or orange in Fig. 6) for both samples measured on IN5 (i.e., nos. 4 and 5 as in Table I). In these spectra, three peaks appear more and more intensely as the sample is being emptied, i.e., their intensity increases with the advance of the evacuation process (see Fig. 6). Taking into account the scientific literature on this subject as well as our previous experience,^15,23 we can interpret the three spectral features as being due to the rattling of a single H₂ inside an LC of the host and, in particular, the presence of the three non-degenerate peaks can be reasonably attributed to the anisotropy of the cage via the H₂-cage potential. Thus, the purpose of the following data reduction procedure will be obtaining for each peak its intensity as a function of Q in order to probe the possible anharmonic effects. The data analysis can be briefly summarized in four steps: first, for each constant-Q spectrum, an appropriate linear background is subtracted exploiting the regions where no spectroscopic signal is expected (i.e., 0.75 $<ℏω<0.95$ and 3.2 $<ℏω<3.5$ meV). Second, the background-corrected INS spectra are fitted using three Gaussian profiles plus an additional linear background. This simple fitting model contains 11 free parameters: intensity I_1,2,3(Q), mean position e_1,2,3(Q), and width w_1,2,3(Q) for each of the three Gaussian curves, plus the intercept p(Q) and the slope m(Q) of an additional linear background. Third, from the fit results, we can easily calculate the average centroid positions $(ē1,2,3)$ and widths $(w̄1,2,3)$ of the three experimental spectral peaks for both samples, as the incoherent nature of the excitations implies that no dispersion in Q is expected. So, one finds $ē1=1.305(7)$ meV, $ē2=1.83(2)$ meV, and $ē3=2.63(4)$ meV. Finally, we re-fit the background-corrected INS spectra using the same model as before, but this time constraining peak positions and widths to their average values, obtaining in this way better estimates of the peak intensities I_1,2,3(Q) (see Fig. 7 for some selected examples). In Fig. 8, we have reported these results as a function of Q for sample nos. 4 and 5, resulting in six bell-shaped curves.

Before studying the Q-behavior of the INS spectral bands observed at ℏω ≈ 1.305, 1.83, and 2.63 meV and plotted in Fig. 8, one needs to discuss their physical origin, show that they arise from the rattling excitation of a single H₂ molecule in the LC, and explain the counterintuitive trend of the spectral intensity, visible in Fig. 6, which grows as the H₂ evacuation process is continued. These three bands correspond to the fundamental rattling transition of the H₂ CM, namely that described by the following quantum numbers: (N₁ = 0, L₁ = 0) → (N₂ = 0, L₂ = 1), of which the meaning is explained in the  Appendix. In order to justify this assignment, we carried out a theoretical calculation by numerically solving the Schrödinger equation (SE) for the H₂-CM states. The potential energy for a single H₂ molecule, as a function of its CM displacement from the center of the hexakaidecahedron cage, has been calculated by summing the contributions due to all the H₂–D₂O pairs (assumed to be given by the standard hydrogen-water potential)³⁸ over 514 D₂O molecules placed around the cage center in a rigid lattice and subsequently averaged over the H₂ molecular orientations, following the procedure described in Ref. 14. The disorder of the deuterons belonging to the host water molecules has been taken into account by performing a second average over several random configurations, all of them complying with the well-known Bernal and Fowler ice rules.³⁹ A further average over the direction of the H₂-CM coordinate, $ŝ$ (see the  Appendix for details), provided the isotropic part of the potential. The solution of the radial Schrödinger equation for an H₂ molecule sitting in such an isotropic potential well determines the first quantum transition, that is, (N₁ = 0, L₁ = 0) → (N₂ = 0, L₂ = 1), characterized by an energy jump equal to 2.21 meV, which compares reasonably well with the observed $ē1,2,3$⁠, whose average value is $ēav=1.92(3)$ meV. Obviously, the splitting of this transition, which is experimentally observed, is not reproduced by this model, but it can be attributed to the small anisotropic components of the real interaction between an LC and a guest H₂ trapped inside it, able to completely lift such a degeneracy if containing both azimuthal (ϕ) and elevation (θ) parts. So, the presence of three bands is very likely due to the anisotropy of the H₂–LC interaction potential. As a matter of fact, in our case, the situation is slightly more complicated since the rotational degrees of freedom of H₂ (described by j and m) have to be properly taken into account. Now, in our experiment, for neutron cross section reasons,⁴⁰ the measured signal is dominated by the transition j₁ = 1 → j₂ = 1, so, in the case of an m degeneracy, no relevant rotational contribution to the energy jump is expected. However, we have to remark that, because of the interplay between the anisotropic components of the H₂–LC interaction and the H₂ intrinsic quadrupole, one may suppose that the j_1,2 = 1 rotational level involved is also split into three sublevels with m_1,2 = −1, 0, 1, so that, in principle, the j₁ = 1 → j₂ = 1 transition could exhibit several components, corresponding to the combinations of the three initial and the three final states (with three “purely elastic” terms having m₁ = m₂ actually coinciding). These terms of the rotational manifold that is always associated with each CM rattling transition are weighted by the rotational form factors a²(Q, j₁, m₁, j₂, m₂) calculated in Table II (see also the  Appendix for further details), and one can simply verify that in our experimental INS conditions, due to the Q values explored, only the three purely elastic terms [namely, (j₁ = 1, m₁ = −1) → (j₂ = 1, m₂ = −1), (j₁ = 1, m₁ = 0) → (j₂ = 1, m₂ = 0), and (j₁ = 1, m₁ = 1) → (j₂ = 1, m₂ = 1)] provide sizable contributions. The mathematical reason for this is given by the interatomic distance in H₂, r₀ ≃ 0.7414 Å, so that one obtains $j02Qr02≫j22Qr02$ for Q ≤ 4 Å⁻¹, which implies that in the present INS experiment, all the transitions with m₁ ≠ m₂ have a negligible intensity. Thus, even if we do not know the magnitude of the possible rotational splitting, we have shown that the observed INS triplet cannot have a rotational origin and, consequently, has to be ascribed to the rattling dynamics in a slightly aspherical potential well.

As for the other issue, the increase in the band intensity during the progress of the H₂ evacuation in sample no. 4, as a matter of fact, the detailed excited states of two (or more) ortho-H₂ molecules in a LC are still unknown, since it is already quite demanding the calculation of the ground state of such a system,⁴¹ where at least ten degrees of freedom have to be tackled. However, noticing that the increase in the mentioned bands in the 1 $<ℏω<$ 3 meV-interval proceeds without any relevant distortion, we can propose a reasonable explanation based on two simple assumptions: (A) the INS spectral features related to H₂, when N_LC > 1, lie outside the studied energy transfer interval, even though at the moment we cannot clearly locate them. This fact is not surprising, since in the extremely crude model of a quantum particle in a box, one observes that energy eigenvalues are proportional to V^−2/3, where V is the box volume available to the particle. Now, if V is halved because N_LC goes from 1 to 2 and then is even further reduced by the existence of an H₂–H₂ short-distance repulsive potential, it is reasonable to assume that the corresponding ℏω values of our INS spectral features get more than doubled. (B) During the evacuation process, the number of singly occupied LCs grows, and, so, the probability of a scattering event occurring from an H₂ trapped in a singly occupied LC becomes also larger. This provides a simple explanation of what is visible in Fig. 6 concerning the step-by-step H₂ removal from sample no. 4. As a matter of fact, we have also considered the possibility of an alternative explanation for the band growth in the 1 meV $<ℏω<$ 3 meV-interval. Hypothetically, one might wonder if the heating and cooling cycles operated on the hydrogen-filled clathrate hydrate sample during the evacuation process could increase the ortho–para H₂ ratio, inducing the mentioned band growth. For this reason, after each evacuation step, we have monitored the j₁ = 1 → j₂ = 0 rotational transition (which is proportional to the ortho-H₂ content) present in the anti-Stokes part of the INS spectra, by integrating their signal in the following zone of the kinematic plane: −16 meV $<ℏω<$ −11 meV and 1.6 Å⁻¹ < Q < 3.28 Å⁻¹. In this way, we have obtained five measurements that exhibit a trend that is opposite to that observed for the 1–3 meV spectral intensity. In other words, while the j₁ = 1 → j₂ = 0 rotational line decreases in parallel with the H₂ evacuation process, the 1–3 meV spectral intensity counterintuitively increases.

A more quantitative study of the INS data plotted in Fig. 8 can be carried out by establishing a comparison between them and the calculated form factor, I_0,0;0,1(Q), reported in Eq. (A14). However, given the spherical symmetry of the potential used for the evaluation of the latter quantity, it is appropriate to take the sum [dubbed I_sum(Q)] of the three curves included in the former dataset, since their splitting, as we have already seen, is caused by the aspherical component of the real potential felt by an H₂ in an LC. As for I_0,0;0,1(Q) in Eq. (A14), this CM form factor has been complemented with the appropriate rotational terms, a²(Q, 1, m₁, 1, m₂), which, as shown in Table II, simply reduce to $j02Qr02+2j22Qr02$ in the case of a sum over m₁ and m₂. As explained above, in our Q range, such a rotational term is largely dominated by $j02Qr02$⁠, which actually coincides with the spectral modulation caused by the rotation of a classical molecular dumbbell. The product between the CM and the rotational term, named I_CM+R(Q), is plotted in Fig. 9 together with the two “toy models” presented in the  Appendix, namely H₂ inside a spherical box (SB) [see Eq. (A17)] and the three-dimensional isotropic harmonic (IH) oscillator [see Eq. (A19)]. It is worth noting that these two analytical calculations have been carried out (see the  Appendix for details), imposing an energy jump ΔE_0,0;0,1 between the respective first two levels equal to the exact experimental value (namely 1.92 meV), rather than the calculated value derived from the eigenvalues of isotropized potential by Alavi et al.³⁸ (namely, 2.21 meV). Even a simple visual inspection of Fig. 9 is sufficient to verify that I_sum(Q) is well reproduced by the form factor I_CM+R(Q) obtained from the exact solutions of the Schrödinger equation (SE) for the aforementioned isotropized potential, while the two analytical models are clearly unable to predict the Q value corresponding to the maximum of the INS dataset, providing only a qualitative description of I_sum(Q). In this respect, the SB model seems slightly more performing than the IH one.

This result might appear as slightly surprising, since one would have expected IH (parabolic potential well) and SB (infinite potential well) to be the two extreme cases with the anharmonic SE lying in between the two. For this reason, it might be of some use to pay attention to the radial density distributions of the initial, (N₁ = 0, L₁ = 0), and final, (N₂ = 0, L₂ = 1), states for the three aforementioned potential schemes. Using the symbols of the  Appendix, one can easily verify that these physical quantities would amount to $s2R00(s)2$ and $s2R01(s)2$⁠, respectively, where s is the distance between the H₂ CM and the LC center. However, since INS, differently from neutron Compton scattering,⁴² does not probe the state itself, but rather the transition (N₁ = 0, L₁ = 0) → (N₂ = 0, L₂ = 1), it is worthwhile to focus on the overlaps between these two states, namely $s2R00*(s)R01(s)$⁠. They are shown in Fig. 10 for IH, SB, and SE. One can see that the IH and SB overlaps appear similar to each other, with the exception of the region where the radial value is larger than s₀ = 2.36 Å, which is forbidden for the SB wavefunctions, but can be accessed by the tail of the IH ones. The reason for such a difference is due to the fact that the SB model, differently from the IH one, is evaluated assuming an infinitely steep potential wall located at s₀ (see the inset of Fig. 10). Here, it is important to note that the s₀ radial value is fixed once the energy jump ΔE_0,0;0,1 is determined, and as shown in Fig. 10, in our case, s₀ turns out to be almost exactly one half of the crystallographic mean value of the LC radius (i.e., 4.73 Ų²). On the contrary, IH is based on a broad parabolic potential (see again the inset of Fig. 10), of which the curvature cannot be adjusted since it determines ΔE_0,0;0,1, giving rise to a rather unrealistic description of the H₂ probability distribution inside the LC. Finally, as for the SE overlap, this quantity is more concentrated in the low-s range and can hardly exceed radial values of about 1.8–1.9 Å. This figure is also very sensible since it turns out to be close enough to 1.775 Å, which is simply the difference between the mean crystallographic LC radius (i.e., 4.73 Å) and the sum of half of the Lennard-Jones particle size for H₂ plus the van der Waals radius for O (namely, σ[H₂] = 2.87 Å⁴³ and r_vw[O] = 1.52 Å,⁴⁴ respectively). This fact proves that SE is by no means an intermediate anharmonic scenario in between a purely harmonic well and an infinitely steep potential wall, but it must be considered as a separate case as clearly shown in the inset of Fig. 10. As a final point, we want to remind that the SB model is actually able to provide a fairly good description of the present intensities of the three low-frequency INS peaks, if ΔE_0,0;0,1 is artificially increased to 3.08–3.09 meV, so that s₀ gets shrunk to about 1.86 Å. However, such an energy jump would be severely incompatible with its experimental value, $ēav=1.92(3)$ meV, reported above and, for this reason, must be rejected.

In the present work, we have experimentally investigated both the structure and the dynamics of a hydrogen-containing sII clathrate hydrate during the evacuation of its guest gas molecules, reaching a mean occupation level near unity (1.0) for the LCs. In particular, we have shown the results of two neutron scattering experiments on such a system: a diffraction measurement on the structure of such a compound and a spectroscopic study on the microscopic dynamics of the H₂ molecules trapped inside the clathrate LCs, with both investigations carried out by varying the gas content in the sample. In the former experiment, by means of a Rietveld refinement of the neutron diffraction patterns, we have been able to estimate the mean occupation level of the small and large cages by guest H₂/D₂ molecules as a function of time. We have also shown that the water host lattice is slightly expanded for a large hydrogen content. In the latter experiment, we have demonstrated that for a small hydrogen content, since in this case the probability of single occupation of the LCs becomes larger, three spectral features can be observed at low energy transfer (namely, around 1–3 meV as in Ref. 23) which can be assigned to the rattling motion of a single H₂ molecule in such a cage type. This model simply explains a remarkable experimental finding, i.e., the fact that increasing the hydrogen filling, these three INS bands gradually lose their intensity and finally completely disappear from the measured spectra.

Focusing on the case of single occupation of the LCs, the spectral intensity of the mentioned H₂ rattling bands has been studied as a function of Q and compared with three distinct quantum models, all of them being able to reproduce the correct energy jump obtained averaging the energy transfer values pertaining to these bands, namely (i) the exact solution of the Schrödinger equation for a semiempirical force field, (ii) a particle trapped inside an infinitely rigid sphere, and (iii) an isotropic three-dimensional harmonic oscillator. The first model provides good agreement between calculations and experimental INS data, while the last two, which are quite simple approaches, reproduce only qualitative trend of the band intensities with Q. Subsequently, the radial wavefunctions of the three aforementioned models, as well as their potential surfaces, have been presented and discussed, shedding some light on the reason why models (ii) and (iii) fail to accurately describe our INS results.

Finally, one can conclude that the present study has proved that, once a hydrogen-containing clathrate hydrate is fully characterized via ND, INS represents a sensitive probe for the H₂ microscopic dynamics in the LCs, similarly to what has been already found for the SCs.¹⁸ In particular, if significant single H₂ occupation of the LCs is obtained in the sample, intense inelastic bands in the measured spectra are due to such a dynamics. This can be interpreted as the quantum rattling motion of an individual particle in an isotropic potential well exhibiting a radial “trough,” which is 7.2 meV deep and is placed at about 1.3 Å from the cage center. However, we have also to admit that a complete description of the INS spectral band splitting cannot be derived in such an isotropic framework, and further work on the aspherical potential components would be needed in order to fully clarify this spectroscopic scenario.

We acknowledge the Institut Laue-Langevin for provision of beam time at the instruments D20 and IN5. We thank Werner F. Kuhs, Andrzej Falenty, and Dirk Wallacher for their help during the sample preparation at Helmholtz-Zentrum Berlin. Some of the authors (L.d.R., M.C., and L.U.) are grateful to Mr. Andrea Donati (IFAC) for his fundamental technical support devoted to the setup of the various high-pressure apparatuses. U.R. thanks Mikhail Kuzovnikov for the useful discussion. L.E.B. acknowledges the financial support by the European Union - NextGenerationEU (PRIN N. 2022NRBLPT), the ANR-23-CE30-0034 EXOTIC-ICE, and the Swiss National Fund (FNS) under Grant No. 212889.

The authors have no conflicts to disclose.

Umbertoluca Ranieri: Data curation (lead); Formal analysis (equal); Investigation (equal); Writing – original draft (equal); Writing – review & editing (equal). Leonardo del Rosso: Data curation (equal); Formal analysis (equal); Writing – original draft (equal); Writing – review & editing (equal). Livia Eleonora Bove: Funding acquisition (equal); Investigation (equal); Project administration (lead); Resources (equal); Writing – original draft (equal). Milva Celli: Data curation (equal); Formal analysis (equal); Software (equal); Writing – review & editing (equal). Daniele Colognesi: Formal analysis (equal); Writing – original draft (lead); Writing – review & editing (equal). Richard Gaal: Data curation (equal); Formal analysis (equal); Supervision (equal); Writing – review & editing (equal). Thomas C. Hansen: Data curation (equal); Formal analysis (equal); Investigation (equal); Supervision (equal). Michael Marek Koza: Data curation (equal); Formal analysis (equal); Investigation (equal); Writing – review & editing (equal). Lorenzo Ulivi: Data curation (equal); Formal analysis (equal); Writing – original draft (equal).

Raw data were generated at the ILL (Grenoble, France) large scale facility and will be publicly available from the ILL depository at DOI: 10.5291/ILL-DATA.6-07-42 after an embargo time. Derived data supporting the findings of this study are available from the corresponding author upon reasonable request.