© B. M. Giuliano et al. 2019

1. Introduction

One of the main problems in unraveling the chemical and physical properties of molecular clouds, in which the star and planet formation process takes place, is to estimate correctly the amount of gas contained. The difficulties in the direct observation of molecular hydrogen constrain the possibility to calculate the total mass of a cloud. An easy alternative could be to use carbon monoxide as a tracer of molecular gas, but in dense and cold regions of the interstellar medium and protoplanetary discs, CO is not a good tracer of gas mass because CO molecules preferentially reside on dust grains, forming thick icy mantles (e.g. Dutrey et al. 1998; Caselli et al. 1999). Alternatively, the dust continuum emission is the best available tool to compute a molecular cloud mass, if dust opacities are known.

The advent of Atacama Large Millimeter/submillimeter Array (ALMA) and Northern Extended Millimeter Array (NOEMA) facilities offers the possibility to observe the dust continuum emission in the millimetre and submillimetre part of the electromagnetic spectrum with very high angular resolution and sensitivity. However, to model the dust continuum emission properly it is necessary to have information about its grain size distribution and its chemical composition, since the dust opacity depends directly on these parameters. If we take into account that the dust grains can be covered by ice mantles at the centre of prestellar cores or in protoplanetary disc mid-planes, we also need to investigate how the presence of ices is changing the dust opacities.

Unfortunately, no experimental data are available for these cases, and the interpretation of the dust continuum emission measurements relies on calculated opacity values, such as those tabulated in Ossenkopf & Henning (1994). The goal of our study is to provide laboratory data on the optical properties of CO ice and utilize these data to calculate the opacities of dust grains covered by CO ice mantles. We compare the opacity values obtained by our study to those available in the literature.

The presently available set of data focusses mainly on the determination of the optical constants in the visible and mid-infrared (MIR) range (Hudgins et al. 1993; Ehrenfreund et al. 1997; Baratta & Palumbo 1998; Loeffler et al. 2005; Dartois 2006; Palumbo et al. 2006; Warren & Brandt 2008; Mastrapa et al. 2009).

Far-infrared (FIR) studies on spectral properties of molecular solids, without deriving optical constants, started early, to deepen the understanding of the infrared-active lattice vibrations of simple species. Anderson & Leroi (1966) studied frequencies of CO and N₂ in the range 40–100 cm⁻¹, and Ron & Schnepp (1967) complemented the available information with CO, N₂ and CO₂ intensity studies in the same frequency range. In 1994 Moore and Hudson published a comprehensive study of FIR spectra of cosmic type ices, including ice mixtures. These data included the analysis of amorphous and crystalline phases of the pure molecular ices, and the authors discussed the implications of the results on the identification based on astronomical observations. An estimation of the band strengths in the FIR region for pure ices and ice mixtures relevant for astrophysical environments can be found in Giuliano et al. (2014, 2016).

Recently, the investigation of the terahertz (THz) spectroscopic properties of ice mantles analogues has gained considerable interest (Allodi et al. 2014; Ioppolo et al. 2014; McGuire et al. 2016). This technique allows direct measurement of the intermolecular vibrations in the ice samples related to the lattice structure, which can be connected to their large-scale structural changes and finally to their thermal history. On the contrary, spectroscopic features measured in the MIR frequency range are indicative of the intramolecular vibrations of the sample, which can provide a wealth of information on the molecular identification and chemical reactivity. A comparison of our THz experimental data with that observed in the MIR range could help us to reveal intra- and inter-molecular vibrations. However, this study is beyond the scope of this paper and will be addressed in our future investigations.

Nowadays, numerous spectroscopic methods are extensively used for dielectric measurements at THz frequencies (Lee 2009); these include the following: Fourier transform infrared (FTIR) spectroscopy (Griffiths & de Haseth 1986), Backward-wave oscillator (BWO) spectroscopy (Komandin et al. 2013), spectroscopy based on photomixing (Preu et al. 2011) or parametric conversion (Kawase et al. 1996; Kiessling et al. 2013), and, finally, THz time-domain spectroscopy (THz-TDS; Auston 1975; Van Exter et al. 1989). These methods exploit either continuous-wave or broadband sources, operate in different spectral ranges, and are characterized with different sensitivity and performance. Among these methods, THz-TDS seems to be the most appropriate for studying laboratory analogues of circumstellar and interstellar ices. In contrast to other approaches, THz-TDS yields detection of both amplitude and phase of sub-picosecond THz pulses in a wide spectral range as a result of a single measurement; thus, the reconstruction of the dielectric response of a sample might be performed without using the Kramers–Kronig relations (Martin 1967) and involving additional physical assumptions. Furthermore, THz-TDS yields analysis of separate wavelets forming the time-domain response of a sample; thus, it is a powerful method for the characterization of multilayered samples. Thereby, we selected THz-TDS as a spectroscopic technique for our experiments.

We aim at the extension of the laboratory data in the FIR/THz region, and we show how the employment of the THz-TDS is able to provide direct measurement of the real and imaginary part of the refractive index of the ice sample. The experimental and theoretical methods employed are explained in Sect. 2, the results obtained and how these data are relevant for astrophysical application are presented in sections Sects. 3.3 and 4, respectively; the conclusions are illustrated in Sect. 5.

2. Experimental and theoretical methods

For this series of experiments a dedicated set-up has been designed and developed in the laboratories of the Center for Astrochemical Studies (CAS) located at the Max Planck Institute for Extraterrestrial Physics in Garching (Germany). The set-up is composed of a closed-cycle He cryocooler coupled to a THz time-domain spectrometer. The cryocooler vacuum chamber is small enough to be hosted in the sample compartment of the THz spectrometer, and it is mounted on a motor controlled translational stage, which ensures the tuning of the cryostat position with respect to the THz beam. The details of the main components of the apparatus and the ice growing procedure are given in the following subsections.

2.1. Cryogenic set-up

The high power cryocooler was purchased from Advanced Research Systems. The model chosen is designed to handle high heat loads thereby ensuring a fast cooling. This instrument is equipped with a special interface capable of reducing the vibration transmitted from the cold head to the sample holder at the nanometer level. This requirement is important in case of spectroscopic measurements in the THz frequency region, where the induced vibration of the sample can cause the increase of the noise level of the recorded spectra.

