JWST detections of amorphous and crystalline HDO ice toward massive protostars
Katerina Slavicinska Ewine F. van Dishoeck Łukasz Tychoniec Pooneh Nazari Adam E. Rubinstein Robert Gutermuth Himanshu Tyagi Yuan Chen Nashanty G. C. Brunken Will R. M. Rocha P. Manoj Mayank Narang S. Thomas Megeath Yao-Lun Yang Leslie W. Looney John J. Tobin Henrik Beuther Tyler L. Bourke Harold Linnartz Samuel Federman Dan M. Watson Hendrik Linz
 




   



   

    affiliations can be found after the references

   







  
   Tracing the origin and evolution of interstellar water is key to understanding many of the physical and chemical processes involved in star and planet formation. Deuterium fractionation offers a window into the physicochemical history of water due to its sensitivity to local conditions.
  
   This work aims to utilize the increased sensitivity and resolution of the James Webb Space Telescope (JWST) to quantify the HDO/H_2O ratio in ices toward young stellar objects (YSOs) and to determine if the HDO/H_2O ratios measured in the gas phase toward massive YSOs (MYSOs) are representative of the ratios in their ice envelopes.
  
   Two protostars observed in the Investigating Protostellar Accretion (IPA) program using JWST NIRSpec were analyzed: HOPS 370, an intermediate-mass YSO (IMYSO), and IRAS 20126+4104, a massive YSO (MYSO). The HDO ice toward these sources was quantified via its 4.1 μm band. The contributions from the CH_3OH combination modes to the observed optical depth in this spectral region were constrained via the CH_3OH 3.53 μm band to ensure that the integrated optical depth of the HDO feature was not overestimated. H_2O ice was quantified via its 3 μm band. New laboratory IR spectra of ice mixtures containing HDO, H_2O, CH_3OH, and CO were collected to aid in the fitting and chemical interpretation of the observed spectra.
  
   HDO ice is detected above the 3σ level in both sources. It requires a minimum of two components, one amorphous and one crystalline, to obtain satisfactory fits. The H_2O ice band at 3 μm similarly requires both amorphous and crystalline components. The observed peak positions of the crystalline HDO component are consistent with those of annealed laboratory ices, which could be evidence for heating and subsequent re-cooling of the ice envelope (i.e., thermal cycling). The CH_3OH 3.53 μm band is fit best with two cold components, one consisting of pure CH_3OH and the other of CH_3OH in an H_2O-rich mixture. From these fits, ice HDO/H_2O abundance ratios of 4.6±1.8×10^-3 and 2.6±1.2×10^-3 are obtained for HOPS 370 and IRAS 20126+4104, respectively.
  
   The simultaneous detections of both crystalline HDO and crystalline H_2O corroborate the assignment of the observed feature at 4.1 μm to HDO ice. The ice HDO/H_2O ratios are similar to the highest reported gas HDO/H_2O ratios measured toward MYSOs as well as the hot inner regions of isolated low-mass protostars, suggesting that at least some of the gas HDO/H_2O ratios measured toward massive hot cores are representative of the HDO/H_2O ratios in ices. The need for an H_2O-rich CH_3OH component in the CH_3OH ice analysis supports recent experimental and observational results that indicate that some CH_3OH ice may form prior to the CO freeze-out stage in H_2O-rich ice layers.

   
             
     
   JWST detections of amorphous and crystalline HDO ice toward massive protostars
    Katerina Slavicinska1,2
            Ewine F. van Dishoeck2,3
            Łukasz Tychoniec2
            Pooneh Nazari2,4
            Adam E. Rubinstein5
            Robert Gutermuth6
            Himanshu Tyagi7
            Yuan Chen2
            Nashanty G. C. Brunken2
            Will R. M. Rocha2
            P. Manoj7
            Mayank Narang8
            S. Thomas Megeath9
            Yao-Lun Yang10
            Leslie W. Looney11, 12
            John J. Tobin12
            Henrik Beuther13
            Tyler L. Bourke14
            Harold Linnartz1
            Samuel Federman9
            Dan M. Watson5
            Hendrik Linz13,17
            
    Received 28 February 2024 / Accepted 23 April 2024
===================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================






§ INTRODUCTION


