© The Authors 2024

1 Introduction

The ingredients constituting interstellar material, stars, and eventually planets, originate in dense molecular clouds where atoms and small molecules are accreted from the gas phase onto cold dust grains to form icy mantles (Herbst & van Dishoeck 2009; Caselli & Ceccarelli 2012; Boogert et al. 2015). The abundance and composition of these ice mantles are therefore a direct probe of the chemical complexity of the molecular cloud at the time of ice formation. Susceptible to changes in these chemical conditions are isotope ratios that have been well studied both in the gas-phase and solid state in multiple celestial bodies (Nomura et al. 2023, and references therein). The variability of the isotope abundance across different astronomical environments is therefore indicative of the chemical evolution of these systems and of the physicochemical processes to which the pristine material has been subjected.

The carbon isotope ratio in particular is an attractive contender for isotope chemistry studies for numerous reasons. First, carbonaceous ices such as ¹²CO and ¹²CO₂ are abundant and readily detected in infrared observations (de Graauw et al. 1996; Whittet et al. 1998; Gerakines et al. 1999; Keane et al. 2001; Pontoppidan et al. 2003; Bergin et al. 2005; Pontoppidan et al. 2008; Öberg et al. 2011; McClure et al. 2023), and if the sensitivity allows, their weaker isotopolog bands ¹³CO and ¹³CO₂ are also observable in the same spectral range. This enables a carbon isotope analysis in the solid state (Boogert et al. 2000, 2002b; Pontoppidan et al. 2003; Brunken et al. 2024). Additionally, CO and CO₂ ice have multiple vibrational modes across the infrared spectrum that facilitate comparison between ratios derived from the different absorption bands.

CO is also the second most abundant gas-phase molecule after H₂ in the interstellar medium (ISM), and its isotopologs ¹²C¹⁶O, ¹³C¹⁶O, and ¹²C¹⁸O therefore have strong gas-phase molecular lines in the submillimeter and infrared, making them ideal for gas-phase studies in various environments (Smith et al. 2015; Nomura et al. 2023), including exoplanet atmospheres (Zhang et al. 2021a,b; Line et al. 2021; Gandhi et al. 2023; de Regt et al. 2024). This plethora of observational possibilities enables us to draw comparisons between gas-phase and solid-state abundances and to further elucidate the origin of interstellar material and its evolution from stellar nurseries to the time of planet formation.

Interstellar material is initially enriched with ¹³C during the CNO-cycle of stellar nucleosynthesis (M_★ > 1.3 M_⊙) when ¹³N atoms are converted into ¹³C atoms. A fraction of these ¹³C atoms remains in the star and is later released into the ISM at the end of its life cycle (Milam et al. 2005, 2009; Prantzos 2008). Once formed, the main fractionation processes that affect the regional carbon isotope ratios are isotope-exchange reactions, isotope-selective photodissociation, and possibly, gas-ice partitioning (Watson 1976; Langer et al. 1984; Langer & Penzias 1993; van Dishoeck & Black 1988; Visser et al. 2009; Smith et al. 2015).

The yardsticks against which all observations are usually measured are the ratios derived for the local ISM by Wilson & Rood (1994, ¹²C/¹³C ~ 77 ± 7) from CO studies and by Milam et al. (2005, ¹²C/¹³C ~ 68 ± 15) from CN studies, as well as the solar abundance ratio (Wilson & Rood 1994, ¹²C/¹³C ~ 89). In diffuse clouds, Ritchey et al. (2009) derived ratios from optical line absorption studies of CH+ that were consistent with the ISM standard (¹²CH+/¹³CH+ ~ 74). Conversely, the effect of selective photodissociation was observed in other diffuse clouds, where higher ratios were extracted from CO measurements (Lambert et al. (1994, ¹²CO/¹³CO ~ 167) and Federman et al. (2003, ¹²CO/¹³CO ~ 125 and 117)).

In dense star-forming regions, Jørgensen et al. (2016, 2018) mainly found sub-ISM ratios from gas-phase studies of complex organic molecules (COMs), with considerable deviations between the different species (¹²C/¹³C ~ 27–67). Large discrepancies were also observed in studies of protoplanetary disks, where the values varied drastically depending on the line of sight and on the region of the disk in which the ratio was measured (Hily-Blant et al. 2019; Booth et al. 2019, 2024; Bergin et al. 2024). Finally, cometary ratios and ratios derived from chon-drites were found to be similar to the standard solar value (~89) (Alexander et al. 2007; Bockelée-Morvan et al. 2015; Hässig et al. 2017; Nomura et al. 2023).

In the solid state, Boogert et al. (2000) examined a number of massive young stellar objects (MYSOs), three low-mass young stellar objects (LYSOs), and two clouds and found ¹²CO₂/¹³CO₂ ~ 52–113, consistent with the standard ISM ratio. Additional studies of CO ice for one MYSO (Boogert et al. 2002a) and one LYSO (Pontoppidan et al. 2003) yielded ¹²CO/¹³CO ~ 71 and 69, respectively. The sample of solar-mass protostars remained limited, however, since previous observational facilities lacked the sensitivity required to observe most of these low-mass objects. Moreover, the partial spectral coverage of these observatories often also meant that CO and CO₂ ice could not be examined simultaneously.

With the exquisite sensitivity of the James Webb Space Telescope (JWST) (Gardner et al. 2023; Wright et al. 2023; Argyriou et al. 2023; Jakobsen et al. 2022; Böker et al. 2022), we are now able to study the various vibrational modes of CO and CO₂, including their weaker isotopolog absorption features, for the first time in solar-mass protostars. Moreover, the high spectral resolution of the JWST offers a unique opportunity to perform a detailed profile analysis of these strong and weak ice features. McClure et al. (2023), for instance, used this unprecedented sensitivity to quantify solid-state ¹²CO₂/¹³CO₂ and ¹²CO/¹³CO ratios in the Cha I dark molecular cloud (¹²CO₂/¹³CO₂ ~ 69– 87 and ¹²CO/¹³CO ~ 99–169). Brunken et al. (2024) also used the high spectral resolution of the JWST to perform a detailed profile analysis of the 4.39 µm band of ¹³CO₂ and showed that the band can be decomposed into five principal components (Pontoppidan et al. 2008) that are representative of the chemical and thermal environment of the CO₂ ice.

In this paper, we present high-resolution JWST NIRSpec and MIRI observations of the ¹²CO₂ v₁ + v₂ and the 2v₁ + v₂ combination modes (2.70 and 2.77 µm), the ¹²CO₂ v₃ stretching mode (4.27 µm), the ¹³CO₂ v₃ stretching mode (4.39 µm), the ¹²CO₂ v₂ bending mode (15.2 µm), the ¹²CO 2-0 overtone mode (2.35 µm), the ¹²CO 1-0 stretching mode (4.67 µm), and the ¹³CO 1-0 stretching mode (4.78 µm) for 17 low-mass Class 0/I protostars that were observed as part of the JWST Observations of Young protoStars (JOYS+) program (van Dishoeck et al. 2023; Beuther et al. 2023, Ressler et al. in prep.). We derive column densities of ¹²CO₂, ¹³CO₂, ¹²CO, and ¹³CO for each vibrational mode and determine the carbon isotope ratio for the sources in this sample. We build on the work of (Boogert et al. 2000) by expanding the sample of solar-mass objects and examining their isotope reservoir.

This paper is structured as follows. In Section 2 we present our data and describe the data reduction process and the methods for analyzing the spectra. In Section 3 we derive column densities and isotope ratios for the bands as a whole, and we discuss the results in Section 4. Finally, in Section 5 we provide a summary and the concluding remarks.

2 Observations and methods

2.1 Observations

Our observations were taken as part of the JWST Observations of Young protoStars (JOYS+) Cycle 1 guaranteed time program (GTO) NIRSpec (PI: E.F. van Dishoeck, ID: 1960) and MIRI (PI: M. E. Ressler, ID: 1236). The data consist of NIRSpec and MIRI spectra of 17 Class 0/I sources, including binaries observed using the G235M, G235H, G395M, and G395H modes (R = λ /Δ λ = 1000 and 2700 for the NIRSpec observations).

The NIRSpec datacubes were reprocessed using the JWST pipeline (version 1.13.4) and the corresponding CRDS context jwst_1231.pmap (CRDS_VER = ‘11.17.19’). Additional corrections outside the JWST pipeline were applied to enhance the data quality. The calwebb_detector1 step was executed with standard parameters. Subsequently, NSClean (Rauscher 2024) was used to address faint vertical banding and picture frame noise in the rate files. Following this, calwebb_spec2 was run to generate the cal files. A systematic search was conducted in the cal files to identify warm pixels for flagging, as detected in the MAST final cubes. Warm pixels labeled UNRELIABLE_FLAT and NO_SAT_CHECK were flagged to prevent their usage in the final product production. The remaining not-flagged warm pixels were identified using sigma clipping with high sigma values to prevent the inadvertent clipping of real data. Subsequently, the DQ extension of the cal files was updated to incorporate these newly flagged pixels. Finally, calwebb_spec3 was executed, configuring the outlier detection in the JWST outlier_detection step with a threshold_percent = 99.9 and a kernel_size = 7 7.