The cryostat is placed in a 15 cm diameter vacuum chamber, equipped with four ports for optical access and for the gas inlet. The optical arrangement is designed to work in the transmission configuration. The optical windows chosen for the measurements at the desired frequency range are made of high-resistivity float-zone silicon (HRFZ-Si), purchased by Tydex. This material features a high refractive index of n_Si = 3.415, negligible dispersion, and impressive transparency in the desired frequency range. The same material was chosen as a substrate for ice growing. In order to suppress the Fabry-Perot resonances in the THz spectra, caused by multiple reflection of the THz pulse within the windows, the thickness of the three Si windows must be different from each other. For this purpose we chose to use the 1 mm thick Si window as substrate for the ice deposition, and we used two windows of 2 mm and 3 mm thickness as optical access to the THz beam.

We performed measurements with a slightly focussed THz beam configuration in order to mitigate vignetting and diffraction of the beam at the metal components of the vacuum chamber. A schematic overview of the chamber arrangement is sketched in Fig. 1. The pumping station is composed of a turbomolecular pump (84 ls⁻¹ nitrogen pumping speed) combined with a backing rotary pump (5 m³ h⁻¹ pumping speed), thus providing a base pressure of about 10⁻⁷ mbar. The minimum temperature measured at the sample holder in normal operation mode is 5 K.

Fig. 1. Sketch of the vacuum chamber of the cryostat coupled to the THz beam at the CAS.

2.2. Terahertz time-domain spectrometer

The THz-TDS set-up used for this work was purchased from the company Batop GmbH. The model chosen is the TDS-1008, which has a customized sample compartment in order to allocate the cryostat. This model is based on two photoconductive antennas made of low-temperature-grown gallium arsenide (LT-GaAs), which constitute the emitter and detector of the THz pulse (Lee 2009). The antennas are triggered by a femtosecond laser (TOPTICA, 95 fs, 780 nm) with a pulse repetition rate of 100 MHz and an average input power of 65 mW. Further details on the set-up are provided in Appendix A.

Figure 2a shows the optical path of the laser beam into the optical bench of the spectrometer. Panel b in Fig. 2 shows standard broadband Fourier spectra. The THz pulse registered with the beam path empty is then converted in the blue spectrum spanning the frequency range from 0.05 THz to 3.5 THz with maximal spectral amplitude centred at about 1.0 THz. In TDS, the time-domain THz waveform is converted to the frequency domain using the direct Fourier transform, which yields the frequency-dependent amplitude and phase of the THz wavelet. Since the frequency-domain data are calculated via the direct Fourier transform, the spectral resolution of measurements is determined as Δν = 1/ΔT, where ΔT is a size of the time-domain apodization filter, chosen to avoid the edge effects (i.e. the Gibbs effect) in the frequency domain. In our experiments we used the 35-ps Tukey apodization filter (see Appendix B), which yields the spectral resolution of about 0.03 THz.

Fig. 2. Set-up of THz-TDS for the spectroscopy of ices. Panel a: chematic of the set-up, where FS laser stands for the femtosecond laser, M stands for the optical mirrors, ATT stands for the attenuator of the laser beam intensity, BS stands for the optical beams splitter, DS stands for the mechanical double-pass delay stage, PCAE and PCAD stand for the photoconductive antenna-emitter and antenna-detector, respectively, L stands for the TPX lenses, and S and W stand for the HRFZ-Si substrate and windows, respectively. Panel b: spectra of THz waveforms E(ν) transmitted through the empty THz beam path or the THz beam path with the cryostat; the shaded region below ≈0.3 THz indicates the spectral range in which distortions from the THz beam diffraction on the aperture of the substrate are expected.

The green spectrum was recorded with the cryostat placed in the sample compartment. It is converted from a waveform which contains both a first THz pulse (i.e. a ballistic one) and a train of satellite pulses originating from the multiple THz wave reflections within the windows. The ballistic THz pulse is delayed in the input/output windows and the substrate of the vacuum chamber. The spectrum of this waveform is slightly suppressed owing to the Fresnel losses and modulated due to the interference of the ballistic pulse and the satellites.

In Fig. 2b, we show the shaded area at lower frequencies, where we expect growing distortions of the experimental data caused by the THz beam diffraction on the aperture of the substrate that is 20 mm in diameter. Assuming that the THz beam spot formed at the substrate is diffraction-limited, the lateral intensity distribution in the spot is defined by the Bessel function of the first kind (Born & Wolf 1980). The resulting width of the first intensity peak is approximately (3.8/π)(f/D)(c/ν), where D = 25 mm and f = 67 mm stand for the diameter and the back focal distance of the focussing lens, respectively. From this model, we deduce the critical frequency of 0.3 THz, below which less than 95% of the beam energy passes through the substrate aperture. Thereby, considering both the spectral sensitivity of our THz-TDS set-up and the diffraction limits, the spectral operation range of our experimental set-up is approximately limited within 0.3–2.0 THz.

The THz-TDS housing is kept under purging with cold nitrogen gas during the entire experiment to mitigate the absorption features due to the presence of atmospheric water in the THz beam path. The residual humidity measured at the sensor was less then 10⁻³%.

2.3. Ice preparation

The ices are prepared using a standard technique in which the molecular sample in its gaseous form is allowed to enter the vacuum chamber through a stainless steel 6 mm pipe. The gas flux is controlled by a metering valve. Once the gas is expanding inside the vacuum chamber it condensates on the substrate.

For this set of experiments an ice thickness of the order of millimetre is required to fulfil the sensitivity characteristic of the THz-TDS set-up. This value is orders of magnitude higher than the usual ice thickness reached using this deposition technique, which is of the order of μm. To deposit such a thick ice in a reasonable amount of time, we chose fast deposition conditions, in which a considerable amount of gas is introduced in the vacuum chamber. In these conditions it is very difficult to obtain an ice sample which is homogeneous enough to determine its optical properties. To overcome this problem the gas inlet characteristics must be set accordingly.

