© The Authors 2024

1 Introduction

Complex organic molecules (COMs) are molecules with six atoms or more, including at least one atom of carbon (Herbst & van Dishoeck 2009). They are intrinsically important to comprehend the chemical complexity developed in star-forming regions since these materials are the feedstock for future exoplanetary systems. When made available in primitive planetary systems, this material can potentially promote the habitability of planets. An important question for delivering organic material to new solar systems is whether the molecules are in the gas phase or in ices as part of icy dust grains. It is only in the latter case that water and organic molecules would be expected to be effectively delivered to terrestrial planets, as discussed by Morbidelli et al. (2012, 2018); van Dishoeck et al. (2014, 2021); O’Brien et al. (2018), coupled with the lower UV photodestruction cross-section of those molecules in the solid-phase (e.g., Öberg 2016). Yet, very little information is available about these organic molecules in ices. The James Webb Space Telescope (JWST) provides a tool to amend this situation by observing ice features during the early protostellar phases with a higher resolution and sensitivity than before.

Gas-phase observations have been exceptionally successful in probing the chemical complexity toward low-mass and massive young stellar objects (LYSOs and MYSOs, respectively) as shown in the literature (e.g., Blake et al. 1987; Cazaux et al. 2003; Bergner et al. 2017; Manigand et al. 2020; Belloche et al. 2020; van Gelder et al. 2020; Jørgensen et al. 2020; Nazari et al. 2021; Gieser et al. 2021). In addition, the observations and abundances of many COMs are summarized in the review by Jørgensen et al. (2020). The current consensus is that COMs are efficiently formed in the solid phase and desorb into the gas phase by thermal and non-thermal mechanisms. Nevertheless, methanol (CH₃OH) is the only COM securely identified in the solid phase based on the infrared (IR) spectra from the United Kingdom Infrared Telescope (UKIRT; Grim et al. 1991; Skinner et al. 1992; Dartois et al. 1999), Infrared Space Observatory (ISO, e.g., Gibb et al. 2004), Very Large Telescope (Pontoppidan et al. 2004; Dartois et al. 2003; Thi et al. 2006), Infrared Telescopy Facility (Chu et al. 2020), Spitzer Space Telescope (Boogert et al. 2008; Bottinelli et al. 2010), AKARI (Shimonishi et al. 2010; Perotti et al. 2021), and, more recently, with JWST (e.g., Yang et al. 2022; McClure et al. 2023). Larger COMs compared to CH₃OH have been tentatively identified based on only a single infrared (IR) vibrational mode or proposed as possible carriers for yet unidentified features, such as ethanol, CH₃CH₂OH, and acetaldehyde, CH₃CHO (Schutte et al. 1999; Öberg et al. 2011; Terwisscha van Scheltinga et al. 2018). These two COMs were also tentatively identified in the JWST high spectral resolution and sensitivity spectrum of IRAS 15398–3359 (Yang et al. 2022) and toward two background stars (McClure et al. 2023), albeit at lower spectral resolution. Other upper limits have been determined for methyl formate (CH₃OCHO; Terwisscha van Scheltinga et al. 2021), methylamine (CH₃NH₂; Rachid et al. 2021), methyl cyanide (CH₃CN; Rachid et al. 2022), and for-mamide (NH₂CHO; Schutte et al. 1999; Slavicinska et al. 2023). Boudin et al. (1998) also estimated upper limits for molecules such as ethane (C₂H₆), acetylene (C₂H₂), hydrazine (N₂H₄), hydrogen peroxide (H₂O₂) and the hydrozonium ion $\left({\text{N}}_2{\text{H}}_5^{+}\right)$.

In parallel, the gas phase formation of COMs through ion-molecule reactions has been proposed as well (e.g., Balucani et al. 2015; Skouteris et al. 2018; Vazart et al. 2022). The recent detection in a protoplanetary disk of gas phase ${\text{CH}}_3^{+}$ (Berné et al. 2023), an important intermediate in such reactions, supports the idea that COM formation is not exclusively realized on icy grains. This study focuses on the latter.

In the mid-IR spectral range, icy COMs have multiple absorption features, at for example, ~5.8 µm, 6.8–7.0 µm, 7.2 µm, 7.4 µm, ~8 µm, and ~9.5–9.8 µm. However, due to the spectral overlap with the features of other molecules (including dust grain species), the bands at 7.2 and 7.4 µm bands have been considered the main fingerprint of COMs (e.g., Schutte et al. 1999; Öberg et al. 2011; Terwisscha van Scheltinga et al. 2018), which is also emphasized in this work. For example, Schutte et al. (1999) proposed that the 7.2 µm band can be due to HCOOH (formic acid) ice, notably, the O–H bending mode of the carboxylic functional group and the C–H bending mode. In particular, Bisschop et al. (2007) finds a good match between the 7.2 µm band of the high-mass protostar W33A and the tertiary mixture HCOOH:CH₃OH:H₂O. The assumption that HCOOH would be the only carrier of the 7.2 µm band would lead to an unrealistic high solid-state abundance which is not in line with gas-phase observations (Bisschop et al. 2007). Nevertheless, this is not an argument for excluding HCOOH as a component of the interstellar ice, because it is also visible at other wavelengths (e.g., 5.8 µm). Instead, this indicates that HCOOH is not the main carrier of the 7.2 µm band. As an alternative, Boudin et al. (1998) and Öberg et al. (2011) propose that the deformation mode in the methyl functional group (CH₃) of ethanol can be the carrier for the 7.2 µm band. In the case of the 7.4 µm absorption feature, Schutte et al. (1999) propose that the formate ion (HCOO⁻) and acetaldehyde can be the carriers of this band. In all cases, these tentative assignments need a convincing profile fit with laboratory data to confirm these chemical species as carriers of those bands. These fits must also consider a larger spectral range around the 7–8 µm range where other typical vibrational modes of COMs are also detectable.

In addition to COMs, ions such as OCN⁻ (cyanate ion), HCOO⁻ and ${\text{NH}}_4^{+}$ have been proposed to be present in interstellar ices mostly as part of salts produced by acid-base reactions (e.g., Geballe 1984; Grim & Greenberg 1987; Schutte et al. 1999; Schutte & Khanna 2003; Gálvez et al. 2010; Maté et al. 2012; Moreno et al. 2013; Bergner et al. 2016; Kruczkiewicz et al. 2021). The presence of OCN⁻ in ice mantles can be considered a secure detection based on comprehensive laboratory experiments and extensive analysis toward several protostars (e.g., Geballe 1984; van Broekhuizen et al. 2004) and toward background stars (McClure et al. 2023), and the good correlation between gas-phase abundances of HNCO and OCN⁻ ice (Öberg et al. 2009a). In the case of HCOO⁻ and ${\text{NH}}_4^{+}$, Boogert et al. (2015) argue that a convincing profile fit is still needed to firmly confirm the presence of these ions in ice mantles. If present in interstellar ices, these ions are likely part of refractory salts such as ammonium formate $\left({\text{NH}}_4^{+}{\text{HCOO}}^{-}\right)$, which have been studied in the laboratory and found in comet 67P/G-C (e.g., Poch et al. 2020). In addition, the presence of OCN⁻ and HNCO in ices has a strong astrobiological appeal: HNCO participates as a peptide bond between two single amino acids (e.g., Fedoseev et al. 2015; Quénard et al. 2018; Colzi et al. 2021; Ligterink et al. 2022).

In this work, we address the presence of simple molecules, COMs, and ions in protostellar ices using newly observed JWST spectra of two protostars with the Medium Resolution Spectro-graph (MRS) of the Mid-Infrared Instrument (MIRI) under the JOYS+ program¹ (JWST Observations of Young protoStars; van Dishoeck et al. 2023). The first source is a high-mass (~220 M_⊙; d = 4.9 kpc; Molinari et al. 1998, 2008) star-forming region called IRAS 23885+6053 (hereafter IRAS 23385). This cluster is highly embedded in its natal molecular cloud and shows maser emissions of H₂O (Casoli et al. 1986) and CH₃OH (Kurtz et al. 2004), characteristic of shocks. Spitzer observations of this source reveal extended emission of polycyclic aromatic hydrocarbons (PAHs) that are excited in this region by surrounding sources with spectral types between B1.5 and B5 (Molinari et al. 2008). NOEMA (NOrthern Extended Millimeter Array) observations of this source show a variety of gas-phase molecules, such as OCS, H₂CO, HNCO, CH₃OH and CH₃CN (Cesaroni et al. 2019). No other gas-phase COM have been detected in this region apart from CH₃OH and CH₃CN (Gieser et al. 2021). Recently, Beuther et al. (2023) present the rich MIRI-MRS spectrum of IRAS 23885, with focus on the outflow tracers H₂(0–0) S(7), [Fe II] (⁴F_9/2–⁶D_9/2) and [Ne II] (²P_1/2–²P_3/2), and an accretion tracer, the Humphreys α HI(7–6) emission line detected at a 3–4σ level. A multiwavelength analysis using MIRI and NOEMA data provides information about the hot, warm and cold molecular components in IRAS 23385 (Gieser et al. 2023). In addition, Francis et al. (2024) analyse the compact gaseous molecular emission in IRAS 23385.

The second source targeted in this work is the low-mass Class 0 protostar NGC 1333 IRAS 2A (hereafter IRAS 2A) (e.g., Jørgensen et al. 2005; Brinch et al. 2009), a well-studied hot-corino located in the Perseus complex, specifically in the NGC 1333 molecular cloud, at a distance of 299 ± 3 pc (Ortiz-León et al. 2018; Zucker et al. 2018). This is a protobinary system (Looney et al. 2000; Reipurth et al. 2002), hosting two collimated jets (Sandell et al. 1994; Tobin et al. 2015). IRAS 2A is also a source where several gas-phase COMs are identified (e.g., Bottinelli et al. 2007), including glycolaldehyde (HCOCH₂OH; Coutens et al. 2015; Taquet et al. 2015; De Simone et al. 2017), an important sugar molecule that participates in the formation of ribose, a component of the ribonucleic acid (RNA).

Compared to the Spitzer Space Telescope, which had enough sensitivity to observe the brighter low-mass protostars but low spectral resolution, the JWST (MIRI/MRS) observations are qualitatively superior. The MIRI-MRS resolving power (R) around the 7–8 µm region, where COMs signatures are present, is R = 3500–4000 (Labiano et al. 2021), whereas the Spitzer Infrared Spectrograph (IRS) offered a resolving power of only R = 60. Yang et al. (2022) demonstrate for IRAS 15398–3359 that the ice absorption features around the 7–8 µm are significantly stronger with much higher S/N profiles due to better resolution and sensitivity in the MIRI/MrS compared with Spitzer/IRS.

An important disclaimer to make is that although this paper presents JWST/MIRI spectra of a high- and low-mass proto-star, this is not a comparative work of ices in those objects. Instead, we show data from two sources that were observed first in our programs. A suitable comparison of ice features in different sources needs a more extended list of protostars and will be performed in future work. This paper is structured as follows: Sect. 2 introduces the JWST/MIRI observations of IRAS 23385 and IRAS 2A, the method for the data reduction, and the background subtraction. In Sect. 3 we describe how the ice bands are isolated, the fitting procedure, the method for quantifying the fitting degeneracy, and which criteria we use for a firm detection. We focus in this paper on the 6.8–8.6 µm range, and therefore molecules with vibrational modes outside this interval are not discussed (e.g., H₂O, CO₂, ${\text{NH}}_4^{+}$). The results are shown in Sect. 4, which includes the spectral fits, degeneracy analysis and ice column densities. The discussions of these results are presented in Sect. 5, outlining the important insights from the JWST observations, and how the ice COMs abundances correlate between the low- and high-mass protostar, as well as with the comet 67P/G-C. The conclusions are summarized in Sect. 6. For a selective reading of this paper the reader can use the following sections: (i) JWST data treatment and subtractions performed on the data – Sects. 2.1, 2.2, 3.1, 3.2, and 3.3. Appendices A and B; (ii) fitting methodology and results – Sects. 3.4 and 4.1. Appendices C and G; (iii) statistical analysis and degeneracy – Sects. 3.5 and 4.2. Appendices H, I, and J; (iv) identification criteria and suggestions for future works – Sects. 3.6 and 3.7. Appendices E and F; (v) ice column densities, abundances, and correlations – Sects. 4.3, 4.4, 4.5, 5.1.3, and 5.3. Appendices D and K; (vi) ice chemical complexity discussions – Sects. 4.1, 5.1 and 5.2. Appendix L.

2 Observations

2.1 Data reduction

IRAS 23385+6053 (RA 23^h40^m54.5^s, Dec +61^d10^m28^s) and NGC 1333 IRAS 2A (RA 03^h28^m55.57^s, Dec +31^d14^m36.97^s) were observed with JWST as part of the guaranteed observation time (GTO) 1290 (P.I. E. F. van Dishoeck) and 1236 (P.I. M. Ressler), respectively.

All data presented in this paper were taken with the Mid-InfraRed Instrument (MIRI; Rieke et al. 2015; Wright et al. 2015, 2023) in medium resolution mode (MRS; Wells et al. 2015; Labiano et al. 2021; Argyriou et al. 2023). For IRAS 23385, the observation was done in two-point dither mode in a four-point mosaic surrounding the central protostellar cluster. For IRAS2A, the observation consists of a single pointing on the protostellar binary in two-point dither mode. For both targets, also dedicated background observations were performed (no dither for IRAS 23385, two-point dither for IRAS 2A). For both targets, all three gratings (A, B, C) were used, providing the full wavelength coverage of MIRI (4.9–28 µm). All data were taken using the FASTR1 readout mode. The integration time in each grating was 200 s and 111 s for IRAS 23385 (per pointing in the mosaic) and IRAS 2A, respectively.

The data were processed through all three stages of the JWST calibration pipeline (Bushouse et al. 2022), using the reference context jwst_1017.pmap (IRAS 23385) and jwst_0994.pmap (IRAS2A) of the JWST Calibration Reference Data System (CRDS; Greenfield & Miller 2016). The raw uncal data were first processed through the Detector1Pipeline of the pipeline, followed by the Spec2Pipeline. In the latter step, the data are corrected for fringes using the fringe flat for extended sources (Mueller et al., in prep.), followed by applying a residual fringe correction (Kavanagh et al., in prep.). Moreover, in the case of IRAS2A, the telescope background was subtracted in this step using the dedicated background observation. However, the background observation of IRAS 23385 includes significant astronomical background emission across all wavelengths, as well as the strong emission of PAHs; thus, it could not be used to subtract the telescope background. The background estimation and subtraction procedure for IRAS 23385 is further discussed in Sect. 2.2. As mentioned in Beuther et al. (2023), an astromet-ric correction was applied for IRAS 23385, namely, 1.6077" in right ascension (RA) and 0.3485" in declination (Dec) based on identified Gαiα-DR3 stars in the parallel images. No such correction was necessary for IRAS2A. The data were further processed with the Spec3Pipeline of the pipeline, which produces cubes of all 12 subbands. In this step, both the master background and outlier rejection routines were switched off.