Water is one of the most abundantly detected interstellar molecules in both the gas and solid phase, 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., ). 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 <cit.>. 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_2, NH_3, and CH_4, resulting in the growth of ice layers that can be observed via infrared spectroscopy ( and references therein). There is mounting evidence that these water-rich prestellar ice mantles also serve as the birthplace and reservoir of many other more complex molecules, some of which could have prebiotic relevance <cit.>.

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 <cit.>. The extent to which such processing of the prestellar water ice occurs within the hot, radiation-dominated regions of protostellar envelopes and protoplanetary disks has important implications for not only 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, ).

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 prestellar 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_2D^+ <cit.>:


H3+ + HD ⇌ H2D+ + H2 + 232 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 <cit.>:


H2D+ + e- → H2 + 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^-5, corresponding to an HDO/H_2O ratio of 4×10^-5 if there is no fractionation, ). The relative abundance of H_2D^+ is further enhanced in the absence of neutral species like O and CO that destroy H_3^+ 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 <cit.>. The relatively high brightness of the hot (>100 K) inner cores of such massive young stellar objects (MYSOs) allows multiple lines of HDO to be detected relatively easily with millimeter spectroscopy. The detection of these lines, along with the lines of H_2^18O, a water isotopolog that is usually optically thin, has enabled the determination of the gas HDO/H_2O ratios toward many of these objects in the last few decades. These ratios typically range between 10^-4 to a few 10^-3 <cit.>, an enhancement of up to two orders of magnitude relative to the cosmic standard D/H. Highly sensitive millimeter interferometers also enable reliable measurement of gas HDO/H_2O ratios in the hot inner regions of the much dimmer low-mass young stellar objects (LYSOs), whose ratios see a similar enhancement <cit.>.

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_2O 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 <cit.>.

While studying HDO/H_2O 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, a glaring question remains: are the HDO/H_2O ratios measured in protostellar hot gas actually representative of the protostellar HDO/H_2O ice ratios, or has the detected water gas already experienced substantial alteration by gas-phase reactions? To answer this question, the HDO/H_2O ratio in ices must be probed directly.

Despite H_2O being 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 a couple orders of magnitude below that of molecules like CH_3OH and H_2CO thought to form only in the later colder stages of cloud formation <cit.>, 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^-1), lies next to the short-wavelength wing of the strong asymmetric C=O stretching mode of CO_2 ice, which sometimes warps the continuum due to grain shape effects and complicates the detection and quantification of spectrally adjacent species <cit.>. In the amorphous state, the O-D stretch of HDO is also broad and relatively weak <cit.>.

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 <cit.>, but the detection was soon called into question by <cit.>, who pointed out, among several other concerns, that the observed absorption could also be attributed to a weak CH_3OH combination mode whose wavelength overlaps with that of the HDO O-D stretch. They went on to report HDO upper limits ranging from 0.16-1.1% with respect to H_2O toward four intermediate-mass and massive young stellar objects (IMYSOs and MYSOs) with low CH_3OH column densities. <cit.> expanded the search for HDO ice to low-mass young stellar objects (LYSOs) and also derived upper limits, ranging from 0.5-2% with respect to H_2O. Almost a decade later, <cit.> reported tentative detections and measured very high HDO/H_2O 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_3OH ice to the HDO ice spectral region. Furthermore, given the recent characterization of how grain growth affects the continuum around the CO_2 ice feature <cit.>, it is possible that some of the observed "broad absorption" that was attributed to HDO ice by <cit.> 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; ). Care is taken to ensure that the observed HDO ice detections cannot be attributed to the CH_3OH combination mode at 4 μm by establishing a CH_3OH column density via fitting high-resolution laboratory spectra to the 3.53 μm band. Additional high-resolution laboratory spectra of HDO are collected to enable 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, ). Comparisons are made between the derived HDO ice abundances to 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.




§ METHODS




 §.§ 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 <cit.>. Out of these sources, two are presently analyzed due to their high S/N and because their IFUs contain lines of sight that have a non-saturated 3 μm O-H stretching feature of H_2O ice, enabling quantification of its column density and characterization of its morphology. Table <ref> 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 <cit.>. It is driving a powerful bipolar outflow observed in the far IR, millimeter, and centimeter wavelength ranges <cit.>.

IRAS 20126 is a deeply embedded MYSO, the most massive and luminous protostar of the sample. The source contains a hot and massive disk <cit.> as well as jet-driven outflows <cit.>. Although initial observations suggested that its environment is unusually isolated for a MYSO <cit.>, later observations revealed the presence of a young cluster within a 1.2 pc vicinity <cit.>. Two likely companion YSOs are present within the 6"x6" FOV in the NIRSpec IFU of this source <cit.>.



 §.§ Observations and data reduction


