© The Authors 2024

1 Introduction

Water is one of the most abundantly detected interstellar molecules in both the gas and solid phases, and it plays countless essential roles in physical and chemical processes involved in the birth and evolution of stars, planets, and life as we know it (e.g., van Dishoeck et al. 2013, 2021). Despite this significance and ubiquity, its formation and evolution in the star formation cycle are still not fully understood.

Detections of abundant water ice toward molecular clouds and dense prestellar cores show that significant quantities of water can be produced in the solid state even before the earliest stages of star formation (Tanaka et al. 1990; Knez et al. 2005; Boogert et al. 2011, 2013; McClure et al. 2023). In these environments, water forms on the surfaces of silicate and carbonaceous interstellar grains via cold atom addition chemistry along with other simple species like CO₂, NH₃, and CH₄, resulting in the growth of ice layers that can be observed via infrared spectroscopy (Boogert et al. 2015 and references therein). There is mounting evidence that these water-rich prestellar ice mantles also serve as the birthplaces and reservoirs of many other more complex molecules, some of which could have prebiotic relevance (Herbst & van Dishoeck 2009; Jørgensen et al. 2020; Nazari et al. 2022, 2024; McClure et al. 2023; Chen et al. 2023; Rocha et al. 2024).

However, it is uncertain how much of this primordial water ice and its associated chemical complexity survive the subsequent intense thermal and radiative processes that occur upon the formation of protostars and protoplanetary disks. Once a protostellar source forms from a prestellar core and begins to generate radiation, it can cause its surrounding water ice to sublimate into the gas phase via photo- or thermal desorption, where it may be destroyed and subsequently reformed via hot gas-phase neutral-neutral reactions (Bergin et al. 1998; Harada et al. 2010; Visser et al. 2011; Furuya et al. 2013; van Dishoeck et al. 2014; Owen & Jacquet 2015). The extent to which such processing of the prestellar water ice occurs within the hot radiation-dominated regions of protostellar envelopes and pro-toplanetary disks has important implications not only for the physicochemical processes that drive their evolution, but also for the chemical inventory available to the planetesimals that eventually begin to form and grow in these environments (i.e., the question of inheritance versus reset, van Dishoeck 2017; Öberg & Bergin 2021).

Quantifying the deuterium fractions (D/H ratios) of water offers a promising path to understanding the extent of this processing because deuterium fractionation reactions are extremely sensitive to the physicochemical conditions in which they occur. Therefore, a molecule’s D/H ratios reflect the environment in which it formed and existed. In molecular clouds and prestel-lar cores, molecules are thought to be deuterium enriched due to an increase in the relative abundance of available atomic D (with respect to atomic H) caused by the enhanced gas-phase formation of H₂D⁺ (Roberts et al. 2003):

$${\text{H}}_3{}^{+}+\text{HD}\rightleftharpoons {\text{H}}_2{\text{D}}^{+}+{\text{H}}_2+232 \text{K}.$$

At the low temperatures relevant to ice formation (T < 50 K), this reaction mostly proceeds in the forward direction due to its energy barrier. The following dissociative recombination reaction can then free D atoms, increasing deuterium atoms that can land on cold dust grains and react to form a higher abundance of deuterated molecules (Tielens 1983):

$${\text{H}}_2{\text{D}}^{+}+{\text{e}}^{-}\to {\text{H}}_2+\text{D}.$$

Water formed in cold prestellar environments is therefore expected to be deuterium enriched relative to bulk interstellar and solar D/H values (2.0 × 10⁻⁵, corresponding to an HDO/H₂O ratio of 4 × 10⁻⁵ if there is no fractionation, Prodanović et al. 2010; Geiss & Gloeckler 2003). The relative abundance of H₂D⁺ is further enhanced in the absence of neutral species like O and CO that destroy H₃⁺ and its isotopologs, so deuterium enrichment is expected to be even greater in the molecules formed during the coldest (<20 K) and densest stages of prestellar cores, when O and CO are depleted from the gas phase as they freeze out on the grains.

Interstellar HDO was first detected in the gas phase in the massive star-forming region Orion KL (Turner et al. 1975). The relatively high brightness of the hot (>100 K) inner cores of such MYSOs allows multiple lines of HDO to be detected relatively easily with millimeter spectroscopy. The detection of these lines, along with the lines of ${\text{H}}_2^{18}\text{O}$, a water isotopolog that is usually optically thin, has allowed the determination of the gas HDO/H₂O ratios toward many of these objects in the last few decades. These ratios typically range between 10⁻⁴ and a few 10⁻³ (Helmich et al. 1996; van der Tak et al. 2006; Neill et al. 2013; Emprechtinger et al. 2013; Coutens et al. 2014b), an enhancement of up to two orders of magnitude relative to the cosmic standard D/H. Highly sensitive millimeter interferometers also provide reliable measurement of gas HDO/H₂O ratios in the hot inner regions of the much dimmer low-mass young stellar objects (LYSOs), whose ratios see a similar enhancement (Persson et al. 2013, 2014; Jensen et al. 2019, 2021).

This enhancement has been interpreted as evidence that at least some of the water gas detected in protostars has its origins in prestellar ices, in line with the assumption that many of the gas-phase species detected in protostellar hot cores and corinos come directly from ices via thermal desorption. Additionally, the gas HDO/H₂O ratios of clustered LYSOs are remarkably similar to those measured in comets, suggesting that prestellar water ice and its associated chemical complexity may be largely inherited by the planetesimals that grow in the outer protoplanetary disk (Persson et al. 2014).

While studying HDO/H₂O ratios in the gas phase has clearly improved our understanding of the chemical history of water throughout star formation across a range of stellar masses, the glaring question remains of whether the HDO/H₂O ratios measured in protostellar hot gas are actually representative of the protostellar HDO/H₂O ice ratios, or if the detected water gas has already experienced substantial alteration by gas-phase reactions. To answer this question, the HDO/H₂O ratio in ices must be probed directly.

Although H₂O is the most abundant observable hydrogen-bearing ice by a large margin, a secure detection of its deuterated isotopologs is an exceptional observational challenge. Because it begins to form in the earliest stages of cloud formation, its D/H ratio in the gas phase is one to two orders of magnitude below that of molecules like CH₃OH and H₂CO that are thought to form only in the later colder stages of cloud formation (Cazaux et al. 2011; Taquet et al. 2012; Caselli & Ceccarelli 2012; Ceccarelli et al. 2014; Furuya et al. 2016), so the column density of HDO ice is expected to be very low. Furthermore, the strongest IR band of HDO, its O–D stretching mode at ~4.1 µm (~2440 cm⁻¹), lies next to the short-wavelength wing of the strong asymmetric C=O stretching mode of CO₂ ice, which sometimes warps the continuum due to grain shape effects and complicates the detection and quantification of spectrally adjacent species (Dartois et al. 2022, 2024). In the amorphous state, the O–D stretch of HDO is also broad and relatively weak (Dartois et al. 2003; Gálvez et al. 2011).

The first tentative HDO ice detection at 4.1 µm was reported in data from the Infrared Space Observatory (ISO) more than two decades ago by Teixeira et al. (1999), but the detection was soon called into question by Dartois et al. (2003), who pointed out, among several other concerns, that the observed absorption could also be attributed to a weak CH₃OH combination mode whose wavelength overlaps with that of the HDO O–D stretch. They went on to report HDO upper limits in the range 0.16–1.1% with respect to H₂O toward four intermediate-mass young stellar objects (IMYSOs) and MYSOs with low CH₃OH column densities. Parise et al. (2003) expanded the search for HDO ice to low-mass young stellar objects (LYSOs) and also derived upper limits, ranging from 0.5 to 2% with respect to H₂O. Almost a decade later, Aikawa et al. (2012) reported tentative detections and measured very high HDO/H₂O ratios (2–22%) toward four LYSOs, but they cautioned that the HDO ice models did not provide a perfect fit to the observations, and their analysis did not take into account contribution from CH₃OH ice to the HDO ice spectral region. Furthermore, given the recent characterization of how grain growth affects the continuum around the CO₂ ice feature (Dartois et al. 2022), it is possible that some of the observed “broad absorption” that was attributed to HDO ice by Aikawa et al. (2012) was caused by a warping of the local continuum via grain shape effects, particularly in the case of IRAS 04302+2247.

It is evident that the infrared telescopes capable of observing ices in the 4 µm range thus far have lacked sufficient sensitivity to conclusively detect solid-state HDO. The unprecedented sensitivity of the Near Infrared Spectrograph (NIRSpec) instrument on the recently launched James Webb Space Telescope (JWST) provides a prime opportunity to revisit the search for interstellar HDO ice.

