Using inelastic neutron scattering, we have measured the VDOS of water at 120 K (solid state) softly confined⁴⁷ between graphene oxide membranes (GOM) with different levels of hydration. The detailed information on the sample preparation, the experimental methods and the data analysis are provided in the Methods.

As a reference, in Fig. 2a, we present the inelastic neutron scattering data of the bulk crystalline ice measured in this work and those for bulk high-density amorphous ice (HDA) and bulk low-density amorphous ice (LDA) ice taken from ref. ⁴⁸. At low frequency, approximately below 4 meV, we observe a leveling-off of the reduced density of states, g(ω)/ω², and a well-defined Debye level, in agreement with the standard Debye model expectation, g(ω) ~ ω². This is further supported by the VDOS of bulk crystalline and disordered ice samples derived from MD simulations (see next section). Hence, the Debye model is a valid description of the low-energy vibrational dynamics in bulk solids, as expected. In sharp contrast, in the sample with lowest hydration level (0.1 g water per gram GOM), where the distance between the neighboring GOM layers is only 7Å (Fig. 2c), the reduced density of states displays a neat low-frequency linear scaling, thus the original VDOS scales as g(ω) ~ ω³ (Fig. 2b). The same feature is shown in logarithmic scale in the inset of Fig. 2b. This cubic regime appears robust in a large range of frequencies from ≈1.6 meV to ≈5 meV. Below 1.6 meV, the measurements are affected by the resolution function of the instrument (0.35 meV upwards to 1 meV or so ref. ⁴⁹) and thus the corresponding data are not presented. In order to make sure that this anomalous scaling does not stem from the GOM structure itself, we measured the VDOS of the dry GOM without any water in between (shown in Fig. 3a). The resulting reduced VDOS, presented in Fig. 3a, slowly decreases with frequency and does not present any relevant features. Based on the dependence of the reduced VDOS on the frequency ω and the absolute amplitude of the dry sample, one can deduce that the observed anomalous ω³ scaling in hydrated GOM must result from the nano-confined water. Further X-ray diffraction measurements reveal that the water confined in GOM at this hydration level is amorphously packed as no characteristic Bragg peak of the crystalline ice is present (see Fig. 3c). At such low temperature, the translational motion of water molecules is expected to be strongly suppressed like in solids. Thus, we name these confined supercooled water molecules without translational freedom as quasi-2D confined amorphous ice.Fig. 2

a The Debye normalized VDOS g(ω)/ω² for high density amorphous (HDA) ice, low density amorphous ice (LDA) obtained from ref. ⁴⁸ and bulk crystalline ice measured in this work at 120K. The horizontal dashed lines indicate the Debye level for each sample. The curves have been manually shifted vertically for better presentation and the vertical axes has arbitrary units. b The experimental normalized VDOS g(ω)/ω² for the hydrated GOM sample at h = 0.1 g water/gram GOM. The dashed line indicates the low frequency ~ω³ scaling due to confinement. The inset shows the original VDOS data in log–log scale. The y-axes are presented in arbitrary units. c The XRD data of GOM sample with different hydration levels. d The scanning electron microscope (SEM) image of the GOM sample which highlights the fuzzy boundaries of the graphene oxide membranes.