The investigated spectra were collected with the NIRSpec IFU using the G395M grating (R=λ / Δλ∼700-1300, ), providing data from 2.9-5.1 μm (3450-1960 cm^-1) 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 <cit.> and <cit.>.

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_2O 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_2O feature is contaminated by an extended weak emission that may originate from a photodissociation 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 <ref> and are shown on the IFU continuum maps in Figure <ref>. 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 <cit.>.




 §.§ Gas emission subtraction

Both spectra contain gas emission features, including those of H_2, H I, [Fe II], and CO <cit.>. 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 over-plotted in color.




 §.§ Laboratory data collection


High-resolution (0.5 cm^-1) spectra of various laboratory ice mixtures containing HDO, H_2O, CO, and CH_3OH 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 set-up can be found in <cit.>, and recent upgrades to the set-up are described in <cit.>. The method used to create ice mixtures with specific ratios is similar to that described in <cit.>. The chemical compositions and mixing ratios of the laboratory ice mixtures are presented in Table <ref>. Additional experimental details can be found in Appendix <ref>.





All of the spectra are made available for public use and enjoyment via LIDA <cit.>. 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_3OH 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_2O) 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 <ref>. 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_2O mixtures for both the amorphous and crystalline phases and are presented in Table <ref> and Figures <ref> and <ref>. As discussed in detail in Section <ref>, 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_3OH 3.53 μm band from 15 to 150 K in all of the CH_3OH-containing spectra are presented in Tables <ref> and <ref> and Figure <ref>.




§ RESULTS




The full extracted spectra of both sources are presented in Figure <ref>. At first glance, it is apparent that in both spectra, the 3 μm O-H stretching modes of H_2O ice have sharp peaks at 3.1 μm, characteristic of a crystalline component <cit.>. 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_2 components in the ^13CO_2 bands <cit.>, enhanced OCN^- abundances <cit.>, and the relatively low absorptions of the CO ice features at 4.67 μm relative to the CO_2 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_3OH ice column densities and physicochemical environments were determined from the ∼3.53 μm (∼2830 cm^-1) C-H stretching mode (Section <ref>). This was done prior to the analysis of the HDO ice because CH_3OH has several weak combination modes between ∼3.7-4.2 μm (∼2700-2400 cm^-1), one of which overlaps with the HDO O-D stretching mode at 4.1 μm; therefore, pre-constraining the expected contribution of CH_3OH ice to this region via a column density derived from a stronger band in a different region is necessary to ensure that any CH_3OH absorption at 4 μm is not misattributed to HDO ice (Section <ref>), a concern previously raised by <cit.>. The HDO ice column densities and morphologies were then determined via fitting the observed absorptions in the ∼3.7-4.2 μm region (Section <ref>). Finally, the 3 μm (∼3300 cm^-1) H_2O O-H stretching mode was fit to corroborate the morphologies found for HDO ice and to determine the ice HDO/H_2O ratios toward these sources (Section <ref>).




 §.§ CH_3OH ice





Both sources present a clear CH_3OH C-H stretching mode at 3.53 μm. The CH_3OH ice column density can be obtained from this mode via a procedure outlined in <cit.> and <cit.>. 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^-1) and a full-width half maximum (FWHM) of 0.11 μm (91.4 cm^-1) along with laboratory CH_3OH 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_3OH fits, along with the local continuum subtractions used to extract the CH_3OH features into the optical depth scale, are shown in Figure <ref>. 

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_3OH ice (Figures <ref> and <ref>). The 3.53 μm bands of laboratory CH_3OH ices that is warm, crystalline, or diluted with any of the three most abundant observable interstellar ice components (H_2O, CO, and CO_2) are too blue-shifted to fit the observed absorption. However, fitting the entire 3.5 μm region with only pure amorphous CH_3OH ice results in the spectrum being overfit at 3.4 μm by the blended CH_3OH overtones and C-H asymmetric stretching modes.




This issue of pure CH_3OH 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_3OH mixed with H_2O <cit.>, 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_2-rich CH_3OH mixtures or in crystalline CH_3OH ices. Using a combination of two laboratory amorphous CH_3OH spectra, one of pure CH_3OH ice and an H_2O-rich CH_3OH ice mixture (1:5 CH_3OH:H_2O), in the fitting procedure results in satisfactory fits of both the 3.53 μm band and the 3.4 μm region (Figure <ref>). A similar combination of pure CH_3OH and CH_3OH:H_2O laboratory spectra was used by <cit.> to fit the 3.53 and 3.9 μm features in the spectra observed toward MYSOs W33A and RAFGL 7009s.

