The spectrum of Ih ice using terahertz time-domain spectroscopy

Terahertz Time-Domain Spectroscopy (THz-TDS) is a potentially powerful tool for studying biomolecules, as its frequency range of 0.1–6 THz coincides with the vibrational frequencies of various biomolecular structures.^1,2 The coupling of the vibrations of such structures to terahertz radiation may allow monitoring and perhaps even modification of these structures in vivo.^3,4 The biomedical research community has used THz-TDS to image freshly excised biological tissue and in vivo to identify disease regions, mainly based on differences in water content.^5–8 It has also been used to study the solvation process of biomolecules and protein dynamics, which revealed non-linear changes in terahertz absorption relative to the concentration of biomolecules that are indicative of the interaction between hydration shells.^9,10 However, the absorption of terahertz by water is high (240 cm⁻¹ @ 1 THz, at room temperature), which limits the penetration depth in imaging and spectroscopy applications and dominates the absorption spectrum of samples in the aqueous phase, particularly in the range of 0.1–6 THz, resulting in featureless spectra, i.e., no resonant peaks.^11,12

Some research groups have tried to overcome the high absorption of terahertz radiation by liquid water by measuring histologically processed tissues, which are essentially dehydrated samples.^13,14 One complication of dehydration is that the dynamics and structures of biomolecules within the samples are altered, and most biological processes do not occur without water.^15,16 For this reason, other research groups have attempted to freeze biological samples, as it is well known that the absorption of ice is ∼60 times smaller than that of liquid water at 1 THz.^17–20 Freezing the samples improved the pathological contrast but still produced largely featureless spectra. It has, therefore, been suggested that the inherent spectra of solvated biomolecules, even when frozen, must be featureless due to a continuum of density of vibrational states associated with a lack of long-range order.¹ 

However, it has been claimed that once solutions of certain biomolecules related to DNA are frozen, specific terahertz resonances are revealed.^4,21,22 The claim was first made in 2016, which reported a spectral feature observed in a frozen solution of methylated cytidine at 1.6 THz.²¹ The observation requires subtracting the absorption spectrum of ice from the absorption spectrum of the sample. The authors in the study²¹ assumed ice’s absorption to be a smooth Gaussian function between 0.2 and 2 THz. This assumption is worth verifying, considering ice, specifically the hydrogen bonds (HBs) within the crystal lattice, can vibrate collectively in the terahertz range.^23–27 Various ice phenomena have been studied with terahertz because of the sensitivity of terahertz to HB vibrations.^28–30 Vibrational modes within ice crystals, if they exist in the 0.2 to 2 THz range, would need to be considered in the spectral analysis of frozen samples, particularly if subtle features are of interest.

In our recent work, we showed that our measurements of the absorption spectrum of ice were not featureless but had a subtle feature around 1.6 THz.³² Here, we expand upon our previous study by performing a comprehensive set of measurements using three commercial systems and further spectral analysis. Early reports of ice’s far-IR spectrum (1–100 THz) showed that ice has prominent absorption peaks at 5, 7, 24, and 96 THz.^35,39 These peaks are assigned to translational lattice vibrations and phonon transitions.^39,40 Still, the frequencies between 0.2 and 2 THz remained largely inaccessible due to the limitations of the available techniques, e.g., Fourier Transform Infrared (FTIR) spectroscopy. Later spectroscopic studies reported ice’s absorption below 2 THz with a higher spectral resolution, and the spectrum in this frequency range follows a profile that consists of a segment with slowly increasing absorption from 0.2 to 1.2 THz and another segment with more rapidly increasing absorption from 1.2 to 2 THz (see Fig. 1). The broad spectral profile is likely the shoulder of a larger spectral peak around 5 THz,³⁶ and its shape can be approximated using a Gaussian function. However, the data from the literature are inconsistent; for example, there is an offset between the spectra reported by Bertie et al.³⁵ and Whalley and Labbé,³³ even though the measurement setups were similar. Spectra acquired with THz-TDS (dashed lines in Fig. 1) are similar. However, these datasets are unsuitable for detailed analysis due to unexplained oscillations, which could be residual Fabry–Pérot etalons (Ashworth et al.³⁸ and Nagai et al.³⁷) or do not cover the spectral range of interest (Takeya et al.³⁴). The data from Funke et al.³⁶ showed a significant deviation from the others. This measurement was taken with an FTIR spectrometer, for which the signal-to-noise ratio is known to rapidly decrease at smaller frequencies, particularly below 2 THz, making them better suited to measurements at higher terahertz frequencies (>3 THz).⁴¹ Note that the temperatures at which the datasets were collected are different, and we would expect to see a trend of the overall absorption increasing with temperature. However, even after correcting for temperature, the variability between the datasets cannot be eliminated, and no trend is seen. Given the amount of variation, verifying a feature as subtle as the one observed in our previous work³² (blue diamonds in Fig. 1) in other published data is impossible. The comparison also highlights the need for reproducibility of published data, particularly those acquired via THz-TDS.⁴² Here, we validate the 1.6 THz feature previously observed in the absorption spectrum of ice³² and investigate its origin. We focus on reproducing the spectrum, particularly in the 1 to 2 THz range, with different THz-TDS systems to verify the 1.6 THz feature beyond system-specific artifacts and noise considerations. We analyze the spectrum using baseline subtraction to highlight the feature centered at 1.6 THz, and we characterize the feature and the underlying broad absorption with two Gaussian functions. We confirmed, via x-ray diffraction, the crystal structure of the ice in our samples as ice Ih, in which each oxygen atom is hydrogen bonded to two hydrogen atoms to form a hexagon with five other oxygen atoms. Ice Ih is modeled using density functional theory (DFT) simulations, and its terahertz absorption is calculated and compared with experimental data. The results improve the understanding of the terahertz spectra of crystalline ice, which could potentially benefit many scientific endeavors, including 6 G capability in below-freezing weather conditions,⁴³ alternative methods for the remote sensing of sea ice in climate change monitoring,⁴⁴ and ice crystal detection in high-pressure cryopreservation.⁴⁵ 