To further explore this ω³ scaling found in the experimental data, we measured a series of samples with different levels of hydration. The results of the VDOS measurements on the various samples are shown in Fig. 3a, b. A lower hydration level corresponds to a smaller confinement size in the vertical z direction, where the inter-layer distance inside the membrane can be varied from 7Å to a few nanometers when increasing h from 0.1 to 0.7. We have fitted the experimental data with a single power-law function using as the upper value for the fitting window the point at which the VDOS of the bulk ice sample deviates from the Debye law and the shoulder of the higher energy vibrational peak appears, ≈4 meV. As shown in Fig. 3b–d, the power law of the low-energy VDOS changes gradually from 3 to 2 by varying h = 0.1 to 0.7. Such an outcome can be rationalized using the following logic. When the hydration level is low, the water molecules are mostly confined between the GOM layers and form a disordered solid phase; increasing the hydration level h, part of the water migrates away from being nano-confined between the layers to some large voids or even to the external surface of the GOM to form bulk crystalline ice, namely ice segregation^50,51. This picture is supported by the small-angle X-ray scattering results (see Supplementary Fig. 1). As seen in Fig. 3c, when the hydration is low, (h = 0.1), no Bragg peak is observed. However, the Bragg peaks appear when h is above 0.4. We note that these Bragg peaks are rather narrow. This indicates that they result from bulk crystals instead of nanometer-sized crystals, where they would significantly broaden as seen in refs. ^52,53. This observation is further supported by the DSC measurement, where no first-order transition is observed at h = 0.3, but it appears at higher hydration levels, h = 0.4 or above. Moreover, as evident by SAXS data (see Supplementary Fig. 1), at high hydration levels (e.g. h = 0.7), the GOM interlayer-distance is suddenly reduced when cooling down to 253 K at which the Bragg peaks appear. Hence, at high hydration, e.g. h = 0.4 and h = 0.7, the crystalline ice is not formed confined between the layers, but it may exist within some large voids or on the external surface of the GOM, otherwise the inter-layer distance should increase rather than decrease as observed. Therefore, at any given hydration level, the system contains two separate components, nano-confined amorphous ice and large bulk ice crystals, which contribute with different scaling laws to the total density of states. The relative ratio between the two can be estimated using the DSC data (see Supplementary Note) and the results are presented in Fig. 3d. The nano-confined component obeys the ω³ scaling (assumed to be the same as the one discovered at h = 0.1) while the bulk crystal follows Debye’s law ω². Therefore, by fitting to a single power law, we obtain a smooth interpolation of the power from 3 to 2 by increasing h from the strictly confined case towards the one dominated by the bulk crystal component (see Fig. 3b–d). The vanishing of the Debye contribution coming from the bulk crystal component is even clearer in the reduced VDOS presented in Fig. 3a. This agreement supports our previous argument that the density of states of the confined amorphous ice obeys a power-law ~ ω³. Such a cubic power law has also been recently reported in ref. ⁵⁴ for a 2D model glass system. Despite the tempting similarities, our setup is different with respect to that of ref. ⁵⁴ in various aspects. First, our is a quasi-2D system, rather than a small 2D system, with excitations also in the vertical z direction. Second, as we will show with simulations in the next section, the ω³ power law is not an exclusive property of amorphous systems but it appears also in crystalline ones under strong confinement. All in all, the nature of the scaling discussed in this manuscript appears to be profoundly different with respect to the one reported in ref. ⁵⁴.

Before moving on, let us comment on another interesting outcome of the experimental analysis. The bulk sample at 120 K displays a sharp resonance peak in the reduced VDOS at ≈7 meV (Fig. 3a). Despite the similarities with the boson peak anomaly in amorphous systems, most of the literature links its presence to a specific optical mode of the crystalline ice structure^55–57. By decreasing the level of hydration, and therefore increasing the strength of confinement, we observe a shift of this peak towards lower frequencies and a broadening of its linewidth. This observation suggests that the effects of confinement decrease the lifetime of this mode and renormalize the energy of the optical modes similarly to what disorder or anharmonicity would do. Indeed, one could argue that confinement is itself a source of anharmonicity. Interestingly, the observed red-shift of the spectrum induced by the soft confinement is consistent with the results of ref. ⁴².

In order to confirm the universality of this power-law ~ ω³ for nano-confined solids and to prove that, as expected, it is not a peculiar property of the confined amorphous state, we perform all-atom molecular dynamics (MD) simulations on bulk ice, slab ice, supercooled water and slab supercooled water at 120 K, where the thickness of the slab along the z direction is fixed to 3 nm. Supercooled water at such low temperature could be considered as amorphous ice. The different structure between the hexagonal crystalline ice and the amorphous ice is highlighted in panels (a) and (b) of Fig. 4. The geometry of the slab sample is shown in panel (c) of Fig. 4. The details of the simulations can be found in the Methods and Supplementary Note. The results of the simulations are shown in Fig. 4d. Both the amorphous and crystalline ices in the bulk phase show a very clear leveling-off in the reduced density of states at low frequency, signaling the presence of a strong Debye scaling law, g(ω) ~ ω². The Debye level is reasonably lower in the crystalline sample as compared to the bulk amorphous case because of a much larger value of the speed of sound. In the slab (confined) samples (Methods and Supplementary Note), the VDOS from the simulations shows clear differences in the low-frequency regime as compared to the bulk sample: the Debye leveling-off is absent and a ω³ scaling shows up, consistently with the neutron scattering experimental findings. In summary, the analysis of the MD simulations data confirms the outcomes of the experiments and it also provides additional evidence for the universality of the phenomenon which does attain to any confined solid regardless of the crystalline or amorphous structure.Fig. 4