The resulting CH_3OH column densities (reported in Table <ref>) are obtained using the band strength of 4.86×10^-18 cm molec^-1 for the C-H stretching mode of pure amorphous CH_3OH at 20 K <cit.>. This band strength was chosen because the profile of the 3.53 μm band indicates that the majority of the detected CH_3OH ice along these lines of sight is pure, amorphous, and cold (>60 K). The ratio of pure CH_3OH ice to CH_3OH ice mixed with H_2O that results in the best fits is highly dependent on the dilution factor of the CH_3OH:H_2O mixture used, so separate column densities for each of the two CH_3OH components were not calculated. Any effects of dilution with H_2O on the band strength of this feature are expected to be well within the 20% uncertainty reported for the pure band strength <cit.>. The reported uncertainties are explained in Appendix <ref>. The implications of the non-detection of warm or crystalline CH_3OH ice toward these lines of sight are further discussed in Section <ref>.






 §.§ 4 μm region: CH_3OH, S-bearing species, and HDO



A local continuum in the 3.7-4.2 μm (2700-2380 cm^-1) region was adopted so that the optical depth of the CH_3OH combination modes predicted by the CH_3OH 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 over-subtraction 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^-1). It is important to note that, as can be seen in the right panel of Figure <ref>, the warping of the continuum by the scattering wings of the strong 4.27 μm CO_2 feature is very minimal in both spectra in comparison to some other icy lines of sight (e.g., ), 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 Figure <ref>. 

Following this procedure, both lines of sight have excess absorption on top of the CH_3OH feature at ∼3.9 μm (∼2550 cm^-1). <cit.> noted a similar excess absorption in the spectrum of the MYSO W33A and used it to derive an H_2S ice upper limit. However, <cit.> recently pointed out that many other S-H bond-containing ices, like CH_3SH and CH_3CH_2SH, 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^-1), in addition to a sharp peak toward IRAS 20126 at ∼4.07 μm (∼2459 cm^-1). We initially suspected these features could be attributed to crystalline CH_3OH ice, which has two sharp peaks at 3.83 and 4.07 μm (2612 and 2455 cm^-1) in the pure form. However, the entire CH_3OH combination mode complexes over the full 3.7-4.2 μm range could not be satisfactorily fit with any available crystalline CH_3OH ice spectra, pure or mixed. The observed spectra are missing an additional sharp and distinct feature at 3.94 μm (2536 cm^-1) that is present in the crystalline CH_3OH laboratory spectra. Annealed crystalline CH_3OH spectra also did not provide a good fit. Furthermore, the 3.53 μm feature lacks any significant spectral contribution from crystalline CH_3OH 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_3OH as a candidate carrier.

Finally, in both lines of sight, excess absorption emerges from the continuum between ∼4.04-4.17 μm (∼2475-2400 cm^-1) which cannot be accounted for by the CH_3OH combination mode without inducing unrealistic inflections in the local continua or overfitting the CH_3OH combination modes between 3.8-3.9 μm. There are no currently known S-H bond-bearing species whose absorptions peak at these wavelengths <cit.>, and the features present in crystalline CH_3OH ice spectra in this spectral region are all too blue-shifted from the observed peak positions at ∼4.14 μm. We therefore assign these features to HDO ice in various morphological states.



 §.§ HDO ice


We performed fits to the local continuum-subtracted 4.04-4.17 μm region using a combination of CH_3OH 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_2O combination modes in the 0.4% HDO:H_2O ice spectra (see Figure <ref>) and subsequently smoothed with Gaussian fits to eliminate experimental noise. The optical depths of the CH_3OH combination modes (as constrained by the CH_3OH column density derived from the 3.53 μm band in Section <ref>) and the Gaussian fit to the 3.9 μm excess were not varied in the fitting procedure. The same combinations of pure and H_2O-rich CH_3OH 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_2O bands toward both sources also have profiles characteristic of the lines of sight containing both amorphous and crystalline H_2O ice (Section <ref>). This correlation between the detected presence of crystalline H_2O 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_2O ice, then any HDO ice in those regions should have crystallized as well.





In total, 84 amorphous HDO ice spectra between 15-140 K and 22 crystalline HDO ice spectra between 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 χ^2 value for both lines of sight are presented in Figure <ref>. The relationships between the fit χ^2 values and the temperatures of the laboratory spectra are shown in Figure <ref>.