Here we report HDO ice detections toward two protostars, HOPS 370, and IRAS 20126+4104 (hereafter IRAS 20126), observed with JWST’s NIRSpec integral field unit (IFUs) as part of the GO program 1802, Investigating Protostellar Accretion (IPA) across the mass spectrum (Program ID 1802, PI Tom Megeath; Megeath et al. 2021). Care is taken to ensure that the observed HDO ice detections cannot be attributed to the CH₃OH combination mode at 4 µm by establishing a CH₃OH column density via fitting high-resolution laboratory spectra to the 3.53 µm band. Additional high-resolution laboratory spectra of HDO are collected to allow the fitting of the observational spectra and subsequently derive HDO column densities. The laboratory spectra measured and used in this work are made publicly available via the Leiden Ice Database (LIDA, Rocha et al. 2022). Comparisons are made between the derived HDO ice abundances and those measured in the gas phase in the interstellar medium. The resulting implications for the chemical evolution of water throughout the star formation process are then discussed.

2 Methods

2.1 Investigated sources

The IPA program observed five YSOs in their primary accretion phase, with their bolometric luminosities spanning approximately six orders of magnitude, using the NIRSpec IFU (Federman et al. 2024). Of these sources, two are presently analyzed due to their high S/N and because their IFUs contain lines of sight that have a nonsaturated 3 µm O–H stretching feature of H₂O ice, allowing quantification of its column density and characterization of its morphology. Table 1 presents an overview of these sources.

HOPS 370 is an IMYSO with a spectral energy distribution (SED) near the Class 0/Class I border located within the very clustered OMC-2 star-forming region in the Orion A molecular cloud. It is likely currently experiencing a high accretion rate, and its compact complex organic molecule (COM) emission is similar to those of hot corinos (Tobin et al. 2020b). It is driving a powerful bipolar outflow observed in the far IR, millimeter, and centimeter wavelength ranges (Shimajiri et al. 2008; González-García et al. 2016; Osorio et al. 2017; Sato et al. 2023).

IRAS 20126 is a deeply embedded MYSO, the most massive and luminous protostar of the sample. The source contains a hot and massive disk (Chen et al. 2016) as well as jet-driven outflows (Caratti o Garatti et al. 2008). Although initial observations suggested that its environment is unusually isolated for a MYSO (Qiu et al. 2008), later observations revealed the presence of a young cluster within a 1.2 pc vicinity (Montes et al. 2015). Two likely companion YSOs are present within the 6″ × 6″ FOV in the NIRSpec IFU of this source (Federman et al. 2024).

2.2 Observations and data reduction

The investigated spectra were collected with the NIRSpec IFU using the G395M grating (R = λ/∆λ ~700–1300, Rubinstein et al. 2023), providing data from 2.9-5.1 µm (3450–1960 cm⁻¹) with no detector gaps, unlike the higher-resolution G395H grating, which has a detector gap in the IFU at 4.1 µm, where the strongest HDO IR band is found. A 4-point dither pattern was used to obtain a 6″×6″ mosaic with a 0.2″ spatial resolution. Detailed descriptions of the data reduction are provided in Federman et al. (2024) and Narang et al. (2024).

The regions in the IFUs from which spectra were extracted for analysis were chosen by searching through the datacube for the highest S/N lines of sight in which the H₂O 3 µm feature is not saturated or extincted. In both cases, the ice features in the selected positions are from the cold envelopes seen in projection in front of bright scattered light emission from the protostars. The HOPS 370 spectrum was extracted near the source center over a knot of shocked gas. Spectra extracted near the central position of IRAS 20126 cannot be used for the analysis in this work because the flux at 3 µm at these coordinates is so low that the profile of the H₂O feature is contaminated by an extended weak emission that may originate from a photodisso-ciation region along the line of sight. Extended weak emission is present toward all of the IPA sources and is seen most distinctly at the edges of the datacubes, where the source continuum fluxes are minimal and a weak PAH emission feature at 3.3 µm is clearly observed. The analyzed spectrum of IRAS 20126 is therefore extracted from a bright position toward an outflow cavity, where the PAH emission flux is negligible relative to the overall observed flux.

The coordinates of the apertures used to extract the spectra are provided in Table 1 and are shown on the IFU continuum maps in Fig. 1. A wavelength-dependent extraction aperture of 7 times the diffraction-limited radius was used, which results in an aperture that is >3 times the full width at half maximum (FWHM) of the NIRSpec IFU point spread function over the entire spectrum. Because the investigated sources have small radial velocities relative to the spectral resolution, the final extracted spectra are not corrected for source velocity (Brunken et al. 2024).

2.3 Gas emission subtraction

Both spectra contain gas emission features, including those of H₂, H I, [Fe II], and CO (Federman et al. 2024; Narang et al. 2024; Rubinstein et al. 2023). The strongest of these features were subtracted out of the spectral regions of interest (3.33– 4.18 µm) via Gaussian fits to prevent their interference with the fitting functions used to analyze the ice features. Most of the gas emission lines are well resolved from other neighboring gas lines and thus required only single Gaussian profiles to be subtracted out. The few emission lines that are blended were subtracted out using multiple Gaussian profiles.

In all subsequent spectral figures, the original spectra before gas emission subtraction are plotted in gray, and the emission-subtracted spectra are overplotted in color.

Fig. 1 5 µm continuum maps with the aperture used to extract the analyzed spectra at 5 µm shown in white (7 × diffraction-limited radius at 5 µm = 0.678″).

2.4 Laboratory data collection

High-resolution (0.5 cm⁻¹) spectra of various laboratory ice mixtures containing HDO, H₂O, CO, and CH₃OH ranging from temperatures of 15–150 K were collected on the InfraRed Absorption Setup for Ice Spectroscopy (IRASIS) in the Leiden Laboratory for Astrophysics for use in fitting the observed spectra. A schematic of the setup can be found in Rachid et al. (2021), and recent upgrades to the setup are described in Slavicinska et al. (2023). The method used to create ice mixtures with specific ratios is similar to that described in Yarnall & Hudson (2022). The chemical compositions and mixing ratios of the laboratory ice mixtures are presented in Table 2. Additional experimental details can be found in Appendix A.

All of the spectra are made available for public use and enjoyment via LIDA (Rocha et al. 2022). As much of this work focuses on fitting very weak features (i.e., the O–D stretch of highly diluted HDO at 4.1 µm and the combination modes of CH₃OH at 3.9 µm), the strongest peak of the primary ice matrix component (i.e., the C≡O stretching mode of CO or the O–H stretching mode of H₂O) was allowed to saturate in three of the experiments so that the weaker ice features have sufficient S/N for accurate fitting. The spectra with such a saturated peak are indicated with an asterisk (*) in Table 2. We strongly warn against using the saturated peaks in these spectra to fit observational data.

The peak positions, FWHMs, and apparent band strengths of the O–D stretching mode from 15 to 150 K were extracted from the HDO:H₂O mixtures for both the amorphous and crystalline phases and are presented in Table A.2 and Figs. A.2 and A.3. As discussed in detail in Sect. 3.3, the variation in amorphous HDO band strength with respect to temperature has particularly important ramifications on the derived amorphous HDO abundances. The peak positions of the CH₃OH 3.53 µm band from 15 to 150 K in all of the CH₃OH-containing spectra are presented in Tables A.3 and A.4 and Fig. A.4.

3 Results

The full extracted spectra of both sources are presented in Fig. 2. At first glance, it is apparent that in both spectra, the 3 µm O–H stretching modes of H₂O ice have sharp peaks at 3.1 µm, characteristic of a crystalline component (Hagen et al. 1981; Smith et al. 1989). This indicates that the ices along these lines of sight have experienced enough thermal processing to crystallize some of the water ice. The presence of pure CO₂ components in the ¹³CO₂ bands (Brunken et al. 2024), enhanced OCN⁻ abundances (Nazari et al. 2024), and the relatively low absorptions of the CO ice features at 4.67 µm relative to the CO₂ features at 4.27 µm further corroborate that a significant portion of the ices toward these sources have experienced a substantial degree of thermal processing.

The ice species of interest to this work were analyzed as follows: first, the CH₃OH ice column densities and physico-chemical environments were determined from the ~3.53 µm (~2830 cm⁻¹) C–H stretching mode (Sect. 3.1). This was done prior to the analysis of the HDO ice because CH₃OH has several weak combination modes in the range ~3.7–4.2 µm (~2700– 2400 cm⁻¹), one of which overlaps with the HDO O–D stretching mode at 4.1 µm; therefore, pre-constraining the expected contribution of CH₃OH ice to this region via a column density derived from a stronger band in a different region is necessary to ensure that any CH₃OH absorption at 4 µm is not misattributed to HDO ice (Sect. 3.2), a concern previously raised by Dartois et al. (2003). The HDO ice column densities and morphologies were then determined via fitting the observed absorptions in the ~3.7-4.2 µm region (Sect. 3.3). Finally, the 3 µm (~3300 cm⁻¹) H₂O O–H stretching mode was fit to corroborate the morphologies found for HDO ice and to determine the ice HDO/H₂O ratios toward these sources (Sect. 3.4).

Fig. 2 Extracted spectra of the investigated sources. Left: full NIRSpec spectra with major ice features labeled. The original spectra before gas emission subtraction are plotted in gray. The spectral region of the strongest HDO absorption, the O–D stretching mode, is indicated in magenta. The spectral region plotted on the right is indicated in gray. Right: zoomed-in view of the spectral region of interest for HDO ice.

3.1 CH₃OH ice