The observation of IRAS 23385+6053 reveals two mid-infrared continuum sources (Beuther et al. 2023) that are resolved at shorter wavelengths (channels 1 and 2, λ < 12 µm), but which start to become marginally resolved at longer wavelengths (channel 3, 12 µm < λ < 17 µm) and are completely unresolved at the longest wavelengths (channel 4, λ > l7 µm). The spectrum is therefore extracted from a large aperture of 2.5" in diameter (which does not increase with wavelength) centred in between the two sources (RA (J2000) 23^h40^m54.49^s, Dec (J2000) 61^d10^m27.40^s) to encompass the flux of both sources across the full wavelength range. The estimated 1σ rms increases from about 0.4 mJy below 20 µm to a few mJy at the longest wavelengths.

In the observation of IRAS2A, only continuum emission related to the primary component of the binary is detected. The spectrum is therefore extracted from the peak of the continuum emission at 5.5 µm (RA (J2000) 03^h28^m55.57^s, Dec (J2000) 31 ^d14^m36.76^s ). We assume that any contribution of the secondary component to the spectrum is negligible. The diameter of the aperture was set to 4 × l.22λ/D to capture as much of the source flux as possible without including too much noise. The estimated 1σ rms is about 0.4 mJy below 15 µm and increases to a few mJy at 19 µm and >10 mJy longwards of 22 µm.

Figure 1 (top) shows the MIRI/MRS spectra of IRAS 23385+6053 (Beuther et al. 2023) and Fig. 1 (bottom) shows IRAS 2A covering the range between 5 and 28 µm, and with a resolving power of 4000–1500 (Labiano et al. 2021). The spectral absorption features are associated with different ice molecules where the principal molecules are labelled in this figure, and those in the cyan region will be further discussed in Sect. 4. We highlight the broad feature of H₂O covering the range between 5.5 and 8 µm (bending mode) and between 10 and 20 µm (libration mode). HCOOH shows prominent features at 5.8 and 8.2 µm which can be distinguished in these sources. Small features attributed to NH₃ (ammonia), CH₄ (methane) and ${\text{NH}}_4^{+}$ (ammonium) are also seen in this spectrum. Silicates are the other main solid-state species contributing to the absorption bands around 9.8 and 18 µm. In addition to the absorption features, the spectra of IRAS 23385 and IRAS 2A have various narrow emission lines which are masked in this work since it is focused on ice absorption features. We also point out that the IRAS 23385 spectrum is binned by a factor of four between 8.6 and 12 µm because of the saturated silicate band, and IRAS 2A spectrum is binned by a factor of two in the entire MIRI-MRS range.

Fig. 1 MIRI MRS spectrum of the high-mass protostar IRAS 23385+6053 (top) and the low-mass protostar NGC 1333 IRAS 2A (bottom). Strong gas-phase emission lines are masked in both spectra. The assignments for the absorption bands are given and differentiated by the colours. The shaded cyan area indicates the “COMs region” that is studied in this work.

2.2 Background subtraction for IRAS 23385+6053