It is important to emphasize that the temperatures resulting in the lowest χ^2 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 <cit.>. 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 Figure <ref>, 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_2O ice, which could further alter the HDO band profile. For example, models from <cit.> and <cit.> suggest that the majority of HDO ice is co-produced with CH_3OH 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_3OH combination modes. However, such chemical environments would not be expected to affect the crystalline component because CH_3OH 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 (Figure <ref>).

Nonetheless, such an approximation is still useful for evaluating the thermal structure of the ice envelope qualitatively <cit.>. 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 suffers from much less uncertainties due to its sharper profile, matches best with annealed laboratory spectra, hinting at a possibility of heated ices being re-cooled within the ice envelope (see further discussion in Section <ref>).

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 Figure <ref>). Between 15 and 135 K, the amorphous HDO apparent band strength increases by a factor ∼1.7 from 4.2×10^-17 to 7.1×10^-17 cm molec^-1. This uncertainty is much smaller for the crystalline HDO component because its apparent band strength only varies by ∼8% between 15-150 K. When calculating the HDO column densities from the fits, we used band strengths of 6.0×10^-17 and 6.7×10^-17 cm molec^-1 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 <ref>). The resulting HDO ice column densities are listed in Table <ref>, and the reported uncertainties are explained in Appendix <ref>.



Within the uncertainties, the ratio of amorphous:crystalline HDO ice column densities are calculated to be between ∼1.2-4 for HOPS 370 and between ∼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.







 §.§ H_2O ice



To estimate the ice HDO/H_2O ratios along these lines of sight, we quantified the H_2O 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_2O ice spectra. The global continua and least-squares fits of the 3 μm feature are presented in Figure <ref>.

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_2O profile extractions much more reliable, and their H_2O profiles are similar to those in the spectra investigated in this work. The best-matching H_2O profiles found for HOPS 370 and IRAS 20126 were those of Orion from ISO (shown in ) and G034.7123-00.5946 from IRTF (shown in ), 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 Figure <ref>).

As done in previous studies <cit.>, the H_2O 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_2O lab spectrum was determined by calculating the χ^2 values of the least-squares fits of every possible combination. As stressed in Section <ref>, the temperatures of the lab H_2O 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 <cit.>, grain shape effects were not modeled in these fits, as a detailed characterization of the H_2O 3 μm band is outside our primary scope (quantifying H_2O to obtain HDO/H_2O 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 <cit.>.

The contribution of the O-H stretching mode of CH_3OH ice to the observed 3 μm feature was also approximated similarly to the HDO fitting procedure, where the column density of CH_3OH was set as a constant value from the fits to the 3.53 μm band. However, a spectrum of only pure CH_3OH was used to model the contribution of CH_3OH to the 3 μm feature, as it is difficult to accurately deconvolve the overlapping O-H stretching modes of CH_3OH and H_2O in the CH_3OH:H_2O 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^-1) 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 <ref>).

Similar to the band strength of the HDO O-D stretch, the band strength of the H_2O O-H stretch varies with temperature. The reported magnitude of this variation differs in the literature <cit.>. 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_2O, while there is only a factor of 1.16 difference between the integrated absorbance at 150 and 15 K in crystalline H_2O. These variations are again treated as a source of uncertainty in the reported H_2O column densities (see Appendix <ref>).

To quantify H_2O, all of the fit H_2O components were integrated, and a band strength of 1.9×10^-16 cm molec^-1 <cit.> was used for the amorphous component. For the crystalline component, a band strength of 2.5×10^-16 cm molec^-1 was used, which is the apparent band strength we calculated for the crystalline H_2O O-H stretch below 120 K in our laboratory spectra. The resulting values are reported in Table <ref>. Within the error margins, the ratios of amorphous:crystalline H_2O 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_2O ice for reasons discussed in Section <ref>.



§ DISCUSSION




 §.§ Inferred ice HDO/H_2O ratios





The ice HDO/H_2O ratios measured toward HOPS 370 and IRAS 20126, 4.6±2.2×10^-3 and 2.6±1.4×10^-3, 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_2O 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_2O ratios.

While the ice HDO/H_2O 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 <cit.> and are particularly promoted between water molecules following the phase transition from amorphous to crystalline <cit.>.

It may therefore be tempting to compare the ice HDO/H_2O 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_2O 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_2O), 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_2O 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 <cit.>.




 §.§ Comparing ice and gas HDO/H_2O ratios

