Survey of Ices toward Massive Young Stellar Objects: I. OCS, CO, OCN$^-$, and CH$_3$OH

An important tracer of the origin and evolution of cometary ices is the comparison with ices found in dense clouds and towards Young Stellar Objects (YSOs). We present a survey of ices in the 2-5 micron spectra of 23 massive YSOs, taken with the NASA InfraRed Telescope Facility SpeX spectrometer. The 4.90 micron absorption band of OCS ice is detected in 20 sight-lines, more than five times the previously known detections. The absorption profile shows little variation and is consistent with OCS being embedded in CH3OH-rich ices, or proton-irradiated H$_2$S or SO$_2$-containing ices. The OCS column densities correlate well with those of CH$_3$OH and OCN$^-$, but not with H$_2$O and apolar CO ice. This association of OCS with CH$_3$OH and OCN$^-$ firmly establishes their formation location deep inside dense clouds or protostellar envelopes. The median composition of this ice phase towards massive YSOs, as a percentage of H$_2$O, is CO:CH$_3$OH:OCN$^-$:OCS=24:20:1.53:0.15. CS, due to its low abundance, is likely not the main precursor to OCS. Sulfurization of CO is likely needed, although the source of this sulfur is not well constrained. Compared to massive YSOs, low mass YSOs and dense clouds have similar CO and CH$_3$OH ice abundances, but less OCN$^-$ and more apolar CO, while OCS awaits detection. Comets tend to be under-abundant in carbon-bearing species, but this does not appear to be the case for OCS, perhaps signalling OCS production in protoplanetary disks.