Both sources present a clear CH₃OH C–H stretching mode at 3.53 µm. The CH₃OH ice column density can be obtained from this mode via a procedure outlined in Brooke et al. (1999) and Boogert et al. (2022). Briefly, this involves subtracting a local continuum using data points ~3.25–3.3 and 3.6–3.7 µm and then fitting a Gaussian representing ammonia hydrates with a peak position of 3.47 µm (2881.844 cm⁻¹) and a full-width half maximum (FWHM) of 0.11 µm (91.4 cm⁻¹) along with laboratory CH₃OH ice spectra to the extracted feature. An identical local continuum subtraction is performed on the laboratory data as well prior to fitting the observed data to isolate the C–H stretching modes from the strong and broad O–H stretching feature with which they overlap. The resulting CH₃OH fits, along with the local continuum subtractions used to extract the CH₃OH features into the optical depth scale, are shown in Fig. 3.

Toward both sources, the peak position of the 3.53 µm band is most consistent with laboratory data of cold (<60 K), amorphous, relatively pure CH₃OH ice (Figs. 4 and A.4). The 3.53 µm bands of laboratory CH₃OH ices that is warm, crystalline, or diluted with any of the three most abundant observable interstellar ice components (H₂O, CO, and CO₂) are too blueshifted to fit the observed absorption. However, fitting the entire 3.5 µm region with only pure amorphous CH₃OH ice results in the spectrum being overfit at 3.4 µm by the blended CH₃OH overtones and C–H asymmetric stretching modes.

This issue of pure CH₃OH ice spectra causing overfitting at 3.4 µm has previously been documented toward some LYSOs and MYSOs and was solved by using spectra of CH₃OH mixed with H₂O (Brooke et al. 1999; Pontoppidan et al. 2003), which decreases the optical depth of the 3.4 µm features relative to that of the 3.53 µm feature. Such a decrease is not observed in CO- or CO₂-rich CH₃OH mixtures or in crystalline CH₃OH ices. Using a combination of two laboratory amorphous CH₃OH spectra, one of pure CH₃OH ice and an H₂O-rich CH₃OH ice mixture (1:5 CH₃OH:H₂O), in the fitting procedure results in satisfactory fits of both the 3.53 µm band and the 3.4 µm region (Fig. 3). A similar combination of pure CH₃OH and CH₃OH:H₂O laboratory spectra was used by Dartois et al. (1999) to fit the 3.53 and 3.9 µm features in the spectra observed toward MYSOs W33A and RAFGL 7009s.

The resulting CH₃OH column densities (reported in Table 3) are obtained using the band strength of 4.86 × 10⁻¹⁸ cm mol⁻¹ for the C–H stretching mode of pure amorphous CH₃OH at 20 K (Luna et al. 2018). This band strength was chosen because the profile of the 3.53 µm band indicates that the majority of the detected CH₃OH ice along these lines of sight is pure, amorphous, and cold (>60 K). The ratio of pure CH₃OH ice to CH₃OH ice mixed with H₂O that results in the best fits is highly dependent on the dilution factor of the CH₃OH:H₂O mixture used, so separate column densities for each of the two CH₃OH components were not calculated. Any effects of dilution with H₂O on the band strength of this feature are expected to be well within the 20% uncertainty reported for the pure band strength (Kerkhof et al. 1999; Luna et al. 2018). The reported uncertainties are explained in Appendix C.1. The implications of the nondetection of warm or crystalline CH₃OH ice toward these lines of sight are further discussed in Sect. 4.3.2.

Fig. 3 Extraction and fitting of the CH3OH 3.53 µm band. Left: adopted local continua. Right: fits to the extracted optical depth spectra. Some of the substructures present at 3.4 µm in the high S/N spectrum of IRAS 20126 are well reproduced by the more prominent CH3OH overtones and C-H asymmetric stretching modes that emerge in H2O-rich CH3OH spectra.

Fig. 4 Observed CH3OH C–H symmetric stretching modes compared to select laboratory data. The peaks were extracted from the 3.5 µm absorption complexes via a local continuum fit using baseline points in the range ~3.47–3.50 and 3.60–3.62 µm (Fig. 3). It is clear from these comparisons that cold, amorphous, and relatively pure CH3OH is the major contributor to this feature. Left: observed peaks vs. laboratory spectra of CH3OH ice in various chemical environments. All spectra were collected at 15 K except the CH3OH:CO2 spectrum, which was collected at 10 K (Ehrenfreund et al. 1999). Right: observed peaks vs. laboratory spectra of pure CH3OH ice at various temperatures.

Fig. 5 Extraction and fitting of the 4 µm region. Left: adopted local continua, along with CH3OH spectra determined from the 3.53 µm band fitting used to aid the continuum placement. The spectral region of the HDO O–D stretching mode is indicated in magenta. Right: fits to the extracted optical depth spectra (am. = amorphous; cr. = crystalline). The shaded gray regions indicate 3σ uncertainty levels. The same combination of pure and H2O-rich CH3OH spectra as used to fit the 3.53 µm feature is used in this fit. The best-fitting temperatures of the amorphous and crystalline HDO ice components are 132 K and 47 K for HOPS 370 and 138 K and 90 K for IRAS 20126, respectively (see Sect. 3.3 and Appendix B for more details).

3.2 4 µm region: CH₃OH, S-bearing species, and HDO

A local continuum in the 3.7–4.2 µm (2700–2380 cm⁻¹) region was adopted so that the optical depth of the CH₃OH combination modes predicted by the CH₃OH column density derived from the 3.53 µm band did not exceed the optical depth of the continuum-subtracted spectrum, providing a constraint to prevent oversubtraction of features due to continuum choice in this region. The RMS error in optical depth was calculated for this region using data from 3.70–3.77 µm (2700–2650 cm⁻¹). It is important to note that, as can be seen in the right panel of Fig. 2, the warping of the continuum by the scattering wings of the strong 4.27 µm CO₂ feature is very minimal in both spectra in comparison to some other icy lines of sight (e.g., Dartois et al. 2022, 2024), greatly reducing the uncertainty of the local continuum placement in this region. The adopted local continua, the subtracted spectra in optical depth units, and the 3σ uncertainty levels from the RMS error are presented in Fig. 5.

Following this procedure, both lines of sight have excess absorption on top of the CH₃OH feature at ~3.9 µm (~2550 cm⁻¹). Jiménez-Escobar & Muñoz Caro (2011) noted a similar excess absorption in the spectrum of the MYSO W33A and used it to derive an H2S ice upper limit. However, Hudson & Gerakines (2018) recently pointed out that many other S–H bond-containing ices, like CH₃SH and CH₃CH₂SH, have strong absorption features in this spectral region with nearly identical peak profiles, making it difficult to conclusively assign the broad excess absorption in this region to a specific molecule. Additionally, there is a lack of publicly available laboratory ice spectra of many of these sulfur-bearing species in astrophysically relevant mixtures. As sulfur-bearing ices are not the focus of this paper and the fitting of this feature has little effect on the analysis of the HDO ice feature, we simply fit this excess with a Gaussian and leave further analysis to a future work.

A small, narrow excess absorption is also present in both lines of sight at ~3.83 µm (~2611 cm⁻¹), in addition to a sharp peak toward IRAS 20126 at ~4.07 µm (~2459 cm⁻¹). We initially suspected these features could be attributed to crystalline CH₃OH ice, which has two sharp peaks at 3.83 and 4.07 µm (2612 and 2455 cm⁻¹) in the pure form. However, the entire CH₃OH combination mode complexes over the full 3.7– 4.2 µm range could not be satisfactorily fit with any available crystalline CH₃OH ice spectra, pure or mixed. The observed spectra are missing an additional sharp and distinct feature at 3.94 µm (2536 cm⁻¹) that is present in the crystalline CH₃OH laboratory spectra. Annealed crystalline CH₃OH spectra also did not provide a good fit. Furthermore, the 3.53 µm feature lacks any significant spectral contribution from crystalline CH₃OH ice. Because no other ice spectra available to us provided a convincing match, we refrain from assigning these features, although we do not completely exclude crystalline CH₃OH as a candidate carrier.

Finally, in both lines of sight, excess absorption emerges from the continuum in the range ~4.04–4.17 µm (~2475– 2400 cm⁻¹) which cannot be accounted for by the CH₃OH combination mode without inducing unrealistic inflections in the local continua or overfitting the CH₃OH combination modes in the range 3.8–3.9 µm. There are no currently known S–H bond-bearing species whose absorptions peak at these wavelengths (Hudson & Gerakines 2018), and the features present in crystalline CH₃OH ice spectra in this spectral region are all too blueshifted from the observed peak positions at ~4.14 µm. We therefore assign these features to HDO ice in various morphological states.

3.3 HDO ice

We performed fits to the local continuum-subtracted 4.04– 4.17 µm region using a combination of CH₃OH ice laboratory spectra, a Gaussian representing S–H bond-bearing species, and HDO ice laboratory spectra. Prior to the fitting, the profiles of the HDO O–D stretching features were isolated from the underlying broad H₂O combination modes in the 0.4% HDO:H₂O ice spectra (see Fig. A.1) and subsequently smoothed with Gaussian fits to eliminate experimental noise. The optical depths of the CH₃OH combination modes (as constrained by the CH₃OH column density derived from the 3.53 µm band in Sect. 3.1) and the Gaussian fit to the 3.9 µm excess were not varied in the fitting procedure. The same combinations of pure and H₂O-rich CH₃OH ice spectra that were used to fit the 3.53 µm feature were also used here.