Samples were measured with three commercial THz-TDS systems: TeraPulse 4000 (TeraView Ltd., UK), TeraSmart (Menlo Systems GmbH, Germany), and TeraFlash Pro (TOPTICA GmbH, Germany). Specifications for each system are listed in Table I.

Measurements were made in transmission mode and in a nitrogen-purged enclosure with a relative humidity of <0.1%. A variable temperature cell, VTC (GS21525, Specac Ltd., UK), allows for the measurement of solid samples at cryogenic temperatures, which in this study were from 160 to 260 K. No windows were used on the VTC, which ensured minimal loss of signal while maintaining low humidity in the beam path. Due to the samples being continuously exposed to room-temperature N₂ purge gas, the sample temperature may not have matched the temperature recorded by the cryostat. For measurements relating to the temperature-dependent effect, where the absolute temperature is crucial, an independent thermocouple (EL-USB-TC-LCD, EasyLog, UK) was embedded near the center of a second sample under the same experimental conditions (e.g., sample position, cryostat setpoints) to record more accurate temperature readings. In the following sections, the temperature associated with each measurement is that recorded by the cryostat at the time of spectral acquisition, unless stated otherwise.

To determine sample thickness, we used the method proposed by Withayachumnankul et al.,⁴⁶ which provides the optimal sample thickness for terahertz transmission measurements that would yield the lowest standard deviation in the measured absorption coefficient. We followed the method and calculated the average optimal sample thickness of ice over the range of 1–2 THz. We selected the Teflon spacer with a thickness that would give us a sample thickness closest to this optimal value. A circular, 32 mm diameter, z-cut quartz⁴⁷ window (Q-W-32-3, ISP Optics) was wrapped with laboratory film made from polyolefin and wax to make a hydrophobic surface. A Teflon spacer (outer diameter 32 mm, inner diameter 26 mm, thickness 1 mm) was then placed on the window, and 640 μl of distilled water was pipetted onto the center of the window (Milli-Q^®, a.k.a. reversed osmosis purified and filtered water with a resistivity of 18.2 MΩ/cm at 25 °C). A second window, also wrapped with laboratory film, was placed on top of the Teflon spacer, forming a space that holds the liquid water. The weight of the second window kept the liquid water sample flat while allowing for volume expansion during the freezing process. The assembly was then placed in a standard laboratory freezer, set to −15 °C, for at least 30 min. Once frozen, the windows were separated, and the ice sample was extracted and placed into the VTC’s sample holder. The resulting ice sample was slightly thicker than the spacer due to the decrease in density with the phase change. For example, a 1 mm thick volume of liquid water would freeze to ∼1.3 mm thick. Note that the sample holder and the VTC were pre-cooled to (at least) 200 K to prevent the sample from melting. This method used to prepare the ice samples resulted in some air bubbles within the samples. However, the air bubbles did not affect the terahertz spectrum, as confirmed by measuring ice made using degassed water. We also confirmed the effectiveness of degassing under a cryogenic scanning electron microscope. The ice sample prepared with degassed water had no detectable air bubbles and produced the same spectrum as those prepared with non-degassed water. For more details, see  Appendix A.

