Hydrogen clathrate hydrates are ice-like crystalline substances in which hydrogen molecules are trapped inside polyhedral cages formed by the water molecules. Small cages can host only a single H_{2} molecule, while each large cage can be occupied by up to four H_{2} molecules. Here, we present a neutron scattering study on the structure of the *sII* hydrogen clathrate hydrate and on the low-temperature dynamics of the hydrogen molecules trapped in its large cages, as a function of the gas content in the samples. We observe spectral features at low energy transfer (between 1 and 3 meV), and we show that they can be successfully assigned to the rattling motion of a single hydrogen molecule occupying a large water cage. These inelastic bands remarkably lose their intensity with increasing the hydrogen filling, consistently with the fact that the probability of single occupation (as opposed to multiple occupation) increases as the hydrogen content in the sample gets lower. The spectral intensity of the H_{2} rattling bands is studied as a function of the momentum transfer for partially emptied samples and compared with three distinct quantum models for a single H_{2} molecule in a large cage: (i) the exact solution of the Schrödinger equation for a well-assessed semiempirical force field, (ii) a particle trapped in a rigid sphere, and (iii) an isotropic three-dimensional harmonic oscillator. The first model provides good agreement between calculations and experimental data, while the last two only reproduce their qualitative trend. Finally, the radial wavefunctions of the three aforementioned models, as well as their potential surfaces, are presented and discussed.

## I. INTRODUCTION

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.^{1} The solid H_{2}O–H_{2} 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.^{10} In this structure, the H_{2}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_{2}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.^{10} 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 $\u2248100$ MPa, with the H_{2} 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_{2} 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_{2} 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_{2} in the LCs has been probed by quasi-elastic neutron scattering^{22} and INS^{23} 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.^{23} 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_{2} dynamics and showing some counterintuitive results, such as an increase in intensity of the mentioned rattling bands as the average H_{2} 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^{23} 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_{2} 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.

## II. EXPERIMENTAL METHODS

As illustrated in Table I, the samples investigated in the present study, i.e., *sII* clathrate hydrates containing H_{2} and D_{2}, 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_{2}O ice I*h* powder, previously prepared by finely grinding a block of ice in an inert N_{2} 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_{2}O molecules in the spectral range 180–310 cm^{−1} together with strong and well distinct bands originated from the rotational transitions *S*_{0}(*j*) (i.e., with *ν* = 0 and Δ*j* = 2) of D_{2} (three bands) and H_{2} (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_{2} and H_{2} (right panels of Fig. 2), one can observe the roto-vibrational transitions *Q*_{1}(*j*) (i.e., with *ν* = 1 and Δ*j* = 0) of D_{2} and H_{2} 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 I*h* spheres, obtained by spraying ultra-pure D_{2}O water into a liquid N_{2} bath within a N_{2}-filled glovebox to prevent moisture condensation, were exposed to an H_{2} 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_{2} or D_{2} molecules on average in their LCs,^{23} while for the samples prepared at IFAC, a mean occupation level of the LCs close to 3 is systematically obtained in both H_{2} and D_{2} 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_{2} 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).

Sample no. . | Host . | Guest . | Batch no. . | Protocol . | Technique . | Instrument . | λ (Å)
. | Raman N_{LC} (T = 77 K)
. | ND N_{LC} (T = 130 K)
. |
---|---|---|---|---|---|---|---|---|---|

1 | D_{2}O | D_{2} | I | IFAC | ND | D20 | 2.419 | 2.96 | 2.30 |

2 | D_{2}O | H_{2} | II | HZB | ND | D20 | 2.419 | ⋯ | 1.25 |

3 | D_{2}O | H_{2} | III | IFAC | ND | D20 | 2.419 | 2.91 | 1.90 |

4 | D_{2}O | H_{2} | III | IFAC | INS | IN5 | 3.50 | 2.91 | ⋯ |

5 | D_{2}O | H_{2} | II | HZB | INS | IN5 | 3.50 | ⋯ | ⋯ |