Toward both sources, a minimum of two HDO profiles, one amorphous and one crystalline, is needed to fit the HDO feature due to its asymmetric profile. The 3 µm H₂O bands toward both sources also have profiles characteristic of the lines of sight containing both amorphous and crystalline H₂O ice (Sect. 3.4). This correlation between the detected presence of crystalline H₂O and crystalline HDO provides additional confidence in the peak assignment, because if some regions of the observed ice envelopes have been thermally processed enough to crystallize the H₂O ice, then any HDO ice in those regions should have crystallized as well.

In total, 84 amorphous HDO ice spectra in the range 15–140 K and 22 crystalline HDO ice spectra in the range 15–150 K were collected. Every possible combination of one amorphous profile and one crystalline profile out of these spectra were fit to the observed spectra using a least-squares fitting procedure. The least-squares fits that result in the lowest χ² value for both lines of sight are presented in Fig. 5. The relationships between the fit χ² values and the temperatures of the laboratory spectra are shown in Fig. B.1.

It is important to emphasize that the temperatures resulting in the lowest χ² values should not be interpreted as the actual temperatures of the observed ices for four reasons. First, laboratory ice temperatures do not correspond directly to interstellar ice temperatures due to the faster heating rates and higher pressures experienced by laboratory ices (Redhead 1962; Minissale et al. 2022). This means that, in the laboratory, ices crystallize and desorb at higher temperatures than in the interstellar medium. Second, these fits are not guaranteed to be unique, especially given that the broad, weak profile of the amorphous HDO could be considered degenerate at a wide range of temperatures when accounting for fitting errors caused by factors such as the observed spectral noise and uncertainties in the local continuum choice. This is particularly clear from Fig. B.1, which shows that a wide range of HDO temperatures produce fits that are degenerate within 3σ. Third, it is possible that much of the HDO ice along these lines of sight is located in different chemical environments than just pure H₂O ice, which could further alter the HDO band profile. For example, models from Taquet et al. (2014) and Furuya et al. (2016) suggest that the majority of HDO ice is co-produced with CH₃OH ice after the heavy CO freeze-out stage. Changes in the amorphous HDO feature caused by such chemical environments are difficult to unambiguously characterize in the laboratory due to the spectral overlap between the main HDO band and the CH₃OH combination modes. However, such chemical environments would not be expected to affect the crystalline component because CH₃OH and CO both desorb at temperatures lower than the water ice crystallization temperature, so by the time HDO ice crystallizes, it is expected to have been effectively distilled by heating. And finally, fitting two laboratory spectra at discrete temperatures to the observations is merely an approximation of the expected continuous thermal gradient experienced by the observed ices along each line of sight (Fig. 6).

Nonetheless, such an approximation is still useful for evaluating the thermal structure of the ice envelope qualitatively (Smith et al. 1989). In particular, the presence of a crystalline HDO component indicates that at least some regions within the ice envelopes must have experienced temperatures high enough to crystallize water ice, while the presence of an amorphous component indicates that other regions of the ice envelopes must have remained colder. Furthermore, the peak position of the crystalline HDO component, which is affected by much less uncertainties due to its sharper profile, matches best with annealed laboratory spectra, hinting at a possibility of heated ices being recooled within the ice envelope (see further discussion in Sect. 4.3.1).

The uncertainty in the temperature of the observed ices translates to an uncertainty in any HDO ice column density calculations because the apparent band strength of the O–D stretching feature in HDO measured in the laboratory varies with temperature, particularly for the amorphous phase (see Fig. A.2). Between 15 and 135 K, the amorphous HDO apparent band strength increases by a factor ~1.7 from 4.2 × 10⁻¹⁷ to 7.1 × 10⁻¹⁷ cm mol⁻¹. This uncertainty is much smaller for the crystalline HDO component because its apparent band strength only varies by ~8% in the range 15–150 K. When calculating the HDO column densities from the fits, we used band strengths of 6.0 × 10⁻¹⁷ and 6.7 × 10⁻¹⁷ cm molec⁻¹ for the amorphous and crystalline HDO ice features, respectively, which correspond to intermediate laboratory ice temperatures of around 80 K for the amorphous component and 90 K for the crystalline component (see Table A.2). The resulting HDO ice column densities are listed in Table 4, and the reported uncertainties are explained in Appendix C.2.

Within the uncertainties, the ratio of amorphous:crystalline HDO ice column densities are calculated to be in the range ~ 1.2– 4 for HOPS 370 and in the range ~ 1.2–6 for IRAS 20126 (i.e., there is more amorphous than crystalline HDO ice toward both sources). However, it must be emphasized once more that the uncertainties of the assignment, fits, and calculated column densities are greater for the amorphous HDO component because the amorphous profile is significantly weaker and broader and its apparent band strengths vary to a much greater extent with temperature, so more precise ratios of amorphous:crystalline HDO ice cannot be determined.

Fig. 6 Simplified diagram showing the difference between the expected thermal profile of the observed ice columns and the two-component fit approximation used to fit the observations in this work.

3.4 H₂O ice

To estimate the ice HDO/H₂O ratios along these lines of sight, we quantified the H₂O ice column density from the 3 µm band by defining global continua over the full NIRSpec G395M range and subsequently fitting the 3 µm band with laboratory H₂O ice spectra. The global continua and least-squares fits of the 3 µm feature are presented in Fig. 7.

An important caveat in the analysis of this feature is the lack of collected continuum data points below 2.9 µm, resulting in greater uncertainties in the slope of the blue wing and optical depth of the 3 µm feature in the global continuum-subtracted spectra. We attempted to mitigate this issue by using spectra of other MYSO ice envelopes as templates to guide our global continuum choice. These spectra contain continuum data points < 2.9 µ m, making their global continuum determinations and subsequent H₂O profile extractions much more reliable, and their H₂O profiles are similar to those in the spectra investigated in this work. The best-matching H₂O profiles found for HOPS 370 and IRAS 20126 were those of Orion from ISO (shown in Dartois & d’Hendecourt 2001) and G034.7123-00.5946 from IRTF (shown in Boogert et al. 2022), respectively. The continuum points on the blue side of the 3 µm band were then defined so that the profile of the 3 µm feature in the JWST spectra replicated the profile in the template spectra as closely as possible (see Fig. 7).

As done in previous studies (Smith et al. 1989; Dartois & d’Hendecourt 2001), the H₂O band in the optical depth scale was then fit with multiple components as an approximation to the thermal gradient expected along the line of sight. The fitting procedure was analogous to that performed on the 4 µm region in that the best-fitting combination of one amorphous and one crystalline H₂O lab spectrum was determined by calculating the χ² values of the least-squares fits of every possible combination. As stressed in Sect. 3.3, the temperatures of the lab H₂O ices fit to the observations should not be taken as indications of the real physical temperature of the observed ices, as these two-component fits are only simple approximations of the thermal gradient that is expected along these lines of sight, and they are not guaranteed to be unique; however, it is certain that some combination of cold and warm ices is needed to sufficiently fit the feature toward both sources. Furthermore, similarly to the analysis performed by the ICE AGE team in McClure et al. (2023), grain shape effects were not modeled in these fits, as a detailed characterization of the H₂O 3 µm band is outside our primary scope (quantifying H₂O to obtain HDO/H₂O ratios). Indeed, recently published follow-up analyses of grain shape effects in the ICE AGE spectra show that their derived relative ice abundances are largely the same regardless of whether their spectral fitting procedure includes grain shape effects (Dartois et al. 2024).

The contribution of the O–H stretching mode of CH₃OH ice to the observed 3 µm feature was also approximated similarly to the HDO fitting procedure, where the column density of CH₃OH was set as a constant value from the fits to the 3.53 µm band. However, a spectrum of only pure CH₃OH was used to model the contribution of CH₃OH to the 3 µm feature, as it is difficult to accurately deconvolve the overlapping O–H stretching modes of CH₃OH and H₂O in the CH₃OH:H₂O lab spectra. Contributions to the red wing by ammonia hydrates were neglected. For this reason (as well as the neglect of grain shape effects), only the region from 2.88–3.13 µm (3470–3191 cm⁻¹) was considered in the fitting procedure, and the long-wavelength wing of the 3 µm feature remains underfit. This underfitting is taken into account in our reported uncertainties (see Appendix C.3).

Similar to the band strength of the HDO O–D stretch, the band strength of the H₂O O–H stretch varies with temperature. The reported magnitude of this variation differs in the literature (Hagen et al. 1981; Gerakines et al. 1995; Mastrapa et al. 2009). In our experiments, the integrated absorbance of the O–H stretch increases by a factor of up to 1.26 between 15 and 135 K in amorphous H₂O, while there is only a factor of 1.16 difference between the integrated absorbance at 150 and 15 K in crystalline H₂O. These variations are again treated as a source of uncertainty in the reported H₂O column densities (see Appendix C.3).