We decided to remove any directionality from the gas inlet, keeping the pipe end cut at the vacuum chamber wall at a distance of approximately 7 cm from the substrate (see Fig. 1). This configuration creates an ice layer on each side of the substrate. During the deposition, the pressure measured inside the chamber is approximately 10⁻² mbar. The ice deposition was divided in steps of 4, 5, and 6 min duration, up to a total of 30 min deposition time, in three different experiments. The final temperature was up to 28.5, 31.2, and 33.1 K for each step at 4, 5, and 6 min deposition time, respectively. This increase is due to the condensation of the gas onto the cold surfaces of the cryostat, which is producing a heating rate too fast to be dynamically removed from the cooling system during the deposition. After each step, the THz spectrum has been recorded. This procedure was performed to rule out possible effects on the ice structure due to different deposition temperatures. As interstellar ices can be commonly found at temperatures as low as 10 K, the temperature of the cold substrate has been kept at 14 K, which is the lowest temperature achievable in the set-up in this configuration, because the radiation shield of the sample holder must be removed to ensure that no directionality of the molecular beam is present. Before moving to the next deposition step, the system was allowed to thermalize and spectra recorded after each deposition step have been taken at a temperature of 14 K.

As stated in Urso et al. (2016) the analysis of the Raman and infrared spectra of experiments performed at increasing tempertaures from 17 to 32 K show no profile variation in the band at 2140 cm⁻¹, which could be ascribed to a structural change in the ice morphology.

The waveform recorded in the time domain is compared, as reference for the measurements, with the waveform recorded for the substrate without ice, kept at a reference temperature of 14 K as well. After the deposition was completed, we measured the spectra in different regions of the sample, to ensure that the ice morphology is spatially homogeneous. The results obtained from the spectra measured on a grid of 11 points spaced by 2 mm are in agreement within 10%, indicating a uniform ice formation over the substrate.

3. Derivation of the optical constants

In order to determine the optical constants, the ice thickness must be known. The laser interference technique is a well-established method to estimate the thickness of an ice sample deposited on a substrate as a function of the time. The absolute accuracy of this method is approximately within 5%, but the maximum CO ice thickness that we can measure with this technique is limited to 5 μm before the reflected laser signal becomes too weak to be detected, due to scattering losses occurring both in the bulk and on surface of the film. Thus, this technique is well suited for studying thin layers, when the ice thickness in total is below 10 μm (5 μm on each side of the substrate), but is not appropriate for experiments on thick ice samples. In turn, for the thickness estimation of the millimeter-size sample of ice featuring rather low THz absorption, we developed a model to perform an initial estimation of the optical properties and the ice thickness directly from the recorded THz spectra, as described in the following subsections.

3.1. Ice parameters modelling

Our model aims to reconstruct the optical properties of ices, which are defined as follows:

(1)

where n′ and n″ are the real and imaginary parts of the complex refractive index n, c is the light speed, and α is the amplitude absorption coefficient, which is defined as half of the value of the power absorption coefficient.

Equivalently, we can write

(2)

where ε′ and ε″ are the real and imaginary parts of the complex dielectric permittivity ε.

The model describes the THz wave propagation through the substrate with the ice deposited on both surfaces. The reconstruction of the ice parameters proceeds following three main steps.

The first task is modelling the reference and sample waveforms. Because of the focussed arrangement of the THz beams, the electromagnetic wave is assumed to be planar and to interact with the sample interfaces at the normal angle of incidence. This is a common and conventional assumption widely applied in dielectric spectroscopy (Pupeza et al. 2007; Zaytsev et al. 2014). It allows us to describe all the features of the THz pulse interaction with the multilayered sample using the Fresnel formulas, which define the THz wave amplitude reflection at (and transmission through) the interface between the media and the Bouguer-Lambert-Beer law, which defines the absorption and phase delay of the THz wave in a bulk medium. Further details on these assumptions are given in Appendix B.

Figure 3 represents the THz wave propagation through the three layers structure: the first ice film, the HRFZ-Si substrate, and the second ice film, where the symbols 0–3 and N correspond to different components of the plane wave passing through the multilayered structure. As shown in Fig. 3, for the sample spectrum, we took into account the contribution of the ballistic THz pulse (1) and the satellite pulses (2 and 3), caused by the multiple THz wave reflection in the ice films. The mathematical description of the wave propagation can be found in Appendix B.

Fig. 3. Time–distance diagram illustrating the THz wave propagation through the HRFZ-Si substrate with ice films deposited on its surfaces. Lines 0–3 illustrate the ballistic pulse and satellite pulses transmitted in the direction of the antenna-detector; N stands for unaccounted satellites with larger time delays. The solid lines represent the pulses used for the analysis, while dotted lines correspond to neglected pulses.

The second step consists in estimating the initial thickness l_CO, I, l_CO, II and the initial complex refractive index n_init of the two ice layers as shown in Fig. 4. The thickness estimation can be derived from the time delay δt₀₁ between the ballistic pulse of the reference and sample waveforms (0 and 1 in Fig. 4), the first satellite pulse (2) and the ballistic pulse (1) of the sample waveform δt₁₂, and the second satellite pulse (3) and the ballistic pulse (1) of the sample waveform δt₁₃, in Fig. 4a. Since the HRFZ-Si has a very high refractive index, we consider that the refractive index of ice n satisfies the inequality n₀ = 1.0 <  n <  n_Si = 3.415, where n₀ is the refractive index in vacuum.

Fig. 4. Calculation of the thicknesses lCO, I, lCO, II and the complex refractive index n of the CO ice films (see Fig. 3). Panel a: time delays between the ballistic THz pulses of the reference (0) and sample (1) waveforms, δt01; the first satellite pulse (2) and the ballistic pulse (1) of the sample waveform, δt12; and the second satellite pulse (3) and the ballistic pulse (1) of the sample waveform, δt13 (the pulses are delineated as in Fig. 3). Panels b and c: estimates for the thicknesses of the two ice films as a function of the total deposition time tdep for the different deposition intervals of Δtdep = 4, 5, and 6 min; the first assumption for the real part of the refractive index of ice is between and 1.255 ± 0.035; for the imaginary part we first set .