INTRODUCTION
Ices are the dominant phase of molecules such as H 2 O, CO, CO 2 , and CH 3 OH in the densest regions of molecular clouds (e.g., Chiar et al. 1995;Whittet et al. 1998;Boogert et al. 2011Boogert et al. , 2015. This paper is structured as follows. In §2, the selection of MYSO targets from the Red MSX Survey is described. The IRTF/SpeX observations and data reduction process are discussed in §3. The continuum is determined in §4.1, the removal of contaminating gas phase CO lines is addressed in §4.2, and the ice band absorption profiles are fitted and column densities are derived in §4.3. Laboratory spectroscopy of solid state OCS, taken from the literature, is summarized in §4.4. In §4.5 we search for correlations between the column densities of H 2 O, OCN − , OCS, CO, and CH 3 OH. The ice abundances are compared across MYSOs, low mass YSOs. dense cloud background stars, and comets in §4.6. In the discussion, we address the chemical origin of OCS and OCN − ( §5 and §5.1), and the origin of cometary ices ( §5.2). The conclusions and future work are summarized in §6.

SOURCE SELECTION
The embedded MYSOs for our 2-5 µm spectroscopic survey were selected from the Red MSX Survey (RMS; Lumsden et al. 2013, Pomohaci et al. 2017. The RMS survey contains more than 700 YSOs and candidate YSOs across the Galaxy. We selected those that are visible from Maunakea at airmass values of < 3 (Declination > −50 o ). They must also be MYSOs in an early evolutionary stage, enhancing the likelihood that strong ice absorption features from their massive envelopes are detectable, much like previously studied MYSOs (Gibb et al. 2004). For this, we selected the most luminous objects (L >3000 L ⊙ ) that are radio quiet (<0.5 Jy at 5 GHz). As such, they likely have not yet started to ionize their surroundings and produce an HII region. This sample was further reduced by selecting only those targets with K−band magnitudes brighter than 13, which facilitates on-slit guiding with the IRTF. This yielded a list of 215 targets. For about half of these, 1-5 µm spectra were obtained with IRTF/SpeX. This final sub-selection was based on the visibility during the observing runs, and on the J − K color. Preference was given to targets with J − K > 5, i.e., with the steepest spectral energy distributions. This further enhanced the likelihood of high quality spectra in the M−band, for which the IRTF/SpeX sensitivity is strongly background-limited. This increases the bias towards the most embedded targets with the deepest ice absorption features. The full sample will be presented in an upcoming paper (K. Emerson, in preparation). Here, 23 targets (Table 1) are presented for which the 4.90 µm ice absorption band of OCS was detected, or hints of it were seen at the 2-3σ level.

OBSERVATIONS AND DATA REDUCTION
All targets were observed with the SpeX spectrometer  at the NASA Infrared Telescope Facility (NASA/IRTF) in a number of nights in the years 2020 and 2021 (Table 1). K-band images were taken with the SpeX guider camera and then compared with the target designations in the RMS database in order to identify the MYSO. In several cases, the field was rather confused, and we indicate the selected target in the last column of Table 1. Subsequently, on-slit K-band guiding was enabled, and spectra were taken with the SpeX LXD Long mode, using the 0.3 or 0.5 arcsec wide slits (Table 1). This yields resolving powers of R = λ/∆λ =2,500 and 1,500, respectively. The instantaneous spectral coverage is 1.95-5.36 µm.
b Date of the IRTF/SpeX observations.
c Telluric standard star.
d The IRTF/SpeX 0.5 ′′ slit is used unless noted otherwise.
usually not nearby. The standard stars we selected include HR 6378 (A2IV-V), HR 6629 (A1Vn), HR 7236 (B8.5V), HR 7235 (A0IV-Vnn), HR 8585 (A1V), HR 8371 (B8Ib), HR 8028 (A0IIIn), and HR 2421 (A1.5IV), where the spectral types were obtained from SIMBAD 2 . The IRTF/SpeX spectra were reduced using Spextool version 4.1 (Cushing et al. 2004). Flat fielding was done using the images obtained with SpeX's calibration unit. The wavelength calibration procedure uses lamp lines at the shortest wavelengths and sky emission lines in much of the L and M-bands. The telluric correction was done using the Xtellcor program (Vacca et al. 2003). This uses a model of Vega to divide out the stellar photosphere. Although the spectral types of the standard stars were not exactly of type A0V, no photospheric residuals were observed that could affect our data analysis. The final continuum signal-to-noise values near the location of the 4.90 µm absorption feature is on average 45, determined from the scatter in the data after the correction for gas phase CO absorption ( §4.2). The lowest signal-to-noise value is 15 (G111.2552) and the highest 105 (G029.8620).

RESULTS
The observed MYSOs have steeply rising spectra in the 1.95-5.36 µm wavelength range (Figs. 1-3). All have deep absorption features centered on 3.0 µm, covering the wavelength range of 2.8 to ∼3.7 µm. At the short wavelength side, the band is cut off by the opaque Earth's atmosphere. The 3 µm band is mostly attributed to H 2 O ice, with minor contributions from NH 3 , NH 3 hydrates, CH 3 OH, and likely the O-H and C-H stretch mode of other, yet unidentified, species.
The M-band spectra show relatively narrow absortion features due to CO (4.67 µm), OCN − (4.62 µm), and OCS (4.90 µm). Most MYSOs also show the presence of the shallow, broad absorption by the combination mode of H 2 O ice. In addition, a number of targets show a weak change in slope at ∼ 5.0 µm. This seems to be due to a new, previously unreported absorption feature centered at ∼5.25 µm. Figure 4 illustrates the presence of all these spectral features by putting them in a wider context using an Infrared Space Observatory/Short Wavelength Spectrometer (ISO/SWS) spectrum (Gibb et al. 2004) for comparison. Many targets also show the unresolved absorption lines of the ro-vibrational transitions of gas phase CO. For a few targets, these lines are in emission.
Finally, a number of MYSOs show unresolved emission lines due to hydrogen Brγ, Brα, and Pfβ. These are not analyzed here, and their wavelength ranges are excluded from the continuum and ice absorption feature analysis.

Continuum Determination
To be able to analyze the profile and depth of the ice absorption features, a global baseline defining the continuum emission was determined over the full observed wavelength range of 1.95-5.36 µm. This baseline was determined by fitting a simple model of two reddened blackbody functions to the data. The first blackbody is that from a hot central source. Its temperature is assumed to be at least  A change in slope is visible near ∼5.0 µm in both the IRTF/SpeX and ISO/SWS spectra, due to what appears to be a shallow absorption feature centered on ∼5.25 µm. The change in slope at ∼ 5.5 µm is due to the onset of the prominent 6.0 µm absorption band. The apparent lack of absorption by gas phase CO lines in the ISO/SWS spectrum is due to its low spectral resolving power (R = 800 versus R =1,500). several tens of thousand Kelvin, so that it is represented by a temperature-independent Rayleigh-Jeans tail in the observed wavelength range. The second blackbody represents the emission by warm dust along the line of sight, perhaps from a circumstellar disk or the inner envelope. Both blackbodies are reddened by the same foreground extinction due to cooler envelope and foreground cloud material. The extinction curve used for this is that derived for the J, H, and K-bands by Indebetouw et al. (2005), and at longer wavelengths it is assumed to be flat. Furthermore, the 3.0 µm ice band was also included in the model, because H 2 O absorption covers such a large portion of the observed wavelength range. For this, we used template 3.0 µm band spectra derived from previously well studied MYSOs (Gibb et al. 2004), covering a range of profiles (AFGL 2136, MonR2 IRS3, AFGL 989, NGC 7538 IRS9, S140 IRS1). This model thus contains four variables: the temperature of the warm dust emission, the foreground extinction A V , the peak optical depth of the 3.0 µm ice band τ 3.0 , and a relative flux scaling factor for the hot (Rayleigh-Jeans) and warm blackbody emission. This model is then normalized to the observed spectrum near 5 µm and fitted to selected wavelength ranges. Those ranges avoid the long-wavelength wing of the 3.0 µm ice band (∼ 3.2-3.9 µm), the OCN − , CO, and OCS ice features (4.5-4.9 µm), as well as regions with poor atmospheric transmission. The global baseline model generally fits the observations well ( Figs. 1-3). The spectra were subsequently converted to optical depth scale. At a detailed level in the M-band optical depth spectra, small offsets from the zero optical depth level were observed, however, and these were removed by subtracting an additional local straight line fit, shown in Figs. 1-3.
Because it is our primary goal to determine the baseline for the ice features, we do not present or further interpret the parameters and their uncertainties for the blackbody fit parameters. Due to the small width of the ice features in the M-band, and the subtraction of an additional local baseline, the uncertainty due to the baseline is small. An error of 0.01 is included in the uncertainties derived for the ice band optical depths. Uncertainties for the 3.0 µm H 2 O ice band optical depths are larger (∼ 10%). This was checked by using a different method to determine the baseline for the 3.0 µm band, applying a spline function (K. Emerson et al., in preparation). These uncertainties are included in the derived H 2 O ice column densities ( §4.3.5).

Gas Phase CO
Many targets show absorption, and some emission, by the ro-vibrational lines of gas phase 12 CO, and sometimes 13 CO. These lines are unresolved at the IRTF/SpeX resolving power of R = 1, 500−2, 500. At the level of a few percent of the continuum, these lines somewhat affect the profiles of the ice absorption features. We therefore attempt to remove them using a rotation diagram analysis (Fig. 5). Integrated absorption line depths were measured, and converted to lower level column densities assuming optically thin absorption in Local Thermodynamic Equilibrium (LTE). The constructed rotation diagrams typically show a steep linear relation at low energy levels, and a shallower linear relation at higher energy levels. Straight lines were fitted to each of these, providing integrated line depths for all rotational levels. These were then used to divide out the gas phase CO lines. The detected 12 CO lines are most likely optically thick, and we do not further use the column densities and excitation temperatures, derived from the rotation diagrams in the optically thin limit, in this paper.

Ice Band Profiles and Column Densities
The absorption features detected at ∼4.62, 4.67, and 4.90 µm are shown in Figures 6-8 on optical depth scale. They are ascribed to OCN − , CO, and OCS, respectively, and are fitted with Gaussian and Lorentzian functions to characterize their profiles and depths. To put the results in an astrochemical context, the 3.0 and 3.53 µm absorption bands of H 2 O and CH 3 OH ice are discussed as well. However, for a more detailed analysis of these two bands, in a much larger sample of MYSOs, we refer to K. Emerson et al. (in preparation).

CO
The 4.67 µm absorption band of CO ice shows significant profile variations (Figs. 6-8). For some targets the profile is much narrower (e.g., G063.1140) than for others (e.g., G025.4118A). Such variations were previously observed for both high and low mass YSOs, and were ascribed to the presence of CO in different molecular environments (Tielens et al. 1991;Boogert et al. 2002;Pontoppidan et al. 2003): a central, narrow absorption due to pure CO (and traces of other low dipole moment molecules), a broad red-shifted absorption due to CO mixed with high dipole moment species (H 2 O, CH 3 OH), and a minor blue-shifted absorption, likely due to CO mixed wth CO 2 ice. By fitting constant profiles for CO apolar , CO polar , and CO blue , and only varying the peak optical depths, good fits are obtained across a wide range of low and high mass YSOs as well as background stars (Pontoppidan et al. 2003). Figure 5. Illustration of the removal of gas phase 12 CO ro-vibrational absorption lines for the MYSO G108.7575−00.9863. Panel a is the IRTF/SpeX spectrum before removal, and panel b after. Panel c shows the rotation diagram used in this process, with E low , N low , and J low the energy level, column density, and rotational quantum number for the lower energy levels, respectively. Two slope regimes are discerned, which are fitted with two straight lines, corresponding to excitation temperatures of 90 and 895 K and column densities of 7.5×10 16 and 2.5×10 17 cm −2 , respectively. These temperatures are likely highly inaccurate, and the column densities strongly underestimated due to the lines being optically thick.
The same approach is adopted for the MYSOs in this work. CO apolar was fitted using a Gaussian function (peak wavelength 4.6731 µm, FWHM=0.0076 µm), CO polar with a Lorentzian (4.6806 and 0.0232 µm), and CO blue with a Gaussian (4.6648 and 0.0065 µm). The fitting process was done simultaneously with the 4.62 µm band of OCN − ( §4.3.2), due to the wavelength overlap. Also, all components were convolved with the Gaussian SpeX instrument profile before fitting them to the data. Satisfactory fits are obtained for all bands (Figs. 6-8). The total CO ice column densities and those of the different components are listed in Table 2. An integrated band strength of 1.1×10 −17 cm/molecule, as measured in the laboratory (Jiang et al. 1975), was used for this. Figure 6. M -band spectra of the MYSO sample on optical depth scale. A zoom in on the OCS ice band is shown in the right panels. The OCS band is fitted with a single Gaussian. To this are added the Gaussian and Lorentzian fits to the CO ice bands (green and black, respectively) and Gausian fits to the 2165.7 and 2175.4 cm −1 components of the OCN − band (solid purple and dashed purple, respectively). A dashed red line in the right panels indicates that OCS was not detected and the Gaussian depth should be regarded as an upper limit.

OCN −
The 4.62 µm band of OCN − was detected toward all targets. In a study of low and high mass YSOs, it was found that this band consists of two components, centered on 2165.7 and 2175.4 cm −1 (4.6174 and 4.5969 µm; van Broekhuizen et al. 2005). The 2165.7 cm −1 component was found to be consistent with solid OCN − in polar environments, while the other component relates perhaps to CO bonding to the grain surfaces or to OCN − in apolar environments (Öberg et al. 2011). We adopt the same decomposition procedure, by fitting two Gaussians, one at 2165.7 cm −1 (F W HM =26 cm −1 ) and one at 2175.4 cm −1 (F W HM =15 cm −1 ). The fitting was done simultaneously with those of the CO ice feature ( §4.3.1), since they overlap. Good quality fits were obtained (Figs. 6-8). The contribution by the 2175.4 cm −1 component is usually at the less than 2σ confidence level, except for three targets at 2-3σ, and three targets at > 3σ. The OCN − column densities are listed in Table 3. An integrated band strength of 1.3×10 −16 cm/molecule, as measured in the laboratory (van Broekhuizen et al. 2004), was used for this.

OCS
The absorption feature at 4.90 µm, attributed to the C-O stretch mode of solid OCS, is well fitted with a single Gaussian for all . Upper limits are derived for three of the targets. The fit parameters are listed in Table 4. Solid OCS column densities were derived using an integrated band strength of 1.5×10 −16 cm/molecule, as measured in the laboratory (Palumbo et al. 1997). In order to highlight similarities and differences among the targets, each target is compared to the same Gaussian in Figure 9. Clearly, any profile variations are small.

CH 3 OH
To put the OCS, OCN − , and CO observations in an astrochemical context, CH 3 OH ice column densities were derived as well. The 3.53 µm C-H stretch mode of solid CH 3 OH is detected towards 15 of the 23 MYSOs. As shown in Figures 10-12, the CH 3 OH band overlaps with the 3.47 µm absorption feature, which is possibly due to NH 3 hydrates (Dartois & d'Hendecourt 2001). These absorption features are located on the wing of the deep 3.0 µm band of H 2 O ice. A low order polynomial was used to define and subtract a local baseline. The absorption features were then disentangled by fitting the combination of a Gaussian and a laboratory spectrum. The laboratory spectrum used is the mixture H 2 O:CH 3 OH:CO:NH 3 =100:50:1:1 at a temperature of 10 K (Hudgins et al. 1993). To provide more stable solutions, the Gaussian peak position and F W HM were kept fixed at 3.470 µm and 0.11 µm, respectively. CH 3 OH column densities were derived by using an integrated band strength of 5.6×10 −18 cm/molecule (Kerkhof et al. 1999) for the C-H stretch mode, a F W HM of 0.04 µm, and the peak optical depth provided by the fits. The uncertainties were estimated to be ∼20% due to the baseline and ∼10% due to the F W HM. The column densities are listed in Table 3. Figure 9. The OCS ice band towards the sample of MYSOs compared to the same Gaussian, scaled to the best-fit peak optical depth of Table 4 in order to highlight any variations in the sample. The Gaussian represents the best fit to the target G012.9090 (W33A; λ =4.9012 µm, FWHM=0.0441 µm).

H 2 O
To be able to derive relative ice abundances, column densities for the most abundant ice species, H 2 O, were derived as well. This was done as part of the continuum fitting process ( §4.1). An integrated band strength of 2.0×10 −16 cm/molecule (Hagen et al. 1981) was used. The column densities reflect those of the small grains portion of the 3.0 µm band, excluding the long-wavelength wing, which may have a varied and uncertain origin. This approach was also adopted for low mass YSOs and background stars in previous work (e.g., Boogert et al. 2008). For a more in depth analysis of the profile of the 3.0 µm band of this and a larger sample of MYSOs, we refer to K. Emerson et al. (in preparation). The H 2 O column densities are listed in Table 2.

OCS Laboratory Spectroscopy
The spectra of solid OCS in a range of astrophysically relevant ice mixtures and temperatures were measured by Hudgins et al. (1993), Palumbo et al. (1995), Ferrante et al. (2008), and Garozzo et al. (2010). Peak positions and FWHM values of the C-O stretch mode of OCS for non-irradiated mixtures are displayed in Figure 13a. The peak position depends strongly on the ice composition: mixtures with apolar species peak near 4.87-4.88 µm, mixtures with H 2 O peak at 4.88-4.89 µm, and mixtures with CH 3 OH near 4.90 µm. The profiles of the apolar mixtures become much broader (factor of 2-3) upon heating, while the polar mixtures are very broad at temperatures of 10 K, and then dramatically narrow at higher temperatures.
The proton-irradiated ices studied in Ferrante et al. (2008), generally exhibit a different behavior (Figure 13b). This is confirmed by the proton-irradiation experiments of Garozzo et al. (2010). Upon heating, the peak position shifts to much longer wavelength (up to 4.915 µm, the longest of all experiments), while its width does not change much, except for a narrowing for the H 2 O:CO:SO 2 ice mixture. The strong shift must at least in part be related to the formation of new species in the ice, such as CO 3 , CS 2 , and CO 2 , changing the dipole interactions. Figure 10. Portion of the L−band spectra of the MYSO sample. All targets show the 3.47 µm absorption band, likely caused by NH 3 hydrates, and some show signs of hydrocarbon absorption as well near 3.4 µm (e.g., G033.5237). The feature at 3.53 µm is attributed to CH 3 OH ice. For each target, the left panel shows the adopted local polynomial baseline, and the right panel the baseline subtracted spectrum. In the right panels, the solid green line is a laboratory spectrum of CH 3 OH ice at a temperature of 10 K (H 2 O:CH 3 OH:CO:NH 3 =100:50:1:1; Hudgins et al. 1993) and the dashed green line a Gaussian. The sum of these fits is shown as a solid red line. For targets without a CH 3 OH ice detection, a dashed red line is shown instead.   Palumbo et al. (1995), andFerrante et al. (2008). Different ice mixtures are distinguished by different symbols and colors. Mixing ratios are labeled relative to OCS (e.g., "CH 3 OH-10" means a mixing ratio of CH 3 OH:OCS=10:1). For a giving ice mixture, different temperatures are connected with a line. For each temperature series, the lowest temperature mixture is indicated with an open symbol. The unlabeled dark blue squares indicate a mixture of H 2 O:CH 3 OH:OCS=3:10:1, and the dark blue triangles of H 2 O:CH 3 OH:OCS=100:50:1. The peak position changes typically by the amount indicated by the black arrow if the OCS is present at concentrations > 10% on grains with shapes of a Continuous Distribution of Ellipsoids (CDE), rather than on a flat surface as in the laboratory transmission spectra. The black bullets with error bars indicate the peak positions and widths observed toward the sample of MYSOs. Panel (b): Similar to panel (a), but here the observations (black bullets) are compared to the proton-irradiated mixtures of Ferrante et al. (2008). Figure 13 also shows a comparison to the observations. Of the non-irradiated mixtures, ices rich in CH 3 OH typically give the best match of peak position and FWHM. Of the mixtures without CH 3 OH, only those with relatively high OCS/H 2 O (∼2) provide a match to some of the observed MYSOs, but this deteriorates when grain shape effects are taken into account, which typically shift the absorption band to shorter wavelengths. Grain shape effects are significant if OCS has a relatively high concentration (≥10%). The best match in peak position and FWHM is in fact provided by the proton-irradiated CO:H 2 S ices, as those are the only ones that cover the longest observed peak positions near 4.905 µm. This might point to a relatively complex molecular environment for OCS.

Ice Correlations
To further constrain the molecular environment and formation history of OCS, its column density is plotted versus that of other species (Figure 14), and correlation coefficients are derived. A positive Spearman's Rank coefficient (ρ) indicates the probability that there is a monotonically increasing relation, with a maximum of 1.0. This value is indicated in each panel, with the second number the significance of there being no correlation. The best correlations are found for OCS with the OCN − and CH 3 OH column densities (ρ ∼0.72), followed by CO polar (ρ ∼0.65). The correlations are rather poor for OCS with the H 2 O column (ρ ∼0.28) and there is no correlation with CO apolar .
We investigated whether the rather poor correlation of OCS with H 2 O could be an artifact related to underestimated H 2 O column densities (and uncertainties) due to saturated 3.0 µm absorption bands. Inspection of Figures 1-3 shows that this could be the case for five MYSOs. But of those, G020.7617 was not included in the OCS version H 2 O corelation due to the non-detection of OCS, G111.2552 was not included due to severe saturation of the 3.0 µm band, and G012.9090 is a well studied MYSO (W33A) for which N(H 2 O) is quite reliable (e.g., Gibb et al. 2004). Excluding the remaining two targets (G105.5072 and G059.7831) from the correlation plot reduces the Spearman's Rank coefficient from 0.28 to 0.21. Including those two targets, and artificially increasing their column density with a factor of 1.6 to force an improved correlation, improves the Spearman's Rank coefficient somewhat (0.34), but it is still much poorer than the correlations of OCS with CH 3 OH, OCN − , and CO polar .
The correlations among CH 3 OH, OCN − , and CO polar are also investigated (Fig. 15). Indeed, as expected, the best correlation is that between CH 3 OH and OCN − (ρ ∼0.64). The correlations of CO polar with CH 3 OH (ρ ∼0.43) and with OCN − (ρ ∼0.57) are poorer than that of OCS with CO polar (ρ ∼0.65).
In summary, OCS, CH 3 OH, and OCN − correlate best with eachother, and somewhat poorer with CO polar . There is no correlation of these species with H 2 O and CO apolar . Despite the reasonably good correlations found for OCS with the OCN − and CH 3 OH column densities, a linear model does not describe the observations well. The scatter in the datapoints is much larger than the uncertainties. At any given CH 3 OH or OCN − column density, the OCS column density varies by a factor of 2-3.

Abundance Distributions
The distribution of the ice abundances for the MYSOs and other target categories is shown in histograms (Figs. 16-18). The corresponding median abundances and upper and lower quartiles are listed in Table 5, following the format in (Öberg et al. 2011) and (Boogert et al. 2015). The distribution of the OCS ice abundances relative to CH 3 OH is much narrower compared to abundances relative to H 2 O ice (Fig. 16), consistent with the much better Spearman's Rank coefficient (Fig. 14). It also shows the large improvement in sample size relative to previous work. No detections or significant upper limits of OCS ice abundances were made towards low mass YSOs and dense cloud background stars (Palumbo et al. 1997;Boogert et al. 2000;Gibb et al. 2004). Towards comets, OCS abundances relative to H 2 O where published by Saki et al. (2020) and Figure 16 shows a distribution quite similar to that for MYSOs.
The abundance of CH 3 OH relative to H 2 O ice towards MYSOs appears to be fairly similar, or perhaps a little larger, compared to low mass YSOs and dense cloud background stars (Fig. 17). For OCN − relative to H 2 O ice, a fairly broad distribution is observed for MYSOs, with a peak at higher abundances compared to low mass YSOs. Towards dense cloud background stars, OCN − has not yet been detected, with upper limits comparable to the detections towards low mass YSOs (Whittet et al. 2001b;Chu et al. 2020).
For CO ice, the abundances relative to H 2 O ice are quite similar for MYSOs, low mass YSOs, and dense cloud background stars (Fig. 18). But while the abundances of the polar component of CO relative to H 2 O ice are also quite similar across the different classes of objects, this is not the case for the apolar CO phase. Previous work on MYSOs (red hatched histograms in Fig. 18) showed a deficiency of the total CO and polar CO ices relative to low mass YSOs and background stars. In the Figure 16. Distribution of the OCS ice abundance relative to H 2 O ice towards MYSOs (panel (a)). The solid black histogram is for MYSO detections presented in this work. The unfilled histogram includes abundance upper limits. The hatched red histogram is for all OCS ice detections of previous work. Panel (b) shows the OCS abundance in comets (Saki et al. 2020), with the open histogram including upper limits. Panel (c) shows the OCS ice abundances relative to CH 3 OH ice for the MYSOs presented in this (black) and previous (red) work. The solid black histogram is for MYSO detections presented in this work. The hatched red histogram is for all CH 3 OH ice detections of previous work. Panel (b) shows the CH 3 OH ice abundances toward low mass YSOs (black) and comets (blue hatched). Panel (c) shows CH 3 OH ice abundances relative to H 2 O ice for lines of sight tracing quiescent dense cloud material. In all panels the unfilled histogram includes abundance upper limits. Panels (d), (e), and (f ) show the abundances of OCN − ice relative to H 2 O ice for MYSOs, low mass YSOs, and background stars, respectively. current sample, this is no longer the case, their abundances are comparable. However, the deficiency of the apolar CO ice component in the current MYSO sample is consistent with previous work. The solid black histogram is for MYSO detections presented in this work. The hatched red histogram is for all CO ice detections of previous work. Panel (b) shows the CO ice abundances toward low mass YSOs (black) and comets (blue hatched). Panel (c) shows CO ice abundances relative to H 2 O ice for lines of sight tracing quiescent dense cloud material. In all panels the unfilled histogram includes abundance upper limits. Panels (d), (e), and (f ) show the abundances of CO ice in the polar phase relative to H 2 O ice for MYSOs, low mass YSOs, and background stars, respectively. Panels (g), (h), and (i) show the abundances of CO ice in the apolar phase towards these same target categories. · · · · · · 0.14 0. b Saki et al. (2020) c Disanti & Mumma (2008) d Disanti & Mumma (2008) and Saki et al. (2020)

DISCUSSION
The presented infrared spectroscopic survey of solid OCS towards 23 MYSOs provides observational constraints to the formation environment of sulfur-bearing species in interstellar and circumstellar ices. When compared with un-irradiated laboratory ice mixtures, the observed 4.90 µm C-O stretch-ing mode band profiles of OCS are most consistent with CH 3 OH-rich OCS ice mixtures, and do not favor H 2 O, CO 2 , and CO-rich environments. This confirms the early work by Palumbo et al. (1997) on a sample of three MYSOs. In our much larger sample, the identification of OCS in CH 3 OH-rich mixtures is supported by a good correlation of the OCS and CH 3 OH column densities, indeed much better than that with the H 2 O and CO column densities.
Observations and models have shown that CH 3 OH formation is strongly enhanced deep into dense clouds (A V > 9, densities ≥ 10 5 cm −3 ), once much of the CO freezes out (Pontoppidan et al. 2004;Cuppen et al. 2009;Boogert et al. 2011;Chu et al. 2020). The good correlation of the OCS and CH 3 OH column densities would then indicate that OCS is preferably formed in this later, dense stage of molecular clouds, or the envelopes of YSOs, where similar conditions occur. A direct accretion of OCS from the gas phase can be excluded, considering that gas phase OCS abundances are ∼ 1 × 10 −5 relative to CO, orders of magnitude smaller than observed in the ices ( Table 5). The route to solid state OCS might then be the oxidation of CS on the grains. At the same time, the oxidation of CO would likely produce CO 2 , and indeed, the presence of CO 2 heavily diluted in CH 3 OH ice is indicated by a long-wavelength shoulder on the 15 µm bending mode of solid CO 2 ice (Dartois et al. 1999;Gerakines et al. 1999). The abundance of this CO 2 relative to the CH 3 OH column density varies between 5-30% , indicating that such oxidation reactions in the CO-rich ice phase are significant.
CS and CO are formed early in the cloud evolution, by gas phase ion-molecule reactions, when at low densities, C, S, and O are mostly in atomic form (e.g., Palumbo et al. 1997). Later gas phase chemistry does not modify CS, but converts S to SO. On the grain surfaces, early oxidation of CO and CS would make CO 2 and OCS and oxidation of S and SO makes SO 2 . Evidently, at the same time, abundant H 2 O ices are produced, and some CH 3 OH formation is expected as well. Therefore, if OCS formation takes place early in the cloud's lifetime, at relatively low densities (∼ 10 3 cm −3 ), it is expected to be mixed with H 2 O, CO 2 , and CH 3 OH. However, while the CO 2 abundance correlates very well with that of H 2 O ice (N(CO 2 )/N(H 2 O)∼ 20-30%, e.g., Gerakines et al. 1999), CH 3 OH and OCS ice have not been detected at the level of < 1% in these early ices. Our observations are consistent with the formation of the bulk of these species much deeper in the cloud.
A key question is if there is sufficient CS to produce the OCS abundance observed for the MYSO sample. In diffuse clouds and dense clouds, the gas phase CS/CO abundance ratio is ∼ 5 × 10 −5 (Laas & Caselli 2019 and references therein). The observed median abundance of OCS ice relative to CO ice is two orders of magnitude larger at 6×10 −3 ( Table 5). Much of the CO was converted to CH 3 OH and CO 2 , however. The ice abundances of CO, CH 3 OH, and CO 2 are comparable (Table 5, Gerakines et al. 1999), resulting in a OCS/(CO+CH 3 OH+CO 2 ) ice ratio of 2×10 −3 . This is still a factor of ∼ 40 larger than the CS/CO ratio produced by gas phase chemistry. In addition, much like CO ice is converted to H 2 CO and CH 3 OH, CS would be converted to H 2 CS and CH 3 SH as well (Lamberts 2018). Observations toward a comet (Calmonte et al. 2016;Altwegg et al. 2019) and the envelopes of Class 0 protostar (Drozdovskaya et al. 2019) show N(CH 3 SH)/N(CS)∼ 1. The order of magnitude enhancement of CS due to grain surface chemistry seen in the models of Laas & Caselli (2019), would not fully remedy this. Overall, CS appears to be only a minor precursor to OCS ice.
Sulfurization of CO ice is a more likely route to OCS. Laas & Caselli (2019) find that S atom addition reactions are dominant on the grains. Sulfur allotropes, such as suggested for the formation of abundant gas phase SO 2 in a hot core (Dungee et al. 2018), could also play a role. H 2 S seems a reasonable candidate as well. The upper limits on the H 2 S ice abundance (<0.3-1% relative to H 2 O; Smith 1991) are larger than the OCS ice abundances and thus the observations currently do not constrain this possibility.
The production of OCS by energetic process may play a role as well. The peak positions and widths of the 4.90 µm band toward many of the observed targets are consistent with the proton-irradiated CO:H 2 S=5:1 and CO:SO 2 =5:1 ices, with or without H 2 O, after a modest amount of heating to temperatures of <50 K ( Fig. 13; Ferrante et al. 2008; see also Garozzo et al. 2010). Such irradiation would produce much CO 2 , and minor amounts of C 3 O 2 , C 2 O, HCO, CS 2 , and for the SO 2 -containing ices also O 3 . Other than CO 2 , which is also abundantly produced by grain surface chemistry, the latter species are not detected in the MYSO spectra, but not with significant upper limits. Indeed, in the irradiation experiments, the absorption strengths of the other species observed is much less than that of the 4.90 µm OCS band. A deep search with JWST is needed to further constrain the irradiation history of the interstellar ices.
The starting mixture in these irradiation experiments is CO:H 2 S=5:1 and CO:SO 2 =5:1. The presence of H 2 S in CO-rich ices locates the ices at lower densities, where H 2 S would be formed via hydrogenation reactions on the grain surface, while the good correlation of OCS with CH 3 OH indicates the formation during a high density phase. SO 2 , on the other hand, was suggested to be present in a CH 3 OH-rich environment, following a tentative detection by Boogert et al. (1997). At an SO 2 abundance of 0.3% relative to H 2 O towards the MYSO W33A (G012.9090-00.2607), the abundance relative to CO ice is 3.5%. This is within an order of magnitude of the composition of the initial laboratory ice, and thus irradiation of SO 2 -containing ices could be the source of the observed OCS ices. A more quantitative assessment is beyond the scope of this work, and would also benefit from SO 2 ice measurements in a wider variety of sight-lines. Also, other S-bearing species in the CO-rich ice, such as CS and S allotropes may produce OCS in the irradiation process as well. Indeed, free S atoms readily react with CO to form OCS (Ferrante et al. 2008).
Finally, we note that despite the correlation of OCS column densities with the species mentioned above, the scatter in the correlations is of the order of a factor of 2. This implies that the OCS production is quite strongly dependent on specific local conditions (radiation fields, ice temperatures). This is reminiscent of conclusions that were drawn for the production of OCN − (van Broekhuizen et al. 2005).

OCN −
The correlation of OCS with OCN − is as good as that with CH 3 OH (Fig. 14), and indeed OCN − and CH 3 OH also correlate well with eachother (Fig. 15). In a study of low mass YSOs, van Broekhuizen et al. (2005) found that OCN − is present in both CO-rich and polar environments, claiming that OCN − is formed in a CO-rich environment (short-wavelength component) and survives in a polar-CO environment (long-wavelength component; see alsoÖberg et al. 2011). This is confirmed by a correlation of the short-wavelength component with CO in an apolar environment, and of the long-wavelength component with polar CO ices. For our MYSO sample, the short wavelength component is detected in just a few cases. The long-wavelength component indeed correlates much better with the polar than with the apolar ices. Thus, overall the idea of OCN − being produced in a CO-rich environment is confirmed. For our MYSO sample, much of the CO has already been converted to other species.
While it was shown that OCN − is easily produced at low temperatures by acid base chemistry of HNCO and NH 3 mixtures (Raunier et al. 2004), the apolar ices are expected to be NH 3 -poor, just like they are H 2 O-poor. Perhaps another base enables this reaction. On the other hand, it was shown that the presence of N in CO-rich ices leads (non-energetically) to the formation of both HNCO and NH 3 (Fedoseev et al. 2015) as a result of hydrogen addition and abstraction reactions.

Comparing YSOs, Dense Clouds, and Cometary Ices
In a statistical comparison of YSO and cometary ice abundances, Öberg et al. (2011) concluded that comets are carbon and nitrogen-poor. This could hint to the effects of ice processing in protoplanetary disks, or that compared to the distribution of ice abundances of low mass YSOs, the early Sun resembled the YSOs with low C and N ice abundances. One possible cause of low C and N abundances could be the presence of nearby massive stars, that sublimate or photodesorb the most volatile (CO, N 2 ) ices. The Sun was indeed most likely formed in a cluster (Adams 2010).
Re-assessing the comparison between YSOs and comets, using the vastly increased MYSO sample, we find that, in contrast to the lower abundances in previous work, MYSOs have CO ice abundances that are comparable to those of low mass YSOs and background stars (Fig. 18a-c). About 50% of the comets have CO ice abundances below those of MYSOs and low mass YSOs. The large heterogeneity of the cometary CO ice abundances is also apparent, probably reflecting their processing history. Comets are also particularly poor in CH 3 OH ices (Fig. 17)), especially when compared to the new sample of MYSOs.
In contrast to carbon-bearing molecules, comets do not show a deficiency in the OCS abundance (Fig. 16)). This is surprising, considering that we find that OCS is mixed with CH 3 OH-rich ices. This could point to the additional production of OCS in the protoplantary disk phase. Studies of OCS in low mass YSOs will be enabled by JWST. A tentative detection by the AKARI space telescope towards the edge-on disk IRC L1041-2 indicates a very high OCS ice abundance of 1.1% relative to H 2 O (Aikawa et al. 2012), which is 50% more than the MYSO with the largest OCS abundance (G010.8856+00.1221; Table 4) As for the comparison between MYSOs, low mass YSOs and dense cloud background stars, MYSOs are CH 3 OH and OCN − -rich. The total and polar CO ice abundances are quite similar, but MYSOs do appear to have significantly less apolar ices. This is consistent with a more efficient conversion of CO to CH 3 OH and OCN − in the ices. One should keep in mind, however, that the low mass YSO and dense cloud background star samples are presently relatively small.

CONCLUSIONS AND FUTURE WORK
Infrared 2-5 µm spectra of a sample of twenty-three deeply embedded MYSOs were studied with the primary goal to study the abundance of solid state OCS via the C-O stretch mode at 4.90 µm. The following conclusions are drawn: 1. The OCS band at 4.90 µm was detected in twenty MYSO lines of sight, increasing the sample by more than a factor of 5.
2. The absorption profile of the OCS band shows little variation from source to source and is consistent with mixtures of OCS embedded in CH 3 OH-rich ices, or mildly heated protonirradiated CO:H 2 S ices.
3. All twenty-three MYSOs show CO and H 2 O ice absorption, and most show absorption by OCN − and CH 3 OH.
4. The OCS column density correlates well with that of CH 3 OH and OCN − , and polar CO ice, but not with H 2 O and apolar CO ice. This firmly establishes that OCS ice is formed together with CH 3 OH and OCN − , deep inside dense clouds or protostellar envelopes.
5. The median ice abundance towards the MYSO sample, as a percentage of H 2 O ice is CO:CH 3 OH:OCN − :OCS=24:20:1.53:0.15 6. It is deduced that CS is not the main precursor to OCS, because its abundance is too low by at least a factor of ∼40. Sulfurization of CO is likely needed. The source of this sulfur is unclear. Perhaps sulfur atoms or allotropes. The H 2 S ice abundance is insufficiently constrained.
7. The OCS abundances relative to H 2 O ice towards MYSOs and comets are remarkably similar. This contrasts with CH 3 OH, which shows a deficit for comets. Perhaps additional OCS is formed in protoplanetary disk environments.
8. Compared to low mass YSOs, MYSOs show an excess of OCN − ice in a significant fraction of the targets. The CO and CH 3 OH ice abundances are comparable to those observed towards low mass YSOs and dense cloud background stars. It should be noted that the low mass YSO and background star samples are fairly small, but this will be improved upon with JWST observations in the near future. Evidently, the current lack of detections of OCS ice towards low mass YSOs and background stars needs to be improved upon to form the trail of sulfur to comets. Indeed, the presented IRTF/SpeX 2-5 µm spectra of this sample of MYSO, most of which would saturate the JWST/NIRSpec and NIRCam WFSS detectors, will serve as a reference for the future JWST observations.