In the case of IRAS 23385+6053, the telescope and other backgrounds could not be subtracted in either the Spec2Pipeline or Spec3Pipeline since this results in negative fluxes due to significant astronomical emission in the dedicated background observation. It is, however, crucial to remove the telescope background from our observations to derive accurate ice column densities. The background was therefore estimated by extracting a spectrum from the science observation off-source from the main infrared continuum sources at the position within the IFU (Integral Field Unit) where the background flux was the lowest (RA (J2000) 23^h40^m54.15^s, Dec (J2000) 61^d10^m26.96^s) using the same aperture size as used for extracting the science data (2.5"). The background subtraction also results in the subtraction of the 8.6 µm and 11.3 µm PAH features, the emission of which was about equally strong in the background-position as at the source position. However, a possible under- or over-subtraction of PAH emission does not alter the results in this work since the PAH bands are broader than the ice bands at the wavelengths targeted in this paper. The background subtracted spectrum of IRAS 23385+6053 is shown in the top panel of Fig. 1. The numerous strong gas-phase emission lines have been masked in this version of the spectrum. The unmasked version is available in Beuther et al. (2023) and van Gelder et al. (2024) for IRAS 23385 and IRAS 2A, respectively.

3 Methodology

In this section, we provide information about the methods used to trace and subtract the dust continuum profile, the procedure to remove the silicate bands, and finally, the technique used to fit and identify the ice absorption bands. The focus of this paper is on the 6.8–8.6 µm range, as indicated in Fig. 1, but the entire MIRI spectrum has to be taken into account for the continuum fitting. For clarity, we summarize the methodology steps in a flowchart shown in Fig. 2.

Fig. 2 Summary of the methodology used in this paper.

3.1 Continuum subtraction and silicate removal

The spectral energy distributions (SEDs) of IRAS 23385+6053 and IRAS 2A show an increasing slope toward long wavelengths which is typical of embedded protostars. Such SED shapes are also observed toward many high-mass protostars (e.g., Orion BN, Orion IRc2; Gibb et al. 2004) and low-mass protostars (e.g., CrA IRAS32, IRAS 03301+311 and L1448 IRS1; Boogert et al. 2008). The shorter wavelength SED is composed of contributions by warm dust at a range of temperatures, whereas the increase in flux beyond 20 µm is due to the coldest envelope material (Adams et al. 1987). Determining the continuum of these protostars in the mid-IR range is not trivial because of the broad absorption bands in this spectral region. Often, a guided polynomial function is used (e.g., Boogert et al. 2008) to trace the continuum SED. In this work, a third-order polynomial function is used to fit selected points between 5.3–5.5 and 27–28.5 µm where there is little or no overlap exists with ice absorption bands. An additional point was added at 7.5 µm slightly above the observational data to avoid unrealistic inflexions in the low-order polynomial that would lead to an unphysical continuum. The reason for using this point is to account for known broad absorption features in this range, for example, the blue wing of the 9.8 µm silicate feature, the H₂O ice broad bending mode, part of the C5 component proposed in Boogert et al. (2008), and some of the complex molecules targeted in this paper. In this case, the observed flux itself at 7.5 µm is not suitable to be used as an anchor point. Uncertainties in the position of this guiding point do not affect the conclusions of this work. Figures 3a and b, show the polynomial fits used for IRAS 23385 and IRAS 2A, respectively, and the emission lines are masked. We note that there is significant absorption with respect to this continuum over the entire wavelength range in both sources.

Once the continuum SED is determined, we convert the MIRI-MRS spectra of the two protostars to an optical depth scale, as shown in Figs. 3c and d by using the equation below: $${\tau }_{\lambda }=-\ln \left(\frac{F_{\lambda }^{\text{source }}}{F_{\lambda }^{\text{cont }}}\right)\text{, }$$(1)

where $F_{\lambda }^{\text{source }}$ is the source spectrum and $F_{\lambda }^{\text{cont}}$ is the continuum SED.

Among the absorption bands seen in both spectra, silicates significantly contribute to the bands around 9.8 µm and 18 µm. Since the profile at 9.8 µm is positioned on top of a few ice bands, we perform a silicate removal to analyse the absorption features related to icy molecules. Often, the silicate profile observed toward the galactic centre source, GCS 3, is used to remove the silicate profile of other protostars (Boogert et al. 2008; Bottinelli et al. 2010). However, this method is not used in this work because the shape of the silicate features toward IRAS 23385+6053 and IRAS 2A are broader than GCS 3 as shown in Figs. A.1 and A.2, which would lead to a spurious absorption profile around 8.6 µm after removing the silicate. As an alternative method, we combined two types of silicate to match the band at 9.8 µm following the approach described in Boogert et al. (2011); Poteet et al. (2015); Do-Duy et al. (2020) and McClure et al. (2023). Specifically, we used a mixture of amorphous pyroxene (Mg_0.5Fe_0.5SiO₃) and olivine (MgFeSiO₄) from Dorschner et al. (1995).

We used the optool code (Dominik et al. 2021) to generate optical depth spectra of the two silicate types. Briefly, we assumed a mixture of silicate and carbon, typical chemical species of interstellar grains, with volume fractions of 87% and 13%, for IRAS 23385 and 82% and 18% for IRAS 2A, as typically used in protostar models. The variation in the fraction of carbon in the models allows us to fit the 9.8 µm and 18 µm bands simultaneously because a higher volume fraction of amorphous carbon reduces the intensities of the silicate bands. Different carbon fractions have been used in the literature to create dust models, for example, 30% (Weingartner & Draine 2001), and 15% (Pontoppidan et al. 2005; Woitke et al. 2016). For the dust models, we assumed a power-law size distribution with an exponent of −3.5 and grain sizes ranging from 0.1 to 1 µm. We adopted a distribution of hollow spheres (DHS; Min et al. 2005) to model the silicate band, as this approach mimics irregular geometries of the dust grains. The generated silicates are combined linearly to match the 18 µm band without exceeding the absorption at 9.8 µm and the blue silicate wing around 8.5 µm as shown in Fig. 3c and d.

The optical depth ratio between the 9.8 and 18 µm silicate features is equal to 2.7 for IRAS 23385 and 2.95 for IRAS 2A. These values are higher compared to the silicates found in the diffuse interstellar medium, which ranges from 1.4 to 2.0 (Chiar & Tielens 2006). It is likely that grain growth plays a role in this case, but a detailed study of this process is beyond the scope of this work. The silicate-removed spectra of IRAS 23385 and IRAS 2A are presented in Figs. 3e and f. By removing the silicate bands, the H₂O ice libration band is revealed around 12 µm. One can also see the ammonia umbrella mode at 9.0 µm and the methanol C-O stretch mode at 9.8 µm in IRAS 2A (Fig. 3f), but not in IRAS 23385 because of the saturation due to the silicate feature, as well as the higher noise level.

For completeness, we show in Appendix B, the silicate subtraction using different silicates for IRAS 2A, the source with high signal-to-noise ratio (S/N). Those spectra are scaled to an optical depth of τ = 4.9, which is the same as in the synthetic silicate fit (top panel of Fig. B.1). The subtracted spectra are shown in the bottom panel. It is clear that the synthetic silicate and the MgSiO₃ model taken from Poteet et al. (2015) show very similar profiles. On the other hand, the use of the GCS 3 silicate leads to an unrealistic absorption excess. Finally, it is worth mentioning that silica (SiO₂) has a blue shoulder at 8.6 µm and a peak at 9 µm. However, the presence of silica is associated with other materials such as enstatite and forsterite whose spectral features are not seen in the sources addressed in this paper. For the different silicates considered in this paper, there are no relevant differences in the spectral shapes between 7.8 and 8.5 µm.

Fig. 3 Continuum (a, b) and silicate subtraction (c, d) steps in IRAS 23385 (left) and IRAS 2A (right). The spectrum of IRAS 23385 is binned by a factor of two between 9 and 11 µm due to the saturated silicate profile. A third-order polynomial function is used to trace the continuum, and the silicate profile is a combination of two laboratory silicate spectra (olivine and pyroxene). Panels e and f show the silicate subtracted spectra of both protostars, with major features labelled.

3.2 SO₂ gas emission subtraction

In the spectrum of IRAS 2A, clear molecular emission is superimposed on the ice absorption features between ~7.25 µm and ~7.45 µm, see Fig. 4. This emission was recently found to originate from warm gas-phase SO₂ (van Gelder et al. 2024). In order to accurately determine the contribution of ices in this wavelength range, the gas-phase emission lines of SO₂ (ν₃) have to be subtracted from the spectrum. This was achieved by subtracting the best-fit gas-phase emission line model of SO₂, which was recently derived by van Gelder et al. (2024). This model is very well constrained by the R branch lines at 7.3 µm.

The SO₂ emission subtracted spectrum is also presented in Fig. 4 as the orange line and clearly reveals an ice absorption feature around 7.4 µm that was hidden by the SO₂ emission. The SO₂ emission is slightly over-subtracted around 7.35 µm (i.e., at the Q-branch) but this does not hamper the analysis of the ice bands since this residual is much narrower than typical ice absorption bands.

Fig. 4 Observed spectrum of IRAS 2A (blue) surrounding the 7.2 µm and 7.4 µm ice absorption features with SO2 (ν3) emission. The best-fit SO2 model of van Gelder et al. (2024) is presented in green and the SO2 subtracted spectrum of IRAS2A is shown in orange. Subtracting the SO2 emission clearly reveals the 7.4 µm absorption feature.

3.3 Isolating ice features between 6.85 and 8.6 µm

Weak absorption features have been measured in the laboratory covering the range between 6.85 and 8.6 µm (e.g., Lacy et al. 1991; Schutte et al. 1999; Öberg et al. 2011), and Terwisscha van Scheltinga et al. (2018), which are compared with observational spectra of massive protostars, by Boogert et al. (2008); Öberg et al. (2011), low-mass protostars (e.g., Zasowski et al. 2009; Yang et al. 2022) and background stars (McClure et al. 2023). Most notable are the absorption features around 7.2 and 7.4 µm. To isolate these bands, previous works used a polynomial fit on the flux scale data to trace a local continuum starting around 7.07.14 µm and finishing around 7.8–8.0 µm (e.g., Schutte et al. 1999). This approach isolates the 7.2 and 7.4 µm features but excludes potential absorption features around 7 µm and at wavelengths long-wards of 7.8 µm where C–H and C–O absorption features of many possible molecules contribute.

In this work, we isolate the 7.2 and 7.4 µm features using a third-order polynomial function, and following a slightly different approach. First, instead of using the spectrum on the flux scale, we perform the polynomial fit on the optical depth scale after removing the silicate absorption. The strong silicate band makes it difficult to observe small features such as those due to COMs. Second, we use guiding points fitted by a third-order polynomial function, as shown in the left panels of Fig. 5. This additional continuum represents blended absorption profiles from the broad H₂O ice bending mode, the red wing of the ${\text{NH}}_4^{+}$ cation, and the C5 component proposed by Boogert et al. (2008). Another small contribution from the O-H bending mode of CH₃OH ice is also considered in this step. The points used for IRAS 23385 are at 6.8, 7.2, 7.7 and 8.5 µm and at 6.8, 7.3, 7.5, 9.4 and 10 µm for IRAS 2A. The positions of the guiding points are distinct because of the differences in the absorption profiles of the two sources. In order to account for possible C–H absorption bands, the first point at 6.8 µm is chosen to be marginally above the wing of the strong 6.85 µm feature. The points long-wards of 7.8 µm are selected where we expect no or weak ice absorption, to account for the C–O features. In the case of IRAS 2A, we use points at 9.4 and 10.0 µm because of the clear absorption profiles at 9 and 9.8 µm.

After subtracting the local continuum fit, we isolate the ice features in the range between 6.8 and 8.5 µm as shown in the right panels of Fig. 5. In these figures, we can see absorption bands at 7, 7.2, 7.4, 7.5–7.8, and 8.2 µm that are fitted and discussed in Sects. 4 and 5, respectively. Notably, the 8.2 µm band of IRAS 23385 seems broad and asymmetric, whereas in IRAS 2A the same band seems broad, weaker and symmetric. The band around 7.7 µm has a prominent blue shoulder in both sources, reflecting the interaction with other neighbouring chemical species.

3.4 Fits and laboratory data

Since molecules have multiple functional groups, performing simultaneous fits at different wavelengths supports secure detection. Additionally, some molecules are expected to absorb in the wavelength range targeted in this paper, such as CH₄, SO₂, HCOOH, CH₃CHO and CH₃CH₂OH (e.g., Schutte et al. 1999; Bisschop et al. 2007; Öberg et al. 2008), which are used as starting point in the MIRI spectral ice decomposition. We fit the region between 6.8 and 8.6 µm of IRAS 2A and IRAS 23385 using the ENIIGMA fitting tool (Rocha et al. 2021). This code searches for the global minimum solution that fits the observations by performing a linear combination of laboratory ice data. The genetic algorithm approach benefits from a simple fitness function (Baeck et al. 2000), and here we use the root-mean-square error (RMSE) given by the equation below: $$\text{RMSE}=\sqrt{\frac{1}{n}{\displaystyle \sum _{i=0}^{n-1}{\left({\tau }_{v,i}^{\text{obs}}-{\displaystyle \sum _{j=0}^{m-1}{w}_j}{\tau }_{v,j}^{\text{lab}}\right)}^2}},$$(2)

where both experimental $\left({\tau }_{v,\text{j}}^{\text{lab}}\right)$ and observational $\left({\tau }_{v,\text{i}}^{\text{obs}}\right)$ spectrum are converted to wavenumber space (v), w_j is the scale factor, and m and n are the mth and nth data point. The absorbance laboratory data (Abs) are converted to an optical depth scale by the equation ${\tau }_v^{\text{lab}}=2.3Ab{s}_v$. In the degeneracy analysis shown in Sect. 3.5, the error of the data is taken into account.

In a nutshell, ENIIGMA uses genetic modelling algorithms for searching the optimal coefficients of the linear combination (w). Genetic algorithms are robust optimization techniques based on the processes of natural selection that aim to find the global minimum solution for complex problems (Holland 1975; Koza 1992). Once the best fit is found, ENIIGMA calculates the ice column density of each component using the following equation: $$N_{\text{ice }}=\frac{1}{A}{\displaystyle {\int }_{v_1}^{v_2}{\tau }_v^{\text{lab}}}dv,$$(3)

where A is the vibrational mode band strength of the molecule, which is listed in Table 1. The band strengths of molecules change depending on the chemical environment. For this reason, we adopt values of corrected band strengths, when available, to derive the ice column densities. The derivation of band strengths is not straightforward because it depends on the ice density. The typical band strength uncertainties are around 15% and 30% for pure and mixed ices, respectively (Rachid et al. 2022; Slavicinska et al. 2023).

The laboratory data considered in this paper are listed in Appendix C. These data were taken mainly from the Leiden Ice Database for Astrochemistry² (LIDA; Rocha et al. 2022) and from the Goddard NASA database³. The methodology used by ENIIGMA to test the data available is detailed in Rocha et al. (2021). Here, we provide a brief description of the method. In the first stage ENIIGMA combines IR spectra of pure ice at low temperature with pure ice at high temperature. The best group of solutions is passed to the second step. In the second stage, ENIIGMA combines the previous best solution with pure ice, with all ice mixtures in the ENIIGMA database in a sequential way. At this stage, all data that has an IR feature in the range fitted was tested. Finally, ENIIGMA passes the best-ranked groups of solutions to a final stage, where species from one group of solutions are mixed and combined with species from another group. This allows the code to diversify the number of combinations and increases the possibility for the code to find the global minimum solution not only among the coefficients but also among the laboratory data available. In total, the code tested 3173 different combinations, where the best solution is the best-ranked group of components based on the RMSE value.

We note that for the purposes of this paper, we subtracted the absorption profiles of H₂O and CH₃OH around 6.8 µm from the ice mixtures using a local subtraction with a polynomial function (see Appendix E). This enables fitting the observational data after local subtraction. As previously mentioned, the contribution of these two molecules is taken into account in the polynomial fit used to trace the local continuum between 6.8 and 8.6 µm. Most of the COM laboratory data have a spectral resolving power of R = 5000, which is degraded to the nominal spectral resolution of the two sources (R ~ 3500) around 7–8 µm. In the case of HCOO⁻ and OCN⁻, we isolated the ice bands of these ions using a local baseline subtraction. These two species are formed from molecules engaged in an acid-base reaction, also known as Bronsted-Lowry acid-base theory (Brönsted 1923; Lowry 1923), and also seen in interstellar ices (e.g., Grim & Greenberg 1987; Schutte & Khanna 2003; van Broekhuizen et al. 2004). The amount of ions formed depends on the initial abundances of the parent molecules. By isolating these bands, we can mimic (at first order) different initial conditions of parent species. It is also worth mentioning that the baselines of HCOOH ice mixture are checked before the analysis. Öberg et al. (2011) has noted that we should be careful when deriving the formic acid and formate ion column densities because of baseline issues in some experimental data measured further in the past. In addition, it is important to note that HCOO⁻ and HCOOH share a band at 7.2 µm. For that reason, we would recommend using the 7.4 µm to quantify HCOO⁻.

Fig. 5 Procedure for isolating weaker absorption features in the spectra of IRAS 23385 and IRAS 2A. Left: local continuum in the 6.8–8.6 spectral range using a fourth-order polynomial function. Right: isolated 6.8-8.6 optical depth spectrum in both protostars.

3.5 Degeneracy analysis

The ENIIGMA fitting tool performs a degeneracy analysis of the coefficients in the linear combination that results in the best fit. Briefly, the code performs a Gaussian variation around each coefficient by using the numpy function numpy.random.normal (Harris et al. 2020), and is given by: $$p(x)=\frac{1}{\sqrt{2\pi {\sigma }^2}}\exp \left[-\frac{{(x-\mu )}^2}{2{\sigma }^2}\right]$$(4)

where µ values are the optimal coefficients, and σ is the standard deviation around µ. This analysis allows us to calculate χ² values for each new linear combination and derive confidence intervals based on a Δχ² map (Avni & Bahcall 1980), which is formulated as: $${\chi }^2=\frac{1}{\text{ d}.\text{o}.\text{f}.}{\displaystyle \sum _{i=0}^{n-1}{\left(\frac{{\tau }_{v,i}^{\text{obs}}-{\displaystyle \sum _{j=0}^{m-1}{w}_j}{\tau }_{v,j}^{\text{lab}}}{{\gamma }_{v,i}^{\text{obs}}}\right)}^2},$$(5a) $$\text{Δ}{\chi }^2(\alpha ,ϵ)={\chi }^2-{\chi }_{\min }^2,$$(5b)

where d.o.f. is the number of degrees of freedom, γ is the error in the observational optical depth spectrum propagated from the flux error assumed to be 10%, α and e are the statistical significance and the number of free parameters, respectively. ${\chi }_{\min }^2$ corresponds to the goodness-of-fit in the global minimum solution.

ENIIGMA also quantifies the statistical significance of a given IR spectrum based on its recurrence, which is defined as: $$ℛ=\frac{f_i}{S},$$(6)

where f_i is the absolute frequency of sample i (i.e., how many times a specific laboratory data participates in the fit) and S is the total number of solutions. If ℛ = 100%, the chemical species cannot be excluded from the fit. On the other hand, lower percentages mean that the spectrum can be replaced by another one with a similar spectral shape without going outside of the confidence intervals. While the confidence interval analysis evaluates the degeneracy among the components in the best fit, the recurrence analysis quantifies the degeneracy among different spectral data.

3.6 Criteria for firm COM ice detections

Boogert et al. (2015) classify three types of detection in ices:

• 1.

  Secure: multiple bands or bands of isotopologues are detected in high-quality spectra;

• 2.

  Likely: A single band is detected and the profile matches the laboratory spectra;

• 3.

  Possibly: A single band is detected and there is no exact match between the profile and the laboratory spectra.

COMs bands are naturally weak, and therefore the isotopologues criteria in (1) will be hardly satisfied for solid-phase detections. In addition, to these three criteria, Jørgensen et al. (2020) discuss gauges for the detection of exotic chemical species in the context of gas-phase observations. They mention that a firm identification needs a complete spectral survey with a synthetic spectrum that accounts for all the bands of the identified molecules instead of using independent analytical functions (e.g., Gaussian fits) of individual lines. These criteria can also be applied in the context of COMs in ice, and it is strongly recommended to use IR laboratory data for comparison instead of analytical functions (e.g., Gaussian, Lorentzian).

We note that different from gas-phase observations, where the emission profiles are narrow and isolated, solid-phase absorption profiles are often blended because of the common functional groups and broader spectral features. In addition, because of the ice matrix in which the molecules are embedded, the shape of the absorption bands change, with mixing ratio and temperature (e.g., Öberg et al. 2007; Bouwman et al. 2007). While ice COMs show peculiar band shapes at high temperatures (>70 K), they are very similar at lower temperatures (see Terwisscha van Scheltinga et al. 2018, 2021; Rachid et al. 2020, 2021, 2022; Slavicinska et al. 2023). In this sense, we add to the criteria presented by Boogert et al. (2015) and Jørgensen et al. (2020) that a degeneracy analysis of molecules sharing similar functional groups is needed in order to claim a firm ice COM detection.

3.7 Remarks for future works

In this section, we point out some aspects important to guide future works on the analysis of COMs fingerprints in protostellar ices:

• The intensities of COM bands in observational data are weaker compared to the major ice components. In this case, special attention must be paid to the baseline correction of COM IR laboratory spectra. Any minor inflexion can mimic a spurious feature and lead to misinterpretation of the observational spectrum. An example is presented in Appendix F. The example is given for H₂O:CH₃CH₂OH ice spectrum where a 7th-order polynomial is used to baseline correct the spectrum. We show that if a lower number of data points is considered around 8 µm, the polynomial function can fluctuate and originate spurious features in the final data. Instead cubic spline functions can be used to mitigate those fluctuations. Another issue, not shown here, is that some weak bands can be removed unintentionally if the user is not familiar with a particular dataset and source details.

• Tracing the local continuum on the observed spectrum is as critical as the spectral fitting. In the COMs fingerprint region, where broad bands are present, one should be careful when attributing zero absorption for tracing the local continuum between 6.8–8.6 µm.

• Analysis of observational spectra in the range between 6.8 and 8.6 µm using laboratory data of H₂O- or CH₃OH-containing ices must have the features of these two molecules subtracted to allow for a direct comparison with observational data after local continuum removal. Appendix E shows an example of how the CH₃CH₂OH bands in the mixture with H₂O and CH₃OH were isolated. While a single fourth-order polynomial was used in the case of H₂O:CH₃CH₂OH, three polynomial functions with different orders were needed in the case of CH₃OH:CH₃CH₂OH, which increases the risk of creating spurious features.

4 Results

In this section, we show the fitting results of the protostars IRAS 23385+6053 and IRAS 2A in the range between 6.8 and 8.4 µm, as well as the confidence interval analysis.

4.1 Spectral decomposition and feature analysis

The fits of the IRAS23385 and IRAS 2A spectra are shown in Fig. 6 top and bottom, respectively (see Appendix G for an incremental version of these figures, following the individual fitting steps). These two MIRI/MRS spectra are decomposed using nine laboratory spectra which were selected by ENIIGMA and they provide the global minimum solution. Among the COMs fitted in this work are, CH₃CHO, CH₃CH₂OH, CH₃OCHO (methyl formate) and CH₃COOH (acetic acid), with the first three robustly detected (see Sects 3.6 and 4.2). The simple molecules identified are CH₄, SO₂, and HCOOH. Additionally, we found a good match of two ions, HCOO⁻ and OCN⁻. Because of the contribution of HCOO⁻ with similar intensities at 7.2 and 7.4 µm, we stress that these two bands are not only associated with COMs.

The band around 7.68 µm was clearly visible in previous ice observations (e.g., Gibb et al. 2004; Öberg et al. 2008), and is now seen in JWST spectra with unprecedented S/N and spectral resolution. CH₄ ice is the main carrier of this feature, and Öberg et al. (2008) suggests that SO₂ can contribute to the blue wing of this band based on a Gaussian decomposition of the 7.68 µm feature. In the present work, this band is decomposed into four components, with CH₄ mixed with H₂O the dominant carrier. This is in line with laboratory experiments suggesting a common formation pathway for CH₄ and H₂O ices (Qasim et al. 2020). The blue wing has contributions of SO₂ mixed with CH₃OH, and the negative cyanate ion (OCN⁻) for both sources. The red wing of this band can be fitted with CH₃COOH in the case of IRAS 2A but is not present in IRAS 23385. A good correlation in ice column density between sulfur-bearing molecules and methanol was observed before by Boogert et al. (1997, 2022). Here, OCN⁻ is one of the ions formed from molecules engaged in acid-base reactions, and it has a band at 7.63 µm, in addition to its well-known feature at 4.61 µm, commonly seen toward protostars. In fact, recent observations of IRAS 2A with the JWST/Near-Infrared Spectrometer (NIRSpec) detect the 4.61 µm band, thus confirming the presence of OCN⁻ in the ices toward IRAS 2A as part of the JOYS program. The analysis of the NIRSpec data of this source will be presented in a future paper, but we calculate the ice column density of OCN⁻ at 4.61 µm to check the consistency of the fit in the mid-IR (see Sect. 4.3). As a result, the blue wing of the 7.68 µm feature is composed both by SO₂ and OCN⁻. We also note that OCN⁻ in both sources has a similar intensity, whereas SO₂ is stronger in IRAS 2A.

Protostar observations in the mid-IR show the 7.2 and 7.4 µm features with similar strengths across a number of sources, which suggests the dominance of a single species. In fact, we find that the formate ion, initially proposed by Schutte et al. (1999), matches well the 7.2 and 7.4 µm bands of IRAS 2A and IRAS 23385. This ion is formed by the acid-base reaction of H₂O:NH₃:HCOOH (100:2.6:2) ice mixture at 14 K (Gálvez et al. 2010). The intensity of the formate ion in IRAS 2A is a factor of three stronger than in IRAS 23385. In Appendix H, we compare the 7.2 and 7.4 µm band with the formate ion at other temperatures. The formate ion spectrum at 150 K has a broader profile at 7.4 µm, whereas laboratory data at 210 K shows a strong feature at 7.3 µm that is not observed in either source. Thus, the 7.2 and 7.4 µm features are a signature of ices dominated by ions in cold regions. Our results also show that despite HCOO⁻ being the main carrier of the 7.2 and 7.4 µm bands, other components can contribute to these two features separately as discussed below.

The presence of HCOO⁻ in the ice is supported by the detection of HCOOH at 8.2 µm, in addition to its band at 5.8 µm. The best fit is found when formic acid is mixed with CH₃OH and H₂O as suggested by Bisschop et al. (2007). The lower intensity of formic acid in IRAS 23385 can be related to the nature of the source. At warmer temperatures (>50 K), HCOOH is more efficiently destroyed via an acid-base reaction (e.g., Schutte et al. 1999; Gálvez et al. 2010). Another possibility is that HCOOH was not efficiently formed in this high-mass source because of the lower amount of CO ice available in a high-mass star-forming region.

CH₃CH₂OH is detected in both protostars through the absorption of four vibrational modes at 6.8–7.15 µm (C–H stretch), 7.25 µm (CH₃ s-deformation), 7.4–7.7 µm (OH deformation) and at 7.85 µm (CH₂ torsion). Among the ethanol data tested, the mixture with H₂O provides the best fit. Ethanol is stronger and contributes significantly to the 7.2 µm band in IRAS 23385, whereas it is less prominent in IRAS 2A.

CH₃CHO shows weak absorption around 7.03 µm (C–H stretch), and at 7.41 µm (CH₃ s-deformation/CH wagging). These bands are stronger in IRAS 2A than in IRAS 23385. Experimental characterization of the CH₃ s-deformation/CH wagging shows a strong dependence on the chemical environment and temperature (Terwisscha van Scheltinga et al. 2018). The closer match with both protostars is for acetaldehyde mixed with H₂O ice, which has a peak at 7.41 µm. The observed band, however, has a peak at 7.38 µm, indicating other carriers for this band (e.g., HCOO⁻).

CH₃OCHO is observed in both sources around 8.2 and 8.6 µm, but with different spectral shapes and intensities. In IRAS 23385, the methyl formate band is better fitted by a mixture with methanol at 15 K, although the data at 30 and 80 K are also statistically likely (see Sect. 4.2.3). Other methyl formate mixtures do not exhibit such an asymmetric profile, which makes them unlikely to be the carrier of this band (see Sect. 5.1.1). While the shape of the 8.2 µm is well defined, this is not the case for the 8.6 µm band, which is affected by saturation due to the strong silicate feature. In IRAS 2A the bands at 8.2 and 8.6 µm are consistent with a mixture containing CO, H₂CO, and CH₃ OH. This ice mixture of CH₃OCHO also fits well the 8.0 µm band associated with H₂CO.

CH₃COOH is found to fit well to the IRAS 2A spectrum, but not that of IRAS 23385. In addition, a better match with the observations is found for acetic acid mixed with H₂O ice. The two identified bands are located at 7.3 and 7.7 µm (salmon colour in Fig. 6, bottom). Specifically, the 7.7 µm band in IRAS 2A has a broad profile that is more than just CH₄ and it is necessary to fit the absorption excess around this region without any strong overlap between 7.7 and 7.85 µm. The 7.3 µm is slightly over-estimated in the fits, because of the amount of absorption needed to fit the 7.7 µm band. This can be because of the SO₂ emission lines subtraction around 7.3 µm or the uncertainties involved in the local continuum choice around the 7.7 µm band. The robustness of these detections based on the difference in local continuum choice is discussed in Sect. 4.2.4 for IRAS 2A, the higher S/N source.

Fig. 6 ENIIGMA fits of IRAS 23385 (top) and IRAS 2A (bottom). Gas-phase lines are masked. The figure labels show the ice mixture used in the fits, and a simplified version with the names of the chemical species names is shown close to the bands. The corresponding temperature of these laboratory spectra ranges between 10 K and 15 K.

4.2 Statistical analysis

4.2.1 Confidence intervals

We derive the confidence intervals (see Sect. 3.5) of the fits for IRAS 23385 and IRAS 2A in three different ranges separately: (i) 6.85–7.5 µm, (ii) 7.5–7.8 µm and (iii) 7.8–8.6 µm.

The components in these three ranges are relatively isolated, and therefore their contribution in one given interval is kept constant when analysing other ranges.

In Fig. 7, we show the confidence intervals for the fit of IRAS 2A in the range between 7.5 and 7.8 µm. The yellow and red contours indicate 2σ and 3σ confidence intervals. Based on these contours, we may note that all the components, apart from CH₃CH₂OH, are required to fit IRAS 2A. In this particular spectral range (7.5 and 7.8 µm), only a small portion of CH₃CH₂OH spectrum at 7.5 µm contributes to the absorption. Therefore, the confidence interval analysis shows that CH₃CH₂OH is not crucial to fit the 7.5 and 7.8 µm. The contribution of CH₃CH₂OH in IRAS 2A is better evaluated using the range between 6.86 and 7.5 µm, which is shown in Appendix I. In Fig. I.1 (top), the CH₃CH₂OH:H₂O is not zero, which reinforces the idea that CH₃CH₂OH is robustly found in IRAS 2A. Another important result from Fig. 7 is that SO₂ (w1) and OCN⁻ (w3) are both required to fit IRAS 2A since their coefficients cannot be zero. Finally, CH₃COOH:H₂O cannot be excluded as a solution based on this statistical analysis. Further analysis of the 6.86– 7.5 µm range allows us to conclude that CH₃CHO, HCOO⁻, and HCOOH are also robust detections. In addition to these chemical species, Fig. I.1 (bottom) also shows that CH₃COOH cannot be excluded from the fit obtained with one specific continuum subtraction. More details for CH₃COOH are presented in Sect. 4.2.4. The analysis of the spectral range between 7.8 and 8.6 µm, shows that another COM, CH₃OCHO is also a robust detection in IRAS 2A spectrum.

A similar analysis is performed for IRAS 23385 (Figs. I.2 and I.3). For the range between 7.5 and 7.8 µm, CH₄ and OCN⁻ are essential, whereas SO₂ can be statistically not required due to the low S/N in IRAS 23385 spectrum. Likewise, the analysis of the 6.86–7.5 µm range shows that CH₃CH₂OH, CH₃CHO, and HCOOH are robust detections, whereas the formate ion (HCOO⁻) absorption may be explained by CH₃CH₂OH or HCOOH ice features. For the 7.8-8.6 µm interval, both HCOOH and methyl formate (CH₃OCHO) are robust detections.

Fig. 7 Corner plots showing the IRAS 2A coefficient confidence intervals for the range between 7.5 and 7.8 µm. The grey scale colour is the Δχ2 map derived from a total of 5000 values. Yellow and red contours represent 2 and 3σ significance, respectively.

Fig. 8 Bar plot showing the recurrence of the difference chemical species to the fit in IRAS 2A (top) and IRAS 23385 (bottom). Recurrences above 50% are considered essential to the fit and robust detection. CH3COOH in IRAS 23385 does not occur in any of the solutions. In the case of IRAS 2A (top panel) we present hatched bars that show the recurrence of the solutions if the errors in optical depth are increased by a factor of three.

4.2.2 Recurrence of the ice components

A complementary statistical analysis is performed on the recurrence of all solutions within a given confidence interval. The difference in this method is that we do not vary the coefficient values of each solution, but the laboratory data instead. Figure 8 (top) shows the recurrence of the IRAS 2A fit of 15 chemical species inside 3σ confidence interval and considering several solutions. This analysis indicates that the COMs providing the best fits (see Fig. 6), are the most recurrent in the bar chart (87.5% ≤ R ≤ 100%), and therefore cannot be excluded as a solution. Other chemical species, such as CH₃NH₂, CH₃OCH₃, NH₂CHO, HCOCH₂OH and CH₃COCH₃, some of which have been suggested previously to contribute in this range, have a recurrence lower than 50%. This is not sufficient to claim a firm detection and at best upper limits can be derived.

The reason that formamide and acetone are not part of the global fit is that the formate ion band shape dominates the absorption profile at 7.2 and 7.4 µm. In addition, acetic acid and ethanol also contribute to these two absorption profiles. In the cases of methylamine and dimethyl ether, the fits indicate that methyl formate contributes more to the 8.5–8.6 µm range.

The same analysis is done for IRAS 23385, shown in Fig. 8 (bottom). Since the MIRI data for this source have a lower S/N compared to IRAS 2A, there exist slightly more variations among the recurrence values. The components found in the global fit have a recurrence above 85%, whereas other COMs are recurrent by ~40% or less. Interestingly, CH₃COOH does not participate in any of the solutions tested because OCN⁻ and CH₃CH₂OH account for all the absorption around 7.7 µm.

4.2.3 Temperature degeneracy

Despite the well-known changes in the band profiles of astro-physical ices with temperature, there are some absorption features that barely vary with temperature. Consequently, the fitting routine is not able to distinguish between these data. Rocha et al. (2021) show that when data that is known to provide a good fit is arbitrarily removed from the database, the ENIIGMA fitting tool uses another data of the same species, but with similar temperature when possible. This slightly increases the fitness function value, but the fit is still good within the confidence intervals. In this section, we perform a statistical analysis to evaluate which temperature ranges are degenerate and provide a good fit, from those temperature ranges that can be excluded as a solution.

Figure 9 shows the recurrence plots for CH₃CH₂OH, CH₃CHO and CH₃OCHO in IRAS 2A and CH₃OCHO in IRAS 23385 for different temperatures settings available. The analysis of CH₃CH₂OH and CH₃CHO is not performed for IRAS 23385 because of the low S/N in the spectral range considered for the fits. To calculate the recurrence plot for these COMs, we selected the solutions that are ranked inside a 3σ confidence interval. For example, in the case of CH₃CH₂OH, we found 352 solutions where ethanol is present. The ethanol mixture at 15 K is present in all of these solutions, and therefore it has a recurrence of 100%. There are solutions that combine the low temperature (15%) ethanol mixture with other temperatures (e.g., 30, 70, 100 K). Because of the similarity of the ethanol bands in the fitted range, the recurrence of these additional ethanol data is slightly reduced, but still high. For the data at 150 and 160 K, there are around 70 solutions using one of these data, which gives a recurrence of ~20%. Another way to perform the same analysis is by using the best global solution as the initial guess and running new fits by replacing the IR spectra of specific COMs at different temperatures. This forces ENIIGMA to use only one CH₃CH₂OH data at a time, and prevents overlaps with data that feature similar spectral shapes. For ethanol, we obtain a total of nine solutions if the presence of the other components is fixed. In this case, the temperature range between 15 and 100 K has a recurrence of 100%, and the higher temperatures, have a recurrence of 0%. In the case of ethanol, IR spectra with temperatures below 100 K are degenerate, and all CH₃CH₂OH:H₂O mixtures fit IRAS 2A in that temperature range. Above 120 K, the ethanol ice mixture no longer fits the IRAS 2A spectrum well because of band broadening and shift.

Repeating the same procedure for the other species, we notice a similar behaviour. The fits with acetaldehyde in IRAS 2A are degenerate below 30 K. Above this temperature, CH₃CHO:H₂O ice mixtures do not offer likely solutions mostly because of band shift. The fits with methyl formate are degenerated below 50 K in IRAS 2A. At higher temperatures, some substructures arise in the CH₃OCHO spectrum that deviates from the observational data. Similarly, in IRAS 23385, methyl formate is degenerated below 50 K, whereas solutions at 80 K and above are less recurrent because of changes in the band profile. In conclusion, it is likely that most of the ices toward IRAS 2A and IRAS 23385 are located in regions with temperatures below 50 K.

Fig. 9 Bar plot showing the recurrence of the same ice mixture, but at different temperatures. Plots inside the blue and green polygons refer to IRAS 2A and IRAS 23385, respectively.

4.2.4 Robustness of detection based on different local continuum choices

Other local continuum fits for IRAS 2A (higher S/N) in the 6.8– 8.6 µm region are also investigated and presented in Appendix J. The first three panels of Fig. J.1 display different continuum profiles, where the top panel is the version adopted for the analysis in this paper that traced a third-order polynomial to the guiding points. The second panel displays the fourth-order polynomial, where the red dot is added to the guiding points. In this case, the continuum is slightly elevated at shorter wavelengths to accommodate the fit to the extra point at 8.5 µm. The third panel presents the continuum when two extra points are added (λ = 7.8 µm and 8.5 µm), and a sixth-order polynomial is used. All subtracted spectra using these three approaches are shown in the bottom panel of Fig. J.1. The major difference is seen in the last case (orange continuum), which completely removes any absorption excess at 7.8 µm, thus excluding any band at this wavelength.

Given the variability of the optical depth spectra of IRAS 2A between 6.8–8.6 µm with the choice of the local continuum, it is worth accessing the robustness of the detections reported in the previous sections considering other continuum profiles. This is shown in Fig. J.2. Figure J.2 top presents the new fit assuming the red continuum profile (inset panel) from Fig. J.1. This spectral fit remains good and all components found in the best fit are still present. The only issue is found at 8.5 µm (added guiding point) where the CH₃OCHO band is slightly over-predicted. This situation changes when the orange continuum is used to isolate the absorption features in IRAS 2A. The bottom panel of Fig. J.2 depicts a fit where CH₃COOH is no longer needed because of the guiding point at 7.8 µm to trace the local continuum. In this version of the fit, the HCOO⁻ and OCN⁻ components become much stronger than in the other cases, and the COMs bands are reduced by a factor between 1.5–2.0. A clear mismatch is seen around 7 µm because of the reduction of the COM bands. Also, there is more excess around 7.3 µm that is not fitted with the other COMs tentatively detected in this work (see Sect. 4.4).

This analysis reinforces that the detections of CH₃CHO, CH₃CH₂OH, and CH₃OCHO, as well as the ions and the simple molecules, are robust and do not depend on the local continuum choice. Clearly, the only exception is for CH₃COOH. However, a valid consideration is what other chemical species would create an absorption profile similar to the orange continuum. There is no trivial answer to this question. A tentative explanation comes from the five components decomposition from Boogert et al. (2008), in which only the C5 broad feature contributes to the 7.8 µm. The nature of the C5 component is not well understood and can be related to a flat profile of high-temperature H₂O ice bending mode, the overlap of other negative ions $({\text{HCO}}_3^{-},{\text{NO}}_3^{-},{\text{NO}}_2^{-})$ or organic refractory residue produced by energetic processing. The flatter bending mode would have little absorption at 7.8 µm and therefore is less likely, but the other two options are strong candidates.

4.3 Ice column densities and abundances

The column densities of the molecules fitting the 6.8–8.6 µm range are presented in Table 2. To derive their abundances with respect to H₂O ice, we use the libration band around 12 µm to calculate the water ice column density as shown in Appendix K. These abundances are compared to literature values for low- and high-mass protostars. In Sect. 5.3 we also show abundances with respect to CH₃OH ice, in which the method to derive the column density is demonstrated in Appendix K.

The ice abundances of CH₄, SO₂, OCN⁻, HCOO⁻ and HCOOH are within or close to the range expected for LYSOs and MYSOs. From the analysis in this paper, CH₄, SO₂, OCN⁻ compose the band around 7.67 µm, whereas previous works assigned this band to only CH₄, with a possible contribution of SO₂ (Öberg et al. 2008). Additionally, instead of using Gaussian profiles to derive the column densities, this paper uses laboratory data of CH₄:H₂O and SO₂:CH₃OH. Gaussian profiles provide first-order approach analysis of the profile of the ice bands, but they can under or overestimate the FWHM of real CH₄ and SO₂ bands.

In the case of OCN⁻, we derived ice column densities and abundances from the band at 7.62 µm for IRAS 2A. For consistency, we also calculated the ice column from the NIRSpec feature at 4.61 µm, which is 2.5 × 10¹⁷ cm⁻². Both values are consistent within the errors, and they result in OCN⁻ abundances in agreement with the range estimated in the literature for low-mass protostars (van Broekhuizen et al. 2005). For IRAS 23385, the OCN⁻ abundance calculated from the MIRI data is within the range found in MYSOs taken from Boogert et al. (2022).

Finally, the formic acid (HCOOH) and the formate ion (HCOO⁻) abundances are closely aligned with the literature values for both LYSOs and MYSOs. In the case of the formate ion in LYSOs, there is only one estimate in the literature for the low-mass protostar HH 46 IRS, which is similar to the abundance calculated for IRAS 2A. The ratios between HCOO⁻ and HCOOH are 0.46 and 0.16 for IRAS 2A and IRAS23385, respectively.

For the COMs, we derive ice column densities and abundances for CH₃CHO, CH₃CH₂OH, CH₃OCHO and CH₃COOH. In terms of ice abundance with respect to H₂O ice, IRAS 2A is more abundant in CH₃CHO by a factor of 1.75 compared to IRAS 23385. CH₃COOH is not detected in IRAS 23385, but has an estimated abundance of 0.3% in IRAS 2A. CH₃CH₂OH and CH₃OCHO are more abundant in IRAS 23385 by a factor of 1.5 and 7, respectively, compared to IRAS 2A. We note that the CH₃CH₂OH abundance for IRAS 23385 is in line with the upper limit derived from the high-mass protostar W33A (Terwisscha van Scheltinga et al. 2018). The other ice abundances are consistent with the upper limits derived in the literature, illustrating that JWST can now probe deeper than previous instruments.

With respect to CH₃OH ice, the abundances are between 9.3–14.6% for CH₃CHO, 16.1–24.6% for CH₃CH₂OH, 0.9–1.3% for CH₃OCHO and 3.9–6.0% for CH₃COOH, depending on the lower and higher CH₃OH listed in Table 2. The values are much lower than the upper limits derived in the literature (Terwisscha van Scheltinga et al. 2018) for CH₃CHO (52%) and for CH₃CH₂OH (42%).

4.4 Tentative detections and upper limits on column densities

A few COMs were not part of the global solution shown in Fig. 6. In this case, we perform a separate comparison of these data to check for tentative detections and derive upper limit column densities in IRAS 2A. This procedure is not applied to IRAS 23385 due to the low signal-to-noise ratio. For this step, we scale laboratory spectra of COMs at specific wavelengths to the MIRI spectrum. This allows us to take into account both the intensity and width of the IR ice band. We performed separate scaling for COMs that exhibit overlaps in terms of the IR features. The molecules used in this step are CH₃OCH₃, CH₃COCH₃, HCOCH₂OH, NH₂CHO, CH₃NH₂, and CH₃CN mixed with H₂O ice. Figure 10 shows the COMs spectra superposed to the IRAS 2A MIRI data, with upper limit column densities also listed in Table 2. To scale CH₃OCH₃, we use the CH₃-rock mode at 8.6 µm as a reference. This band has also contributions of ammonia and methyl formate. CH₃COCH₃ and CH₃CN have a CH₃ symmetric deformation mode at ~7.3 µm, and HCOCH₂OH a CH₂ deformation mode at the same position. Both CH₃COCH₃ and CH₃OCH₃ may contribute to the band at 8.1 µm due to the CCC asymmetric stretch. In the case of CH₃NH₂, we use the CH₃-rock mode at 8.5 µm as a reference, whereas the C–H bend at 7.2 µm is considered for NH₂CHO.

Finally, we also check the potential presence of solid-phase glycine, the simplest amino acid, in the spectrum of IRAS 2A. In Ioppolo et al. (2021), the solid-state formation of this simplest amino acid was proven. Glycine has many absorptions features in the IR, and the ωCH₂ mode at 7.46 µm is the strongest band of β-glycine. Depending on the glycine structure (e.g., neutral -NH₂CH₂COOH or zwitterion – ${\text{NH}}_3^{+}{\text{CH}}_2{\text{COO}}^{-})$, the position of this band may shift (Potapov et al. 2022) or be narrowed if the molecule is mixed in an argon sample (Ehrenfreund et al. 2001). In this paper, we use the zwitterionic form (${\text{NH}}_3^{+}{\text{CH}}_2{\text{COO}}^{-}$; Pilling et al. 2011), which is available in the UNIVAP database⁴. In Fig. 11, we compare the local subtracted spectrum around the 7.46 µm band with the β-glycine data at three different column densities. We use the band strength from Holtom et al. (2005) of the 7.46 µm band, calculated as 1.16 × 10⁻¹⁷ cm molecule⁻¹. This range of column densities is compatible with Gibb et al. (2004) who estimated an upper limit glycine column density of 3 × 10¹⁶ cm⁻², based on the absorption feature around 5.8 µm band. With respect to H₂O ice, we find an upper limit of between 0.03–0.2% (see sect. 4.3 and Table 2). This is close to the upper limit derived for W33A (<0.3; Gibb et al. 2004), and in agreement with the theoretical models (0.03–0.7%; Ioppolo et al. 2021).

Fig. 10 IR spectra of COMs not contributing to the global fit, but with absorption features in the range between 6.8 and 8.6 µm. These spectra are scaled to the IRAS 2A spectrum (black) to derive upper limit column densities. The vertical dotted lines are colour-coded and indicate the features used to derive the ice column densities. A small horizontal shift in the vertical lines is used to distinguish the bands at 7.3 µm.

4.5 Testing more complex alcohols and other molecules

In this section, we compare the IRAS 2A spectrum with two other alcohols more complex than CH₃CH₂OH, for instance, propanal (CH(O)CH₂CH₃) and 1-propanol (CH₃CH₂CH₂OH). The goal is to check for similarities and differences between the functional groups of these alcohols with ethanol since these more complex alcohols are expected to have features at similar locations. It is worth mentioning that propanal and 1-propanol have been synthesized in experiments with ice analogues (Qasim et al. 2019) using atom addition reactions, and in the case of propanol, via CH₃OH ice UV irradiation (Tenelanda-Osorio et al. 2022). Both propanal (starless core TMC-1; Agúndez et al. 2023) and propanol (Galactic Centre; Belloche et al. 2022; Jiménez-Serra et al. 2022) have been securely detected in the gas phase. Figure 12 shows that these three alcohols have absorption features around 6.8, 7.2 and 7.5 µm. However, the relative intensities among these bands for these three molecules are different. For example, ethanol has similar intensities at these three bands, whereas 1-propanol and propanal have different intensities by factors of 3–4. In the case of 1-propanol, the fit of the 7.2 µm band would require twice more absorption at 6.8 µm than observed in IRAS 2A. On the other hand, propanal could contribute to the blue wing of the 7.2 µm band. We highlight, however, that in addition to HCOO⁻, NH₂CHO can also contribute to the same absorption feature (see Sect. 4.4 and Fig. 10). Another caveat is that both 1-propanol and propanal IR spectra correspond to pure molecules, and therefore spectral differences because of the chemical environment are not perceived in this analysis.

Other than alcohols, hydrocarbons may also contribute to the 7.2 and 7.4 µm bands. In Appendix L, we show a comparison of IRAS 2A spectrum with pure C₂H₂, C₂H₄ and C₂H₆ around the 7 µm bands and beyond 11 µm. As a result, it can be seen that C₂H₂ and C₂H₆ could contribute to the red wing of the 7.2 µm band. However, their absorption bands longwards of 11 µm exceed the absorption profile in IRAS 2A. Based on these comparisons, the statistical analysis, and different choices for the silicate (Fig. B.1) and local continuum (Fig. J.1) we can conclude that CH₃CH₂OH is the alcohol that contributes most to the 7.2 µm band in addition to HCOO⁻. The systematic analysis of additional high S/N MIRI data is necessary to obtain robust constraints of larger COMs.

Fig. 11 Comparison of the ωCH mode of glycine IR spectrum (red) with local subtracted spectrum of IRAS 2A (black). The intensities of the glycine band are given for three ice column densities.

5 Discussion

In this section, we discuss the implications of our results from the analysis of the 6.8–8.6 µm region in a low- and high-mass source. This discussion is focused on the presence of simple and neutral species (CH₄, SO₂, HCOOH, H₂CO), simple and ionic species (OCN⁻, HCOO⁻), and complex organic molecules (CH₃CHO, CH₃CH₂OH, CH₃OCHO, CH₃COOH). This section ends with a direct comparison between the ice abundances in the low-mass protostar, IRAS 2A, and the ice bulk abundances in comet 67P/G-C. The goal of this comparison is to evaluate the COM inheritance scenario in solar-type protostellar environments. Because of this, we do not perform the same comparison with IRAS 23385.

Fig. 12 Experimental IR spectra of ethanol, propanal and propanol compared to IRAS 2A spectrum in the range between 6.8–8.6 µm. The grey area is not considered in the fits, but it is shown here to highlight the C–H stretching mode of these molecules. For better readability of this figure, the gas-phase emission lines between 6.8–7.2 µm are masked.

5.1 Chemical complexity of protostellar ices

In this work, we expand the COMs ice inventory by reporting the detection of at least two vibrational modes of COMs frozen in protostellar ices. We consider secure detections in the cases of CH₃CHO, CH₃CH₂OH, and CH₃OCHO. CH₃COOH also has two vibrational modes and a high recurrence in the degeneracy analysis of IRAS 2A. However, its presence depends on silicate subtraction and local continuum determination and therefore it is considered a tentative detection. This series of molecules are chemically related and strongly supports laboratory experiments and computational simulations that suggest COMs formation in the solid phase. We also note that all these COMs are commonly detected in the gas phase in hot core sources and are among the most abundant ones (e.g., Chen et al. 2023). A detailed comparison between gas and ice abundances is left to a future paper (Chen et al., in prep.). Below we discuss particular aspects of the COMs solid phase detections reported in this work.

5.1.1 A polar ice environment

In both protostars, the CH₃CHO, CH₃CH₂OH, and CH₃COOH (continuum-dependent) molecules are diluted in a polar environment (i.e., molecules with high dipole moment), in particular, dominated by H₂O ice. In the general case where H₂O is the major ice component, this means that COM features are under the strong influence of the water ice polarity. In Fig. 13, we compare the MIRI observations with the band shape of three COMs in the best fit, the same COMs in an apolar ice matrix, for instance, mixed with CO ice, and finally, in a CH₃OH-rich environment. CH₃CH₂OH mixed with H₂O band shape at 7.25 µm agrees better with the data than the mixture with CO that has narrow and separate bands compared with the observations of IRAS 2A and IRAS 23385. In the case of the IRAS 2A observations, which has a higher signal-to-noise ratio, CH₃CH₂OH:CO could, potentially, contribute to faint features at 7.16 and 7.29 µm. However, this particular ice mixture is not part of the possible solutions selected by the ENIIGMA code. If present, it would have an ice column density one order of magnitude lower than ethanol mixed with water ice. The other mixtures with CH₃OH are less recurrent because of the narrow 7.2 µm and the stronger 7.5 µm that in the global fit makes the 7.2 µm less prominent. In the case of CH₃CHO, the apolar mixture (CH₃CHO:CO) makes the CH₃ deformation mode narrow by a factor of 2 and red-shifted by 0.02 µm. In the case of CH₃OH-rich mixtures, the peak is red-shifted by 0.03 µm and the feature is broader.

This result does not necessarily contradict laboratory experiments that show that COMs are formed via hydrogenation ofCO molecules (e.g., Fuchs et al. 2009) and C atoms (Fedoseev et al. 2022), as well as that a fraction of CH₃OH ice is mixed with CO (Cuppen et al. 2011). Instead, it points toward a strong effect of H₂O ice on the spectral IR bands of COMs. Moreover, results from the JWST-Ice Age program (McClure et al. 2023), suggest that a fraction of CH₃OH, the most abundant COM, coexists with H₂O in the same ice matrix in cold prestellar clouds. This can be linked to another formation scheme, CH₄ + OH, as studied by Qasim et al. (2020). Regardless of whether CH₃CHO and CH₃CH₂OH reside in a H₂O- or CH₃OH-rich ice, the important message is that the IR band shapes that resemble better the observations are induced by a polar environment.

Finally, CH₃OCHO seems to be primarily mixed with another polar environment that includes CH₃OH instead of H₂O. The spectral shape of CH₃OCHO mixed with H₂O has a broader profile centred at 8.1 µm which does not match well the observations. From the degeneracy analysis of IRAS 23385, it is not completely excluded as a solution, but in the global fits, the presence of HCOOH discards this component as part of the fit. The mixture with CO creates a sharp peak at 8.25 µm that deviates from the observations. The differences between the chemical environment of CH₃CH₂OH and CH₃CHO compared to CH₃OCHO could indicate that methyl formate in IRAS 23385 may have experienced a physical process different from those in IRAS 2A. We discuss this possibility below in Sect. 5.1.2.

Fig. 13 Comparisons between the IR spectra of COMs in different ice matrices with the observed bands of IRAS 2A (top) and IRAS 23385 (bottom). The green checkmarks indicate the data providing the best fit. Excluded data are given by the red cross. Data not part of the best fit and with lower recurrence are indicated by the yellow triangle. The term “polar” in the right panels refers to CO:H2CO:CH3OH and it is the original label published in Terwisscha van Scheltinga et al. (2021).

Fig. 14 Ethanol/acetaldehyde ratio in the gas and solid phases toward different sources.

5.1.2 Possible evidence of ice thermal processing

Among the ice features identified in this work, CH₃OCHO has different spectral shapes in IRAS 2A and IRAS 23385, regardless of the local continuum choice. While the fit of the IRAS 2A spectrum contains CH₃OCHO mixed with CO:H₂CO:CH₃OH, in IRAS 23385, the CH₃OCHO mixed with CH₃OH provides the best fit, and excludes other solutions (see right panels in Fig. 13). A possible cause for this difference can be due to thermal processing, in particular ice distillation. In this process, the most volatile species desorbs from the ice while other species with higher sublimation temperatures remain. For example, during the CH₃OH formation via CO hydrogenation (e.g., Watanabe & Kouchi 2002; Fuchs et al. 2009) both CH₃O and HCO are formed (e.g., Chuang et al. 2016; Garrod et al. 2022; Chen et al. 2023), leading to CH₃OCHO. H₂CO is also an intermediate step toward the CH₃OH formation. In a trivial situation, it is expected that at low temperature (<30 K) CH₃OCHO should be mixed with CO, H₂CO and CH₃OH, the exact components of the ice mixture used to fit the IRAS 2A spectrum. When this ice is warmed up to 50–100 K, CO ice is fully desorbed and H₂CO is partially des-orbed, and the ice would be composed mostly of CH₃OCHO and CH₃OH. Because of the low S/N spectrum of IRAS 23385 we cannot definitely conclude if H₂ CO is present in the ices toward this source.

5.1.3 Gas versus ice: CH₃CH₂OH/CH₃CHO ratio

CH₃CH₂OH and CH₃CHO are chemically linked since the double hydrogenation of acetaldehyde in the solid state leads to ethanol (Fedoseev et al. 2022). In this regard, the ratio of these two molecules provides a way to understand the hydrogenation efficiency of solid-phase molecules to form larger species in protostars.

Figure 14 compares the inferred CH₃CH₂OH/CH₃CHO ratios in the solid phase with those found in the gas phase of many protostars. A persistent ratio of CH₃CH₂OH and CH₃CHO above the unity has been observed in the gas phase with sensitive observations toward protostars in both low- and high-mass star-forming regions (e.g., Yang et al. 2021; Jørgensen et al. 2020; van Gelder et al. 2020; Chen et al. 2023). For example, the ratios summarized by Jørgensen et al. (2020) for IRAS 16293B and Sgr(B2) N2 are 1.9 and 4.5, respectively. For low-mass proto-stars, the CH₃CH₂OH/CH₃CHO ratio taken from van Gelder et al. (2020) are 3.2 (B1-c), 3.0 (S68N), and <4.7 (B1-bS). For the first time, a comprehensive analysis of the 6.8–8.6 µm JWST data allows us to observe the same trend in the solid phase. The CH₃CH₂OH/CH₃CHO ratios derived in this work are 1.7 and 5.8 for IRAS 2A and IRAS 23385, respectively.

This overall trend in protostars suggests that the ethanol/acetaldehyde ratio in the gas phase reflects those in the solid phase. This is a strong hint that ices are the birthplace of COMs detected in the gas phase. Additionally, these COMs are inherited from molecular cloud stages where the ice hydrogenation is efficient to form saturated molecules (e.g., Fedoseev et al. 2022). We also mention that a full three-phase chemical modelling of these sources, assuming reactions within the bulk of the ice, on the ice surface, and in the gas phase can help us constrain the ice chemistry from the observed ice column densities, and provide a more complete explanation for the ethanol/acetaldehyde ratio in both gas and solid-phase. In fact, recent three-phase chemical modelling work by Garrod et al. (2022) assuming general physical conditions of molecular clouds during collapse and warm-up as well as with non-diffusive ice chemistry confirms that the solid-and gas-phase ratios of ethanol/acetaldehyde are consistent with the results presented in this paper.

5.2 Ions in icy grain mantles

The presence of ions in interstellar ices has long been proposed in the literature as a result of the energetic processing of ice mantles (e.g., Grim & Greenberg 1987; Strazzulla & Palumbo 1998; Allamandola et al. 1988; Martinez et al. 2014; Pilling et al. 2010) and acid-base reactions (e.g., Novozamsky et al. 2001). Among these ions, OCN^- has been detected in several protostel-lar envelopes (e.g., van Broekhuizen et al. 2005), and recently with JWST in the Chameleon I molecular cloud toward background stars (McClure et al. 2023). Other ions, such as HCOO^-, were proposed as one of the carriers of the 7.2 and 7.4 µm by Schutte et al. (1999) based on the IR peak position.

In the analysis presented in this work, we find that both HCOO⁻ and OCN⁻ are important contributors to signals in the range between 6.8–8.6 µm. In particular, an important aspect to highlight regarding HCOO⁻ is the same relative intensities of the 7.2 and 7.4 µm bands at low temperature (see Appendix H), which is considered the main carrier of the absorption profiles in the two protostars presented in this work. The fact that higher temperature HCOO⁻ profiles do not match well the observations indicates that the acid-base reaction occurs in cold regions of the protostellar envelope without further thermal processing. One could argue that the HCOO⁻ profile at 14 K is degenerate with the data at 150 K as seen in Appendix H. However, the full IR spectrum of HCOO⁻ has an intense peak at 6.3 µm that is around 8 times stronger than the 7.2 and 7.4 µm bands, and therefore can not be the main carrier of those bands. This strongly suggests that acid-base reactions occur in IRAS 2A and IRAS 23385 and that ions cannot be neglected when interpreting the IR observations.

The formation of HCOO⁻ can occur via the acid-base reaction involving H₂O:HCOOH:NH₃. This is explained by theoretical calculations of HCOOH·(H₂O)_n clusters (Park & Woon 2006). In particular, formic acid partially ionizes when mixed with H₂O, forming HCOO⁻ and H₃O⁺ via proton exchange (Theule et al. 2011). Another reaction to form HCOO⁻ was studied by Bergner et al. (2016), in which only HCOOH and NH₃ are present in the ice, without any H₂O.

Similar to HCOO⁻, OCN⁻ is efficiently formed via acid-base reaction, which has been studied in laboratory experiments starting with HNCO and NH₃ (Schutte & Khanna 2003; Raunier et al. 2003; van Broekhuizen et al. 2004). A large fraction of OCN⁻ is produced in these experiments, as well as the counter ion, ${\text{NH}}_4^{+}$ to maintain the electric neutrality. A possible difficulty in this approach is the high abundance of HNCO in the ice needed to synthesize OCN⁻, associated with the non-detection of this molecule in ices yet. However, gas-phase observations (e.g., Hernández-Gómez et al. 2019) show high abundances of HNCO and laboratory experiments in the ice by Fedoseev et al. (2015) and Noble et al. (2015) show that HNCO can be rapidly formed via an exothermic solid-phase reaction between NH and CO. Additionally, the non-detection of HNCO in ices (e.g., McClure et al. 2023) could be related to the low abundance after being converted into OCN⁻. The presence of OCN⁻ and HNCO in ices has a strong astrobiological appeal. HNCO participates as a peptide bond between two single amino acids as shown by Fedoseev et al. (2015). Moreover, irradiation experiments of ice samples containing OCN⁻/HNCO would lead to the formation of amino acids and their anions, as well. Finally, we highlight that although ${\text{NH}}_4^{+}$ is one of the byproducts of this acid-base reaction, it does not contribute to the absorption bands covered in this work, and therefore, it is not discussed here.

Another mechanism for the OCN⁻ formation is via UV irradiation (van Broekhuizen et al. 2005), which is discussed as being less dominant in low-mass protostars. This route requires high UV fluxes and an abundance of around 30% of NH₃ in the ice. The former condition is satisfied for high-mass protostars, but not otherwise. In fact, a recent paper by Onaka et al. (2022) shows a clear correlation in the high-mass-source AFGL 2006 between the OCN⁻ ice column density with the flux intensity of the HI Brα line. This hydrogen recombination line is a tracer of strong UV radiation and supports an OCN⁻ formation induced by UV photons. Other mechanisms, such as UV-induced flux by cosmic rays are not enough to produce significant amounts of OCN⁻, and low-mass stellar UV would not reach the regions where ice is located. The latter condition exceeds the abundances estimated in the literature for both low- and high-mass protostars, which is between 2% and 15% (Bottinelli et al. 2010).

5.3 Similarities and differences between IRAS 2A, IRAS 23385 and the Comet 67P/G-C

Gas-phase comparisons between the abundances of CHO-bearing COMs with respect to CH₃OH in low- and high-mass protostars, and with the bulk composition of the comet 67P/G-C were made by Drozdovskaya et al. (2019) and Jørgensen et al. (2020). The main conclusion is that there is a good correlation between the abundances of gas-phase COMs in high- and low-mass star-forming regions. On the other hand, there are also differences between the low-mass protostar (IRAS 16293B) and the comet 67P/C-G. In particular, the COMs abundances onto the comet 67P/C-G are enhanced by a factor of up to 10. This difference is interpreted as ice inheritance followed by chemical alteration toward later protostellar phases.

In Fig. 15, we show a comparison between the ice CHO-bearing COMs abundances in IRAS 2A and in the comet 67P/G-C, both with respect to methanol ice. We stress that the peak of the methanol band in IRAS 2A is saturated, and the comparisons are made by assuming that the real methanol ice column density is higher by a factor of 2–3 based on the wings of the C–O band at 9.74 µm (see Appendix K). Under this assumption, we see that COMs and the volatile CH₄ correlate well with the cometary abundances within a factor of 5. This agrees with previous conclusions that these COMs are inherited by comets from early protostellar stages. The scatter, however, can be attributed to further chemical alteration at later stages, or slightly different initial composition in the parental molecular cloud. It is also interesting that the volatiles, CH₄ and SO₂ are enriched in the comet 67P/G-C. This could indicate that these molecules are also formed in the gas phase and condensed at later stages in the comet. However, more analysis of other JWST observations is needed to draw strong conclusions from this correlation.

In Fig. 16, we show a comparison between the COMs abundances with respect to H₂O ice in IRAS 2A, IRAS 23385 and the comet 67P/G-C. The abundances in the low- and high-mass pro-tostar are very close to a linear correlation, which agrees with the results from Jørgensen et al. (2020) between IRAS 16293B and Sgr B2(N2). On the other hand, the COMs abundances compared to H₂O ice are lower in the comet 67P/G-C than in IRAS 2A. Such a high abundance of molecules in low-mass protostars with respect to H₂O ice was noticed before for other molecules, such as NH₃ (Kawakita & Mumma 2011), and CO, CH₄ and CH₃OH (Öberg et al. 2011). The reason is unclear, but it can indicate selective ice destruction of these species compared to H₂O ice in the protosolar nebula or that those COMs were formed in a carbon-poor protostellar envelope (Öberg et al. 2011).

Fig. 15 Comparison between ice abundances with respect to CH3OH ice in IRAS 2A and the comet 67P/G-C (Rubin et al. 2019). Full and hollow circles show the abundances considering $N_{{\text{CH}}_3\text{OH}}^{\text{ice}}=1.5\times {10}^{18} {\text{cm}}^{-2}$ and 1.3 × 1018 cm−2, respectively. The solid line indicates the 1:1 abundance relation, and the dashed lines indicate a cometary abundance lower and higher by a factor of 5.

6 Conclusions

This paper presents a comprehensive analysis in the range between 6.8–8.6 µm of JWST spectra observed toward the IRAS 2A and IRAS 23385 protostars. We performed spectral fits exploring a vast number of IR laboratory spectra that cover, simple molecules, ions, and COMs. Our main conclusions are summarized below:

• The 7.7 µm band is fully resolved in the MIRI spectrum of IRAS 2A and IRAS 23385. As found in previous studies, CH₄ ice is the main carrier of this band. In addition to CH₄, our analysis shows that OCN⁻ and SO₂ contribute to the blue wing of the 7.7 µm band. SO₂ was tentatively detected based on ISO and Spitzer observations, and based on a statistical analysis we confirm that SO₂ is present in ices toward IRAS 2A. For IRAS 23385, SO₂ is classified as likely detection because of a lower S/N and statistical arguments. In the case of OCN⁻, this is the first time that this ion is attributed to the 7.7 µm band. Besides the statistical analysis confirming this detection, OCN⁻ is also seen at 4.59 µm in IRAS 2A with NIRSpec observation, thus confirming the feature at 7.7 µm.

• The 7.2 and 7.4 µm bands observed in IRAS 2A and IRAS 23385 are mainly due to the formate ion (HCOO⁻). Since this ion is efficiently formed via acid-base reactions, and given the detection of OCN⁻ at 7.65 µm, these results demonstrate that this type of chemical reaction may be rather common in interstellar ices.

• The earlier suggestions that CH₃CH₂OH (ethanol) and CH₃CHO (acetaldehyde) contribute to the absorption profiles at 7.2 and 7.4 µm, respectively, are confirmed in this work. In addition, we find that CH₃CH₂OH contributes to other bands in the range considered in this work, most notably, at 6.8–7.05 µm and 7.4–7.6 µm. Other possibilities, such as hydrocarbons (C₂H₂, C₂H₄, C₂H₆) and more complex alcohols (1-propanol, propanal) are less likely based on their absorption profiles. In the case of CH₃CHO, it also contributes to the range between 6.9–7.2 µm.

• We find robust evidence that CH₃OCHO is present in the ices toward IRAS 2A and IRAS 23385. At least two strong bands of this molecule are found between 8.1–8.35 µm and 8.45–8.6 µm. The statistical analysis shows that this spectral component cannot be excluded as a solution for the fit. Based on the criteria for a firm identification, CH₃OCHO is classified as a secure detection. We also find that another COM, CH₃COOH (acetic acid), is present in the fits of IRAS 2A and cannot be excluded based on the confidence interval analysis for two out of three choices of continuum. However, since it is not found in IRAS 23385, and the local continuum strongly affects the shape of the bands around 7.8 µm, we classify acetic acid only as a likely detection. More comparisons with other JWST data will elucidate at which level if acetic acid is present in interstellar ices.

• The COMs found in this work are likely mixed in a polar environment. For example, CH₃CH₂OH and CH₃CHO fit better the observations when mixed with H₂O. These molecules mixed with CO have spectral profiles that deviate from the observed protostellar. In the case of CH₃COOH, the mixture with H₂O ice also provides a good fit to the IRAS 2A spectrum. In the case of CH₃OCHO, the two protostars show different spectral shapes. This differentiation could be related to ice distillation in IRAS 23385, the high-mass protostar.

• Ice-gas ratios between CH₃CH₂OH and CH₃CHO show values above unity in both gas and ice phases. This suggests solid-phase reactions for the formation of these molecules. In addition, the high abundance of CH₃CH₂OH shows that saturated molecules efficiently form in ices because of the high amount of hydrogen available.

• For the first time, we compare ice COMs abundances relative to CH₃OH ice in a protostar and the comet 67P/G-C. Our results indicate that the COMs ice abundances in the comet 67P/G-C correlate well with those in the protostar within a factor of 5, which strongly suggests that COMs in comets are significantly inherited from earlier protostellar phases. In the case of CH₄ and SO₂ with respect to CH₃OH ice, we find an enhancement in the comet 67P/G-C. On the other hand, the abundance comparison with respect to H₂O ice shows that the COMs, SO₂ and CH₄ are depleted in the comet.

The results presented in this paper illustrate how JWST, aided by laboratory experiments, is fully capable of probing the chemical complexity in interstellar ices. Future work using more JWST data with high S/N (e.g., IRAS 15398; Yang et al. 2022) will enable us to verify the recurrence of the COMs found in this paper in other sources. Consequently, this will also allow us to assess the robustness of these detections. Another outlook would be the determination of ice-gas ratios of COMs more complex than CH₃OH. Ultimately, an analysis based on a larger sample will help to answer the question of what the full extent of chemical complexity in interstellar ices might be.

Fig. 16 Comparison between ice abundances with respect to H2O ice. The left panel shows the ice abundances in IRAS 2A and IRAS 23385. The solid line indicates the 1:1 abundance relation, whereas the dashed and dot-dashed lines indicate the enhancement and diminution of IRAS 23385 abundance by a factor of 3. The right panel compares the ice abundances between IRAS 2A and the comet 67P/G-C from Rubin et al. (2019). The solid line indicates the same as in the left panel. The dot-dashed and dashed lines indicate a cometary abundance lower by a factor of 10 and 100, respectively.

Acknowledgements

This paper is dedicated to Professor Harold Linnartz who unexpectedly has passed away on December 31st 2023 between the acceptance of this paper and its publication. Professor Linnartz was chair of the Leiden Laboratory for Astrophysics and was well known for his research on “unlocking the chemistry of the heavens”. Many of the IR ice spectra of COMs used in this paper were collected by students under his supervision. He was also otherwise actively involved in the JOYS+ program and its science. Professor Linnartz was thrilled with the results presented in this paper, given the quality of the JWST data and their significance for astrochemistry. As additional and important acknowledgements, the following National and International Funding Agencies funded and supported the MIRI development: NASA; ESA; Belgian Science Policy Office (BELSPO); Centre Nationale d'Études Spa-tiales (CNES); Danish National Space Centre; Deutsches Zentrum fur Luftund Raumfahrt (DLR); Enterprise Ireland; Ministerio De Economiá y Competivi-dad; The Netherlands Research School for Astronomy (NOVA); The Netherlands Organisation for Scientific Research (NWO); Science and Technology Facilities Council; Swiss Space Office; Swedish National Space Agency; and UK Space Agency. We thank the anonymous reviewer for the careful reading of our manuscript and the comments and suggestions that improved the clarity of this work. W.R.M.R., E.v.D., K.S. N.B. M.v.G., L.F. and H.L. acknowledge the funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No. 291141 MOLDISK). W.R.M.R. thanks Niels Ligterink and Maria Droz-dovskaya for useful discussions about the 67P/G-C comet abundances. We are grateful for continuing support through NOVA, the Netherlands Research School for Astronomy, the NWO through its Dutch Astrochemistry Program (DANII). The present work is closely connected to ongoing research within INTERCAT, the Center for Interstellar Catalysis located in Aarhus, Denmark. The work of M.E.R. was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration (80NM0018D0004). L.M. acknowledges the financial support of DAE and DST-SERB research grants (SRG/2021/002116 and MTR/2021/000864) from the Government of India. L.F acknowledges the grant TOP-1 614.001.751 from the Dutch Research Council (NWO). T.R. acknowledges support from ERC grant no. 743029 EASY. T.H. acknowledges support from the ERC Advanced grant no. Origins 83 24 28. H.B. acknowledges support from the Deutsche Forschungs-gemeinschaft in the Collaborative Research Center (SFB 881) “The Milky Way System” (subproject B1). P.J.K. acknowledges financial support from the Science Foundation Ireland/Irish Research Council Pathway programme under Grant Number 21/PATH-S/9360. A.C.G. has been supported by PRIN-INAF MAINSTREAM 2017 “Protoplanetary disks seen through the eyes of new generation instruments” and from PRIN-INAF 2019 “Spectroscopically tracing the disk dispersal evolution (STRADE)”. K.J. acknowledges the support from the Swedish National Space Agency (SNSA).

Appendix A 9.8 µm band compared to ISO and Spitzer sources

Figures A.1 and A.2 show comparisons with the silicate feature toward GCS 3 and with other low- and high-mass pro-tostars observed with Spitzer and ISO. The first note is that both IRAS 2A and IRAS 23385 have a broader 9.8 µm silicate profile compared to GCS 3. This highlights the need of considering different grain compositions when subtracting the silicate absorption band. Secondly, IRAS 2A and IRAS 23385 show a similar blue side of the spectrum with other sources. On the other hand, the red profile has more differences, which are related to the amount of H₂O ice toward the source. In fact, icy-grain models by Ossenkopf & Henning (1994) show that coagulated icy grains have a prominent spectral bump around 12 µm because of the H₂O ice libration band.

Fig. A.1 Comparison among MIRI/JWST spectrum of IRAS 2A, the silicate profile of GCS 3 and Spitzer/IRS spectrum of different low-mass protostars. These spectra are normalized by the optical depth at 9.8 µm.

Fig. A.2 Comparison among MIRI/JWST spectrum of IRAS 23385, the silicate profile of GCS 3 and Spitzer/IRS spectrum of different high-mass protostars. These spectra are normalized by the optical depth at 9.8 µm.

Appendix B Effect of using different silicates for subtraction in IRAS 2A

In Section 3.1, we describe how the silicate feature is removed from IRAS 23385 and IRAS 2A MIRI-MRS spectra using a synthetic silicate profile combining olivine and pyroxene. In this section, we compare subtractions using different profiles to check how this step affects the shape of the 6.8–8.6 µm region. Figure B.1 (top) shows three silicate profiles, scaled to IRAS 2A spectrum, i) the synthetic silicate used in Section 3.1, ii) the silicate profile of GCS 3, and iii) the MgSiO₃ profile taken from Poteet et al. (2015). The peak optical depths are set to match the synthetic silicate profile. We point out that this comparison is focused on the 6.8–8.6 µm region, and therefore mismatches of these silicate spectra at longer wavelengths are not relevant for this specific analysis. We briefly mention SiO₂ (silica), as another possible candidate for the blue wing of the 9.8 µm band. However, silica has not been found in absorption in protostars so far and its relatively narrow profile at 18 µm is not seen in both sources investigated in this paper. Nevertheless, if present, SiO₂ would not affect the COMs bands investigated in this paper because of its broadband at around 8.3 µm. At most, it could slightly reduce the HCOOH ice column density.