a The snapshot of the structure of the bulk ice sample. For visual purposes, only a few layers are kept. The zoom emphasizes its hexagonal structure. b The snapshot of the structure of the bulk amorphous ice sample. Also here, only a few layers are kept. c The amorphous ice slab sample with size 180 Å and confinement length along the z direction 30 Å. Here, the original structure is kept. d The normalized VDOS in function of frequency obtained from the MD simulations for slab ice (blue), bulk crystalline ice (green), slab amorphous (red) and bulk amorphous ice (yellow). The horizontal dashed lines indicate the Debye level for the bulk systems. The diagonal dashed lines guide the eyes towards the g(ω) ~ ω³ low-frequency scaling. The y-axes are presented in arbitrary units.

We provide a concise analytic derivation of the ω³ scaling observed for nano-confined solids both in the experimental data and MD simulations. Let us consider a three-dimensional system of linear size L. If the system were confined by atomically smooth and infinitely rigid boundaries (which is not the case for our experimental setup, see Fig. 2d), hard-wall boundary conditions (BCs) would apply, implying a net zero displacement of the systems’ atoms/molecules near the boundary. In turn, the hard-wall BCs would lead to the usual “quantization” of the wavevector of the acoustic (elastic) waves that can propagate in the system: these are standing waves (eigenmodes of the wave equation) with k = nπ/L, with n being an integer (see Supplementary Note for more details).

In our system, however, the smooth hard-wall BCs do not apply because the confining boundaries are graphene oxide sheets, in which the basal plane is highly screwed and possesses different and spatially random-distributed oxide groups (hydroxyl, expoxy, and carboxy groups) (see the scanning electron microscope (SEM) image in Fig. 2d). An additional proof of this fact is provided in Supplementary Fig. 4 using the MD simulations. A crucial consequence of the lack of hard-wall BCs is the fact that the minimum length of wavevector is not k = π/L and that the wavevector is not discrete but continuous. We have verified this statement using the numerical simulations. There (Supplementary Fig. 4 and Supplementary Note), it is found that the minimum wavevector is ≈14πL\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\approx \frac{1}{4}\frac{\pi }{L}$$\end{document}, indeed smaller than π/L. No sign of the discreteness of the spectrum is observed either. Similar conclusions about the fundamental role of the boundary conditions in nano-confined systems have been reached in^35,42.

Going back to the main discussion, the number of states with wavenumber between k and k + dk, i.e. in a spherical shell in k-space are given by1dn=1Vk4πk2dk.\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{{{{\rm{d}}}}}}n=\frac{1}{{V}_{k}}\,4\pi \,{k}^{2}{{{{{\rm{d}}}}}}k.$$\end{document}where V_k = (2π)³/L³ is the k-state volume occupied by a single wavevector. Assuming a linear low-energy dispersion relation for the vibrational modes ω = v k, which as we will see persists even at high level of confinement (see Supplementary Fig. 4 and Supplementary Note), one finally obtains the density of states as2g(ω)=g(k)dkdω=1Vdndkdkdω~ω2,\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$g(\omega )\,=\,g(k)\,\frac{{{{{{\rm{d}}}}}}k }{{{{{{\rm{d}}}}}}\omega}=\frac{1}{{V}}\,\frac{{{{{{\rm{d}}}}}}n}{{{{{{\rm{d}}}}}}k}\,\frac{{{{{{\rm{d}}}}}}k }{{{{{{\rm{d}}}}}}\omega} \sim {\omega }^{2},$$\end{document}where we omit the numerical prefactor since it is undetermined due to the non-smooth BCs (but importantly frequency independent). The ω² scaling in Eq. (3) is the famous Debye’s result^1,3. Here, we have considered a single sound mode and we have neglected the existence of different polarizations as in realistic solids. Given the additive property of the density of states g(ω), one can easily consider this by simply summing the contributions from the different polarizations as3g(ω)≈∑pω2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$g(\omega )\,\approx \,\mathop{\sum}\limits_{p}\,{\omega }^{2}$$\end{document}where p is the index labelling the polarizations. Importantly, the Debye derivation relies on an integration over all directions in the solid angle, reflected in the factor 4π in the above Eq. (1).