To quantify H₂O, all of the fit H₂O components were integrated, and a band strength of 1.9 × 10⁻¹⁶ cm mol⁻¹ (Mastrapa et al. 2009) was used for the amorphous component. For the crystalline component, a band strength of 2.5 × 10⁻¹⁶ cm mol⁻¹ was used, which is the apparent band strength we calculated for the crystalline H₂O O–H stretch below 120 K in our laboratory spectra. The resulting values are reported in Table 5. Within the error margins, the ratios of amorphous:crystalline H₂O ice column densities range from 3.8–5.8 toward HOPS 370 and 5.7–9.1 toward IRAS 20126. These amorphous:crystalline ratios are higher than those of HDO but remain in agreement when considering uncertainties. We refrain from drawing chemical conclusions from the comparison between amorphous and crystalline column densities of HDO and H₂O ice for reasons discussed in Sect. 4.1.

Fig. 7 Extraction and fitting of the H2O 3 µm band. Left: adopted global continua and the archival MYSO ice spectra used as templates to guide the continuum choice. Right: fits to the extracted optical depth spectra (am. = amorphous; cr. = crystalline).

4 Discussion

4.1 Inferred ice HDO/H₂O ratios

The ice HDO/H₂O ratios measured toward HOPS 370 and IRAS 20126, 4.6 ± 2.2 × 10⁻³ and 2.6 ± 1.4 × 10⁻³, respectively, are both two orders of magnitude greater than what is expected from the cosmic standard D/H abundance. This enhancement far exceeds our reported error margins and is consistent with the deuterated water ice that we detect along these lines of sight having cold, prestellar origins. The ice HDO/H₂O of the IMYSO is higher by almost a factor of 2 compared to the MYSO, but this difference is not significant within the error margins of both values. A larger sample size of ice detections in IMYSOs and MYSOs is needed to draw concrete conclusions regarding any systematic differences between IMYSO and MYSO ice HDO/H₂O ratios.

While the ice HDO/H₂O measured toward these protostars are expected to largely reflect their prestellar conditions and timescales, experimental works have shown that it is also possible for the prestellar D/H ratios of protostellar ices to be altered prior to sublimation via inter-species deuterium exchange reactions during ice thermal processing. Specifically, deuterium exchange reactions proceed very efficiently in warm ices (T≳70– 120 K) between atoms involved in hydrogen bonds in a variety of astrophysically relevant simple molecules (Kawanowa et al. 2004; Ratajczak et al. 2009; Faure et al. 2015; Lamberts et al. 2015) and are particularly promoted between water molecules following the phase transition from amorphous to crystalline (Gálvez et al. 2011).

It may therefore be tempting to compare the ice HDO/H₂O ratios of the separate amorphous and crystalline components to investigate if a difference in deuteration ratio in the unprocessed versus thermally processed ices can be observed. Although the HDO/H₂O ice ratios of the reported crystalline components are higher than those of the amorphous components toward both sources (i.e., the ratios of amorphous:crystalline column densities are lower for HDO than for H₂O), we caution against drawing chemical conclusions from these differences for multiple reasons: 1) the differences are close to the reported error margins and therefore are far from statistically significant; 2) because the two-component fits only approximate the expected thermal gradients, the relative abundances of amorphous to crystalline ices are also only approximate; and 3) grain shape effects were not taken into account when the H₂O bands were fit. This last point is particularly important to remember because, while including grain shape effects in spectral fitting procedures may not greatly affect the derived overall relative molecular abundances of ice species, the warping of the strongest ice bands by grain shape effects do preclude detailed quantitative analysis of ice temperature, morphology, or chemical environment via the profiles of such bands without accounting for grain shape correction (Dartois et al. 2024).

Fig. 8 Comparison of HDO/H2O ratios measured in the gas and ice toward protostars of various masses (adapted from Jensen et al. 2019 and Andreu et al. 2023). Ratios measured in the gas phase are plotted with red diamonds for clustered Class 0 LYSOs, yellow squares for isolated Class 0 LYSOs, blue stars for isolated Class 1 LYSOs, and green pentagons for MYSOs. Ice upper limits derived from pre-JWST spectra are plotted with teal triangles. The JWST ice values from this work are plotted with orange hexagons. The literature values used in this figure are provided in Table D.1.

4.2 Comparing ice and gas HDO/H₂O ratios

Figure 8 provides a comparison of the ice HDO/H₂O ratios measured in this work with previously published protostellar gas phase ratios and ice upper limits from predecessor IR telescopes.

The MYSO gas HDO/H₂O ratios reported in the literature over the past couple decades span an order of magnitude from values as low as 2–3 × 10⁻⁴ (Emprechtinger et al. 2013; Coutens et al. 2014b) to as high as 2–3 × 10⁻³ (van der Tak et al. 2006; Neill et al. 2013). The ice HDO/H₂O ratio measured in this work toward the MYSO IRAS 20126 is in agreement with the highest gas-phase MYSO HDO/H₂O ratios. The agreement supports the frequently made assumption that deuterium abundances detected in the hot inner regions of protostars are representative of deuterium abundances in ices.

An ice HDO/H₂O ratio on the order of 10⁻³ toward a MYSO is also consistent with the conclusions drawn by van der Tak et al. (2006), who observed an inverse correlation between the gas HDO/H₂O ratios and the mass-weighted average envelope temperatures of the four MYSOs they studied. They interpreted this as evidence that the HDO/H₂O ratios on the order of 10⁻³ (measured toward the coldest envelopes) were the most representative of ice HDO/H₂O ratios, while the HDO/H₂O ratios on the order of 10⁻⁴ (measured toward the warmer, more evolved envelopes) had already been partially altered by gas-phase reactions in the hot cores where water ice desorbs. It is interesting to also note that, according to the fits of Boogert et al. (2000), the ¹³CO₂ ice features of the two sources studied by van der Tak et al. (2006) with lower HDO/H₂O ratios show a significant degree of CO2 ice segregation and therefore thermal processing, while the¹³ CO₂ ice feature of the van der Tak et al. (2006) source in which the gas HDO/H₂O ratio was the highest, W33A, is indicative of a much lower degree of thermal processing.

However, it is also possible that MYSOs could vary significantly in their prestellar timescales and conditions, and therefore their ice HDO/H₂O ratios may greatly differ. A larger sample size of MYSO ice HDO/H₂O measurements will be needed to ascertain what causes the observed spread in MYSO gas HDO/H₂O values. In any case, the agreement between the IRAS 20126 ice HDO/H₂O ratio and the gas HDO/H₂O ratios of several other MYSOs suggests that the water gas detected toward at least some hot cores has not been substantially processed by gas-phase reactions.

van der Tak et al. (2006) additionally noted that the gas HDO/H₂O ratios measured toward their MYSOs are one to two orders of magnitude lower than those measured in LYSOs, consistent with the theory that LYSOs experience longer prestellar dense stages than MYSOs. An analogous trend is observed in the deuterium ratios of protostellar gas-phase CH₃OH, which are on average more than an order of magnitude lower in MYSOs than LYSOs; such a difference can be explained with models if MYSOs have shorter and/or warmer prestellar stages relative to LYSOs (van Gelder et al. 2022). However, the most recent gas HDO/H₂O LYSO values that utilize interferometric techniques to target only the hot inner regions of these objects and detect multiple high S/N HDO lines indicate that some of the previously measured LYSO gas HDO/H₂O ratios were likely overestimated by one to two orders of magnitude (Persson et al. 2013). In fact, recent gas phase ratios of isolated LYSOs (which are thought to experience long and cold dense stages) ranging from 1.8–6.3 × 10⁻³ (Jensen et al. 2019, 2021) are very similar to the highest measured MYSO gas ratios as well as the IMYSO and MYSO ice ratios measured here. The clustered LYSO gas ratios are on average a factor of a few lower, ranging from 5.4– 9.2 × 10⁻⁴ (Persson et al. 2014; Jensen et al. 2019), and are more similar to the low-end of measured MYSO ratios.

This similarity in the recently published gas HDO/H₂O ratios of LYSOs and MYSOs could be interpreted in a number of ways. One potential explanation is that the most recently published gas-phase LYSO HDO/H₂O ratios have been altered by gas-phase reactions, meaning that the measured hot inner gas molecular abundances are not fully representative of the molecular abundances of icy species (similar to the interpretation of MYSO gas HDO/H₂O ratios on the order of 10⁻⁴ made by van der Tak et al. 2006). Measurements of ice HDO/H₂O ratios toward LYSOs may prove useful in resolving if this is the case.

Alternatively, some MYSOs may experience prestellar stages with timescales and temperatures similar to LYSOs. Measuring the gas-phase D₂O/HDO ratios in MYSOs could elucidate the credibility of this latter scenario. This is because D₂O/HDO ratios, especially when compared to HDO/H₂O ratios, are also an excellent diagnostic of formation conditions (Furuya et al. 2017). Multiple studies have measured gas D₂O/HDO ratios in LYSOs that are more than a factor of two higher than their gas HDO/H₂O ratios (Coutens et al. 2014a; Jensen et al. 2021), a result that can only be reproduced by models if the majority of deuterated water ice is produced in the cold dense prestellar stage after the heavy CO freeze-out stage alongside CH₃OH ice (Furuya et al. 2016).