Figure B.1 (bottom) shows the silicate subtracted optical depth spectrum of IRAS 2A. Both synthetic and enstatite silicates result in similar spectral profiles between 6.8–8.6 µm. On the other hand, the silicate subtraction using GCS 3 creates an unrealistic absorption excess (also observed in Boogert et al. (2008)) that deviates from the other two profiles.

Fig. B.1 Comparing silicate profiles in IRAS 2A. The top panel shows the synthetic silicate providing the best fit in IRAS 2A (blue) compared to GCS 3 (red) and enstatite (green) silicate model from Poteet et al. (2015). The bottom panel shows the IRAS 2A spectrum after removing the three silicate models.

Appendix C List of laboratory data

A comprehensive list of molecules was used in this paper to search for the best fit of the region between 6.8 and 8.6 µm. This list includes COMs in different mixtures as well as simple molecules. The full list is shown in Table C.1.

Appendix D Acetic acid and OCN⁻ band strengths

Fig. D.1 Optical contacts of acetic acid. Vertical dotted lines indicate the features with band strength calculated in this paper.

Most of the band strengths for the molecules detected in the range addressed in this paper are available in the literature. However, to the best of our knowledge, the band strengths of acetic acid and OCN⁻ between 6.8 and 8.6 µm were not calculated before. In particular, the band strength of acetic acid at 5.8 µm is often assumed to be the same as, or corrected from, the gas-phase acetic acid (e.g., Öberg et al. 2009b; Chuang et al. 2020) based on Marechal (1987).