Let us now turn to a different situation in which one of the three bulk dimensions is confined but the atoms are still free to vibrate along it. This corresponds exactly to the setup considered in the experimental and simulation parts discussed above. We follow closely the framework presented in⁵⁸, and we shall work without hard-wall BCs consistent with the discussion above. In particular, let us consider a geometry which is rotational invariant in the (x, y) plane but in which the z direction is confined in a range z ∈ [0, L]. We assume the size of the sample in the (x, y) direction to be much larger than that in the z direction, L_z ≪ L_x, L_y. We use spherical coordinates, measuring the polar angle θ from the vertical z-axis. A pictorial representation of the geometry considered is shown in the panel a of Fig. 5. At a fixed angle θ, the maximum wavelength allowed is λmax=L/cosθ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\lambda }_{{{{{{\rm{max}}}}}}}=L/\cos \theta$$\end{document} because of the vertical confinement. Going to momentum space, this translates into a minimum wavevector allowed given by 2πcosθ/L\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$2\pi \cos \theta /L$$\end{document}. This implies that below a certain crossover momentum k_× ≡ 2π/L, the phase space is reduced by the effects of confinement. In particular, two spheres of radius π/L centered at (0, 0, ± π/L) have to be excluded, as shown in panel b of Fig. 5. In the non-confined direction, i.e. for θ = π/2, the macroscopic geometry does not provide any large-distance cutoff and the minimum value of k is simply zero (ignoring the cutoff imposed by the size of the sample). Importantly, this minimum condition does not arise from the hard-wall BCs but simply from the geometry of the slab sample. The absence of hard-wall BCs is fundamental in allowing wavevectors below k_×. Also, the effects of confinement are relevant only for wavevectors below the crossover scale k_×. At higher wavevectors, the phase space is not affected since no minimum value nor angular dependence appears (see panel b of Fig. 5). This point will result important in the following.Fig. 5

a The confined geometry considered. The system size in the (x, y) directions is much larger than that in the z direction, L_z ≪ L_x, L_y. The system is assumed to be rotational invariant in the 2-dimensional (x, y) subspace. The angle θ is defined with respect to the vertical z axes. b The wavevector phase space corresponding to the confined geometry of panel (a). c The normalized VDOS of slab amorphous ice with different confinement sizes L (from the bulk sample in yellow color to the strongly confined sample in dark blue color). The black cross symbols show the crossover frequency ω_× defined in the main text. d The crossover scale ω_× as a function of the inverse confinement length 2π/L. The dashed line is a linear fit whose result is reported in Eq. (7). The inset shows the dispersion curve of the transverse acoustic modes in amorphous ice and the fitting curve giving v_TA ≈ 1300 m s⁻¹.

As a result of this geometric constraint on k-space, below the crossover scale k_×, the integration over the polar angles, which gave the 4π factor in Eq. (1), is now replaced by42∫cos−1Lk/2ππ/2sinθdθ∫02πdϕ=2Lk\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$2\int\nolimits_{{\cos }^{-1}\left(L\,k/2\pi \right)}^{\pi /2}\sin \theta {{{{{\rm{d}}}}}}\theta \int\nolimits_{0}^{2\pi }{{{{{\rm{d}}}}}}\phi \,=\,2\,L\,k$$\end{document}where the complete integral in the ϕ angle reflects the SO(2) rotational invariance in the (x, y) plane. We note that the expression in Eq. (4) reduces to the standard Debye result when the wavevector approaches the crossover value k_×. In particular, in that limit, the lower limit of integration tends to zero. Hence, we obtain that, below a certain threasold k = k_×, the number of states with wavenumber in the range ∈ [k, k + dk] is given by5dn~Lk3dk.\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${{{{{\rm{d}}}}}}n \sim \,L\,{k}^{3}{{{{{\rm{d}}}}}}k.$$\end{document}which clearly deviates from the Debye result. In Eq. (5), the numerical pre-factors have been neglected since they are not relevant for the present analysis. Following the same steps as before, one can deduce the VDOS at low frequency as6g(ω)~ω3\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$g(\omega )\, \sim {\omega }^{3}$$\end{document}which is our main result of this section.

Equation (6) predicts a fundamental deviation from the Debye ~ ω² law which arises because of the geometric constraints on the momentum phase space induced by spatial confinement in real space. This analysis is in agreement with the experimental and simulations outcomes of the previous sections and it is able to explain with a simple argument the ω³ universal scaling in solid systems under slab-confinement. Importantly, the same results would not be obtained for geometries with confinement along two spatial directions (cylinder) or three spatial directions (sphere) where, at least at low frequency, Debye’s law is expected to work, as also confirmed in refs. ^{35,37–39,41,43}.