Then, by neglecting imaginary parts in the complex refractive indexes of media, the first assumptions for the real part of the complex refractive index of ice films and both their thicknesses l_CO, I and l_CO, II are described as follows:

(3)

(4)

Equation (3) is obtained from a mathematical model of sample and reference waveform. Further details on the derivation can be found in Appendix B. From Eq. (4) it is possible to obtain information only on the thicknesses of ice films, while the identification of the specific layer (I or II) is not allowed. We observe a linear increase of the ice thickness with the total deposition time. A first assumption for the real part of the complex refractive index of CO ice is in the range to 1.255 ± 0.035 for all the considered deposition intervals; here, the error accounts for an accuracy of the THz pulse peak position estimation. The first assumption for the imaginary part of the complex refractive index of ice has been done considering α_init = 0; thus, .

Finally, from the first estimation of the ice thickness and complex refractive index, it is possible to reconstruct the THz dielectric response of ice. The reconstruction procedure is reported in Appendix B, while the results are summarized in Fig. 4. Panels b and c show the growth of the two CO ice layers in time (t), considering different deposition steps of Δt = 4, 5, and 6 min.

We can compare the results obtained from the THz spectral data to the thickness calculation performed with the well-established laser interference techniques described in Sect. 3.2 to validate the calculation of the ice thickness using this model. The good agreement between the two methodologies confirm the validity of the present analysis.

3.2. Laser interference technique

In the adopted experimental configuration, a He–Ne laser beam (λ = 632.8 nm) is directed towards the sample and reflected at near normal incidence both by the vacuum-sample and sample-substrate interfaces. The reflected beam is detected by an external diode detector. It is possible to follow the accretion of the ice film by looking at the interference curve (intensity vs. time) of the reflected laser beam. Further details on the laser interference technique can be found in Urso et al. (2016)¹.

The results obtained with the two techniques have been compared. The data on the accretion of the ice versus time obtained from the analysis of the interference curve, measured at the early stage of the deposition process, are in good agreement with the data obtained from the THz spectra. In addition, using the laser interference technique, we obtain for the CO ice a refractive index n_CO = 1.27, that is close to the value obtained using the THz technique. The good agreement between the two methodologies confirm the validity of the present analysis.

3.3. Reconstruction of the THz response

The recorded THz waveforms and their Fourier spectra are presented in Fig. 5 for the optical substrate (used as a reference) and CO ice samples at increasing thicknesses. In panel a the waveform E(t) is shown for the reference (green) and five subsequent deposition steps (black to light red) of approximately 0.45 mm total thickness for each step. The thickness reached after the total deposition time is approximately 2.3 mm, split in two ice layers of ≈0.85 mm and ≈1.45 mm on top of each side of the substrate. A small source of inhomogeneity can be ascribed to the position of the pipe connected to the pumping system, which is located in a lateral position of the vacuum chamber with respect to the cold substrate.

Fig. 5. Evolution of the THz pulse and its spectra during the CO ice deposition. Panel a: reference waveform E(t) transmitted through the cryostat with the empty substrate (green), and sample waveforms (black to light red) transmitted through the substrate with the CO ice deposited on its surfaces. Panel b: Fourier spectra |E(ν)| of the reference and sample THz waveforms calculated with the use of Tukey apodization. The waveforms in (a) and spectra in (b) correspond to different values of the total deposition time tdep (indicated); the deposition intervals are Δtdep = 6 min.

In the THz spectrum of CO ice, we observe a Lorenz-like resonant peak centred near 1.5 THz (50 cm⁻¹) and a second blurred feature close to 2.5 THz (83 cm⁻¹), masked by the sharp bands produced by the atmospheric water contamination in the spectrometer sample compartment. The highest frequency accessible in our set-up is presently limited by the strong absorption features of residual water and carbon dioxide in the spectrometer case. We plan to change the current set-up to an evacuated case to get rid of the contamination from the residual atmosphere, and we expect to extend the accessible frequency range up to 4 THz.

The estimated deposition rate for these deposition conditions is ≈0.05 mm min⁻¹ for layer I and ≈0.03 mm min⁻¹ for layer II. These values are in a reasonable agreement with the results obtained employing the laser technique, with which the deposition rate is calculated to be 0.02 mm min⁻¹. This agreement validates our hypothesis that the ice structure of thick ices does not differ significantly from the structure of thin ices, growing homogeneously over time during the deposition. The calculated optical properties are independent from the total thickness of the ice sample, allowing us to relate the laboratory data to the astrophysical ice conditions.

Figure 6 shows the determination of CO ice parameters, which are the refractive index a, the amplitude absorption coefficient b, and both the real c and imaginary d parts of the complex dielectric permittivity; see Eqs. (1) and (2).

Fig. 6. Optical properties of CO ice. (a) Real part of the refractive index, (b) amplitude absorption coefficient, and (c) real and (d) imaginary parts of the dielectric permittivity (see Eqs. (1) and (2)). For all deposition intervals, the dielectric curves demonstrate the existence of a Lorenz-like absorption peak, centred near 1.5 THz and featuring similar bandwidth. Distortions of the results seen at frequencies below 0.3 THz (such as an oscillatory character of n and ε′ and an increase of ε″ with decreasing frequency) are due to diffraction effects (see Sect. 2.2).

4. Discussion

A benefit of the direct reconstruction of the optical properties of the ices, provided by THz-TDS, is the detection of the frequency-dependent amplitude and phase of the waveform in a broad frequency range as a result of a single measurement. These data eliminate the need to use the Kramers–Kronig relations for the reconstruction of the optical properties, excluding additional distortion of the experimental data by edge effects, which frequently appears as a result of the Hilbert integral transformation. This is of particular importance when dealing with broadband spectral kernels which are usually present even when operating at low temperature.

We could also compare the refractive index of CO ice at THz frequencies with that previously calculated in the MIR range, from Hudgins et al. (1993), Ehrenfreund et al. (1997), and Baratta & Palumbo (1998) for a CO ice deposited at 10 K. However, a direct comparison between refrative index values found by different authors is not straightforward. As discussed by Loeffler et al. (2005) and Baratta & Palumbo (2017), the density and in turn the refractive index of an ice sample could strongly depend on the experimental conditions such as temperature, growth angle, and deposition rate. We plan to test the effect of the change in the deposition conditions on the ice structure and check if this change will affect the optical constants. These results will be published in a forthcoming paper.