We derive the absolute (A) and apparent (A') band strengths of acetic acid. The absolute band strengths are derived from the imaginary refractive index, by the following equation: $$A=\frac{m}{\rho N_A}{\displaystyle {\int }_{v_1}^{v_2}4}\pi vk(v)dv$$(D.1)

where m is the molar mass of acetic acid in g mol⁻¹ (60.052), p is the density of acetic acid (0.892 g cm⁻³; Hudson et al. 2020), N_A is the Avogrado’s number, v is the wavenumber and k is the imaginary refractive index. We derive k using the recent version of the NKABS code (Rocha & Pilling 2014). This code calculates the real (n) and imaginary (k) refractive index from the absorbance spectrum (Abs_v). The input spectrum is taken from Hudson & Gerakines (2019), who also estimated the thickness of the ice (d = 2.1 µm) and the refractive index around 700 nm (n₀ = 1.29). The n and k values are shown in Fig. D.1. Finally, the band strengths derived for two acetic acid bands are listed in Table D.1.

The apparent band strength is calculated by: $$A^{\prime }=2.303\frac{m}{\rho N_{A}d}{\displaystyle \int A}b{s}_{v}dv$$(D.2)

As we can see, both A and A^' values for acetic acid are similar.

In the case of OCN⁻, we derive the apparent band strength at 7.62 µm (1312 cm⁻¹). The OCN⁻ band at 4.59 µm (2175 cm⁻¹) is used as a reference since the band strength is known (1.3×10⁻¹⁶ cm molecule⁻¹; van Broekhuizen et al. 2005). We use the following equation to derive the OCN⁻ apparent band strength at 7.62 µm: $$A_{7.62\mu m}^{\prime }=1.3\times {10}^{-16}\frac{{\displaystyle \int A}b{s}_{1312{\text{cm}}^{-1}}dv}{{\displaystyle \int A}b{s}_{2175{\text{cm}}^{-1}}dv},$$(D.3)

where $Ab{s}_{2175{\text{cm}}^{-1}}$ and $Ab{s}_{1312{\text{cm}}^{-1}}$ are the OCN⁻ bands at the given wavelengths.

Appendix E Laboratory spectra: removing H₂O and CH₃OH ice features

In this section, we demonstrate how to remove the H₂O and CH₃OH features of COMs spectra containing these two molecules in the range between 6.5 and 9.0 µm. This process is necessary when analysing the observational spectrum between 6.8 and 8.6 µm using local continuum subtraction. Figure E.1 (top) shows the H₂O:CH₃CH₂OH ice spectrum (Terwisscha van Scheltinga et al. 2018) and the polynomial (4th-order) function used to trace a baseline under the CH₃CH₂OH bands. In this spectrum, the ethanol features overlap with the broad shoulder of the H₂O ice bending mode, which is represented by the polynomial fit. In Figure E.1 (bottom), we show the CH₃OH:CH₃CH₂OH IR spectrum (Terwisscha van Scheltinga et al. 2018). Since both molecules are alcohols they share functional groups, which makes it harder to isolate the ethanol features of methanol. Because of the high dilution factor (20:1), it is not feasible to disentangle the features of the two molecules below 6.8 µm, and therefore, we separate the CH₃CH₂OH bands between 6.9 and 8.6 µm. The step is performed in three stages: (i) high-order polynomial (7th) between 6.9–7.15 µm to remove the CH₃OH shoulder, (ii) 3rd-order polynomial between 7.15– 7.45 to isolate the ethanol band at 7.2 µm and account for the strong curvature in the data, and (iii) 5th-order polynomial between 7.45–8.7 µm to extract the other ethanol bands.

Fig. E.1 Isolating the CH3CH2OH ice bands from H2O and CH3OH ices. Polynomial functions anchored at the points indicated by the grey areas are used. The top panel shows the ethanol mixed with H2O ice, and the bottom panel displays the mixture with CH3OH.

Appendix F Laboratory baseline correction: avoiding spurious features

IR spectra of ices are recorded in the laboratory using Fourier Transform Infrared Spectroscopy (FTIR) and interference effects are corrected using a spline or polynomial function, the so-called baseline correction. In Figure F.1 we show the baseline correction of the H₂O:CH₃CH₂OH IR spectrum at 15 K. Panels a and b show the laboratory spectrum and a 7th-order polynomial used to correct the interference effect. The difference between these two panels is that in panel a, we use as many points as possible to trace the baseline, whereas in panel b, we use only half of the points available in the range of 7.8–8.5 µm. This creates a small fluctuation in the polynomial function. Zoom-ins of both cases are shown in panels c and d, respectively. In panel e, we show the baseline data. The spectrum with no polynomial inflexion contains the real CH₃CH₂OH IR features, whereas the other data has spurious features at 8.0 and 8.3 µm. This demonstrates that care must be taken when correcting IR spectrum baselines to avoid creating spurious features.

Fig. F.1 Baseline correction of the IR ice spectrum of H2O:CH3CH2OH data. Panels a and b show the IR ice spectra and the baselines. The blue dotted squares highlight the zoom-ins shown in panels c and d displaying polynomial baselines without and with little inflexion, respectively. Panel e demonstrates the effect of non-accurate baseline subtraction in the ice spectrum.

Appendix G Incremental version of the fits for IRAS2A and IRAS 23385

Figure G.1 displays the best fits for IRAS 2A (left) and IRAS 23385 (right) by adding one component at a time in each panel. This allows for a better understanding of how each component contributes to the fit.

Fig. G.1 Incremental version from top to bottom of the best fits with the ENIIGMA fitting tool for IRAS 2A (left) and IRAS 23385 (right). The molecule label follows the colour code of the laboratory components in the fit.

Appendix H HCOO⁻ at 14,150 and 210 K

Fig. H.1 Comparison between the 7.2 and 7.4 µm band of IRAS 2A (left) and IRAS 23385 (right) with the HCOO− absorption profiles at 14, 150 and 210 K. These IR ice spectra are taken from Gálvez et al. (2010).

We compare the absorption profiles of HCOO⁻ at three temperatures (14 K, 150 K, 210 K) with the 7.2 and 7.4 µm bands in IRAS 2A and IRAS 23385 (Figure H.1). Our analysis shows that the spectrum at 14 K provides the best fit. The spectrum at 150 K is excluded because the second peak (~ 7.4 µm) is broader and shifted compared to the observations. Likewise, the band shape of the highest temperature data (150 K) does not match with both protostars.

Appendix I Confidence intervals

The top and bottom panels in Figure I.1 show the confidence intervals for IRAS 2A in the ranges of 6.8–7.5 and 7.8–8.6, respectively. The confidence intervals for IRAS 23385 are shown in Figures I.2 and I.3. For IRAS 2A, it can be noted that all components are essential to the fit, and cannot be excluded. On the other hand, for IRAS 23385, the SO₂ band can be excluded as a solution if the OCN⁻ band is slightly intense.

Fig. I.1 The top and bottom corner plots show the IRAS 2A coefficient confidence intervals for the range between 6.86–7.5 (top) and 7.8–8.6 µm (bottom), respectively.

Fig. I.2 Top and bottom corner plots show the IRAS 23385 coefficient confidence intervals for the range between 6.86–7.5 (top) and 7.5–7.8 µm (bottom), respectively.

Fig. I.3 Same as in Figure I.2, but for the range 7.8–8.6 µm.

Appendix J Different local continuum profiles between 6.8 and 8.6 µm – IRAS 2A

The first three panels of Figure J.1 show different continuum profiles between 6.8–8.6 µm in the IRAS 2A spectrum. The top panel is the version adopted for the analysis in this paper that traced a third-order polynomial to the guiding points. The second panel displays the fourth-order polynomial where the red dot is added to the guiding points. In this case, the continuum is slightly elevated at shorter wavelengths to accommodate the fit to the extra point at 8.5 µm. The third panel presents the continuum when two extra points are added (λ = 7.8 µm and 8.5 µm), and a sixth-order polynomial is used. All subtracted spectra using these three approaches are shown in the bottom panel of Figure J.1.

New fits for IRAS 2A are performed on the other two optical depth spectra obtained from different local continuum choices, which are shown in Figure J.2. The panel at the top shows that all components remain needed to reach the best fit. Only the 8.5 µm band of CH₃OCHO exceeds the observations, and this is caused by the guiding point added at 8.5 µm. The bottom panel of Figure J.2 shows another fit, where two extra points are added at 7.8 and 8.5 µm. This version of the fit keeps all components, except CH₃COOH, which is excluded because of the anchor point added at 7.8 µm, where CH₃COOH has a strong feature. Other issues are seen in this fit, for example, around 7 µm, 7.3 µm, and a poor fit between 7.8 and 8.6 µm. Despite all these differences, this analysis shows that CH₃CHO, CH₃CH₂OH and CH₃OCHO are still robust detections, and cannot be excluded from the fits of the IRAS 2A spectrum.

Fig. J.1 Effect of the continuum choice on the local continuum subtracted spectrum of IRAS 2A. The first panel shows the continuum adopted as the best model in this paper. The second and third panels show two other continuum options by adding the red and orange points, respectively. The optical depth subtracted spectra are shown in the fourth panel.

Fig. J.2 Alternative fits of IRAS 2A optical depth spectrum with different local continuum choices (see Figure J.1). The top panel shows the fits with all the components after subtracting the red local continuum with an extra point at 8.5 µm. The bottom panel shows the same as in the top panel, but considering two extra points for the local continuum (7.8 and 8.5 µm). Only CH3COOH is excluded in this fit.

Appendix K Water and methanol ice column densities

The water ice column densities in both protostars are calculated from the libration mode around 12.5 µm. To determine the profile of the band, we combine water ice spectra at different temperatures (Figure K.1). As pointed out by Boogert et al. (2008), this band is sensitive to the grain geometry. For this reason, we assume small spherical water ice grains which is consistent with Boogert et al. (2008) to fit the libration mode best. The optical constants for ices at 15, 75, and 160 K are taken from Rocha et al. (2022). The water libration band is fitted with two components representing different temperatures. It is likely that the libration band is sensitive to a range of temperatures between 15 and 160 K, but addressing this is beyond the scope of this work. Despite this simplification, we may note that solely cold water ice is not enough to fit the libration band of IRAS 23385 and IRAS 2A. In particular, IRAS 2A has a strong blue wing excess that requires H₂O ice at 160 K. The water ice column density is shown in Table 2.

In the case of CH₃OH ice, we use the band at 9.8 µm to derive a column density. We fit a Gaussian profile to the feature at 9.8 µm (solid curve in Figure K.2) and multiply it by a factor of 2 (dashed curve) and 3 (dotted curve). It is unlikely that the solid curve accounts for all CH₃OH absorption, and therefore the CH₃OH ice column densities used in this paper correspond to the other two Gaussian profiles.

Fig. K.1 Fits for the H2O ice libration band for IRAS 23385 (top) and IRAS 2A (bottom). The best fit is found by combining two H2O ice grain-shaped corrected spectra: 15K and 75 K for IRAS 23385 and 15 K and 160 K for IRAS 2A.

Fig. K.2 NH3 and CH3OH features in the bottom of the silicate band and H2O ice subtracted spectra of IRAS 2A. Three Gaussian profiles are scaled to the CH3OH band to indicate different ice column densities.

Appendix L Comparison with hydrocarbons

Because of the degeneracy intrinsic to the ice fittings, particularly regarding COMs that share the same functional groups, we show in this section a comparison of hydrocarbons (C_xH_x) with the IRAS 2A spectrum. These molecules participated in the global fits between 6.8 and 8.6 µm, but are discarded as solutions. Therefore, these comparisons serve as an additional check that these components are not part of the global minimum solution. Figure L.1 shows scaled IR spectra of pure C₂H₂, C₂H₄, C₂H₆ compared to IRAS 2A. The scaling factors are arbitrarily chosen to match the IRAS 2A absorption profile at 7.3 µm. It can be seen that C₂H₄ does not have any contribution to the 7.2 and 7.4 µm. On the other hand, the absorption bands of C₂H₂ and C₂H₆ that could be hidden at 7.3 µm, exceed the IRAS 2A absorption profiles around 12 and 13.5 µm.

Fig. L.1 Comparison between the 7.2 and 7.4 µm band of IRAS 2A (silicate and H2O ice subtracted) and hydrocarbons (C2H2, C2H4 and C2H6).

References

All Tables

All Figures

	Fig. 1 MIRI MRS spectrum of the high-mass protostar IRAS 23385+6053 (top) and the low-mass protostar NGC 1333 IRAS 2A (bottom). Strong gas-phase emission lines are masked in both spectra. The assignments for the absorption bands are given and differentiated by the colours. The shaded cyan area indicates the “COMs region” that is studied in this work.
In the text

	Fig. 2 Summary of the methodology used in this paper.
In the text

Fig. 3 Continuum (a, b) and silicate subtraction (c, d) steps in IRAS 23385 (left) and IRAS 2A (right). The spectrum of IRAS 23385 is binned by a factor of two between 9 and 11 µm due to the saturated silicate profile. A third-order polynomial function is used to trace the continuum, and the silicate profile is a combination of two laboratory silicate spectra (olivine and pyroxene). Panels e and f show the silicate subtracted spectra of both protostars, with major features labelled.

In the text

	Fig. 4 Observed spectrum of IRAS 2A (blue) surrounding the 7.2 µm and 7.4 µm ice absorption features with SO2 (ν3) emission. The best-fit SO2 model of van Gelder et al. (2024) is presented in green and the SO2 subtracted spectrum of IRAS2A is shown in orange. Subtracting the SO2 emission clearly reveals the 7.4 µm absorption feature.
In the text

	Fig. 5 Procedure for isolating weaker absorption features in the spectra of IRAS 23385 and IRAS 2A. Left: local continuum in the 6.8–8.6 spectral range using a fourth-order polynomial function. Right: isolated 6.8-8.6 optical depth spectrum in both protostars.
In the text

	Fig. 6 ENIIGMA fits of IRAS 23385 (top) and IRAS 2A (bottom). Gas-phase lines are masked. The figure labels show the ice mixture used in the fits, and a simplified version with the names of the chemical species names is shown close to the bands. The corresponding temperature of these laboratory spectra ranges between 10 K and 15 K.
In the text

	Fig. 7 Corner plots showing the IRAS 2A coefficient confidence intervals for the range between 7.5 and 7.8 µm. The grey scale colour is the Δχ2 map derived from a total of 5000 values. Yellow and red contours represent 2 and 3σ significance, respectively.
In the text

Fig. 8 Bar plot showing the recurrence of the difference chemical species to the fit in IRAS 2A (top) and IRAS 23385 (bottom). Recurrences above 50% are considered essential to the fit and robust detection. CH3COOH in IRAS 23385 does not occur in any of the solutions. In the case of IRAS 2A (top panel) we present hatched bars that show the recurrence of the solutions if the errors in optical depth are increased by a factor of three.

In the text

	Fig. 9 Bar plot showing the recurrence of the same ice mixture, but at different temperatures. Plots inside the blue and green polygons refer to IRAS 2A and IRAS 23385, respectively.
In the text

Fig. 10 IR spectra of COMs not contributing to the global fit, but with absorption features in the range between 6.8 and 8.6 µm. These spectra are scaled to the IRAS 2A spectrum (black) to derive upper limit column densities. The vertical dotted lines are colour-coded and indicate the features used to derive the ice column densities. A small horizontal shift in the vertical lines is used to distinguish the bands at 7.3 µm.

In the text

	Fig. 11 Comparison of the ωCH mode of glycine IR spectrum (red) with local subtracted spectrum of IRAS 2A (black). The intensities of the glycine band are given for three ice column densities.
In the text

	Fig. 12 Experimental IR spectra of ethanol, propanal and propanol compared to IRAS 2A spectrum in the range between 6.8–8.6 µm. The grey area is not considered in the fits, but it is shown here to highlight the C–H stretching mode of these molecules. For better readability of this figure, the gas-phase emission lines between 6.8–7.2 µm are masked.
In the text

Fig. 13 Comparisons between the IR spectra of COMs in different ice matrices with the observed bands of IRAS 2A (top) and IRAS 23385 (bottom). The green checkmarks indicate the data providing the best fit. Excluded data are given by the red cross. Data not part of the best fit and with lower recurrence are indicated by the yellow triangle. The term “polar” in the right panels refers to CO:H2CO:CH3OH and it is the original label published in Terwisscha van Scheltinga et al. (2021).

In the text

	Fig. 14 Ethanol/acetaldehyde ratio in the gas and solid phases toward different sources.
In the text

Fig. 15 Comparison between ice abundances with respect to CH3OH ice in IRAS 2A and the comet 67P/G-C (Rubin et al. 2019). Full and hollow circles show the abundances considering $N_{{\text{CH}}_3\text{OH}}^{\text{ice}}=1.5\times {10}^{18} {\text{cm}}^{-2}$ and 1.3 × 1018 cm−2, respectively. The solid line indicates the 1:1 abundance relation, and the dashed lines indicate a cometary abundance lower and higher by a factor of 5.

In the text

Fig. 16 Comparison between ice abundances with respect to H2O ice. The left panel shows the ice abundances in IRAS 2A and IRAS 23385. The solid line indicates the 1:1 abundance relation, whereas the dashed and dot-dashed lines indicate the enhancement and diminution of IRAS 23385 abundance by a factor of 3. The right panel compares the ice abundances between IRAS 2A and the comet 67P/G-C from Rubin et al. (2019). The solid line indicates the same as in the left panel. The dot-dashed and dashed lines indicate a cometary abundance lower by a factor of 10 and 100, respectively.

In the text

	Fig. A.1 Comparison among MIRI/JWST spectrum of IRAS 2A, the silicate profile of GCS 3 and Spitzer/IRS spectrum of different low-mass protostars. These spectra are normalized by the optical depth at 9.8 µm.
In the text

	Fig. A.2 Comparison among MIRI/JWST spectrum of IRAS 23385, the silicate profile of GCS 3 and Spitzer/IRS spectrum of different high-mass protostars. These spectra are normalized by the optical depth at 9.8 µm.
In the text

	Fig. B.1 Comparing silicate profiles in IRAS 2A. The top panel shows the synthetic silicate providing the best fit in IRAS 2A (blue) compared to GCS 3 (red) and enstatite (green) silicate model from Poteet et al. (2015). The bottom panel shows the IRAS 2A spectrum after removing the three silicate models.
In the text

	Fig. D.1 Optical contacts of acetic acid. Vertical dotted lines indicate the features with band strength calculated in this paper.
In the text

	Fig. E.1 Isolating the CH3CH2OH ice bands from H2O and CH3OH ices. Polynomial functions anchored at the points indicated by the grey areas are used. The top panel shows the ethanol mixed with H2O ice, and the bottom panel displays the mixture with CH3OH.
In the text

	Fig. F.1 Baseline correction of the IR ice spectrum of H2O:CH3CH2OH data. Panels a and b show the IR ice spectra and the baselines. The blue dotted squares highlight the zoom-ins shown in panels c and d displaying polynomial baselines without and with little inflexion, respectively. Panel e demonstrates the effect of non-accurate baseline subtraction in the ice spectrum.
In the text

	Fig. G.1 Incremental version from top to bottom of the best fits with the ENIIGMA fitting tool for IRAS 2A (left) and IRAS 23385 (right). The molecule label follows the colour code of the laboratory components in the fit.
In the text

	Fig. H.1 Comparison between the 7.2 and 7.4 µm band of IRAS 2A (left) and IRAS 23385 (right) with the HCOO− absorption profiles at 14, 150 and 210 K. These IR ice spectra are taken from Gálvez et al. (2010).
In the text

	Fig. I.1 The top and bottom corner plots show the IRAS 2A coefficient confidence intervals for the range between 6.86–7.5 (top) and 7.8–8.6 µm (bottom), respectively.
In the text

	Fig. I.2 Top and bottom corner plots show the IRAS 23385 coefficient confidence intervals for the range between 6.86–7.5 (top) and 7.5–7.8 µm (bottom), respectively.
In the text

	Fig. I.3 Same as in Figure I.2, but for the range 7.8–8.6 µm.
In the text

	Fig. J.1 Effect of the continuum choice on the local continuum subtracted spectrum of IRAS 2A. The first panel shows the continuum adopted as the best model in this paper. The second and third panels show two other continuum options by adding the red and orange points, respectively. The optical depth subtracted spectra are shown in the fourth panel.
In the text

Fig. J.2 Alternative fits of IRAS 2A optical depth spectrum with different local continuum choices (see Figure J.1). The top panel shows the fits with all the components after subtracting the red local continuum with an extra point at 8.5 µm. The bottom panel shows the same as in the top panel, but considering two extra points for the local continuum (7.8 and 8.5 µm). Only CH3COOH is excluded in this fit.

In the text

	Fig. K.1 Fits for the H2O ice libration band for IRAS 23385 (top) and IRAS 2A (bottom). The best fit is found by combining two H2O ice grain-shaped corrected spectra: 15K and 75 K for IRAS 23385 and 15 K and 160 K for IRAS 2A.
In the text

	Fig. K.2 NH3 and CH3OH features in the bottom of the silicate band and H2O ice subtracted spectra of IRAS 2A. Three Gaussian profiles are scaled to the CH3OH band to indicate different ice column densities.
In the text

	Fig. L.1 Comparison between the 7.2 and 7.4 µm band of IRAS 2A (silicate and H2O ice subtracted) and hydrocarbons (C2H2, C2H4 and C2H6).
In the text