Interestingly, our theoretical framework is also able to predict the crossover scale above which the Debye scaling re-appears. The latter is indeed given by k_× = 2π/L, or alternatively, making use of the low-energy dispersion relation, by ω_× ≡ 2πv/L, with v the characteristic sound speed of the material. Because of the simplicity of the theoretical model, the crossover appears as a sharp transition at which the density of states is continuous but its derivative is not. Obviously, in more realistic situations, we do expect the crossover to be continuous and smoothed out by various effects including thermal fluctuations. In order to test this prediction of the theory, we have performed additional molecular dynamic simulations by dialing the size of the confined direction z from very large (corresponding to a bulk system) to very small (corresponding to a nano-confined quasi-2D system). The results are presented in panel c of Fig. 5 by using a Debye reduced representation g(ω)/ω². In the bulk sample (yellow markers), a clear Debye level-off is visible at low frequencies. By decreasing the size of the confined region L, we observe the appearance of a low-frequency ω³ scaling as reported in the previous sections both in experiments and simulations. Importantly, using the data from simulations we are able to track the frequency at which the scaling of the density of state changes from the cubic scaling to the more standard Debye one (indicated with a × in panel c of Fig. 5). We then plot the position of this crossover scale as a function of the inverse confinement length 2π/L. The data show a clear linear behaviour which is consistent with7ω×=vfit2π/Lwithvfit≈1500ms−1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${\omega }_{\times }={v}_{{{{{{\rm{fit}}}}}}}\,2\pi /L\quad \,{{\mbox{with}}}\,\quad {v}_{{{{{{\rm{fit}}}}}}}\approx 1500\,{{{{{{{\rm{m}}}}}}}}\,{{{{{{{{\rm{s}}}}}}}}}^{-1}$$\end{document}The fitted value for the sound velocity is in good quantitative agreement with the one extracted directly from the dispersion relation, v_TA ≈ 1300 m s⁻¹ within a 14% error(see the inset of Fig. 5d). Here, the reference sound velocity (inset of Fig. 5d) is the transverse mode. If one used the average speed of transverse and longitudinal modes, 3/v¯3=2/vTA3+1/vLA3\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$3/{\bar{v}}^{3}=2/{v}_{{{{{{\rm{TA}}}}}}}^{3}+1/{v}_{{{{{{\rm{LA}}}}}}}^{3}$$\end{document}, with v_LA ≈ 3900 m s⁻¹, an even better agreement with the value derived from Eq. (7) will be obtained, as v¯≈1480ms−1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\bar{v}\approx 1480\,{{{{{{{\rm{m}}}}}}}}\,{{{{{{{{\rm{s}}}}}}}}}^{-1}$$\end{document}.

In summary, the results from simulations confirm our theoretical framework and also prove that, despite the simplicity of the model, even its quantitative predictions are accurate.

The GOM sample was synthesized using the modified Hummers’ method⁶². The GOM sample was first dehydrated by heating it from room temperature to 313 K and then annealed at this temperature for 12 h in a vacuum to the dry condition. The oxidation rate of the GOM sample is 28%, which is determined by XPS. The dehydrated sample was sealed in a desiccator and exposed to the water vapor to allow water molecules to adsorb to the surface and the interlamination of the GOM sheets. The hydration levels were controlled by adjusting the expose time of the sample and the final values of hydration levels were determined by measuring the weight before and after the water adsorption.

Differential scanning calorimetry (DSC) was used to measure the ratio between bulk ice and confined amorphous ice in GOM sample at low temperature. The DSC results of GOM at different hydration levels were performed by the DSC1 (METTLER TOLEDO). The samples were first annealed at 293 K for 5 min, and then cooled down to 213 K at a cooling rate of 2 K/min to obtain the DSC data (see Supplementary Fig. 1). The DSC data were analyzed by the TA Trios software.

The powder X-ray diffraction data for GOM at different hydration levels were collected using a Rigaku Mini Flex600 X-ray diffractometer, with a Cu Kα source (λ = 1.5406 Å) operated at 40 kV and 15 mA at a scan rate of 10^∘/min from 10^∘ to 60^∘. The PXRD data were analyzed by MDI Jade software.