We compare the spectroscopic signature of the CO ice in our experiments with previous studies available in the literature. Data on FIR spectra of solid CO were reported by Anderson & Leroi (1966) and Ron & Schnepp (1967). These studies investigated the absorption spectra of amorphous CO, deposited at 10 K on a crystalline quartz substrate between 250 and 30 cm⁻¹. Two bands are visible at 50 and 83 cm⁻¹ (1.5 and 2.5 THz). The spectral features observed in our experiments are in excellent agreement with these data, even though the 2.5 THz feature is masked by atmospheric water bands in our set-up. The THz-TDS technique has a low sensitivity, requiring very thick ice layers to be detectable. Using our fast deposition rate the MIR vibrational bands of CO were strongly saturated within the first 30 s of deposition. We did some preliminary tests on the ice growing using the vibrational bands in the NIR range and we were able to follow the ice growing up to approximately 160 μm thickness. Also in this case the ice growing is constant and linear in time during deposition, but the error associated with the thickness estimate using the NIR band is large (20–30%).

It is not surprising then that in the data reported by Ioppolo et al. (2014), on THz and MIR spectroscopy of interstellar ice analogues, the CO ice absorption band in the THz region was not observed. The ice thickness in these experiments was estimated in all cases to be less than 10 μm. In our case, the minimum thickness required to be able to observe the 1.5 THz feature is estimated to be of the order of hundreds μm.

The data obtained are then employed to calculate the dust opacity for a given grain size distribution, as reported for example in Ossenkopf & Henning (1994); see Appendix C for additional details on the method we developed to reproduce their results. Following their approach, we report in Fig. 7 the calculated dust opacities assuming different ice coatings and different experimental data for the optical constants of ices. Data from Ossenkopf & Henning (1994) are labelled OH94 and reported for bare grains (green dotted line), thin (blue dotted), and thick ice mantles (orange dotted). The labels V = 0, 0.5, and 4.5 indicate the volume ratio of the core of refractory material to the ice mantle (see Appendix C). Conversely to the present work, Ossenkopf & Henning (1994) assumed a H₂O:CH₃OH:CO:NH₃ = 100:10:1:1 mixture ice mantle composition, i.e. water-based, with a minor amount of methanol, carbon monoxide, and ammonia. As described in detail in Appendix C, following the same procedure as in Ossenkopf & Henning (1994), we extended the real and the imaginary parts of the refractive index from Hudgins et al. (1993) to longer wavelengths and we included spherical carbonaceous impurities.

Fig. 7. Calculated and reference opacities of astrophysical dust with CO ice and ice mixtures as a function of the wavelength. The dotted lines labelled with OH94 refer to bare grains and ice mixtures by Ossenkopf & Henning (1994), the dashed lines with BP98 to CO ice by Baratta & Palumbo (1998), and the solid lines to the CO data by the present work. The value V indicates the volume ratio of refractory core to ice mantle, for which we follow Ossenkopf & Henning (1994), where V = 0 (black) is the bare grain, V = 0.5 (blue) thin ice, and V = 4.5 (orange) thick ice. See text for additional details.

The dielectric functions found by our experiments refer to pure CO ice, which explains the differences in the opacities shown in Fig. 7 calculated in the 100–1000 μm range (solid lines). To have a more relevant comparison, we calculated the opacity using the same grain distribution and refractory materials as in Ossenkopf & Henning (1994). But in this work we extrapolated the optical values of CO ice coating from Baratta & Palumbo (1998, labelled BP98), where the refractive index of CO ice deposited at 12 K is calculated from the spectrum recorded in the 4400–400 cm⁻¹ (2.27–25 μm) infrared spectral range. The opacities of Baratta & Palumbo (Fig. 7, dashed lines) show an agreement with our results, except for the contribution from the presence of a CO ice absorbing feature at approximately 200 μm, which is absent in the extrapolated data, as expected.

When we compare the refractive index, we note that the real part by BP98 (n′ = 1.28) is reasonably close to our data (n′ = 1.24), while for the imaginary part the discrepancy is very large, since n″ is completely determined by the absorption feature at 200 microns that is out of the range investigated by Baratta & Palumbo (1998) (see Fig. 6). It is worth mentioning that Anderson & Leroi (1966) and Ron & Schnepp (1967) reported a CO absorption feature at 2.5 THz, which is not visible in our spectra because it is masked by the atmospheric water contamination. This feature should further decrease the actual value of n′ at higher frequencies (above 2.5 THz). If we could measure the 2.5 THz feature and extend the calculation of the optical properties to this value, we would probably obtain a lower value of n′, slightly increasing the discrepancy with the data by Baratta & Palumbo (1998). The reconstruction of the THz dielectric response of ices without the use of the Kramers–Kronig relations, which is provided by THz-TDS, can provide an independent methodology to determine the optical properties of ice samples and validate the previous studies.

Since the real and imaginary parts of the dielectric constant are employed to compute the opacity, we can infer from the calculated opacity curve that the imaginary part, which shows the biggest difference from the data presented by Baratta & Palumbo (1998), does not play a major role in the determination of the opacity. This conclusion might be different for other absolute values or in a different spectral range.

Values of opacities for some selected wavelengths are also given in Table 1. The motivation for the choice of CO ice arises from the need of investigating the ice mantle properties for sources in which drastic CO depletion is expected, such as prestellar cores or protoplanetary discs mid-planes (e.g. Caselli et al. 1999; Pontoppidan et al. 2003). In these cases, an ice mantle rich in CO can be formed and influence the optical properties of the dust grains. Thus, it is interesting to compare how the opacities change when ice mantles with diverse chemical composition are present. The common molecular species in astrophysical ices, such as H₂O, CO₂, NH₃ and possibly N₂, present absorption features in the FIR range. Therefore, the study of the influence of the spectroscopic features on the opacity of the ice mantles is important and we plan to extend this study to other pure ices and ice mixtures, which could be representative of the different molecular ice compositions in various astrophysical environments.