Due to the extremely low abundance of D₂O, few gas D2O/HDO ratios have been measured in either LYSOs or MYSOs. Currently, one D₂O/HDO ratio has been measured toward the hot core of the Orion KL MYSO (Neill et al. 2013). This value, 1.6 × 10⁻³, is an order of magnitude lower than the gas D₂O/HDO ratios that have been measured in LYSOs (Jensen et al. 2021), which may indicate that less of the deuter-ated water detected toward Orion KL was produced in the cold dense prestellar stages. However, this is again a comparison of very few data points, and more D₂O/HDO measurements in both LYSOs and MYSOs are needed to draw such a conclusion with certainty.

4.3 Ice envelope evolution

In addition to the ice HDO/H₂O ratios, the ice fitting performed in this work can also provide information about the thermal and chemical evolution of the ice envelopes surrounding the investigated objects.

4.3.1 Ice thermal cycling

The peak position of the crystalline components observed toward HOPS 370 and IRAS 20126 match best with that of annealed HDO ice (i.e., HDO that has been heated to its crystallization temperature and then cooled back down), as fits using nonannealed crystalline HDO ice at 150 K result in Δχ² values that lie outside the 3σ confidence intervals (Fig. B.1). Such an observed redshift from the laboratory crystalline HDO ice peak positions at 150 K may be indicative of thermal cycling within the protostellar envelopes, where ice that was heated to high enough temperatures to transition from an amorphous to a crystalline state was either transported back to cooler regions in the cold envelope (e.g., via motions driven by outflows), or the thermal gradient within the cold envelope changed (e.g., due to an accretion burst).

4.3.2 CH₃OH ice formation

The CH₃OH ice abundances of 16% and 27% with respect to H₂O ice toward both sources are on the high end of the values typically measured for both low-mass and massive protostars (Boogert et al. 2015 and references therein). Although unusual, it is not unprecedented for such high values to be measured toward MYSOs (Dartois et al. 1999). These high CH₃OH ice abundances may be due to particularly long cold dense prestel-lar stages relative to the length of the warmer and less dense stages when the bulk of the H₂O ice is produced. As models suggest that the majority of HDO ice is produced during the cold dense prestellar stage with CH₃OH ice (Furuya et al. 2016), this would imply that the ice HDO/H₂O ratios measured here could also be on the high end of most protostars, and other proto-stars with lower CH₃OH ice abundances relative to H₂O ice that had shorter cold dense prestellar stages could then have lower ice HDO/H₂O ratios. However, it is also possible that thermal processing could affect CH₃OH ice abundances. Correlations between ice HDO/H₂O ratios, CH₃OH ice abundances, and the profiles of ice thermal processing tracers like ¹³CO₂ may be interesting to investigate once the sample size of measured ice HDO/H₂O ratios is larger.

The fact that the current fitting procedure of the 3.5 µm absorption complex requires a component of CH₃OH in H₂O-rich mixtures also carries astrochemical significance. Although the earliest fits to CH₃OH ice features toward protostars often utilized H₂O-rich CH₃OH ice mixtures (Allamandola et al. 1992; Dartois et al. 1999; Brooke et al. 1999; Pontoppidan et al. 2003), subsequent studies popularized an ice evolution model where most of the observed CH₃OH ice forms after CO freeze-out in a relatively CO-rich and H₂O-poor ice environment. This model is supported by analysis of the profiles of the 9.75 µm CH₃OH ice C–O stretching mode (Bottinelli et al. 2010), the “red component” of the 4.67 µm CO ice C≡O stretching mode (Cuppen et al. 2011; Penteado et al. 2015), the relationships between CH₃OH ice column densities, CO ice column densities, and visual extinction A_v (Boogert et al. 2015), and the experimentally demonstrated efficiency of solid-state CH₃OH formation via CO hydrogenation (Watanabe & Kouchi 2002; Fuchs et al. 2009).

However, recent experimental works indicate that CH₃OH can also be produced in H₂O-rich ices via a reaction between CH₄ and OH radicals (Qasim et al. 2018), albeit less efficiently than via CO hydrogenation. Recent JWST data toward prestel-lar cores, where the 9.75 µm CH₃OH ice feature is well fit with a combination of CO-rich and H₂O-rich CH₃OH ice mixtures, provide observational evidence supporting this pathway as a secondary means of prestellar CH₃OH ice production (McClure et al. 2023). In our fits, the lack of a CO-rich CH₃OH ice component is likely due to the fact that most of the CO ice toward these lines of sight has thermally desorbed, as indicated by the low CO ice column density and presence of CO gas lines toward both lines of sight (Federman et al. 2024; Rubinstein et al. 2023), leaving behind a layer of nearly pure CH₃OH ice that likely formed via CO hydrogenation, as well as some CH₃OH ice in a water-rich chemical environment that may have formed during the less dense prestellar stages, perhaps via the CH₄ + OH pathway.

Notably, despite the presence of thermal processing tracers in both spectra, there is a lack of strong spectroscopic evidence of heated or crystalline CH₃OH ice along either line of sight. Perhaps this could be explained by the heated ices having mostly experienced temperatures too hot for CH₃OH to remain in the ice (i.e., a steep thermal gradient along the line of sight). However, both the dominant component of the CH₃OH feature being pure CH₃OH ice as well as the low optical depth of the CO ice feature indicate that the majority of the “cold” ices in the outermost envelope must have also experienced some amount of heating (equivalent to at least 30–40 K in the laboratory).

5 Conclusions

Here we present >3σ detections of both amorphous and crystalline HDO ice at 4.1 µm in JWST NIRSpec spectra taken toward the IMYSO HOPS 370 and the MYSO IRAS 20126+4104. The HDO ice column densities were quantified in a fitting procedure that accounted for potential spectral contributions from CH₃OH ice. The H₂O ice column densities were quantified via the 3 µm feature to allow comparisons between HDO/H₂O ratios measured in the ice and those measured in the gas toward protostars. Our main findings are summarized as follows:

1. Toward both sources, the detection of a crystalline HDO ice component at 4.1 µm correlates with the detection of a crystalline H₂O ice component at 3 µm as well as other tracers of ice thermal processing, increasing confidence in our assignment of the 4.1 µm feature to HDO ice;

2. We measure ice HDO/H₂O ratios of 4.6 ± 2.2 × 10⁻³ and 2.6 ± 1.4 × 10⁻³ toward HOPS 370 and IRAS 20126, respectively. These ratios are enriched by approximately two orders of magnitude relative to the cosmic standard deuterium abundance, consistent with the water ice toward these objects having cold prestellar origins;

3. The ice HDO/H₂O ratio measured toward the MYSO IRAS 20126 is consistent with the highest reported gas HDO/H₂O ratios toward MYSOs (a couple × 10⁻³), supporting the assumption that the deuterium abundances of the gas-phase waterobserved toward these objects are representative ofthe deuterium abundances of their surrounding water ice;

4. The ice HDO/H₂O ratios toward these objects are remarkably similar to each other and to the gas HDO/H₂O ratios measured toward isolated LYSOs. A larger sample of ice HDO/H₂O ratios toward LYSOs and MYSOs, as well as more measurements of D₂O/H₂O ratios in the gas phase, could help to explain the cause of this similarity;

5. The peak position of the crystalline HDO ice component toward both sources suggests that the thermally processed ice in these lines of sight experienced recooling to some degree;

6. The CH₃OH ice requires a water-rich CH₃OH ice mixture to obtain a satisfactory fit to the 3.4 µm region, supporting recent experimental and observational evidence that some CH₃OH ice may form prior to the heavy CO freeze-out stage in H₂O-rich environments;

7. As part of this work, several new laboratory spectra of HDO, H₂O, CH₃OH, and CO ice mixtures are now available for public use on the Leiden Ice Database (LIDA).

Acknowledgements

K.S. thanks Brian Ferrari, Adwin Boogert, Jenny Noble, Helen Fraser, Melissa McClure, Sergio Ioppolo, Neal Evans, and Alessio Caratti o Garatti for helpful discussions, Adwin Boogert for providing IRTF spectra, and Julia Santos for experimental support. This work is based on observations made with the NASA/ESA/CSA James Webb Space Telescope. The data were obtained from the Mikulski Archive for Space Telescopes at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-03127 for JWST. These observations are associated with program #1802. All the JWST data used in this paper can be found in MAST: 10.17909/3kky-t040. Astrochemistry at Leiden is supported by funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No. 101019751 MOLDISK), the Netherlands Research School for Astronomy (NOVA), and the Danish National Research Foundation through the Center of Excellence “InterCat” (Grant agreement no.: DNRF150). Support for SF, AER, STM, RG, JJT and DW in program #1802 was provided by NASA through a grant from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-03127. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. Y.-L.Y. acknowledges support from Grant-in-Aid from the Ministry of Education, Culture, Sports, Science, and Technology of Japan (20H05845, 20H05844, 22K20389), and a pioneering project in RIKEN (Evolution of Matter in the Universe).