SAXS characterizations were employed to monitor the interlayer distance evolution in GOM with temperature decreasing and ice freezing. The SAXS measurements were carried out at the BL16B1 beamline of the Shanghai Synchrotron Radiation Facility (SSRF). The wavelength of the X-ray was 1.24 Å. The SAXS patterns were collected by using a Pilatus 2M detector with a resolution of 1475 pixels × 1679 pixels and a pixel size of 172 μm × 172 μm. The data acquisition time for SAXS was set as 10 s for each frame. The sample-to-detector distances of the SAXS is 258 mm, respectively.

SEM images were taken by a MIRA 3 FE-SEM with a 5 kV accelerating voltage.

The dynamic neutron scattering is described in terms of the intermediate scattering function, I(q, t) including incoherent and coherent terms I_inc(q, t) and I_coh(q, t) :8Iinc(q,t)=∑jNbj,inc2<exp[−iq⋅rj(0)]⋅exp[iq⋅rj(t)]>,\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${I}_{{{{{{\rm{inc}}}}}}}(q,t)\,=\,\mathop{\sum }\limits_{j}^{N}\,{b}_{j,{{{{{\rm{inc}}}}}}}^{2} < \exp [-iq\cdot {r}_{j}(0)]\cdot \exp [iq\cdot {r}_{j}(t)] > ,$$\end{document}9Icoh(q,t)=∑jN∑iNbj,cohbi,coh<exp[−iq⋅ri(0)]⋅exp[iq⋅rj(t)]>,\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${I}_{{{{{{\rm{coh}}}}}}}(q,t)\,=\,\mathop{\sum }\limits_{j}^{N}\,\mathop{\sum }\limits_{i}^{N}{b}_{j,{{{{{\rm{coh}}}}}}}\,{b}_{i,{{{{{\rm{coh}}}}}}} < \exp [-iq\cdot {r}_{i}(0)]\cdot \exp [iq\cdot {r}_{j}(t)] > ,$$\end{document}where N is the total number of atoms, b_j,inc and b_j,coh are the incoherent and coherent scattering lengths of a given atom j, r_j is the coordination vector of that atom, the bracket <⋯> denotes an ensemble and orientation average and q is the scattering wavevector. I_inc(q, t) contains the information about self-motions of atoms, and I_coh(q, t) probes mostly interatomic motions. As the incoherent scattering cross section of hydrogen is at least one order of magnitude larger than incoherent and coherent scattering cross sections of other elements, the neutron signals collected on GOM hydrated in H₂O are dominated by incoherent intermediate scattering function, I_inc(q, t) and primarily reflect the self-motion of the water molecules. When t = 0, the coherent intermediate scattering function, I_coh, becomes the static structure factor, I(q), characterizing the atomic structure of the system10I(q)=∑jN∑iNbj,cohbi,coh<exp[−iq⋅ri(0)]⋅exp[iq⋅rj(0)]>.\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$I(q)\,=\,\mathop{\sum }\limits_{j}^{N}\,\mathop{\sum }\limits_{i}^{N}{b}_{j,coh}\,{b}_{i,coh} < \exp [-iq\cdot {r}_{i}(0)]\cdot \exp [iq\cdot {r}_{j}(0)] > .$$\end{document}The measured dynamic structure factor, S(q, ω) corresponds to the time Fourier transform of the intermediate scattering function,11Sinc(q,ω)=∫−∞+∞Iinc(q,t)exp(iωt)dt,\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${S}_{{{{{{\rm{inc}}}}}}}(q,\omega )\,=\,\int\nolimits_{-\infty }^{+\infty }{I}_{{{{{{\rm{inc}}}}}}}(q,t)exp(i\,\omega \,t){{{{{\rm{d}}}}}}t,$$\end{document}12Scoh(q,ω)=∫−∞+∞Icoh(q,t)exp(iωt)dt,\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${S}_{{{{{{\rm{coh}}}}}}}(q,\omega )\,=\,\int\nolimits_{-\infty }^{+\infty }{I}_{{{{{{\rm{coh}}}}}}}(q,t){{{{{\rm{exp}}}}}}(i\,\omega \,t){{{{{\rm{d}}}}}}t,$$\end{document}where ω is the frequency and ℏω = Δ E is the energy transfer between the incident and scattered neutrons. S(q, ω) provides information about the amplitude-weighted distribution of the dynamical modes in the sample with respect to frequency at any given wavevector q.