Sample no. . | Host . | Guest . | Batch no. . | Protocol . | Technique . | Instrument . | λ (Å)
. | Raman N_{LC} (T = 77 K)
. | ND N_{LC} (T = 130 K)
. |
---|---|---|---|---|---|---|---|---|---|

1 | D_{2}O | D_{2} | I | IFAC | ND | D20 | 2.419 | 2.96 | 2.30 |

2 | D_{2}O | H_{2} | II | HZB | ND | D20 | 2.419 | ⋯ | 1.25 |

3 | D_{2}O | H_{2} | III | IFAC | ND | D20 | 2.419 | 2.91 | 1.90 |

4 | D_{2}O | H_{2} | III | IFAC | INS | IN5 | 3.50 | 2.91 | ⋯ |

5 | D_{2}O | H_{2} | II | HZB | INS | IN5 | 3.50 | ⋯ | ⋯ |

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 I*h* and ice I*c* 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 I*c*. 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.^{24}. 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_{2}Ca_{3}Al_{2}F_{14} 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_{2} (which requires 8 h as a maximum to get a mean *N*_{LC} ≤ 1.0) with respect to D_{2}, we chose to focus our INS experiments only on the H_{2}–D_{2}O clathrate hydrate system, namely on sample nos. 4 and 5 in Table I. The INS measurements were performed on the IN5 instrument^{25} 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 Å^{−1}. Making use of a chopper velocity *v*_{0} = 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_{2} guests in an LC) and at different stages of the process of emptying the LCs from the guest H_{2} 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_{2} 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 Å^{−1}, with a selected *Q*-spacing equal to 0.2 Å^{−1}. 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.

## III. NEUTRON DIFFRACTION RESULTS

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_{2} or D_{2} guest molecules. As a function of time, for the D_{2}–D_{2}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.^{26} As for the H_{2}–D_{2}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_{2} 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_{2} molecules with respect to D_{2}. 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\u0304$*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_{2} or D_{2} 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_{2} or D_{2} 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

_{2}molecules leave the LCs slightly faster than D_{2}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,^{32}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

_{2}-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

_{2}and D_{2}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_{2}.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

_{2}and H_{2}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.^{24}This fact might be explained by the hindered diffusion of a guest molecule through the (small) pentagonal faces of an SC.^{33}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.

^{24}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_{2}- and D_{2}-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_{2} or D_{2} 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_{2}. 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_{2} 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.^{11} 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_{2} clathrate hydrate for a high occupation of the SCs. An expansion of the lattice was also reported for the phase-*C*_{0} hydrogen hydrate with increasing H_{2} filling^{36} and similarly for the isostructural Ne- and O_{2}-filled hydrates.^{37} Finally, the empty *sII* clathrate remarkably exhibits a lattice constant larger than those in the N_{2} 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^{24}). 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,^{34} 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.