We used x-ray diffraction to determine the crystal structure of the ice in our samples, which is subsequently used in the DFT simulation. X-ray diffraction measurements were performed on ice samples prepared as described earlier. A PANalytical Empyrean X-ray powder diffractometer using CuKα radiation collected in reflection (Bragg–Brentano) geometry using a programmable spinner stage was used. Data were collected in a continuous scan mode over a range of angles from 20° to 50°. The sample holder was pre-cooled in dry ice (195 K) to ensure the sample would not melt during data collection.

The sample thickness, refractive index, and absorption coefficient between 0.2 and 2 THz were calculated using commercial software, TeraLyzer Software (Menlo Systems GmbH, Germany).⁴⁸ The software uses an iterative algorithm based on the method described in Refs. 49 and 50. The TeraLyzer eliminates the need for sample thickness measurements as it determines the sample refractive index and the sample thickness from the Fabry–Pérot reflections. It also means that a pre-processing window for filtering out Fabry–Pérot pulses in the time trace is not required. Details about TeraLyzer and the algorithm used can be found on its web page.⁴⁸ 

Although the spectral shape of the absorption coefficient of ice over 0.2–2 THz can be approximated using a Gaussian function, any attempt to fit a function tends to be distorted by the feature around 1.6 THz. Therefore, we assume that the broad spectral profile is a baseline, which is calculated using the “Asymmetric Least Square Smoothing (ALSS)” algorithm available in OriginPro 2022b (OriginLab, USA). ALSS aims to produce a smooth baseline by minimizing the second derivative and deviation from the data points.⁵¹ The asymmetric factor for ALSS was set to 0.001, the threshold was 0.1, the smoothing factor was 2, and the number of iterations was set to 10. More information about asymmetric least square smoothing can be found in OriginPro’s manual.⁵² We then fit a Gaussian function to the ALSS baseline, which is then subtracted from the absorption spectrum, leaving a residual that clearly shows the 1.6 THz feature. The residual containing the 1.6 THz feature can also be fitted with a second Gaussian function to characterize its profile.

Baseline subtraction was performed in the OriginPro 2022b built-in “Peak Analyzer” tool (OriginLab, USA). The baseline was generated using the “Asymmetric Least Square Smoothing” algorithm with the asymmetric factor set to 0.001, threshold 0.2, smoothing factor 2, and the number of iterations set to 10.

Asymmetric least square smoothing is a baseline generation algorithm that aims to produce a smooth baseline that is close to the data.⁵¹ Both the second derivative of the baseline and the deviation from the data points are minimized iteratively. More information can be found in OriginPro’s manual.⁵² 

The far-IR spectrum of ice Ih was simulated in the Cambridge Serial Total Energy Package (CASTEP) based on density functional theory (DFT). The ultrasoft pseudopotential was applied in the generalized gradient approximation (GGA) with the parameterized Perdew–Burke–Ernzerh (PBE) hybrid density functional as the exchange-correlation function to determine the electronic bandgap accurately. The specific calculation parameters were set as follows: On the 2 × 3 × 3 K-point grid uniformly generated by the Monkhorst–Pack method, the convergence threshold of the computational geometric optimization is set to 0.01 eV/A. The spatial reciprocal integral used a kinetic energy cutoff wave function of 600 eV. To model the numerous positions the protons can be in (a phenomenon in ice Ih known as proton disorder⁵³), a primitive cell of Ih ice containing four water molecules was used. For each water molecule in the primitive cell, six different orientations of the H-bonding exist, leading to a total of 24 variations of the primitive cell model.⁵⁴ Four of the 24 variations were identified as representative and optimized. One of the four representative models resulted in converging energy (Fig. 2) and was used for terahertz spectrum calculation. For further details, please see  Appendix B.

The x-ray diffraction data confirmed the structure of the ice in our samples to be hexagonal crystalline ice Ih. The presence of vitreous ice was negligible, based on the flat background of the crystalline sample (Fig. 3).