As the incoherent scattering cross-section of hydrogen is at least 1 order of magnitude larger than incoherent and coherent scattering cross-sections of other elements, the neutron signals are dominated by incoherent scattering function, and primarily reflect the self-motion of the water molecules. The experimental vibrational density of states (DOS) g(ω) can be obtained from the dynamic structure factor, S(q, ω) using the function⁶³:13g(ω)=∫ℏωq2S(q,ω)(1−e−ℏωkBT)dq\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$g(\omega )\,=\,\int \frac{\hslash \omega }{{q}^{2}}S(q,\omega )(1-{{{{{{\rm{e}}}}}}}^{-\frac{\hslash \omega }{{k}_{{{{{{\rm{B}}}}}}}T}}){{{{{\rm{d}}}}}}q$$\end{document}where ℏ is the Planck constant, q is the scattering wavevector, ω is the frequency related to the energy transfer, k_B is the Boltzmann constant, and T is the temperature. The experiments for samples of pure water, GOM absorbed by H₂O with h (gram water/gram GOM) of 0.1, 0.4 and 0.7 were conducted by using a time-of-flight (TOF) cold neutron polarization analysis spectrometer PELICAN at ANSTO in Australia with an energy resolution ΔE = 0.35 meV (upwards to 1 meV or so⁴⁹) and the energy ranges up to 24.8 meV with the energy gain mode used in the q range from 0.08 Å⁻¹ to 4.5 Å⁻¹. The incident energy is 14.9 meV with the wavelength of 2.345 Å. The samples were contained inside aluminum foils in a solid form and sealed in aluminum sample cans in a helium atmosphere. The empty can signal was subtracted at each temperature. The detector efficiency in the data was normalized using a vanadium standard. All steps were performed with standard routines within the LAMP software package⁶⁴ and the scripts are available upon request. The experiment for pure dry GOM was conducted by using a high-intensity Fermi-chopper spectrometer 4SEASONS at J-PARC in Japan⁶⁵. The measurement was done with multi-incident energies⁶⁶ and the data with incident energy 27.1 meV was chosen to cover the energy range up to 17.2 meV and the q range from 0.225 Å⁻¹ to 7 Å⁻¹. The energy resolution at the elastic line is ΔE = 0.8 meV. Similar processing steps to other samples data above were performed on Utsusemi⁶⁷ and Mslice software packages.

The MD simulations were performed using LAMMPS⁶⁸ to simulate ice and supercooled water at 120 K. The supercooled water is simply obtained by simulating bulk liquid water at room temperature and then cooling it down to 120 K. A temperature of 120 K is low enough to freeze out the translational degrees of freedom; therefore, one can consider the disorderly packed water at such temperature as amorphous ice. The equilibration of the MD systems was performed in constant temperature and constant pressure ensemble, using the Nosé-Hoover thermostat and Parrinello–Rahman barostat to control the temperature and pressure, and then switched to NVT ensemble to calculate dynamical properties. The timestep is set as 2 fs. The inter-molecule potential of H₂O used is TIP4P/2005⁶⁹, which shows good accuracy for ice and supercooled water⁷⁰. To reduce the finite-size effect, the simulation are performed using 360,000 and 216,000 molecules for the ordered and disordered structures respectively, with the configuration edge sizes of 180 and 200 Å. Both crystal and amorphous structures were equilibrated at 120 K at 1 atm. The slab structures were then cut from the bulk system to a thickness of 30 Å. We freeze the bottom and top layer (~3 Å) of the slab systems to force the structure to remain 2D confined during the whole simulation. In order to find the optimal volume of the ice slab, we performed a simulation for ice slab in the NPT ensemble. The position of the top and bottom layers are rescaled to new positions when the simulation box changes. A snapshot of the slab geometry can be found in Fig. 4c in the main text.

The VDOS is calculated by the Fourier transform of the oxygen velocity autocorrelation function:14Cv(ω)=∫−∞∞Cv(t)exp(−iωt)dt.\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${C}_{v}(\omega )=\int\nolimits_{-\infty }^{\infty }{C}_{v}(t)\exp (-i\omega t){{{{{\rm{d}}}}}}t.$$\end{document}The velocity auto-correlation function (VAF) is defined as15Cv(t)=<v(0)⋅v(t)>\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${C}_{v}(t)= < {{{{{{{\bf{v}}}}}}}}(0)\cdot {{{{{{{\bf{v}}}}}}}}(t) > $$\end{document}where v(0) are the oxygen velocities.