## IV. INELASTIC NEUTRON SCATTERING RESULTS

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_{2} 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_{2}-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 $<\u210f\omega <0.95$ and 3.2 $<\u210f\omega <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 $(e\u03041,2,3)$ and widths $(w\u03041,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 $e\u03041=1.305(7)$ meV, $e\u03042=1.83(2)$ meV, and $e\u03043=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_{2} molecule in the LC, and explain the counterintuitive trend of the spectral intensity, visible in Fig. 6, which grows as the H_{2} evacuation process is continued. These three bands correspond to the fundamental rattling transition of the H_{2} CM, namely that described by the following quantum numbers: (*N*_{1} = 0, *L*_{1} = 0) → (*N*_{2} = 0, *L*_{2} = 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_{2}-CM states. The potential energy for a single H_{2} 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_{2}–D_{2}O pairs (assumed to be given by the standard hydrogen-water potential)^{38} over 514 D_{2}O molecules placed around the cage center in a rigid lattice and subsequently averaged over the H_{2} 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.^{39} A further average over the direction of the H_{2}-CM coordinate, $s\u0302$ (see the Appendix for details), provided the isotropic part of the potential. The solution of the radial Schrödinger equation for an H_{2} molecule sitting in such an isotropic potential well determines the first quantum transition, that is, (*N*_{1} = 0, *L*_{1} = 0) → (*N*_{2} = 0, *L*_{2} = 1), characterized by an energy jump equal to 2.21 meV, which compares reasonably well with the observed $e\u03041,2,3$, whose average value is $e\u0304av=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_{2} 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_{2}–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_{2} (described by *j* and *m*) have to be properly taken into account. Now, in our experiment, for neutron cross section reasons,^{40} the measured signal is dominated by the transition *j*_{1} = 1 → *j*_{2} = 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_{2}–LC interaction and the H_{2} 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} = 1 → *j*_{2} = 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*_{1} = *m*_{2} 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*^{2}(*Q*, *j*_{1}, *m*_{1}, *j*_{2}, *m*_{2}) 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} = 1, *m*_{1} = −1) → (*j*_{2} = 1, *m*_{2} = −1), (*j*_{1} = 1, *m*_{1} = 0) → (*j*_{2} = 1, *m*_{2} = 0), and (*j*_{1} = 1, *m*_{1} = 1) → (*j*_{2} = 1, *m*_{2} = 1)] provide sizable contributions. The mathematical reason for this is given by the interatomic distance in H_{2}, *r*_{0} ≃ 0.7414 Å, so that one obtains $j02Qr02\u226bj22Qr02$ for *Q* ≤ 4 Å^{−1}, which implies that in the present INS experiment, all the transitions with *m*_{1} ≠ *m*_{2} 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.

. | m_{2}
. | ||
---|---|---|---|

m_{1} | −1 | 0 | 1 |

−1 | $13j02Qr02+115j22Qr02$ | $15j22Qr02$ | $25j22Qr02$ |

0 | $15j22Qr02$ | $13j02Qr02+415j22Qr02$ | $15j22Qr02$ |

1 | $25j22Qr02$ | $15j22Qr02$ | $13j02Qr02+115j22Qr02$ |

. | m_{2}
. | ||
---|---|---|---|

m_{1} | −1 | 0 | 1 |

−1 | $13j02Qr02+115j22Qr02$ | $15j22Qr02$ | $25j22Qr02$ |

0 | $15j22Qr02$ | $13j02Qr02+415j22Qr02$ | $15j22Qr02$ |

1 | $25j22Qr02$ | $15j22Qr02$ | $13j02Qr02+115j22Qr02$ |