Figure 4 shows the optical parameters from 0.2 to 2.4 THz calculated using the TeraLyzer software. The top panel [Fig. 4(a)] displays the refractive index profile, the middle panel [Fig. 4(b)] shows the absorption coefficient, and the bottom panel [Fig. 4(c)] shows the residual absorption after the Gaussian baseline has been subtracted from the absorption coefficient spectra. The gray zone above 2 THz indicates where the measurements from the three systems start to deviate, which suggests the onset of noise and system-dependent artifacts. We, therefore, limit the analysis to data below 2 THz. The refractive index, shown in the top panel [Fig. 4(a)], shows a gentle, positive slope from 0.2 THz to around 1.6 THz, at which point the slope turns negative, indicating a clear change in the dispersion. The Kramers–Kronig relations relate the real and imaginary portions of the complex refractive index. Therefore, we expect a corresponding feature around 1.6 THz in the absorption coefficient [middle panel, Fig. 4(b)]; however, such a feature is not immediately obvious. The bottom panel, Fig. 4(c), shows the residual absorption spectra after the Gaussian baseline has been subtracted. The profile is consistent across the systems, with a maximum around 1.6 THz with similar amplitude and full width at half maximum (FWHM). Given such an obvious feature in the residual absorption of ice Ih, we can now confidently be certain of the feature’s existence and that it is not a noise or system-dependent artifact because it is reproduced across three independent systems. Note that there is a $<3$% system-to-system difference between the refractive index spectra that remains relatively constant across frequencies. We attribute this difference to changes in the relative path length of the terahertz beam due to the contraction of metallic components at lower temperatures, the degree of which is determined by the specific component layout of each system. Figure 5 illustrates the reconstruction of the absorption coefficient spectrum using Gaussian functions. The first Gaussian, G1, is fitted to the broad baseline feature of the absorption spectrum shown in Fig. 5(a), while a smaller, narrower Gaussian function, G2, was chosen to fit the residual with the feature around 1.6 THz, shown in Fig. 5(b). Both the Lorentzian and Gaussian models were tested for fitting the residual; however, the Gaussian model provided a marginally better fit. Figure 5(c) displays the combination of the two Gaussian functions (G1 + G2), effectively reconstructing the absorption coefficient across the three systems. This dual Gaussian approach aligns closely with the observed data (correlation >0.99) and accurately represents the feature at 1.6 THz.

The measured spectrum also matched the literature well,^{33–35,37,38} although the data in the literature does not have the same resolution and/or range to allow a more detailed comparison. Figure 1 compares our previously published data (173 K, Tao et al. 2021)³² with other THz-TDS measurements in the literature. Takeya et al.³⁴ produced the closest match with our data, but it only extended to around 1.4 THz. Two other datasets, by Nagai et al.³⁷ and Ashworth et al.,³⁸ had unexplained oscillations that interfered with peak detection. The DFT simulation predicted peaks at 1.8, 5, and 7 THz, which correspond with those given in the literature, verifying the accuracy of the model (see Fig. 11 in  Appendix B). Figure 6 shows the predicted spectral features from the DFT simulation over the range of 0–4 THz. A large peak is shown at 1.8 THz, with a smaller peak on the right-hand shoulder of the main peak at 2 THz. The simulation is run at zero point energy, meaning that any experimentally measured spectra above 0 K will be spectrally broadened. Therefore, we would expect that the 1.8 and 2.0 THz peaks are not resolvable from each other experimentally in this work, given that the lowest temperature measured was 162 K. The presence of the smaller peak at 2 THz and its corresponding blue shift may be verified in future experiments conducted at lower temperatures. The temperature at which the measurements were performed also explains the discrepancy between the simulated peak at 1.8 THz and the position of the feature in the experimental data around 1.6 THz, as we expect a blue shift in peak position with temperature.^56,57 We estimated the peak location at 0 K from temperature-dependent measurements to assess whether the frequency discrepancy is reasonable. As shown in Fig. 7, the blue-shift with decreasing temperature is approximately linear (dashed black arrow), and when extrapolated to 0 K, it matches the frequency of the larger peak at 1.8 THz. Please see  Appendix B for a more detailed comparison between the calculated spectrum and data from the literature.

The HB in the unit cell of ice Ih can, in principle, have many orientations. Our model is based on the most stable orientations tested, which we would expect to be prevalent in actual ice samples. Further verification of the model could be achieved by measuring ice samples at lower temperatures, which could narrow spectral features and resolve the two peaks shown in the simulated spectrum at 1.8 and 2.0 THz. All three THz-TDS systems tested here are widely used commercial systems, which should make reproducing the terahertz spectra of ice below 2 THz more reliable. Future studies could also focus on verifying the result with other THz-TDS systems and modalities, including custom-built systems and in reflection mode, for example. In addition, vitreous ice is arguably more prevalent in biological samples that are flash-frozen,⁵⁸ and its spectrum could be a more accurate baseline for analyzing the spectra of flash-frozen samples. The ice samples prepared in this study contain minimal vitreous ice, as shown in the powder x-ray diffraction pattern (Fig. 3). A natural extension of this research would be to acquire the terahertz spectrum of other forms of ice and compare it to the ice Ih spectrum presented here, which would provide more reliable information for baseline subtraction when studying bio-molecules and other solutes in the frozen state.