As a result, the final cubes exhibit an improved cleanliness and a significantly higher signal-to-noise ratio than those obtained directly from the MAST.

For the MIRI cubes, a two-point dither pattern in the negative direction was performed for the majority of the sources on the source position. One exception is B1-c, for which a four-point dither pattern was used. For each star-forming region (i.e., Taurus, Perseus, and Serpens), one dedicated background was taken in a two-point dither pattern, except for Taurus, where a single dither was used, in order to properly subtract the telescope background and to better remove detector artifacts. All observations were carried out using the FASTR1 readout mode and used all three MIRI-MRS gratings (A, B, and C), providing the full 4.9–28.6 µm wavelength coverage.

The data were processed through all three stages of the JWST calibration pipeline version 1.13.4 (Bushouse et al. 2024) using the same procedure as described by van Gelder et al. (2024). The reference contexts jwst_1188.pmap of the JWST Calibration Reference Data System (CRDS; Greenfield & Miller 2016) was used. In short, the raw data were processed through the Detercor1Pipeline using the default setting. The dedicated backgrounds were subtracted on the detector level in the Spec2Pipeline. This step also included the correction by the fringe flat for extended sources (Crouzet et al. in prep.) and a residual fringe correction (Kavanagh et al., in prep.). In order to remove any remaining warm and bad pixels from the calibrated detector data, an additional bad-pixel routine was applied using the Vortex Imaging Processing (VIP) package (Christiaens et al. 2024). Last, the final datacubes were constructed separately for each band of each channel using the Spec3Pipeline.

The properties of the Perseus objects studied in this work can be found in Tobin et al. (2016). An overview of the ice bands discussed in this work is presented in Figure 1 for the two sources EDJ183-b and Per-emb35.

2.2 Spectral extractions

Spectra were extracted for both the NIRSpec and MIRI range at central position using a cone diffraction extraction method with a 3λ/D cone aperture. We opted for this aperture size to ensure that a maximum amount of flux was included while avoiding overlap between the binary sources as much as possible. Overlap in the long-wavelength MIRI channels was unavoidable due to the increasing size of the point spread function. The identifying names and spectral extraction coordinates are given for all sources in Table C.1. Binary systems are denoted with ‘a’ and ‘b’.

2.3 Continuum subtraction

Prior to the analysis, the spectra were converted from flux scale $\left({\text{F}}_{\lambda }^{\text{obs}}\right)$ into optical depth scale using the Equation (1) $${\tau }_{\lambda }^{\text{obs}}=-ln\left(\frac{F_{\lambda }^{\text{obs}}}{F_{\lambda }^{\text{cont}}}\right),$$(1)

where ${\text{F}}_{\lambda }^{\text{cont}}$ is the flux of the continuum. In the following sections, we discuss the continuum placement for each vibrational mode.

2.3.1 The 2 µm region

The NIRSpec G235M and G235H modes covering the CO and CO₂ bands in the 2 µm region are available for 9 of the 17 sources. Local continuum subtractions were performed for the ¹²CO overtone mode (2.35 µm) and the ¹²CO₂ combination modes (2.70 and 2.77 µm) by fitting a third-order polynomial to the data, as illustrated in Figure B.1. For the overtone mode, we selected points between 2.30–2.34 and between 2.36–2.38 µm while being mindful of the noise spikes and gas-phase lines in this region. Local continua for the 2.70 µm bands were drawn using a similar method as in Keane et al. (2001).

2.3.2 The 4 µm region

The ¹²CO₂ stretching (4.27 µm), the ¹³CO₂ stretching mode (4.39 µm), the ¹²CO stretching mode (4.67 µm), and the ¹³CO stretching mode (4.78 µm) all have absorption features in the 4 µm spectral region. For many sources, this region is also heavily dominated by the CO gas-phase rotation-vibrational lines that in some cases even overlap with the ¹³CO₂ ice bands (Brunken et al. 2024). For this region, we performed a local continuum subtraction by selecting points between 4.30–4.3,4.41–4.60,4.71–4.80, and 4.91–5 µm. We opted to place a local continuum over the 4.30–4.9 µm range since the shape of the continuum could be better traced when we expanded the wavelength range, especially for sources with a strong gas-line forest (Rubinstein et al. 2023; Federman et al. 2023). In Figure B.2 we show the continuum fittings for the two sources B1a-b (strong CO gas-line forest) and L1527 (weak CO gas-line forest).

We performed separate local continuum subtractions for the ¹²CO₂ 4.27 µm stretching mode because this band is sensitive to scattering effects (Dartois et al. 2022, 2024). Continuum points were placed between 4.17–4.19 and 4.31–4.35 µm following the methods described in Dartois et al. (2022). The results are presented in Figure B.3.

2.3.3 The 15 µm region

Continuum subtraction in the 15 µm region is particularly challenging due to several deep absorption and, in some cases, emission features that collectively alter the shape of the continuum. Two significant features, for instance, are the silicate bending and stretching modes, which have long been suspected of distorting the wing of the ¹²CO₂ band. The two strong absorption features are located at 9.7 and 18 µm, with the latter likely pushing down the long-wavelength wing of the 15.2 µm CO₂ band when the silicate is in absorption. This produces a band profile that is eerily similar to the grain-shape effects observed for the CO₂ stretching mode (4.27 µm), where the short-wavelength wing of the band is raised relative to the long-wavelength wing, creating a negative absorption feature (Dartois et al. 2022, 2024).

Because of this striking resemblance, we briefly considered the possibility that the distortion at 15.2 µm might also be the result of scattering effects, but we argue that in order to produce such a strong scattering feature at these long wavelengths, grain sizes of ~10 µm are required, which is an unlikely scenario in dense protostellar envelopes where the grain sizes are ~0.9 µm (Dartois et al. 2024). Furthermore, for sources in which the silicate band is in emission (Figure B.13), we see the opposite effect: The short-wavelength wing of the CO₂ band is lowered with respect to the long-wavelength wing since the silicate band pushes the long-wavelength wing up instead of down. In addition to the silicate absorption (and emission) band, water has a broad liberation mode that peaks at 13.6 µm and extends over the 15–28 µm region, and this might further contribute to the total absorption.

Taking all of these factors into account, we first placed a local continuum by selecting points between 12–14.5 and 25–28 µm and converted the spectra to optical depth scale (Figures 2, B.9, B.10, and B.11, first column). We subsequently modeled the asymmetric shape of the silicate bending mode with two Gaussian profiles centered at 18.3 µm and 22 µm, and subtracted this from the optical depth spectra (Figures 2, B.9, B.10, and B.11, second column). After subtraction, the bands still bear an extended long-wavelength wing (Figures 2, B.9, B.10, and B.11, third column), but the short-wavelength wing is no longer significantly raised with respect to the long-wavelength wing. This method of modeling the silicate band with Gaussian profiles was also used by Gerakines et al. (1999) and Pontoppidan et al. (2008). For the two sources in which the silicate band is in emission, we fit an additional Gaussian profile centered at 17 µm to simulate the shape of the silicate band.

Finally, we investigated the contribution of the water liberation mode, but found that the wing of this water band is almost identical to the local continuum we placed between the 12 and 28 µm, and it therefore does not additionally contribute to the total absorption after we subtracted the continuum from the spectra.

Fig. 1 Full NIRSpec and MIRI spectra for EDJ183-b and Per-emb35. The ice absorption features are shaded in color. The bottom panels show zoomed-in sections of the spectra.

Fig. 2 Continuum and silicate subtraction for the 12CO2 bending mode in the 15 µm region. The left panel shows the fitting of the continuum (purple line). The middle panel shows the two Gaussian profiles (red and pink lines) and the combined final fit (orange line), which simulates the silicate absorption feature at 18 µm. The third panel shows the subtracted spectra after the silicate band was removed. The continuum subtraction for the remaining sources can be found in Appendix B.2.2.

2.4 Column densities

The column densities for the vibrational modes were derived using $$N=\frac{{{\displaystyle \int }}^{\text{}}\tau \text{d}v}{A},$$(2)

where τ_v is the integrated absorbance under the absorption feature, and A is the corresponding band strength of the vibrational mode.