Appendix A Laboratory data

Here we present the methodologies used to collect and analyze the laboratory data in this work.

Appendix A.1 Experimental details

In this work, independent dosing lines were used to deposit the CH₃OH-containing binary ice mixtures using the procedure outlined in Yarnall & Hudson (2022), where each leak valve is calibrated so that its position results in the desired deposition rate of each mixture component in its pure form. After this calibration step, shut valves between the gas or vapor reservoirs and the leak valves are opened so that the mixture components deposit simultaneously upon the substrate. As mentioned in Yarnall & Hudson (2022), the leak valve calibration depends on knowing the refractive index and density of each ice component, as well as the measurement of the laser interference pattern as the ice deposits in the calibration step. The benefits of utilizing such a system of independent dosing lines to create mixtures in situ on the ice substrate rather than using a single dosing line and mixing gases and vapors prior to their introduction into the chamber have been previously discussed in Gerakines et al. (1995), Yarnall & Hudson (2022), and Slavicinska et al. (2023).

A HeNe laser with a wavelength of 632.8 nm was used for the laser interference measurements. The refractive indexes and densities of the ices used for our calibrations are provided in Table A.1.

Because the physical properties and masses of H₂O and HDO are very similar, their flow rates through the leak valves, pumping speeds, and deposition rates are also expected to be very similar, minimizing the expected error of the HDO:H₂O ice mixing ratio when premixing the components prior to deposition. Therefore, the HDO:H₂O ices were created via the traditional premixing method, following the procedure described in Gálvez et al. (2011) to generate HDO via mixing H₂O and D₂O.

During the warm-ups, the ices were heated at a rate of 25 K hr⁻¹, and spectra ranging from 4000–500 cm⁻¹ (2.5–20 µm) were continuously collected and averaged every 128 scans, meaning a new spectrum was obtained every 3.5 min. This results in an uncertainty of approximately ± 1 K in the reported temperature of each amorphous spectrum. During the recooling of the annealed crystalline HDO:H₂O ice mixtures, a cooling rate of −2 K min⁻¹ was used so that the experiment could be performed within a single workday, resulting in a higher temperature uncertainty of approximately ± 4 K for the annealed crystalline HDO spectra. All spectra were baseline-corrected via cubic spline functions prior to their analysis.

Appendix A.2 HDO O-D stretch characterization

The peak position, FWHM, and integrated absorbance of the HDO O-D stretching mode was extracted from the thick HDO:H₂O ice spectra (see Table 2) using a cubic spline to isolate the feature from the broad 4.5 µm H₂O combination mode followed by a Savitzky-Golay filter to smooth out experimental noise (see Figure A.1).

The apparent band strengths of the HDO O-D stretch at different temperatures, $A_i^{\prime }$, were derived by assuming the apparent band strength of the crystalline band at 150 K, $A_{150K}^{\prime }$, is equal to 6.4 × 10⁻¹⁷ cm molec⁻¹ (Dartois et al. 2003; Gálvez et al. 2011) and using the approximation:

$$A_i^{\prime }=\frac{{\displaystyle \int {\mathrm{abs}}_{\text{i}}}(v)\text{d}v}{{\displaystyle \int {\mathrm{abs}}_{150K}}(v)\text{d}v}\times A_{150K}^{\prime },$$(A.1)

where ∫ abs_150K(ν)dν is the integrated absorbance of the crystalline HDO peak at 150 K, and ∫ abs_i(v)dv is the integrated absorbance of the HDO peak at the phase and temperature of interest.

The resulting peak positions, FWHMs, and apparent band strengths are plotted in Figures A.2 and A.3 and reported in Table A.2. It is evident that the peak extractions are much less reliable for the amorphous ice, where the O-D stretch is very weak and broad (see Figure A.1), resulting in a greater experimental scatter between the amorphous data points in both plots. We smoothed this scatter with fourth-order polynomial fits to the individual amorphous data points, which were used to calculate the final peak positions, FWHMs, and apparent band strengths reported for the amorphous ice. The crystalline data was not treated in this way because the profile extraction of the crystalline feature from the laboratory spectra was much more straightforward, resulting in much smaller experimental scatter.

Despite the problematic profile of the amorphous HDO ice band, the apparent band strength we derived for amorphous HDO ice at 15 K, 4.2 × 10⁻¹⁷ cm molec⁻¹, is in excellent agreement (within ~2%) with the band strengths reported by Dartois et al. (2003) and Gálvez et al. (2011) for amorphous HDO ice at similar temperatures.

Appendix A.3 CH₃ OH C-H symmetric stretch characterization

The peak position of the CH₃OH C-H symmetric stretching mode was extracted from the CH₃OH-containing ice spectra using a fourth-order polynomial to isolate the feature from the entire 3.4 µm absorption complex followed by a Savitzky-Golay filter to smooth out experimental noise. These values are given in Tables A.3 and A.4 and are plotted as a function of temperature in Figure A.4, which also shows the observed peak positions of the C-H symmetric stretching mode in the analyzed spectra. Only peak positions are provided here because a full characterization of this peak (i.e., including FWHMs and relative apparent band strengths) is complicated by the fact that this feature is blended with the CH₃OH overtones and C-H asymmetric stretching modes.

Fig. A.1 Example extractions of the profile of the 4.1 µm HDO O-D stretching mode from the H2O combination mode in 0.4% HDO:H2O laboratory spectra (top: amorphous; bottom: crystalline). The laboratory data (already baseline corrected) is plotted in black, the local continuum is plotted in dashed red, and the smoothed and extracted profile is plotted in green. The red shading indicates the wavelengths used to define the local continuum, and the green shading indicates the integrated area used to derive the apparent band strengths.

Fig. A.2 Temperature-dependent variation of the apparent band strengths of amorphous and crystalline HDO ice. The dashed blue line is a fourth-order polynomial fit to the amorphous band strength data used to extract the final reported apparent band strengths reported in Table A.2.

Fig. A.3 Range of peak positions and FWHMs of amorphous and crystalline HDO ice in the heating and annealing experiments. The arrows indicate the direction in which the data changes with increasing temperature. The dashed blue line is a fourth-order polynomial fit to the amorphous peak position and FWHM data used to extract the final reported peak positions and FWHMs reported in Table A.2.

Fig. A.4 Peak positions oſthe3.53 µm CH3OH C-H stretching mode in laboratory CH3OH ices from 15–145 K (points), and peak positions of the 3.53 µm band observed toward HOPS 370 and IRAS 20126 (thick horizontal lines). The observed peak positions were extracted from the observed spectra via Gaussian fits to the local-continuum subtracted 3.53 µm peaks (error margins indicated as shaded areas with dashed borders). The CH3OH:CO2 values were extracted from data from Ehren-freund et al. (1999) available on LIDA.

Appendix B: χ² plots

Least-squares fits using every possible combination of one amorphous and one crystalline HDO lab spectra were performed, and the χ² values were calculated for all of these possible fits using the following equation:

$${\chi }^2=\frac{{\displaystyle \sum {\left({\tau }_{obs}-{\tau }_{fit}\right)}^2}}{{\sigma }^2},$$(B.1)

where τ_obs is the observed optical depth and τ_fit is the optical depth of the fit. The σ value is calculated by propagating the uncertainties in the optical depths that result from the RMS errors and the local continuum placement uncertainties (Figures C.l and C.2). The 1, 2, and 3cr bounds plotted in Figure B.l correspond to the critical ∆χ² (χ²-χ²_min) values at the 68.3, 95.5, and 99.7% confidence intervals, respectively, for fits with four free parameters (in this case, the temperatures and scaling factors of the amorphous and crystalline HDO ice spectra).

Fig. B.1 Contour maps of the ∆χ2 values of the least-squares fits with every possible combination of one amorphous and one crystalline HDO component. The white X marks χ2min, and the white, purple, and yellow contour lines show the ∆χ2 values corresponding to 1, 2, and 3σ, respectively.

Appendix C Estimated continuum errors

Appendix C.1 CH₃OH ice

The 20% band strength uncertainty (Luna et al. 2018) is propagated to determine the uncertainty in the calculated CH₃OH ice column densities. The implications of these fits and column densities on the thermal and chemical structures of the ice envelopes as well as the physicochemical histories of the sources are discussed in Section 4.3.2.

Appendix C.2 HDO ice

The reported errors account for two sources of uncertainty: 1) the uncertainty of the HDO O-D stretching mode band strength, and 2) the uncertainty of the local continuum choice. The uncertainty in the band strength accounts for both the aforementioned variation of the band strength with temperature as well as the reported experimental uncertainties of the amorphous and crystalline HDO band strengths from Gálvez et al. (2011) (~10% and 5%, respectively), which were used to calculate the band strengths of HDO at higher temperatures presented in Table A.2.

The uncertainty of the local continuum is somewhat subjective, as it depends on what one considers a “reasonable” local continuum. We attempted to quantify this uncertainty by identifying the limits of what we consider a reasonable local continuum for each line of sight and then quantifying the variation in the calculated HDO column density that resulted from these continuum shifts. These lower and upper limit continua are shown in Figure C.1, and the resulting lower and upper limits on the optical depths of the 4 µm region are shown in Figure C.2.