Measurements of ice using three different commercially available THz-TDS systems agree and confirm the presence of a feature which blue shifts between 1.55 and 1.65 THz with a decreasing temperature from 260 to 160 K. As expected by Kramers–Kronig relations, the feature is observed in both the absorption coefficient and the refractive index profiles. DFT computer simulations of the ice crystal suggest that the origin of this subtle feature is likely to be a translation in the ice Ih lattice, specifically a mode perpendicular to the 110 lattice plane and another mode perpendicular to the $1̄$10 lattice plane. The feature should be expected in the terahertz spectra of frozen aqueous solutions in general when analyzing for subtle features.

The authors acknowledge the facilities and the scientific and technical assistance of Microscopy Australia at the Center for Microscopy, Characterisation and Analysis, University of Western Australia, a facility funded by the University, State and Commonwealth Governments.

The authors have no conflicts to disclose.

Yu Heng Tao: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Validation (equal); Visualization (equal); Writing – original draft (equal). Xiangyu Dai: Formal analysis (equal); Investigation (equal); Methodology (equal); Validation (equal); Writing – review & editing (supporting). Stephen A. Moggach: Data curation (equal); Formal analysis (equal); Investigation (equal); Visualization (equal); Writing – review & editing (equal). Peta L. Clode: Investigation (equal); Resources (equal); Supervision (supporting); Visualization (equal); Writing – review & editing (equal). Anthony J. Fitzgerald: Formal analysis (equal); Funding acquisition (equal); Methodology (equal); Supervision (equal); Validation (equal); Visualization (equal); Writing – review & editing (equal). Stuart I. Hodgetts: Formal analysis (supporting); Funding acquisition (equal); Supervision (equal); Visualization (supporting); Writing – review & editing (equal). Alan R. Harvey: Conceptualization (equal); Formal analysis (supporting); Funding acquisition (equal); Resources (equal); Supervision (equal); Validation (equal); Visualization (equal); Writing – review & editing (equal). Vincent P. Wallace: Conceptualization (equal); Formal analysis (equal); Funding acquisition (equal); Methodology (equal); Resources (equal); Supervision (equal); Validation (equal); Visualization (equal); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.

Samples were imaged using a JEOL IT800 cryoFESEM. Frozen samples were cut into radial strips and placed onto the sample platform with a cross-section facing up. The sample platform was then transferred into the Leica ACE600 cryo-preparation system for etching and coating. The sample was etched under vacuum (5 × 10⁻⁷ mbar) at 183 K for 2 min and chromium-coated for 6 nm at 0° tilt angle and 5 nm at 27° tilt angle.

CryoSEM images showed entrapped air pockets confined to a single layer in the middle of the sample cross-section [see Fig. 8(a)]. The average diameter of the air pockets is 30 μm. In contrast, ice made with degassed water showed no air pockets [see Fig. 8(b)] but still produced a similar feature at 1.6 THz (see Fig. 9).

The four representative H-bonding orientations of the primitive ice Ih cell are shown in Fig. 10(a) and their corresponding energies in Fig. 10(b). Model “D” had converging energy within 25 optimization steps, while others either oscillated or diverged [see Fig. 10(b)].

The calculated and experimental spectra agree not only in the THz-TDS frequency range (0.2–2 THz) but also from 2 to 12 THz (see Fig. 11). The only mismatch was a slight blue-shifting and an additional peak at around 7 THz, both of which could be due to the temperature difference between the simulation and the experimental data, as demonstrated for ice XI.⁵⁶ Data on ice Ih in the literature already suggest similar effects; comparing the peaks reported by Bertie et al.^35,36 and Funke et al.³⁶ in Fig. 11(b), the former was measured at a lower temperature and showed a slight blue shift with a more prominent second peak at 5 THz compared to the latter.

Figure 12 shows the motion of the atoms in the unit cell for each of the vibrational modes calculated between 0.2 and 12 THz. The modes at 1.77 and 2.01 THz, which are vibrations perpendicular to the 110 and −110 planes, respectively, are likely the source of the 1.6 THz feature observed in the experimental spectrum.