For the band strengths, we used the values determined by Gerakines et al. (1995) and Gerakines et al. (2005) and corrected by Bouilloud et al. (2015). All band strengths used in this work are listed in Table 1 and are further discussed along with the error analysis on the column densities in Appendix A. The column densities derived in this work can be multiplied by the correction factors given in Table 1 to compare them with the values derived using the band strengths reported in Gerakines et al. (1995) and Gerakines et al. (1995). We note that the recent band strengths Gerakines & Hudson (2015) derived for amorphous CO₂ are similar to the corrected values we used in this work. The change in the band strength due to the ice mixture and the temperature are accounted for in the error analysis. Last, the ¹²C/¹³C ratios were determined for the bands as a whole and not for the individual components.

3 Analysis

We determined the column densities and derived the carbon isotope ratio for the ¹²CO₂ v₁ + v₂ combination mode (2.70 µm), the ¹²CO₂ v₃ stretching mode (4.27 µm), the ¹³CO₂ v₃ stretching mode (4.39 µm), the ¹²CO₂ v₂ bending mode (15.2 µm), the ¹²CO 2-0 overtone mode (2.35 µm), the ¹²CO 1-0 stretching mode (4.67 µm), and the ¹³CO 1-0 stretching mode (4.78 µm). We also detected the ¹²CO₂ 2v₁ + v₂ combination mode (2.77 µm) in a number of sources. The list of sources, their luminosities, distances, and the column densities for each individual vibrational mode are presented in Tables 2 and 3 for CO₂ and CO, respectively.

3.1 CO₂

The 2.70, 4.27 and 15.2 µm vibrational modes of ¹²CO₂ were collectively detected in a number of sources. For these sources, we derived the ¹²CO₂/¹³CO₂ ratios and found values that agree within the reported error bars (Table 2). The largest source of uncertainty in the derived ratios stems from the error on the band strengths, the observational errors are negligible in contrast (~3%). For the error on the band strengths, we accounted for the change in the band strength due to the temperature and due to the ice mixture, more specifically, a water-rich ice mixture, because it is known that ~50% of the CO₂ resides in a water-rich matrix (Pontoppidan et al. 2003; Brunken et al. 2024).

The CO₂-H₂O binary ices affect each vibrational mode of CO₂ differently. The band strengths of the 2.70, 4.27, and 4.39 µm bands, for instance, decrease relative to the band strength of pure CO₂ when CO₂ is mixed with water. This subsequently results in higher ¹²CO₂ and ¹³CO₂ column densities (Gerakines et al. 1995). The band strength of the 15.2 µm band, in contrast, increases by a significant amount relative to the band strength of pure CO₂, resulting in lower ¹²CO₂ column densities. This therefore introduces a systematic error when the ratios from the different vibrational modes are calculated because ratios derived from the 15.2 µm band are lower in general than ratios derived from the other ¹²CO₂ vibrational modes when a band strength for a CO₂-H₂O mixture is used. We used the band strength of pure CO₂ and accounted for this systematic error in the error analysis. For a detailed analysis of these uncertainties, we refer to Appendix A.

3.1.1 The ¹²CO₂ 2.70 and 2.77 µm combination modes

The assignment of the 2.70 µm band has been the topic of multiple studies. Keane et al. (2001), for instance, first fit this feature with laboratory spectra of CO-H₂O (100:1) ices, but later assigned this feature to CO₂ because the fraction of CO in the binary CO-H₂O ices was inconsistent with the amount of CO present in the other bands. The band has also been assigned to the OH-dangling mode of porous water ice because water has multiple weak features in this spectral region (Noble et al. 2024). In our sample, we detect this 2.70 µm feature in 9 out of 17 sources on the short-wavelength side of the OH-dangling mode, situated at 2.75 µm (Figure 3), and we assigned this band to the combination mode of CO₂ for the following reasons.

First, we fit several spectra of pure CO₂ ice and CO₂ in mixtures with CO, H₂O, and CH₃OH (Pontoppidan et al. 2008), and we found an absorption feature at 2.70 µm for every mixture containing CO₂, as illustrated in Figure 4. The dangling modes of pure porous water ice and water ice mixed with other ice species do not sufficiently fit the observed profile of this feature. Additionally, we detect the OH-dangling mode of water ice interacting with other species at 2.75 µm in every source in which the combination mode is present, except for Per-emb35, where it disappears (Figure 3). This source shows clear signs of thermal processing with double-peak profiles at 4.39 µm (¹³CO₂) (Brunken et al. 2024) and 15.2 µm (¹²CO₂) (Pontoppidan et al. 2008) (Figures B. 9, B. 10, B. 4 and B. 5).

OH-dangling modes originate from non-H bonded OH-groups of cold amorphous water ice that is either porous or in mixtures with other species that disrupt the H-bonding network (Keane et al. 2001; Noble et al. 2024), and experimental studies have shown that the strength of these bands decreases with increasing temperature (Isokoski et al. 2014). This is because the elevated temperatures cause the molecules to restructure and form bonds, which in turn will lead to a decrease in the number of free OH-groups and the strength of the OH-dangling modes. A gradual decrease of OH trapping sites due to ice heating was also shown in experimental results by Sandford & Allamandola (1988), where the H₂O-CO shoulder that was observed at 4.64 µm disappeared when the amorphous H₂O ice structure was annealed. Since the 4.39 and 15.2 µm bands of Per-emb35 show clear signs of thermal processing, the absence of the 2.75 µm band in this source is therefore likely the result of the amorphous porous water structure disappearing due to the elevated temperatures in the envelope.