5. Conclusions

In this work we have presented a study on the optical properties of solid CO at temperature and pressure conditions significant for astrophysical applications. While previous data in the MIR frequency range are available in the literature, to the best of our knowledge, this study is the first to provide the complex refractive index and complex dielectric permittivity of CO ice in the THz range.

We have shown that the ability of THz-TDS to measure both the amplitude and phase information about the transmitted pulse provides direct reconstruction of the complex dielectric function of ices without the use of the Kramers–Kronig relations. The THz spectral features of ices can have a large bandwidth, such as the CO absorption line at 1.5 THz. In this case, an implementation of the Kramers–Kronig relations (i.e. the Gilbert transform), relying only on the power transmission/reflection spectrum, could lead to edge effects and resulting distortion of the dielectric response. Such distortions are of particular importance when a spectral feature of the sample is located near the border of the spectral operation range. Our results justify that THz-TDS set-up is an appropriate instrument for accurate measurements of dielectric properties of ices at THz frequencies.

These results have been used to calculate the opacities of the dust grains covered by a CO ice layer. Discrepancies with currently available opacities suggest that measurements such as those presented in this work are needed to provide a better interpretation of dust continuum emission, including dust and gas mass estimates. In addition, they will provide further insight into the radiative transfer processes based on ice analogues optical and physical properties.

Acknowledgments

The authors acknowledge the anonymous referee for providing very useful comments which significantly improved the quality of the manuscript. Also, the authors acknowledge Mr. Christian Deysenroth for the very valuable contribution in the designing and development of the experimental set-up. We would like to thank Volker Ossenkopf for fruitful discussion to offering insight into the dust modelling methodology. The work of Arsenii A. Gavdush on alignment of the THz-TDS set-up and digital processing of the THz waveforms was supported by the Russian Foundation for Basic Research (RFBR), Project #18-32-00816. The work of Gennady A. Komandin in solving the inverse problem of THz time-domain spectroscopy was supported by the Russian Science Foundation (RSF), Project #18-12-00328. Tommaso Grassi acknowledges the support by the DFG clusterof excellence “Origin and Structure of the Universe” (http://www.universe-cluster.de/). This work was funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – Ref no. FOR 2634/1 ER685/11-1. This work has been partially supported by the project PRIN-INAF 2016, The Cradle of Life – GENESIS-SKA (General Conditions in Early Planetary Systems for the rise of life with SKA).

References

Appendix A: THz-TDS optics

In this section we provide a detailed description of the configuration of the laser beam inside the optical compartment of the set-up, already shown in Fig. 2. The power of the laser beam is attenuated and successively divided into equivalent channels. Thus, the antenna-emitter is pumped and the antenna-detector is probed with an equal average power of about 20 mW. The optical delay between the pump and probe beams is varied using a double-pass linear mechanical delay stage from Zaber with the travel range of 101.6 mm and the positioning accuracy of < 3 μm. The THz radiation undergoes 10 kHz electrical modulation to detect synchronously the THz amplitude using a lock-in detection principle.

The THz beam emitted by the photoconductive antenna is collimated by an integrated HRFZ-Si hemispherical lens and then focussed on a substrate window using a polymethylpentene (TPX) lens with a focal length of 67 mm and diameter of 25 mm. After passing through the cryostat vacuum chamber, the beam is collimated by an equal TPX lens in the direction of the antenna-detector. Finally, the THz beam is focussed onto the photoconductive gap of the antenna-detector by an equal integrated HRFZ-Si hemispherical lens. In our measurements, during the waveform detection, we used a time-domain stride of 50 fs, which allows for satisfying the Whittaker–Nyquist–Kotelnikov–Shannon sampling theorem (Nyquist 1928), at time-domain window size of 100 ps and an averaging time of 0.1 s at each time-domain step with no waveform averaging.

The signal measured at the antenna-detector is recorded in the time domain E(t) and converted into its Fourier-spectrum E(ν) via

(A.1)

where t and ν stand for time and frequency.

Appendix B: Reconstruction of the terahertz dielectric permittivity

Equations (3) and (4) are obtained as follows. We defined the complex amplitude of electromagnetic wave E₀ = |E₀|exp(iφ₀), which interacts either with a bare substrate (which we use as reference) or with two ice layers (when detecting the sample waveform). The pulses complex amplitudes depend on the initial complex amplitude of electromagnetic wave E₀ and on ices and reference layers, which define the THz-wave reflection and transmission at the interfaces (R-operators and T-operators, respectively) and its absorption and phase delays in a bulk material (P-operators).

Then, we define the amplitudes of the ballistic reference pulse as follows:

(B.1)

where l is the thickness of the medium, and the symbols 0 and Si indicate the vacuum and HRFZ-Si, respectively,

the ballistic sample pulse as follows:

(B.2)

where l_CO, I and l_CO, II are defined as in Sect. 3.1,

and the two satellite sample pulses as follows:

(B.3)

which are clearly observed in reference and sample TDS waveforms in Fig. 4a.

Then, by neglecting the phase changes during the reflection at the interfaces of absorbing media as well as the distortion of optical pulses due to dispersion of material parameters, we calculated the phases of these pulses as

(B.4)

where n₀ and n_Si are defined as in Sect. 3.1 and n_CO is the refractive index of the CO ice. The phase changes are weak since we study rather low-absorbing dielectric materials. The dispersion of HRFZ-Si is very low, while the dispersion in ice is negligible owing to its small thickness.

These phases are used to calculate the time delays between the pulses, which are indicated in Fig. 4a, i.e.

(B.5)

Solving this system of equations yields Eqs. (3) and (4).

Figure 5b shows the Fourier spectra E(ν) of the reference and sample TDS waveforms. In order to filter out the contribution of the satellite THz pulses (caused by the interference in the input and output windows and the substrate) and to improve the analysis of the frequency-domain data, we apply equal apodization procedure (window filtering) to all waveforms,

(B.6)

where E(t) and E_filt(t) stand for the initial and filtered waveforms, H(t) defines the apodization function, and t₀ defines a position of the apodization filter towards the THz waveform, i.e.