Figure <ref> provides a comparison of the ice HDO/H_2O 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_2O ratios reported in the literature over the past couple decades span an order of magnitude from values as low as 2-3×10^-4 <cit.> to as high as 2-3×10^-3 <cit.>. The ice HDO/H_2O ratio measured in this work toward the MYSO IRAS 20126 is in agreement with the highest gas-phase MYSO HDO/H_2O 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_2O ratio on the order of 10^-3 toward a MYSO is also consistent with the conclusions drawn by <cit.>, who observed an inverse correlation between the gas HDO/H_2O ratios and the mass-weighted average envelope temperatures of the four MYSOs they studied. They interpreted this as evidence that the HDO/H_2O ratios on the order of 10^-3 (measured toward the coldest envelopes) were the most representative of ice HDO/H_2O ratios, while the HDO/H_2O ratios on the order of 10^-4 (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 <cit.>, the ^13CO_2 ice features of the two sources studied by <cit.> with lower HDO/H_2O ratios show a significant degree of CO_2 ice segregation and therefore thermal processing, while the ^13CO_2 ice feature of the <cit.> source in which the gas HDO/H_2O 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_2O ratios may greatly differ. A larger sample size of MYSO ice HDO/H_2O measurements will be needed to ascertain what causes the observed spread in MYSO gas HDO/H_2O values. In any case, the agreement between the IRAS 20126 ice HDO/H_2O ratio and the gas HDO/H_2O 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.

<cit.> additionally noted that the gas HDO/H_2O 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_3OH, 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 <cit.>. However, the most recent gas HDO/H_2O 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_2O ratios were likely overestimated by one to two orders of magnitude <cit.>. 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^-3 <cit.> 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^-4 <cit.>, and are more similar to the low-end of measured MYSO ratios.

This similarity in the recently published gas HDO/H_2O 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_2O 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_2O ratios on the order of 10^-4 made by ). Measurements of ice HDO/H_2O 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_2O/HDO ratios in MYSOs could elucidate the credibility of this latter scenario. This is because D_2O/HDO ratios, especially when compared to HDO/H_2O ratios, are also an excellent diagnostic of formation conditions <cit.>. Multiple studies have measured gas D_2O/HDO ratios in LYSOs that are more than a factor of two higher than their gas HDO/H_2O ratios <cit.>, 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_3OH ice <cit.>.

Due to the extremely low abundance of D_2O, few gas D_2O/HDO ratios have been measured in either LYSOs or MYSOs. Currently, one D_2O/HDO ratio has been measured toward the hot core of the Orion KL MYSO <cit.>. This value, 1.6×10^-3, is an order of magnitude lower than the gas D_2O/HDO ratios that have been measured in LYSOs <cit.>, which may indicate that less of the deuterated 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_2O/HDO measurements in both LYSOs and MYSOs are needed to draw such a conclusion with certainty.



 §.§ Ice envelope evolution



In addition to the ice HDO/H_2O 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.



  §.§.§ 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 non-annealed crystalline HDO ice at 150 K result in Δχ^2 values that lie outside the 3σ confidence intervals (Figure <ref>). Such an observed red-shift 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).




  §.§.§ CH_3OH ice formation



The CH_3OH ice abundances of 16% and 27% with respect to H_2O ice toward both sources are on the high end of the values typically measured for both low-mass and massive protostars ( and references therein). Although unusual, it is not unprecedented for such high values to be measured toward MYSOs <cit.>. These high CH_3OH ice abundances may be due to particularly long cold dense prestellar stages relative to the length of the warmer and less dense stages when the bulk of the H_2O ice is produced. As models suggest that the majority of HDO ice is produced during the cold dense prestellar stage with CH_3OH ice <cit.>, this would imply that the ice HDO/H_2O ratios measured here could also be on the high end of most protostars, and other protostars with lower CH_3OH ice abundances relative to H_2O ice that had shorter cold dense prestellar stages could then have lower ice HDO/H_2O ratios. However, it is also possible that thermal processing could affect CH_3OH ice abundances. Correlations between ice HDO/H_2O ratios, CH_3OH ice abundances, and the profiles of ice thermal processing tracers like ^13CO_2 may be interesting to investigate once the sample size of measured ice HDO/H_2O ratios is larger.