If both the 2.70 and the 2.75 µm bands were indeed water OH-dangling modes, however, then the 2.70 µm feature should also disappear when the ices undergo thermal processing. We instead observe a sharpening of the profile and the appearance of a sharp peak on the short-wavelength side of the band, which we successfully reproduced with the spectrum of pure CO₂ and the alternative spectral decomposition for warm ices described in Brunken et al. (2024) (Figure 4). This behavior is similar to the behavior observed at 4.39 and 15.2 µm when pure CO₂ ice begins to segregate from the ice matrix. This strongly indicates that the 2.70 µm feature therefore cannot primarily be the product of OH-dangling modes. We note that the 2.75 µm OH-dangling mode is detected in the extended sources Ser-S68N-N and Ser-SMM3 (Figure B. 12, where the ices are also processed, but this could be due to the fact that we might be tracing a fraction of cold-ice components in the extended envelopes of these sources.

Finally, while laboratory spectra have shown OH-dangling bonds at 2.70 µm for water ice with high porosity, we note that the method of depositing gas-phase water to form very porous water is inconsistent with the mechanism of interstellar ice formation through atom addition, which produces more compact ices. Therefore, we assigned the 2.70 µm feature to the combination mode of CO₂ and quantified the column density of ¹²CO₂ from this band (Figure 5). The ratios extracted from this band are presented in Table 2.

We note that the column densities we derived from the 2.70 µm band are consistent with the column densities extracted from the other two ¹²CO₂ vibrational modes in both our cold and heated sources, except for Ser-SMM3. This is further evidence that the majority of this band is indeed for the most part the ¹²CO₂ combination mode.

We also observed the 2v₁ + v₂ combination mode at 2.77 µm in Per-emb35 and Ser-SMM3 (Figure B. 12). The band strength of this feature is lower by a factor ~3 than that of the 2.70 µm, however, and it is also highly susceptible to changes in the temperature, even more so than its counterpart at 2.70 µm. Because of these high uncertainties, we refrained from quantifying carbon isotope ratios from this weak feature.

Fig. 3 12CO2 combination mode (2.70 µm) and water OH-dangling mode (2.75 µm) for the sources EDJ183-a (not heated) and Per-emb35 (heated). The feature at 2.75 µm is not detected in Per-emb35, while the 2.70 µm feature is detected in both sources.

Fig. 4 Profile decomposition of the 12CO2 2.70 µm combination mode using selected laboratory spectra based on the analysis performed in Brunken et al. (2024). In the left and middle panel, the black line shows the observed spectrum, and the blue line shows the linear combination of all the five different components at low temperature. The yellow line corresponds to the CO2:CH3OH (10 K) component. The pink line shows the contribution of the broad CO2:H2O (10 K) component. The purple line corresponds to the diluted CO2:CO (25 K) component, and the orange line shows the contribution of CO2 and CO mixed in an equal ratio (15 K). Finally, the dark red line corresponds to the pure CO2 (80 K) component. In the right panel, the black line shows the observed spectrum, and the blue line shows the linear combination of all the three different components at high temperatures. The green line corresponds to the hot CO2:CH3OH (105 K) component. The gray line shows the contribution of the hot CO2:H2O (160 K) component. Finally, the dark red line corresponds to the pure CO2 (80 K) component. All laboratory spectra are taken from Ehrenfreund et al. (1997, 1999); van Broekhuizen et al. (2006) and are publicly available in the Leiden Ice Database (LIDA) (Rocha et al. 2022). Further details of the laboratory spectra can be found in Brunken et al. (2024).

3.1.2 The ¹²CO₂ 4.27 µm stretching mode

In Figure 6 we present the 4.27 µm stretching mode for 6 out of 17 sources. While this band is detected in the remaining sources, they are all highly saturated (τ > 5 ). As a test for using the 4.27 µm band for saturated sources, we fit the laboratory spectrum of CO₂:H₂O 1:10 at 10 K (Ehrenfreund et al. 1999) to the wings of the 4.27 µm band in Per-emb35 using the 15.2 µm bending mode as an anchor.

Overall, the band profiles of the selected six sources appear to be very similar, and only the slope of the short-wavelength wing varied depending on the source. Additionally, the strength of the scattering feature on the short-wavelength wing differed per source. The column densities derived from this band agree well with the column densities derived from the 2.70 µm combination mode for each source.

For the fitted band of Per-emb35, we find a ¹²CO₂/¹³CO₂ ratio of ~80 ± 21 (Figure B. 14). This value is consistent within the reported error bars with the ratios derived from the other two vibrational modes in this same source. Due to the ambiguity of the true band profiles and optical depths of these saturated bands, however, we omitted the band of these sources from this study.

3.1.3 The ¹²CO₂ 15.2 µm bending mode

The 15.2 µm bending mode is detected in all 17 sources, each showing distinctive band profiles (Figures 2, B.9, B.10, and B.11). The sources with a double-peak profile are most notable because this is a telltale sign of ice heating (Ehrenfreund et al. 1998; Pontoppidan et al. 2008).

In Section 2.3.3 we discussed the continuum subtraction method and the extended wing-profile after subtracting the silicate model. We investigated the effect of this extended wing in comparison to the wing profile of Pontoppidan et al. (2008), which extends to ~16 µm instead of ~18 µm, and found a difference in column density of ~9%, which is negligible compared to the uncertainties of the band strengths (see also Appendix A). The column densities derived from the bending mode are presented in Table 2.

Fig. 5 Overview of the 12CO2 combination modes.

3.1.4 The ¹³CO₂ 4.39 µm stretching mode

¹³CO₂ column densities were derived from the 4.39 µm band for all 17 sources (Table 2). Consistent with the 15.2 µm band, we observe a double-peak profile at 4.39 µm for every source with a split peak at 15.2 µm. One hindrance is the CO line forest, which bleeds into the ¹³CO₂ band of several sources in this sample, potentially causing a slight underestimation of the ¹³CO₂ column densities. We briefly investigated the effect of the gaseous CO lines in EDJ183-b by fitting a Gaussian profile to the bottom of the ¹³CO₂ band and quantifying the column density from this curve. We find that the ¹²CO₂/¹³CO₂ ratio decreases by 35% when the integrated optical depth is determined from the Gaussian curve instead of directly from the band. In Figure B.7 we show the Gaussian curves that were ultimately used to determine the integrated optical depths for the ¹³CO₂ bands that are most affected by the gaseous CO lines. For the sources where the ¹³CO₂ features were well isolated from the line forest, the integrated optical depths were directly determined from the ice bands.

Fig. 6 Overview of the 12CO2 4.27 µm stretching mode. The features are shown for sources for which the band is not saturated.

Fig. 7 Source Per-emb55b and the 2.35 µm feature. The laboratory spectrum of pure CO is plotted over the v = 2−0 absorption feature (pink).

3.1.5 ¹²CO₂/¹³CO₂

In Table 2 we compare the column densities and ratios derived from the different CO₂ ice features. For the sources in which the 2.70 µm combination mode, the 4.27 µm stretching mode, and the 15.2 µm bending mode are detected, we find consistent ratios, except for TMC1-W, where the ratio derived from the 4.27 µm stretching mode is ~25% lower than the values extracted from both the 2.70 and the 15.2 µm bands. This deviation still lies within our reported error bars, however. We also find a discrepancy of ~50% between the ratio derived from the 2.70 and the 15.2 µm bands in Ser-SMM3.

For the 2.70 µm band, the 4.27 µm band, and the 15.2 µm band, we find mean ¹²CO₂/¹³CO₂ ratios of 85 ± 23, 76 ± 12, and 97 ± 17, respectively. These are similar to the median values we derived from these vibrational modes (see also Section 4.1). We note that the variation in ratios derived for the different sources and the uncertainty margin on the mean ratios could in fact be due to genuine differences that reflect the different chemical conditions of these systems.

Fig. 8 Overview of the 13CO2 (left), 12CO (middle), and 13CO (right) stretching modes in the 4 µm region. The 12CO gas-phase lines are shown superimposed in several sources.

3.2 CO

3.2.1 The ¹²CO 2.35 µm overtone mode

We detected a weak absorption feature at 2.35 µm above 3σ level for 6 out of the 17 sources (TMC1-E, TMC1-W, Per-emb55-b, EJD183-a, EDJ183-b, and Ser-S68N-N). In previous experimental work, this band was assigned to the ν = 2-0 overtone mode of CO, which peaks at this wavelength in the near-infrared (Fink & Sill 1982; Gerakines et al. 2005). We used experimental data of pure CO to analyze the band profiles of these absorption features and found that the peak position of all sources except that of Per-emb55-b were significantly redshifted (Figure 7).

The sources in which the 2.35 µm features are redshifted also display strong gas-phase CO overtone photospheric absorption lines that peak in the 2 µm spectral region. These gaseous CO overtone modes originate from the central protostellar embryo and appear in absorption due to the strong thermal gradient of the gas surrounding the protostar (Le Gouellec et al. 2024). These gas-phase absorption lines include the ν = 4-2 line at 2.3525 µm, which overlaps the 2.35 µm CO ice absorption band and causes the redshifted ice bands we observed in most of these sources.

In order to disentangle and subtract the gaseous CO overtone modes, the 2 µm hot dust emission was first modeled with a blackbody emission. The photosphere was modeled using the BT-Settl grid of photospheric models (Allard et al. 2012). The Starfish modeling framework from Czekala et al. (2015) was used (see also Figure B.15). Further details of these modeling procedures are presented in Le Gouellec et al. (in prep.). The findings show that the models converged well in TMC1-E, TMC1-W, EDJ183-a, EDJ183-b, and Ser-S68N-N, and after subtraction, no residual ice features were observed in these sources. A faint photosphere is detected in Per-emb55-a, and a weak absorption feature is visible at 2.35 µm. This feature, however, is faint and hard to separate from the noise spikes in the spectrum. In Per-emb55-b, no photosphere was detected, and we can conclude that the 2.35 µm feature is indeed the CO overtone ice absorption band.

Because of the gas contamination, we discarded the other five sources in which these gas signatures were observed and only extracted the ¹²CO column density from the band observed in Per-emb55-b (Figure 7). The column density we derived from this band is lower by a factor ~1.6 than the column density derived from the 4.67 µm band in this same source. One possible explanation for the discrepancy between these two vibrational modes is the band strength of this relatively weak absorption feature.

Fig. 9 Zoom-in of the 13CO absorption features. The Gaussian curves (pink) fit to the bands illustrate the bounds we used to determine the integrated optical depths.

3.2.2 The ¹²CO 4.67 µm stretching mode

The bands of the ¹²CO 4.67 µm symmetric stretching mode are presented in Figures 8, B.4, B.5, and B.6, and the column densities are listed in Table 3. We find variations of more than an order of magnitude between the different sources and ascribe this to the volatility of CO. Because of its lower desorption temperature, CO ice is highly sensitive to the temperature structure of the protostellar envelope, which in turn can significantly impact the CO ice budget and cause variations depending on the source.

3.2.3 The ¹³CO 4.78 µm stretching mode

The ¹³CO 4.78 µm stretching mode is observed in 14 out of 17 sources, only 2 of which are well isolated from the CO gasphase rotation-vibrational lines (B1-b and B1-c). The strong line forest that dominates most of the sources in this sample poses a challenge to the study of this weak feature because the actual shape and optical depth of this band are likely affected by the gas lines. We briefly investigated the effect of the emission lines by fitting Gaussian profiles and subtracting the two gaseous CO lines closest to the bands, and we find that the ¹³CO column density increases by a factor 1.4 compared to the nonsubtracted band. This is still within our reported uncertainties. We note, however, that the gas-subtracted ¹³CO ice bands are significantly broader than the nonsubtracted bands and are also much broader than the band of the laboratory spectrum of pure CO, which could be the main contributor of this band (Boogert et al. 2002a; Pontoppidan et al. 2003). On the other hand, a broad band profile is possible when the ¹³CO ice resides in a polar CH₃OH-rich ice matrix (Penteado et al. 2015). A careful modeling and removal of the gaseous CO lines is needed in order to properly isolate the ice feature and quantify the ¹³CO ice column density. This CO gas modeling is beyond the scope of this paper, however. In Figures 9 and B.8, we show the Gaussian curves that were fit to the ¹³CO features in order to determine the integrated optical depths.

3.2.4 ¹²CO/¹³CO

The carbon isotope ratios derived from the CO ice absorption bands are presented in Table 3. In general, we find high ¹²CO/¹³CO ratios. The sources B1-a-S, Per-emb22, Per-emb33, Ser-S68N-N, and Per-emb27 in particular have very high ratios, all of which also have very weak ¹³CO ice absorption features. Moreover, the ¹³CO ice bands of some of these sources are buried in the gas line forest, which could be the reason for some of these extreme ratios. We derived a sub-ISM ratio for source B1-c, but we note that the spectrum of this source was taken at medium resolution unlike the other sources in this sample. Consequently, the narrow ¹³CO band might not be fully resolved. Additionally, the gaseous CO lines in this source are in absorption and could be blended with the unresolved ice feature further contributing to its width.

Finally, we derived a ¹²CO/¹³C ratio of 96 ± 14 from 2.35 µm overtone mode in Per-emb55-a, which is lower by a factor ~1.6 than the ratio of 152 ± 23 derived in this same source from the 4.67 µm band. For the 4.67 µm band, we find a mean ¹²CO/¹³CO ratio of 165 ± 52 for the sample as whole when we exclude the three extreme outliers Per-emb22, Per-emb33, and Per-emb27. Due to the fluctuating values derived from the CO vibrational modes, we recommend referring to the ratios derived from the CO₂ vibrational modes rather than the CO ice bands.

Fig. 10 Solid 12C/13C for various astronomical environments. The filled red, gray, and white stars represent the ratios derived from the 12CO2 15.2, 4.27, and 2.70 µm bands, respectively. We also include molecular clouds (Ice Age) (McClure et al. 2023) and a massive YSO (Infrared Space Observatory). (Boogert et al. 2000). The dashed purple line shows the solar abundance, and the dot-dashed gray line shows the local ISM ratio (Boogert et al. 2000). The small solid error bars of the JOYS+ data points illustrate when the band strengths were excluded from the error analysis. The faint and generally larger error bars of the JOYS+ data points illustrate the uncertainty, including the error on the band strength. Further details of the error analysis can be found in Appendix A.

3.3 Limitations and future work

As shown in the previous sections, we examined the carbon isotope ratios for multiple CO and CO₂ vibrational modes for a large sample of low-mass sources. We found consistent values for the most part, especially between the ratios derived from the CO₂ ice bands. One source of uncertainty are the band strengths, especially for the weaker ¹²CO overtone mode. The 4.40– 4.5 µm spectral region is also heavily dominated by gaseous CO rotation-vibrational lines, which poses a challenge when quantifying ¹³CO column densities from the absorption feature at 4.78 µm. Therefore, a careful modeling and subtraction of these gas lines is needed to properly isolate it in the ¹³CO ice features and determine ¹³CO column densities. Alternatively, spectral extractions at locations where the gas-phase lines are weak or absent could also aid in the analysis of the ¹³CO ice features. Finally, future work should also focus on finding similarities or variations between close binaries, and a spatial ice mapping should be performed since this feat is now possible with the JWST.

4 Discussion

The carbon isotope ratio is a crucial link between the evolutionary phases of star and planet formation because of its sensitivity to the chemical and physical conditions that reign during each stage. In the following section, we examine how the results presented in this work relate to the trends observed in different astronomical environments both in the solid state and in the gas phase.

4.1 Protostars: Solid state

4.1.1 CO₂

The ¹²C/¹³C ratios we derived from solid CO₂ for the JOYS+ sample are presented in Figure 10, along with the isotope ratios derived in the solid state for the Infrared Space Observatory (ISO) study of young stellar objects (Boogert et al. 2000) and two highly extincted lines of sight (McClure et al. 2023) in the Cha I molecular cloud. We note that although previous studies used the uncorrected Gerakines et al. (1995) band strengths, the corrections do no affect the relative band strengths, and consequently, do not affect the ¹²C/¹³C ratios. For this reason, the ratios determined in these previous studies can be directly compared to the ratios derived in this work. The light error bars of the JOYS+ data points illustrate the uncertainties, including the error on the band strengths, and the solid error bars illustrate the uncertainties, excluding the error on the band strengths. The errors in the band strengths account for the change in the band strength due to the temperature and ice mixture. The small error bars when the uncertainty on the band strengths are excluded demonstrate the great capabilities of the JWST and the high quality of these data. Previous studies generally did not include a detailed analysis of the uncertainty in the band strengths in their error bars.

The results show that the ratios obtained from the 15.2 µm bands are clustered slightly above the ISM ratio of 68 (Boogert et al. 2000; Milam et al. 2005). This suggests that the ices in our sources are slightly deficient in ¹³C in comparison to the local ISM. When we detect all three ¹²CO₂ vibrational modes, we find consistent column densities and ratios (Figure 10). This is further evidence that the 2.70 µm feature primarily consists of the ¹²CO₂ combination mode. These consistent values also show that the 4.27 µm stretching mode can be used to determine ¹²CO₂ column densities when the other ¹²CO₂ bands are not available and when the stretching mode is not saturated.

When compared to other solid-state studies of carbon isotope ratios, the ¹²CO₂/¹³CO₂ ratios measured in our sample are fairly consistent with the those found in the ISO-SWS high-mass sources by Boogert et al. (2000) (ranging 52–113). Similar to the JOYS+ sources, the authors reported values that are scattered around the local ISM ratio, but they also reported sub-ISM values as low as 52. These sub-ISM ratios still agree with our values within the reported error bars. Figure 11 illustrates the distribution of the ¹²CO₂/¹³CO₂ ratio for the JOYS+ sample and the ISO high-mass sources (Boogert et al. 2000). While a large fraction of the ISO MYSOs also displays elevated ratios with respect to the ISM standard, the ratios of the JOYS+ protostars are in general higher, with median values of 81, 75, and 100 for the 2.70 µm combination mode, the 4.27 µm stretching mode, and the 15.2 µm bending mode, respectively. These are similar to the mean values of 85 ± 23, 76 ± 12, and 97 ± 17 derived for these bands, respectively. The agreement between the ratios derived from the different CO₂ vibrational modes of the JOYS+ protostars again shows the great capability of the JWST to accurately probe different lines of sight and observe the weak ice features in the 2 µm region, the generally saturated strong stretching mode at 4.27 µm, and the bending mode at 15.2 µm.

Our values are also consistent if slightly higher than the values reported in McClure et al. (2023) for the two lines of sight in the Cha I molecular cloud (¹²CO₂/¹³CO₂ ~ 69–87). We note, however, that the ¹²CO₂ column densities derived for Cha I were extracted from the 4.27 µm bands, which are saturated in both lines of sight. Consequently, the ¹²CO₂ column density could be somewhat underestimated in these sources. Furthermore, the uncertainties reported in McClure et al. (2023) did not include errors on the band strengths. Nonetheless, the consistency between the ratios derived from the protostellar stage and the dense-cloud stage paint the harmonious picture that the CO₂ was likely inherited from its parent molecular cloud.

Fig. 11 Distribution of the 12CO2/13CO2 isotope ratio. The filled orange bars show the distribution of the carbon isotope ratios determined for the JOYS+ sample. These ratios were extracted from the 15.2 µm bending mode. The striped purple bars show the distribution of the carbon isotope ratio derived for the ISO high-mass sources (Boogert et al. 2000).

4.1.2 CO

The ¹²CO/¹³CO ratios derived from solid CO for the JOYS+ sample are shown in Figure 12, along with the isotope ratios derived from various studies both in the gas and solid state. Similar to CO₂, the corrections on the CO band strengths (Gerakines et al. 1995) have no effect on the ¹²C/¹³C ratio and values derived in previous solid-state studies can therefore be directly compared to the findings in this work. Overall, our results show that the isotope ratios obtained from the CO 4.67 µm stretching mode are considerably higher than the local ISM value. The ratios are also higher than the ratios previously derived for one MYSO (Boogert et al. 2002a) and one LYSO (Pontoppidan et al. 2003), both of which are very close to the ISM standard. This discrepancy likely arises because our ¹³CO column density is underestimated due to the weak ¹³CO ice absorption feature and the strong CO gas lineforest that dominates this spectral region. This is particularly the case for the three extreme outliers Per-emb22, Per-emb33, and Per-emb27 (shown only in Table 3), both of which exhibit strong gas-phase lines and weak ¹³CO ice bands (Figure B.5).

The ratio derived from the 2.35 µm feature in Per-emb55-b is also higher than the local ISM standard. When compared to the ratios derived in molecular clouds, however, our values are consistent with those found by McClure et al. (2023), who also reported super-ISM ¹²CO/¹³CO ratios (99–184) (Figure 12).

4.2 CO versus CO₂ isotope ratios

The ratios derived from CO and CO₂ from the JOYS+ sample are presented in Figure 12, along with gas-phase and solid-state ratios derived from various studies. Normally, the carbon isotope ratio of solid CO and CO₂ is expected to be similar because the formation route of CO₂ is tied to that of CO through the equation (3) (Ioppolo et al. 2011) $$\text{CO}+\text{OH}->\text{HO}-\text{CO}->{\text{CO}}_2+\text{H}.$$(3)

In the ices of our sources, however, the ¹²CO/¹³CO ratios are higher in general than the ¹²CO₂/¹³CO₂ ratios derived from the 2.70, 4.27, and 15.2 µm bands, although they still agree within the error bars for a number of our sources. As noted in Section 4.1.2, a reason for this deviation could be that the CO ratios are less accurate because the ¹³CO ice feature is weak. Overall, we consider the CO₂ -derived isotope ratios to be more reliable given the multiple methods we used to quantify this ratio.

4.3 Protostars: Gas phase

In the gas-phase, Smith et al. (2015) reported significant ¹²CO/¹³CO variability for several YSOs, with an overall trend of elevated ratios with respect to the local ISM standard (Figure 12). One possible scenario they proposed is gas-ice partitioning (Section 4.4.3), with an additional mechanism for removing ¹³CO from the solid state through the formation of larger complex organic molecules (COMs) from ¹³C-rich CO ices. Consequently, the COMs will be enriched ¹³C and exhibit low carbon isotope ratios, such as the sub-ISM values found in COMs by Jørgensen et al. (2016, 2018). We note, however, that the gasphase values reported by Smith et al. (2015) are similar to those derived in this study, which suggests that the ices might have sublimated from the grains, delivering the carbon isotope ratio to the gas during the pre-stellar stage.

The quantifying of gas-phase and solid-state isotope ratios in the same source was in fact proposed by Charnley et al. (2004) as a way of tracing the origin of organic material. Therefore, observations of gas-phase CO isotopologs for the JOYS+ protostars will be crucial for enabling this method of posterior isotopic labeling and linking the solid-state chemistry with the gas-phase chemistry in these systems. In the following section, we study the different mechanisms that can affect the carbon isotope ratio in both the solid and gas phase in the early stages of star and planet formation.

Fig. 12 12C/13C for various astronomical environments. The filled red, gray, and white stars represent the ratios derived from the 12CO2 15.2, 4.27, and 2.70 µm bands, respectively. The closed and open yellow stars represent the ratios derived from the 12CO 4.67 and 2.35 µm bands, respectively. The error bars on the JOYS+ data points include the error on the band strengths, which is the main contributor to these uncertainties. The ratios were taken from the following works: Ice age (McClure et al. 2023), (Boogert et al. 2000), KECK-NIRSpec (Boogert et al. 2002a), VLT (Pontoppidan et al. 2003), VLT-CRIRES (Smith et al. 2015), PILS (Jørgensen et al. 2016, 2018), Disks (Yoshida et al. 2022; Bergin et al. 2024), Exo-planets (Line et al. 2021; Zhang et al. 2021a,b; Gandhi et al. 2023; de Regt et al. 2024), Comets (Bockelée-Morvan et al. 2015; Hässig et al. 2017). The dashed purple line shows the solar abundance, and the dot-dashed gray line shows the local ISM ratio (Boogert et al. 2000). Further details of the error analysis can be found in Appendix A.

4.4 Fractionation processes

4.4.1 Isotope-exchange reactions

Isotope-exchange reactions are exothermic reactions that are efficient at low temperatures and that prompt preferential incorporation of ¹³C⁺ in CO (Equation (4)) (Watson 1976; Langer et al. 1984; Langer & Penzias 1993). As a result, molecular ¹³CO is enhanced with respect to ¹²CO in the gas phase, and this low ¹²CO/¹³CO ratio is subsequently conveyed to the solid state upon CO freeze-out. ¹²C⁺, in contrast, is enhanced in the gas phase, leading to high ¹²C⁺/¹³C⁺ gas-phase ratios. $${}^{13}{\text{C}}^{+}{+}^{12}\text{CO}{\leftrightarrow }^{12}{\text{C}}^{+}{+}^{13}\text{CO}+\text{ΔE}=34 \text{K}.$$(4)

4.4.2 Selective photodissociation

Selective photodissociation of ¹³CO is a destruction mechanism that occurs because self-shielding of the less abundant ¹³CO takes place at higher extinctions in comparison to ¹²CO (van Dishoeck & Black 1988; Visser et al. 2009). This leaves a large portion of the ¹³CO gas vulnerable to ultraviolet radiation. As a result, ¹³CO is destroyed, which increases the gaseous ¹²CO/¹³CO ratio and decreases the ¹²C⁺/¹³C⁺ ratios as the molecules dissociate and the gas is enriched with ¹³C⁺ and O atoms. This counters the effect of isotope selective photodissociation (Section 4.4.1). This high molecular ¹²CO/¹³CO ratio translates into the solid state upon CO freeze-out, producing ices that are depleted in ¹³C. Although we expect that under dark cloud conditions, most of the radiation is attenuated by the high-density environment, thus limiting the effect of selective photodissociation, studies by Furuya & Aikawa (2018) have shown that turbulence can transfer selective photodissociated material to the dark inner regions of clouds.

4.4.3 Gas-ice partitioning

A third mass-dependent mechanism was proposed by Smith et al. (2015), who reported that a small difference in the binding energies of ¹²CO and the slightly heavier ¹³CO (∆E_bind~10 K) resulted in ¹²CO desorbing at a slightly lower temperature from the grains. Consequently, the ¹²CO/¹³CO in the gas phase is enhanced while the ices become enriched with ¹³CO, leading to low carbon isotope ratios in the solid state. It is worth noting, however, that gas-ice partitioning is only effective during a very narrow temperature window (Smith et al. 2015).

Of the mechanism discussed above, selective photodissociation could be causing the ¹³C deficiency that is observed in this sample, although its efficiency under these conditions is open to question.

4.5 Evolution of the ¹²C/¹³C ratio

As the stellar system matures, different process begin to dominate. They slowly erode the initial carbon isotope ratio of the pristine material. In the following sections, we discuss the evolution of the carbon isotope ratio during the later epochs of star formation.

4.5.1 Protoplanetary disks

Studies of carbon isotope ratios in disks revealed that at this stage, the conditions in the system are irrevocably changed and that the initial carbon isotope imprint can be effectively erased in certain regions of the disk. Yoshida et al. (2022) and Bergin et al. (2024), for instance, found a dichotomy in the TW Hya disk, where two different carbon isotope ratios were measured for CO and C₂H (21 ± 5 and 65 ± 20, respectively). In the same disk, Hily-Blant et al. (2019) found a ratio of 86 ± 4 from HCN rotational lines.

The variability of the carbon isotope ratio in protoplane- tary disks arises because the ratio is contingent on the location in which the molecules were formed. This is because different processes dominate in different regions of the disk. Photodissociation, for instance, plays a significant role in the upper layers of the disk, where the high ratio derived from C₂H is found, while isotope-exchange reactions are more efficient in the midplane and at larger radii, where the temperatures are low and where the low value for CO is measured (Willacy & Woods 2009; Visser et al. 2018; Bergin et al. 2024).

4.5.2 Exoplanets

The large variability in the carbon isotope abundance is also observed in the atmospheres of exoplanets, where the formation mechanism of the planet plays a crucial role in its carbon isotope reservoir. Recent observational efforts by Zhang et al. (2021a); Line et al. (2021); Gandhi et al. (2023), for instance, revealed a significant enrichment of ¹³CO in the atmospheres of hot Jupiters (31±17, 10.2–42.6 with a confidence of 68%, and 62 ±2, respectively). This ¹³C enrichment is likely because the planet was formed beyond the CO snowline and accreted most of its material from ¹³C-rich ices. These ¹³C-rich ices are likely the result of fractionation processes such as isotope-exchange reactions and possibly gas-ice partitioning, which must have occurred after the protostellar stage.

Zhang et al. (2021b); de Regt et al. (2024), in contrast, measured very high ¹²C/¹³C ratios (97 ± 15 and 184 ± 61, respectively) for brown dwarfs, which indicates a formation route similar to that of stars, where the object forms from cloud fragmentation and gravitational collapse and inherits most of its gaseous material from a parent cloud that is also deficient in ¹³C and ¹³C-ices, as we found here. This shows that the carbon isotope ratio is a potentially powerful tool for tracing the formation pathways of extrasolar planets.

4.5.3 Comets and meteorites

The chemical complexity of the infant planetary system is at last fossilized in the remains of its building blocks, the comets. Observations of C₂, CN, and HCN in numerous comets, for instance, showed that the carbon isotope ratios of Solar System objects are all similar to the solar abundance (~89), with little variation between the ratios (Bockelée-Morvan et al. 2015). Similarly, Alexander et al. (2007) reported ratios consistent with the solar abundance from a large study of 75 chondrites. Finally, Hässig et al. (2017) also derived a ratio of 84 ± 4 from CO₂ measurements in the coma of comet 67P/Churyumov-Gerasimenko. It is worth noting that the carbon isotope ratio found in comets and chondrites resembles the values derived in the protostel- lar stage. This could be an indication that while fractionation processes erode the initial carbon isotope imprint during the later stages, a fraction of the pristine material still survives this journey and is eventually incorporated into the planetary system.

5 Conclusions

We have analyzed JWST NIRSpec and MIRI data of 17 Class 0/1 low-mass protostars and determined the ¹²CO, ¹³CO ¹²CO₂, and ¹³CO₂ column densities and ¹²C/¹³C isotope ratios from the ¹²CO₂ v1 + v₂ and 2v1 + v₂ combination modes (2.70 µm and 2.77 µm), the ¹²CO₂ v₃ stretching mode (4.27 µm), the ¹³CO₂ v₃ stretching mode (4.39 µm), the ¹²CO₂ v₂ bending mode (15.2 µm), the ¹²CO 2-0 overtone mode (2.35 µm), the ¹²CO 1-0 stretching mode (4.67 µm), and the ¹³CO 1-0 stretching mode (4.78 µm). The most significant finding is that the ratios are consistent if slightly higher than the local ISM value. These results show that the ices leave the protostellar stage with high ¹²C/¹³C ratios, after which a series of fractionation processes during the later stages could modify the initial isotope abundance.

• The absorption feature observed at 2.70 µm in 9 out 17 sources is assigned to the combination mode of ¹²CO₂. We also detect the 2.77 µm combination mode of ¹²CO₂ in two sources.

• We find consistent ¹²CO₂/¹³CO₂ ratios between the 4.27 µm, 15.2 µm, and the 2.70 µm bands. Our values are slightly elevated with respect to the standard ISM value, but they are consistent with carbon isotope ratios observed in other protostars and in dark molecular clouds. We observe variations in the ¹²CO₂/¹³CO₂ ratio from one source to the next, which might indicate genuine differences in the chemical complexity of their protostellar envelopes.

• We report a detection of the 2.35 µm ν 2-0 ¹²CO overtone mode in one source.

• The ¹²CO/¹³CO ratios derived from the 4.67 µm band are higher than the local ISM ratio and the ¹²CO₂/¹³CO₂ ratios. The ratio derived from the overtone mode is lower by a factor ~1.6 than the ratio extracted from the 4.67 µm band in the same source.

Our findings show that fractionation processes begin to dominate after the protostellar stage. Future work should be focused on ice spatial mapping of close binary pairs and extended protostellar envelopes, and on obtaining observations of the gas-phase CO isotopologs of the JOYS+ sample in order to compare solid-state and gas-phase carbon isotope ratios.

Data availability

Appendices B and C are available on Zenodo.

Acknowledgements

Astrochemistry in Leiden is supported by the Netherlands Research School for Astronomy (NOVA), by funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No. 101019751 MOLDISK), and by the Dutch Research Council (NWO) grant 618.000.001. Support by the Danish National Research Foundation through the Center of Excellence “InterCat” (Grant agreement no.: DNRF150) is also acknowledged. This work is based on observations made with the NASA/ESA/CSA James Webb Space Telescope. The data were obtained from the Mikulski Archive for Space Telescopes at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5- 03127 for JWST. The following National and International Funding Agencies funded and supported the MIRI development: NASA; ESA; Belgian Science Policy Office (BELSPO); Centre Nationale d’Etudes Spatiales (CNES); Danish National Space Centre; Deutsches Zentrum fur Luft- und Raumfahrt (DLR); Enterprise Ireland; Ministerio De Economiá y Competividad; Netherlands Research School for Astronomy (NOVA); Netherlands Organisation for Scientific Research (NWO); Sci- ence and Technology Facilities Council; Swiss Space Office; Swedish National Space Agency; and UK Space Agency. N.B. thanks K. Pontoppidan for his contribution to the NIRSpec data reduction and helpful discussions. P.J.K. acknowledges support from the Science Foundation Ireland/Irish Research Council Pathway program under grant No. 21/PATH-S/9360. L.M. acknowledges the financial support of DAE and DST-SERB research grant (MTR/2021/000864) of the Government of India.