As for the other issue, the increase in the band intensity during the progress of the H_{2} evacuation in sample no. 4, as a matter of fact, the detailed excited states of two (or more) *ortho*-H_{2} molecules in a LC are still unknown, since it is already quite demanding the calculation of the ground state of such a system,^{41} where at least ten degrees of freedom have to be tackled. However, noticing that the increase in the mentioned bands in the 1 $<\u210f\omega <$ 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_{2}, 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_{2}–H_{2} 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_{2} 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_{2} 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 $<\u210f\omega <$ 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_{2} ratio, inducing the mentioned band growth. For this reason, after each evacuation step, we have monitored the *j*_{1} = 1 → *j*_{2} = 0 rotational transition (which is proportional to the *ortho*-H_{2} 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 $<\u210f\omega <$ −11 meV and 1.6 Å^{−1} < *Q* < 3.28 Å^{−1}. 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} = 1 → *j*_{2} = 0 rotational line decreases in parallel with the H_{2} 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_{2} 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*^{2}(*Q*, 1, *m*_{1}, 1, *m*_{2}), which, as shown in Table II, simply reduce to $j02Qr02+2j22Qr02$ in the case of a sum over *m*_{1} and *m*_{2}. 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_{2} 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.*^{38} (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*_{1} = 0, *L*_{1} = 0), and final, (*N*_{2} = 0, *L*_{2} = 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_{2} CM and the LC center. However, since INS, differently from neutron Compton scattering,^{42} does not probe the state itself, but rather the transition (*N*_{1} = 0, *L*_{1} = 0) → (*N*_{2} = 0, *L*_{2} = 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*_{0} = 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*_{0} (see the inset of Fig. 10). Here, it is important to note that the *s*_{0} 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*_{0} turns out to be almost exactly one half of the crystallographic mean value of the LC radius (i.e., 4.73 Å^{22}). 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_{2} 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_{2} plus the van der Waals radius for O (namely, *σ*[H_{2}] = 2.87 Å^{43} and *r*_{vw}[O] = 1.52 Å,^{44} 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*_{0} gets shrunk to about 1.86 Å. However, such an energy jump would be severely incompatible with its experimental value, $e\u0304av=1.92(3)$ meV, reported above and, for this reason, must be rejected.

## V. CONCLUSIONS

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_{2} 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_{2}/D_{2} 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_{2} 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_{2} 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_{2} microscopic dynamics in the LCs, similarly to what has been already found for the SCs.^{18} In particular, if significant single H_{2} 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.

## ACKNOWLEDGMENTS

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.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**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).

## DATA AVAILABILITY

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.

### APPENDIX: THEORETICAL BACKGROUND AND RELEVANT FORMULAS

^{14}due to the fact that the incoherent neutron scattering length for a proton is by far larger than the coherent one. A second customary approximation that we will consider from the beginning, in order to arrive at approximate formulas for the description of the H

_{2}molecular motion, is the hypothesis of the decoupling of the rotations of the molecule from its CM motion (i.e., the so-called rattling) so that the double differential neutron scattering cross section can be written as

*k*

_{i}and

*k*

_{f}denote the initial and final neutron wavevectors, respectively, and the symbol ⊗ represents a convolution product. The Dirac

*δ*functions are centered at the energies of the rotational transitions

*ℏω*

_{12}of the H

_{2}molecule, and the rotational form factor

*a*

^{2}(

*Q*,

*j*

_{1},

*j*

_{2}) is a function of the momentum transfer

*Q*depending on the rotational transition

*j*

_{1}→

*j*

_{2}of the molecule. The dependence on the

*m*

_{1}and

*m*

_{2}quantum numbers will be considered in the following. The scattering law

*S*

_{self}(

*Q*,

*ω*) is the self-part of the dynamic structure factor for the CM motion.

^{40,45}the energy levels can be identified by the same rotational quantum numbers

*j*and

*m*of the free molecule, plus the translational quantum numbers

*N*,

*L*, and

*M*for the CM rattling motion in the cage, with

*N*= 0, 1, …,

*L*= 0, 1, …, and −

*L*≤

*M*≤

*L*. The general expression of the rotational form factor

*a*

^{2}(

*Q*,

*j*

_{1},

*j*

_{2}) is given in Ref. 40. In principle, the rotational energy levels and wavefunctions of a single, rigid H

_{2}molecule can be calculated [e.g., via density functional theory (DFT)] knowing the interaction of the H

_{2}molecule with the crystal field of the water lattice, which alters the free molecule rotational states, lifting the degeneracy of the

*m*sublevels, even without considering the effect of the other H

_{2}molecules in the nearby cages. It is well known that for a freely rotating molecule, the eigenfunctions of the rotational states are

*spherical harmonics*$Yjm(u\u0302)$, where $u\u0302$ represents the molecular orientation, while the energy of these eigenstates depends only on

*j*and is given, for a rigid molecule, simply by

*E*(

*j*) =

*B*

_{0}

*j*(

*j*+ 1). However, in the presence of an interaction with the cage, the energy levels of this guest molecule do depend on the value of

*m*. In view of the considerations above, we have to take into account the

*m*dependence of the molecular rotational states, so one needs to write the double differential neutron cross section reported in Eq. (A1) in the following, slightly modified, form:

*a*

^{2}(

*Q*,

*j*

_{1},

*m*

_{1},

*j*

_{2},

*m*

_{2}) can be written as

*l*th order, the quantity

*r*

_{0}representing the internuclear distance in the H

_{2}molecule, and

*σ*(

*j*

_{1},

*j*

_{2}) depending only on the parity of

*j*

_{1}and

*j*

_{2}, which will be indicated by the letters “e” or “o” (for even or odd parity, respectively). Written in a more explicit way, we have

^{45}that

*σ*

_{coh}and

*σ*

_{inc}stand for the neutron coherent and incoherent cross section for H, respectively.

*V*(

*s*), where

*s*is the distance of the particle CM from the cage center. In this approximation, the angular part of the wavefunction describing the CM motion can be factorized so that one can write

*Y*

_{LM}(

*θ*,

*ϕ*) must not be confused with the rotational (or orientational) spherical harmonic $Yj,m(u\u0302)$, while

*R*

_{NL}(

*s*) is the radial part of the CM wavefunction, which can be calculated either analytically for few model cases (e.g., three-dimensional harmonic oscillator or infinite spherical well) or numerically, once

*V*(

*s*) is known. It is worth noticing that

*R*

_{NL}(

*s*) complies with the standard normalization condition: $\u222b0\u221e|RNL(s)|2s2ds=1$.

*ψ*

_{1}(with quantum numbers

*N*

_{1},

*L*

_{1},

*M*

_{1}) to a final state

*ψ*

_{2}(with quantum numbers

*N*

_{2},

*L*

_{2},

*M*

_{2}) is proportional to

**orientation, $Q\u0302$, due to the polycrystalline nature of our clathrate hydrate samples. Equation (A9) can be simplified making use of the well-known spherical harmonics expansion of a plane wave, namely,**

*Q***joining the H**

*s*_{2}CM with the cage center and, as we have seen above, is expressed by

*ϕ*(azimuth) and

*θ*(elevation). For practical calculations, one uses the well-known expression for the integral of three spherical harmonics, namely,

*L*(and different

*M*) are degenerate, so the calculation proceeds performing the average over the initial states labeled by

*M*

_{1}and the sum over the final states labeled by

*M*

_{2}. In this way, the result for the transition from the initial state identified by

*N*

_{1},

*L*

_{1}to the final state identified by

*N*

_{2},

*L*

_{2}is simply given by

*N*

_{1}= 0 and

*L*

_{1}= 0, so that the 3

*j*-symbol is zero unless

*L*=

*L*

_{2}, and Eq. (A13) reads

_{2}molecule, one should consider that the intensity of a specific roto-rattling transition also contains the additional form factor

*a*

^{2}(

*Q*,

*j*

_{1},

*m*

_{1},

*j*

_{2},

*m*

_{2}), which takes into account the rotational transition, as clearly shown in Eq. (A3). In the practical case related to the interpretation of the present INS data, one can select

*N*

_{2}= 0,

*L*

_{2}= 1.

*R*

_{00}(

*s*) and

*R*

_{01}(

*s*) to be the eigenstates of a particle confined inside an infinitely rigid sphere of radius

*s*

_{0}, namely,

*κ*

_{0}=

*π*/

*s*

_{0}and

*κ*

_{1}≈ 4.493 41/

*s*

_{0}.

^{46}In this simple model, the energy jump, Δ

*E*

_{0,0;0,1}, for the (

*N*

_{1}= 0,

*L*

_{1}= 0) → (

*N*

_{2}= 0,

*L*

_{2}= 1) transition is given by the difference of the two kinetic energies,

*s*

_{0}. So, one can write an analytic expression for the approximate

*I*

_{0,0;0,1}(

*Q*),

*ω*

_{0}),

^{46}where now

*α*being directly related to the harmonic oscillator energy jump, Δ

*E*

_{0,0;0,1}=

*ℏω*

_{0}, for the (

*N*

_{1}= 0,

*L*

_{1}= 0) → (

*N*

_{2}= 0,

*L*

_{2}= 1) transition:

*α*

^{2}= 2

*M*

_{H}

*ω*

_{0}

*ℏ*

^{−1}. Plugging Eq. (A18) into Eq. (A14), one finally obtains

## REFERENCES

*Inclusion Compounds*

*Clathrate Hydrates of Natural Gases*

*Molecular Theory of Gases and Liquids*