The fact that the current fitting procedure of the 3.5 μm absorption complex requires a component of CH_3OH in H_2O-rich mixtures also carries astrochemical significance. Although the earliest fits to CH_3OH ice features toward protostars often utilized H_2O-rich CH_3OH ice mixtures <cit.>, subsequent studies popularized an ice evolution model where most of the observed CH_3OH ice forms after CO freeze-out in a relatively CO-rich and H_2O-poor ice environment. This model is supported by analysis of the profiles of the 9.75 μm CH_3OH ice C-O stretching mode <cit.>, the "red component" of the 4.67 μm CO ice C≡O stretching mode <cit.>, the relationships between CH_3OH ice column densities, CO ice column densities, and visual extinction A_v <cit.>, and the experimentally demonstrated efficiency of solid-state CH_3OH formation via CO hydrogenation <cit.> 

However, recent experimental works indicate that CH_3OH can also be produced in H_2O-rich ices via a reaction between CH_4 and OH radicals <cit.>, albeit less efficiently than via CO hydrogenation. Recent JWST data toward prestellar cores, where the 9.75 μm CH_3OH ice feature is well-fit with a combination of CO-rich and H_2O-rich CH_3OH ice mixtures, provide observational evidence supporting this pathway as a secondary means of prestellar CH_3OH ice production <cit.>. In our fits, the lack of a CO-rich CH_3OH 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 <cit.>, leaving behind a layer of nearly pure CH_3OH ice that likely formed via CO hydrogenation, as well as some CH_3OH ice in a water-rich chemical environment that may have formed during the less dense prestellar stages, perhaps via the CH_4 + 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_3OH ice along either line of sight. Perhaps this could be explained by the heated ices having mostly experienced temperatures too hot for CH_3OH to remain in the ice (i.e., a steep thermal gradient along the line of sight). However, both the dominant component of the CH_3OH feature being pure CH_3OH 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).




§ 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_3OH ice. The H_2O ice column densities were quantified via the 3 μm feature to enable comparisons between HDO/H_2O ratios measured in the ice with those measured in the gas toward protostars. Our main findings are summarized as follows:




  * Toward both sources, the detection of a crystalline HDO ice component at 4.1 μm correlates with the detection of a crystalline H_2O 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.


  * We measure ice HDO/H_2O ratios of 4.6±2.2×10^-3 and 2.6±1.4×10^-3 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.


  * The ice HDO/H_2O ratio measured toward the MYSO IRAS 20126 is consistent with the highest reported gas HDO/H_2O ratios toward MYSOs (a couple ×10^-3), supporting the assumption that the deuterium abundances of the gas-phase water observed toward these objects are representative of the deuterium abundances of their surrounding water ice.


  * The ice HDO/H_2O ratios toward these objects are remarkably similar to each other as well as to the gas HDO/H_2O ratios measured toward isolated LYSOs. A larger sample of ice HDO/H_2O ratios, both toward LYSOs and MYSOs, as well as more measurements of D_2O/H_2O ratios in the gas phase could help to explain the cause of this similarity.


  * The peak position of the crystalline HDO ice component toward both sources suggests that the thermally processed ice in these lines of sight experienced re-cooling to some degree.


  * The CH_3OH ice requires a water-rich CH_3OH ice mixture to obtain a satisfactory fit to the 3.4 μm region, supporting recent experimental and observational evidence that some CH_3OH ice may form prior to the heavy CO freeze-out stage in an H_2O-rich environment.


  * As part of this work, several new laboratory spectra of HDO, H_2O, CH_3OH, and CO ice mixtures are now available for public use on the Leiden Ice Database (LIDA).





KS 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: https://archive.stsci.edu/doi/resolve/resolve.html?doi=10.17909/3kky-t04010.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).






aa