Appendix A Error analysis

The uncertainties on the derived isotope ratios are mostly sensitive to the uncertainties in the band strengths. Minor contributors are also the error on the placement of the continuum and the noise level. To ultimately set constraints on the ratios we must take all these factors into account when calculating column densities.

The ratio for the CO and CO₂ bands are calculated as follows: $$ratio=\frac{N_{{}^{12}}{}_{C{O}_2}}{N_{13}{}_{C{O}_2}}=\frac{\frac{{\displaystyle \int \tau }d{v}_{12}{}_{C{O}_2}}{A_{12}{}_{C{O}_2}}}{\frac{{\displaystyle \int \tau }d{v}_{13}{}_{C{O}_2}}{A_{13}{}_{C{O}_2}}}=\frac{{\displaystyle \int \tau }d{v}_{12}{}_{C{O}_2}}{{\displaystyle \int \tau }d{v}_{13}{}_{C{O}_2}}\frac{A_{{}^{12}}{}_{C{O}_2}}{A_{{}^{13}C{O}_2}},$$(A.1)

where $N_{12}{}_{C{O}_2}$ and ${\text{N}}_{13}{}_{C{O}_2}$ are the column densities of ¹²CO₂ and ¹³CO₂ respectively, ${{\displaystyle \int }}^{\text{}}\tau d{v}_{12}{}_{C{O}_2}$ and ${{\displaystyle \int }}^{\text{}}\tau d{v}_{13}{}_{C{O}_2}$ are the integrated optical depths and $A_{12}{}_{C{O}_2}$ and $A_{13}{}_{C{O}_2}$ are the band strengths of the ¹²CO₂ and ¹³CO₂ vibrational modes respectively.