(B.7)

The parameter H(t) is a Tukey apodization filter (Tukey et al. 1986) with the width τ, and ω stands for a parameter of the filter smoothness; for ω = 0 the window has a rectangular form and for ω = 1 it is the Hann window Harris (1978). As shown in Fig. 5a, we use the Tukey window centred at the maximum of the reference THz waveform with the smoothness parameter of ω = 0.1 and the width of 40 ps, which yields the frequency-domain resolution of 0.025 THz.

Let us consider the Fresnel formulas (Born & Wolf 1980), defining the THz wave amplitude reflection at (and transmission through) the interface between media m and k as follows:

(B.8)

where R_m, k(ν) and T_m, k(ν) stand for coefficients of the complex amplitude reflection and transmission, respectively, while n_m(ν)+in_k(ν) is the complex refractive index of the media. The relation between complex amplitudes of the THz wave right after the emitter (z = 0), E₀(ν) and at the position z along the beam axis, E(ν,z) is given by

(B.9)

If the thicknesses and refractive indexes of all layers are known, Eqs. (B.8) and (B.9) yield description of all peculiarities of the THz pulse interacting with multilayered structures (Zaytsev et al. 2014, 2015; Gavdush et al. 2019).

We derive the equations describing the complex amplitudes of the reference E_R(ν) and sample E_S(ν) spectra, assuming only wavelets inside the Tukey apodization. For the reference spectrum, we obtain

(B.10)

where the indexes 0 and Si correspond to the free space and the HRFZ-Si medium; l₀ and l_Si are the total length of the THz beam path and the thickness of the HRFZ-Si substrate, respectively; T_0, Si(ν) and T_Si, 0(ν) are the transmission coefficients for the respective interfaces (between the free space and the HRFZ-Si), defined by Eq. (B.8); and P₀(ν,z) and P_Si(ν,z) are operators describing the THz wave propagation in the free space and the HRFZ-Si, respectively, as given by the exponential factor in Eq. (B.9).

As shown in Fig. 3, for the sample spectrum we take into account the contribution of the ballistic THz pulse (1) and the satellite pulses (2 and 3), caused by the multiple THz wave reflection in the ice films. This yields the following equation:

(B.11)

where the summation terms in the brackets correspond to the wavelets 1, 2, and 3; l_CO, I and l_CO, II stand for thicknesses of the first and second ice films; R_CO, Si(ν) and R_CO, 0(ν) are the transmission coefficients for the respective interfaces, see Eq. (B.8); the definition of the remaining factors (P and T) is similar to that in Eq. (B.10). We point out that the complex amplitudes (E) and all factors (T, R, and P) in Eqs. (B.10) and (B.11) are frequency-dependent. Equations (B.10) and (B.11) form a basis for the reconstruction of the THz dielectric response of ices.

The reconstruction is performed via the following minimization procedure:

(B.12)

where argmin is an operator which determines the minimum argument of the vector error functional Φ. The latter is formed from the complex theoretical T_Th and experimental T_Exp transfer functions

(B.13)

where |…| and ϕ[…] are the absolute values and phases of the complex functions, respectively.

We define the theoretical transfer function T_Th as the sample spectrum (Eq. (B.11)) normalized by the reference spectrum (Eq. (B.10))

(B.14)

Considering all reflection and transmission operators in Eq. (B.9), we note that T_Th depends only on the refractive index of the HRFZ-Si substrate n_Si, which is known a priori and on the parameters of the CO ice to be determined; it excludes the contribution of several factors, such as the unknown complex amplitude of the TDS source E₀(ν), the unknown total length of the THz beam path l₀, and, finally, the thickness of the HRFZ-Si substrate l_Si, which is known, too, but can slightly vary owing to angular deviations of the substrate during the vacuum chamber assembling.

The experimental transfer function T_Exp is calculated in a similar manner, relying on the Fourier spectra of the experimental sample E_S and reference E_R waveforms (after applying the Tukey apodization),

(B.15)

We note that all the functions and operators in the theoretical and experimental transfer functions are frequency-dependent, and both transfer functions take into account only the ballistic pulses of the reference and sample waveforms, as well as the first and second ice-related satellite pulses of the sample waveform.

By introducing equal confidence intervals for the refractive index n′ and the amplitude absorption coefficient α, as and [0, α_max] with Δn′=0.25 and α_max = 15 cm⁻¹, we use the non-linear trust region approach (Coleman & Li 1996) to reconstruct the THz dielectric response of the CO ice in the spectral operation range of 0.3–2.0 THz.

Appendix C: Opacity model benchmark

In order to verify the correctness of the machinery employed to calculate the opacity from the dielectric constants, we reproduce the results found by Ossenkopf & Henning (1994) in their Fig. 5, panels a–c, compact grains case. We describe the methodology employed. The results described can be reproduced by running test_04.py from the publicly available code², while Fig. 7 from our paper can be reproduced with test_05.py.

C.1. Dielectric constants

The complex dielectric functions ε′−iε″ of the refractory components are taken from Ossenkopf et al. (1992, cool oxygen-rich silicates, their Fig. 10) and from Preibisch et al. (1993, amorphous carbon, their Table 1), while ice is assumed to be a H₂O:CH₃OH:CO:NH₃ = 100:10:1:1 mixture at 10 K from Hudgins et al. (1993, their Table 2A).

For the aims of this work we need to extrapolate the ice data relative to the H₂O:CH₃OH:CO:NH₃ = 100:10:1:1 mixture to longer wavelengths, as done by Ossenkopf & Henning (1994). First we fit the last 45 data points³ of the imaginary part (i.e. approximately 71 − 194 μm) with a f(λ)∝λ⁻¹ function (Ossenkopf, priv. comm.), and we use this to extrapolate ε″ with 200 linearly spaced wavelength points, in the range 71 − 800 μm. To retrieve the real part of ε at each ω = 2πcλ⁻¹ point, we apply the Kramers–Kronig (e.g., Bohren & Huffman 1983) relations in the discrete form

(C.1)