1. Laboratory for Astrophysics, Leiden Observatory, Leiden University, P.O. Box 9513, 2300 RA Leiden, NL

   slavicinska@strw.leidenuniv.nl

   2. Leiden Observatory, Leiden University, P.O. Box 9513, NL 2300 RA Leiden, NL

   3. Max Planck Institut für Extraterrestrische Physik (MPE), Giessenbachstrasse 1, 85748 Garching, DE

   4. European Southern Observatory, Garching, DE

   5. Department of Physics and Astronomy, Bausch & Lomb Hall, University of Rochester, Rochester, NY 14627, USA

   6. Department of Astronomy, University of Massachusetts Amherst, 710 North Pleasant Street, Amherst, MA 01003, USA

   7. Department of Astronomy & Astrophysics Tata Institute of Fundamental Research, Homi Bhabha Rd, Colaba, Mumbai, Maharashtra, IN

   8. Academia Sinica Institute of Astronomy & Astrophysics, No. 1 Sec. 4 Roosevelt Rd., Taipei 10617, TW, ROC

   9. Ritter Astrophysical Research Center, Department of Physics and Astronomy, University of Toledo, Toledo, OH 43606, USA

   10. Star and Planet Formation Laboratory, RIKEN Cluster for Pioneering Research, Wako, Saitama 351-0198, JP

   11. Department of Astronomy, University of Illinois, 1002 West Green St, Urbana, IL 61801, USA

   12. National Radio Astronomy Observatory, 520 Edgemont Rd., Charlottesville, VA 22903, USA

   13. Max Planck Institute for Astronomy, Heidelberg, Baden-Württemberg, DE

   14. SKA Observatory, Jodrell Bank, Lower Withington, Macclesfield SK11 9FT, UK

   15. Department of Astronomy, University of Texas at Austin, 2515 Speedway, Stop C1400, Austin, TX 78712, USA

   16. INAF Osservatorio Astronomico di Capodimonte, Salita Moiariello 16, I-80131 Napoli, IT

   17. Friedrich-Schiller-Universität, Jena, Thüringen, DE






§ LABORATORY DATA



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



 §.§ Experimental details

In this work, independent dosing lines were used to deposit the CH_3OH-containing binary ice mixtures using the procedure outlined in <cit.>, 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 <cit.>, 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 <cit.>, <cit.>, and <cit.>.

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 <ref>.



Because the physical properties and masses of H_2O 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_2O ice mixing ratio when pre-mixing the components prior to deposition. Therefore, the HDO:H_2O ices were created via the traditional pre-mixing method, following the procedure described in <cit.> to generate HDO via mixing H_2O and D_2O.

During the warm-ups, the ices were heated at a rate of 25 K hr^-1, and spectra ranging from 4000-500 cm^-1 (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 re-cooling of the annealed crystalline HDO:H_2O ice mixtures, a cooling rate of -2 K min^-1 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.



 §.§ 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_2O ice spectra (see Table <ref>) using a cubic spline to isolate the feature from the broad 4.5 μm H_2O combination mode followed by a Savitzky-Golay filter to smooth out experimental noise (see Figure <ref>).

The apparent band strengths of the HDO O-D stretch at different temperatures, A'_i, were derived by assuming the apparent band strength of the crystalline band at 150 K, A'_150K, is equal to 6.4×10^-17 cm molec^-1 <cit.> and using the approximation:


    A'_i = ∫ abs_i(ν)  dν/∫ abs_150K(ν)  dν× A'_150K,


where ∫ abs_150K(ν) dν is the integrated absorbance of the crystalline HDO peak at 150 K, and ∫ abs_i(ν) dν 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 <ref> and <ref> and reported in Table <ref>. 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 <ref>), 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^-17 cm molec^-1, is in excellent agreement (within ∼2%) with the band strengths reported by <cit.> and <cit.> for amorphous HDO ice at similar temperatures.












 §.§ CH_3OH C-H symmetric stretch characterization

The peak position of the CH_3OH C-H symmetric stretching mode was extracted from the CH_3OH-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 <ref> and <ref> and are plotted as a function of temperature in Figure <ref>, 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_3OH overtones and C-H asymmetric stretching modes.













§ Χ^2 PLOTS



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


    χ^2 = ∑(τ_obs - τ_fit)^2/σ^2,


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 <ref> and <ref>). The 1, 2, and 3σ bounds plotted in Figure <ref> correspond to the critical Δχ^2 (χ^2-χ^2_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).







§ ESTIMATED CONTINUUM ERRORS





 §.§ CH_3OH ice


The 20% band strength uncertainty <cit.> is propagated to determine the uncertainty in the calculated CH_3OH 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 <ref>.



 §.§ 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 <cit.> (∼10% and 5%, respectively), which were used to calculate the band strengths of HDO at higher temperatures presented in Table <ref>.

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 <ref>, and the resulting lower and upper limits on the optical depths of the 4 μm region are shown in Figure <ref>.

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.



 §.§ H_2O ice


The uncertainty in the H_2O 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_2O profiles that appeared reasonable. The lower and upper limit continua are shown in Figure <ref>, and the resulting variations in the optical depths and profiles of the 3 μm feature are shown in Figure <ref>.

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 ), and the uncertainties caused by the continuum choice, which were estimated by noting the variation in H_2O ice column density obtained from repeating the fitting procedure on the spectra that were extracted using the lower and upper limit global continua.















§ LITERATURE HDO/H_2O VALUES