From equation A.1, the uncertainty on the ratio is determined following the error propagation rules: $$\frac{{\sigma }_{ratio}}{ratio}=\sqrt{{\left(\frac{{\sigma }_{A_{12}{}_{C{O}_2}}}{A_{{}^{12}C{O}_2}}\right)}^2+{\left(\frac{{\sigma }_{A_{13C{O}_2}}}{A_{{}^{13}C{O}_2}}\right)}^2+{\left(\frac{{\sigma }_{\tau d{v}_{12}{}_{C{O}_2}}}{\tau d{v}_{v_{12C{O}_2}}}\right)}^2+{\left(\frac{{\sigma }_{\tau d{v}_{13C{O}_2}}}{\tau d{v}_{13}{}_{C{O}_2}}\right)}^2},$$(A.2)

where ${\sigma }_{\tau d{v}_{12}{}_{{}_{C{O}_2}}}$and ${\sigma }_{\tau d{v}_{13}{}_{{}_{C{O}_2}}}$are the uncertainties on the integrated optical depth for ¹²CO₂ and ¹³CO₂ respectively.

The uncertainty on the integrated optical depth is influenced by the uncertainty in the placement of the continuum and the noise level and is determined as follows: $${\left(\frac{{\sigma }_{\tau d{v}_{12}{}_{C{O}_2}}}{\tau d{v}_{12C{O}_2}}\right)}_{12C{O}_2}^2=\sqrt{{\left(\frac{{\sigma }_{int.OD,cont}}{int\text{.}OD\text{, }cont}\right)}_{{}^{12}C{O}_2}^2+{\left(\frac{{\sigma }_{noise}}{int\text{.}OD\text{, }noise}\right)}_{{}^{12}C{O}_2}^2},$$(A.3)