with Ω(ω) the discrete integral over the positive frequency ranges using the composite trapezoidal rule with the integrand , and excluding ω, where the denominator of the argument vanishes (i.e. the Cauchy principal value of the corresponding finite integral).

After the extrapolation, the ice dielectric functions are modified by mixing spherical inclusions of amorphous carbon using the Bruggeman effective medium approximation (Bruggeman 1935), with a volume filling fraction of 0.11 and 0.013 for the thin and the thick ice cases, respectively.

C.2. Absorption coefficients

The absorption efficiency Q_abs(λ, a) is function of the wavelength and grain size a, and computed with the routine BHMIE.PY for the bare grains⁴ and with BHCOAT.PY for the coated grains⁵, both from Bohren & Huffman (1983). In Ossenkopf & Henning (1994) there are three cases: bare grains without ice; thin ice, which has a volume ratio of V = 0.5; and thick ice, where V = 4.5. The radius of the refractory core a and the fraction V determine the radius of the mantle a_coat = a(V + 1)^1/3.

C.3. Opacities

With Q_abs it is possible to retrieve the opacities averaged on the grain size distribution φ(a) as

(C.2)

where to is the range in which the size distribution is valid including coating, and

(C.3)

where a_min to a_max is the range where the size distribution is valid, but considering only the refractory material, i.e. silicates or carbonaceous dust.

Ossenkopf & Henning (1994) assume a_min = 5 × 10⁻⁷ cm, a_max = 2.5 × 10⁻⁵ cm, ρ₀ = 2.9 g cm⁻³ (silicates) and ρ₀ = 2 g cm⁻³ (amorphous carbon), and φ(a)=a^−3.5.

The total opacity is κ_tot(λ)=0.678κ_Si(λ)+0.322κ_AC(λ), where the two terms are, respectively, the silicates and carbonaceous opacities including ice coating, and the two coefficients are calculated from the volume ratio discussed in Sect. 3.1 of Ossenkopf & Henning (1994). In particular, Ossenkopf & Henning assume a volume ratio of the refractory components V_AC/V_Si = 0.69 (their Sect. 3.1), which can be converted into the corresponding mass ratio M_AC/M_Si = 0.4758 using the relation M_i = ρ_iV_i, where ρ_i are the bulk densities of the two refractory components (being in Ossenkopf & Henning the opacity defined per unit mass of the refractory material), so that κ_tot = (M_ACκ_AC + M_Siκ_Si)/(M_AC + M_Si).

All Tables

All Figures

	Fig. 1. Sketch of the vacuum chamber of the cryostat coupled to the THz beam at the CAS.
In the text

Fig. 2. Set-up of THz-TDS for the spectroscopy of ices. Panel a: chematic of the set-up, where FS laser stands for the femtosecond laser, M stands for the optical mirrors, ATT stands for the attenuator of the laser beam intensity, BS stands for the optical beams splitter, DS stands for the mechanical double-pass delay stage, PCAE and PCAD stand for the photoconductive antenna-emitter and antenna-detector, respectively, L stands for the TPX lenses, and S and W stand for the HRFZ-Si substrate and windows, respectively. Panel b: spectra of THz waveforms E(ν) transmitted through the empty THz beam path or the THz beam path with the cryostat; the shaded region below ≈0.3 THz indicates the spectral range in which distortions from the THz beam diffraction on the aperture of the substrate are expected.

In the text

Fig. 3. Time–distance diagram illustrating the THz wave propagation through the HRFZ-Si substrate with ice films deposited on its surfaces. Lines 0–3 illustrate the ballistic pulse and satellite pulses transmitted in the direction of the antenna-detector; N stands for unaccounted satellites with larger time delays. The solid lines represent the pulses used for the analysis, while dotted lines correspond to neglected pulses.

In the text

Fig. 4. Calculation of the thicknesses lCO, I, lCO, II and the complex refractive index n of the CO ice films (see Fig. 3). Panel a: time delays between the ballistic THz pulses of the reference (0) and sample (1) waveforms, δt01; the first satellite pulse (2) and the ballistic pulse (1) of the sample waveform, δt12; and the second satellite pulse (3) and the ballistic pulse (1) of the sample waveform, δt13 (the pulses are delineated as in Fig. 3). Panels b and c: estimates for the thicknesses of the two ice films as a function of the total deposition time tdep for the different deposition intervals of Δtdep = 4, 5, and 6 min; the first assumption for the real part of the refractive index of ice is between and 1.255 ± 0.035; for the imaginary part we first set .

In the text

Fig. 5. Evolution of the THz pulse and its spectra during the CO ice deposition. Panel a: reference waveform E(t) transmitted through the cryostat with the empty substrate (green), and sample waveforms (black to light red) transmitted through the substrate with the CO ice deposited on its surfaces. Panel b: Fourier spectra |E(ν)| of the reference and sample THz waveforms calculated with the use of Tukey apodization. The waveforms in (a) and spectra in (b) correspond to different values of the total deposition time tdep (indicated); the deposition intervals are Δtdep = 6 min.

In the text

Fig. 6. Optical properties of CO ice. (a) Real part of the refractive index, (b) amplitude absorption coefficient, and (c) real and (d) imaginary parts of the dielectric permittivity (see Eqs. (1) and (2)). For all deposition intervals, the dielectric curves demonstrate the existence of a Lorenz-like absorption peak, centred near 1.5 THz and featuring similar bandwidth. Distortions of the results seen at frequencies below 0.3 THz (such as an oscillatory character of n and ε′ and an increase of ε″ with decreasing frequency) are due to diffraction effects (see Sect. 2.2).

In the text

Fig. 7. Calculated and reference opacities of astrophysical dust with CO ice and ice mixtures as a function of the wavelength. The dotted lines labelled with OH94 refer to bare grains and ice mixtures by Ossenkopf & Henning (1994), the dashed lines with BP98 to CO ice by Baratta & Palumbo (1998), and the solid lines to the CO data by the present work. The value V indicates the volume ratio of refractory core to ice mantle, for which we follow Ossenkopf & Henning (1994), where V = 0 (black) is the bare grain, V = 0.5 (blue) thin ice, and V = 4.5 (orange) thick ice. See text for additional details.

In the text