The column density uncertainties were determined by adding in quadrature the band strength uncertainties and the uncertainties caused by the continuum choice, which were estimated by noting the variation in HDO ice column density obtained from repeating the fitting procedure on the spectra that were extracted using the lower and upper limit global continua.

Appendix C.3 H₂O ice

The uncertainty in the H₂O column density due to the uncertainty of the continuum choice was quantified via a procedure similar to that used for the 4 µm region by exploring the range of artificial short-wavelength continuum points which resulted in continuum shapes and H₂O profiles that appeared reasonable. The lower and upper limit continua are shown in Figure C.3, and the resulting variations in the optical depths and profiles of the 3 µm feature are shown in Figure C.4.

The column density uncertainties were determined by adding in quadrature the uncertainties caused by the variation in band strengths with temperature, a 30% uncertainty accounting for the underfit red wing caused by neglecting ammonia hydrates and grain shape effects (as estimated by Boogert et al. 2008), and the uncertainties caused by the continuum choice, which were estimated by noting the variation in H₂O ice column density obtained from repeating the fitting procedure on the spectra that were extracted using the lower and upper limit global continua.

Fig. C.1 Alternative local continua used to calculate the uncertainty in HDO column density from local continuum choice. Left: Lower limit continua (i.e., continua that result in the lowest possible optical depth in the region while remaining reasonable). Right: Upper limit continua (i.e., continua that result in the highest possible optical depth in the region while remaining reasonable).

Fig. C.2 Change in the optical depth of the features at 4 µm caused by varying the local continua, as shown in Figure C.1. The spectra used to calculate the column densities in Table 4 are plotted in darker shades, and the spectra resulting from the lower and upper continuum limits used to estimate the effect of continuum uncertainty are plotted in lighter shades.

Fig. C.3 Alternative global continua used to calculate the uncertainty in H2O column density from local continuum choice. Left: Lower limit continua (i.e., continua that result in the lowest possible optical depth in the region while remaining reasonable). Right: Upper limit continua (i.e., continua that result in the highest possible optical depth in the region while remaining reasonable).

Fig. C.4 Change in the optical depth of the 3 µm feature caused by varying the global continua, as shown in Figure C.3. The spectra used to calculate the column densities in Table 5 are plotted in darker shades, and the spectra resulting from the lower and upper continuum limits used to estimate the effect of continuum uncertainty are plotted in lighter shades.

Appendix D Literature HDO/H₂O values

References

All Tables

All Figures

	Fig. 1 5 µm continuum maps with the aperture used to extract the analyzed spectra at 5 µm shown in white (7 × diffraction-limited radius at 5 µm = 0.678″).
In the text

Fig. 2 Extracted spectra of the investigated sources. Left: full NIRSpec spectra with major ice features labeled. The original spectra before gas emission subtraction are plotted in gray. The spectral region of the strongest HDO absorption, the O–D stretching mode, is indicated in magenta. The spectral region plotted on the right is indicated in gray. Right: zoomed-in view of the spectral region of interest for HDO ice.

In the text

	Fig. 3 Extraction and fitting of the CH3OH 3.53 µm band. Left: adopted local continua. Right: fits to the extracted optical depth spectra. Some of the substructures present at 3.4 µm in the high S/N spectrum of IRAS 20126 are well reproduced by the more prominent CH3OH overtones and C-H asymmetric stretching modes that emerge in H2O-rich CH3OH spectra.
In the text

Fig. 4 Observed CH3OH C–H symmetric stretching modes compared to select laboratory data. The peaks were extracted from the 3.5 µm absorption complexes via a local continuum fit using baseline points in the range ~3.47–3.50 and 3.60–3.62 µm (Fig. 3). It is clear from these comparisons that cold, amorphous, and relatively pure CH3OH is the major contributor to this feature. Left: observed peaks vs. laboratory spectra of CH3OH ice in various chemical environments. All spectra were collected at 15 K except the CH3OH:CO2 spectrum, which was collected at 10 K (Ehrenfreund et al. 1999). Right: observed peaks vs. laboratory spectra of pure CH3OH ice at various temperatures.

In the text

Fig. 5 Extraction and fitting of the 4 µm region. Left: adopted local continua, along with CH3OH spectra determined from the 3.53 µm band fitting used to aid the continuum placement. The spectral region of the HDO O–D stretching mode is indicated in magenta. Right: fits to the extracted optical depth spectra (am. = amorphous; cr. = crystalline). The shaded gray regions indicate 3σ uncertainty levels. The same combination of pure and H2O-rich CH3OH spectra as used to fit the 3.53 µm feature is used in this fit. The best-fitting temperatures of the amorphous and crystalline HDO ice components are 132 K and 47 K for HOPS 370 and 138 K and 90 K for IRAS 20126, respectively (see Sect. 3.3 and Appendix B for more details).

In the text

	Fig. 6 Simplified diagram showing the difference between the expected thermal profile of the observed ice columns and the two-component fit approximation used to fit the observations in this work.
In the text

	Fig. 7 Extraction and fitting of the H2O 3 µm band. Left: adopted global continua and the archival MYSO ice spectra used as templates to guide the continuum choice. Right: fits to the extracted optical depth spectra (am. = amorphous; cr. = crystalline).
In the text

Fig. 8 Comparison of HDO/H2O ratios measured in the gas and ice toward protostars of various masses (adapted from Jensen et al. 2019 and Andreu et al. 2023). Ratios measured in the gas phase are plotted with red diamonds for clustered Class 0 LYSOs, yellow squares for isolated Class 0 LYSOs, blue stars for isolated Class 1 LYSOs, and green pentagons for MYSOs. Ice upper limits derived from pre-JWST spectra are plotted with teal triangles. The JWST ice values from this work are plotted with orange hexagons. The literature values used in this figure are provided in Table D.1.

In the text

Fig. A.1 Example extractions of the profile of the 4.1 µm HDO O-D stretching mode from the H2O combination mode in 0.4% HDO:H2O laboratory spectra (top: amorphous; bottom: crystalline). The laboratory data (already baseline corrected) is plotted in black, the local continuum is plotted in dashed red, and the smoothed and extracted profile is plotted in green. The red shading indicates the wavelengths used to define the local continuum, and the green shading indicates the integrated area used to derive the apparent band strengths.

In the text

	Fig. A.2 Temperature-dependent variation of the apparent band strengths of amorphous and crystalline HDO ice. The dashed blue line is a fourth-order polynomial fit to the amorphous band strength data used to extract the final reported apparent band strengths reported in Table A.2.
In the text

	Fig. A.3 Range of peak positions and FWHMs of amorphous and crystalline HDO ice in the heating and annealing experiments. The arrows indicate the direction in which the data changes with increasing temperature. The dashed blue line is a fourth-order polynomial fit to the amorphous peak position and FWHM data used to extract the final reported peak positions and FWHMs reported in Table A.2.
In the text

Fig. A.4 Peak positions oſthe3.53 µm CH3OH C-H stretching mode in laboratory CH3OH ices from 15–145 K (points), and peak positions of the 3.53 µm band observed toward HOPS 370 and IRAS 20126 (thick horizontal lines). The observed peak positions were extracted from the observed spectra via Gaussian fits to the local-continuum subtracted 3.53 µm peaks (error margins indicated as shaded areas with dashed borders). The CH3OH:CO2 values were extracted from data from Ehren-freund et al. (1999) available on LIDA.

In the text

	Fig. B.1 Contour maps of the ∆χ2 values of the least-squares fits with every possible combination of one amorphous and one crystalline HDO component. The white X marks χ2min, and the white, purple, and yellow contour lines show the ∆χ2 values corresponding to 1, 2, and 3σ, respectively.
In the text

	Fig. C.1 Alternative local continua used to calculate the uncertainty in HDO column density from local continuum choice. Left: Lower limit continua (i.e., continua that result in the lowest possible optical depth in the region while remaining reasonable). Right: Upper limit continua (i.e., continua that result in the highest possible optical depth in the region while remaining reasonable).
In the text

	Fig. C.2 Change in the optical depth of the features at 4 µm caused by varying the local continua, as shown in Figure C.1. The spectra used to calculate the column densities in Table 4 are plotted in darker shades, and the spectra resulting from the lower and upper continuum limits used to estimate the effect of continuum uncertainty are plotted in lighter shades.
In the text

	Fig. C.3 Alternative global continua used to calculate the uncertainty in H2O column density from local continuum choice. Left: Lower limit continua (i.e., continua that result in the lowest possible optical depth in the region while remaining reasonable). Right: Upper limit continua (i.e., continua that result in the highest possible optical depth in the region while remaining reasonable).
In the text

	Fig. C.4 Change in the optical depth of the 3 µm feature caused by varying the global continua, as shown in Figure C.3. The spectra used to calculate the column densities in Table 5 are plotted in darker shades, and the spectra resulting from the lower and upper continuum limits used to estimate the effect of continuum uncertainty are plotted in lighter shades.
In the text