where σ_int.OD,cont is the uncertainty due to the placement of the continuum, int.OD,cont is the integrated optical depth with the continuum of choice, σ_noise is the uncertainty due to the noise level and int.OD,noise is the integrated optical depth. These parameters are calculated in a similar manner for ¹³CO₂.

A.1 Continuum placement

We determine the uncertainty on the placement of the continuum by fitting the bands with two different baselines and using the difference in integrated optical depth as the uncertainty.

A.2 Noise level

The uncertainty due to the noise level is determined by calculating the rms in line free regions on flux scale (σ_{f, λ}) and propagating this to optical depth scale (σ_ƒ,τ) following the error propagation rules for logarithms: $${\sigma }_{\tau ,lambda}=\frac{{\sigma }_{f,\lambda }}{f_{cont,\lambda }},$$(A.4)

since the observed fluxes $F_{\lambda }^{obs}$ are converted to optical depth scale following equation (A.5), $${\tau }_{\lambda }^{obs}=-ln\left(\frac{F_{\lambda }^{obs}}{F_{\lambda }^{cont}}\right),$$(A.5)

where $F_{\lambda }^{cont}$ is the flux of the continuum.

We take F_cont,λ as constant since the uncertainty for that parameter is determined separately.

We use σ_ƒ,λ to generate a noise level following a normal distribution and resample the spectrum with this added noise. We calculate the integrated optical depth for these resampled spectra which produces a sampling distribution that we can fit a Gaussian curve to (Figure A.1). From the Gaussian curve we can extract σ_noise and int.OD,noise. This approach is similar to bootstrapping statistics and the contribution of this uncertainty is very small.

Fig. A.1 Sampling distrubition of the integrated optical depth. The dashed shows the Gaussian curve fitted to the distribution from which we extract σnoise and int.OD,noise.

A.3 Band strengths

The column density for each vibrational mode is determined as follows: $$N=\frac{{\displaystyle \int \tau }dv}{A},$$(A.6)

where A is the band strengths for of the corresponding vibrational mode.

As mentioned, we use the Gerakines et al. (1995) and Gerakines et al. (2005) band strength corrected by Bouilloud et al. (2015) to calculate the column densities.

The contributing factors to the uncertainty on the band strength for $A_{C{O}_2}$ are the placement of the continua (since we are considering relative errors), the uncertainty due to changes in the temperature of the ice and the uncertainty due to the ice mixtures. Since we know that the water component is a has a large contribution to the total band we took the uncertainties determined for the ice mixture CO₂:H₂O 1:24 binary ices. All values used in this work are reported in Table 1.

For the error due to the continuum we assume an uncertainty of 5% since Gerakines et al. (1995) only reports uncertainties >5% and no uncertainty is reported for this band strength. The uncertainty due to change in temperature is ~ 10% (Gerakines et al. 1995). Uncertainties for $A_{C{O}_2-H_{2}O}$ are reported in Gerakines et al. (1995).

Taking all these assumption into consideration we determine the uncertainties on the band strengths as follows: $$\sqrt{{\sigma }_{A,cont}^2+{\sigma }_{A,temp}^2+{\sigma }_{A,H_{2}O-C{O}_2}^2}$$(A.7)

The uncertainties on the CO column densities are determined in a similar manner with following exceptions:

1. We do not have mixture component for the error on the band strengths since we expect CO to be mainly in its pure form Pontoppidan et al. (2003)

2. We do not have a temperature component for the error on the band strengths since CO disorbs above 40 K.

3. We take a 10% error on the band strength for both the ¹²CO bands and the ¹³CO to the placement of the continuum in laboratory spectra.

4. The uncertainty on the integrated optical depth of the ¹³CO ice feature between gas-subtracted spectra and non-gas subtracted spectra is ~ 50 %. This is not included in the error bars of the figures.

References

All Tables

All Figures

	Fig. 1 Full NIRSpec and MIRI spectra for EDJ183-b and Per-emb35. The ice absorption features are shaded in color. The bottom panels show zoomed-in sections of the spectra.
In the text

Fig. 2 Continuum and silicate subtraction for the 12CO2 bending mode in the 15 µm region. The left panel shows the fitting of the continuum (purple line). The middle panel shows the two Gaussian profiles (red and pink lines) and the combined final fit (orange line), which simulates the silicate absorption feature at 18 µm. The third panel shows the subtracted spectra after the silicate band was removed. The continuum subtraction for the remaining sources can be found in Appendix B.2.2.

In the text

	Fig. 3 12CO2 combination mode (2.70 µm) and water OH-dangling mode (2.75 µm) for the sources EDJ183-a (not heated) and Per-emb35 (heated). The feature at 2.75 µm is not detected in Per-emb35, while the 2.70 µm feature is detected in both sources.
In the text

Fig. 4 Profile decomposition of the 12CO2 2.70 µm combination mode using selected laboratory spectra based on the analysis performed in Brunken et al. (2024). In the left and middle panel, the black line shows the observed spectrum, and the blue line shows the linear combination of all the five different components at low temperature. The yellow line corresponds to the CO2:CH3OH (10 K) component. The pink line shows the contribution of the broad CO2:H2O (10 K) component. The purple line corresponds to the diluted CO2:CO (25 K) component, and the orange line shows the contribution of CO2 and CO mixed in an equal ratio (15 K). Finally, the dark red line corresponds to the pure CO2 (80 K) component. In the right panel, the black line shows the observed spectrum, and the blue line shows the linear combination of all the three different components at high temperatures. The green line corresponds to the hot CO2:CH3OH (105 K) component. The gray line shows the contribution of the hot CO2:H2O (160 K) component. Finally, the dark red line corresponds to the pure CO2 (80 K) component. All laboratory spectra are taken from Ehrenfreund et al. (1997, 1999); van Broekhuizen et al. (2006) and are publicly available in the Leiden Ice Database (LIDA) (Rocha et al. 2022). Further details of the laboratory spectra can be found in Brunken et al. (2024).

In the text

	Fig. 5 Overview of the 12CO2 combination modes.
In the text

	Fig. 6 Overview of the 12CO2 4.27 µm stretching mode. The features are shown for sources for which the band is not saturated.
In the text

	Fig. 7 Source Per-emb55b and the 2.35 µm feature. The laboratory spectrum of pure CO is plotted over the v = 2−0 absorption feature (pink).
In the text

	Fig. 8 Overview of the 13CO2 (left), 12CO (middle), and 13CO (right) stretching modes in the 4 µm region. The 12CO gas-phase lines are shown superimposed in several sources.
In the text

	Fig. 9 Zoom-in of the 13CO absorption features. The Gaussian curves (pink) fit to the bands illustrate the bounds we used to determine the integrated optical depths.
In the text

Fig. 10 Solid 12C/13C for various astronomical environments. The filled red, gray, and white stars represent the ratios derived from the 12CO2 15.2, 4.27, and 2.70 µm bands, respectively. We also include molecular clouds (Ice Age) (McClure et al. 2023) and a massive YSO (Infrared Space Observatory). (Boogert et al. 2000). The dashed purple line shows the solar abundance, and the dot-dashed gray line shows the local ISM ratio (Boogert et al. 2000). The small solid error bars of the JOYS+ data points illustrate when the band strengths were excluded from the error analysis. The faint and generally larger error bars of the JOYS+ data points illustrate the uncertainty, including the error on the band strength. Further details of the error analysis can be found in Appendix A.

In the text

	Fig. 11 Distribution of the 12CO2/13CO2 isotope ratio. The filled orange bars show the distribution of the carbon isotope ratios determined for the JOYS+ sample. These ratios were extracted from the 15.2 µm bending mode. The striped purple bars show the distribution of the carbon isotope ratio derived for the ISO high-mass sources (Boogert et al. 2000).
In the text

Fig. 12 12C/13C for various astronomical environments. The filled red, gray, and white stars represent the ratios derived from the 12CO2 15.2, 4.27, and 2.70 µm bands, respectively. The closed and open yellow stars represent the ratios derived from the 12CO 4.67 and 2.35 µm bands, respectively. The error bars on the JOYS+ data points include the error on the band strengths, which is the main contributor to these uncertainties. The ratios were taken from the following works: Ice age (McClure et al. 2023), (Boogert et al. 2000), KECK-NIRSpec (Boogert et al. 2002a), VLT (Pontoppidan et al. 2003), VLT-CRIRES (Smith et al. 2015), PILS (Jørgensen et al. 2016, 2018), Disks (Yoshida et al. 2022; Bergin et al. 2024), Exo-planets (Line et al. 2021; Zhang et al. 2021a,b; Gandhi et al. 2023; de Regt et al. 2024), Comets (Bockelée-Morvan et al. 2015; Hässig et al. 2017). The dashed purple line shows the solar abundance, and the dot-dashed gray line shows the local ISM ratio (Boogert et al. 2000). Further details of the error analysis can be found in Appendix A.

In the text

	Fig. A.1 Sampling distrubition of the integrated optical depth. The dashed shows the Gaussian curve fitted to the distribution from which we extract σnoise and int.OD,noise.
In the text