© The Authors 2025

1 Introduction

The ¹²C/¹³C isotopic ratio is one of the most extensively studied in the interstellar medium (ISM). The main ¹²C and the ¹³ C isotopes are produced in the interior of stars via nucleosynthesis, where ¹²C is a mostly primary product and ¹³C is a primary and secondary product (e.g. Romano 2022). Therefore, the elemental ¹²C/¹³C ratio is set by the local stellar activity and changes over the lifetime of the galaxy. In the Milky Way, the ¹²C/¹³C ratio increases with galactocentric distance, where values of 25 ± 7 are expected for the Galactic centre region and 68 ± 15 for the solar neighbourhood ISM (Milam et al. 2005, based on observation of CN, CO, and H₂CO).

Similar to the main isotope, ¹³C can react to form molecules. The ¹²C/¹³C isotopologue ratios derived in the Solar System are quite constant, while there can be significant deviations from the local ISM value of 68 (see Fig. 1) between nearby star forming regions or between different organic molecules (e.g. Nomura et al. 2023). The ¹³C gets incorporated into CO via the following exchange reaction (e.g. Colzi et al. 2020): $${}^{13}{\mathrm{C}}^{+}+{}^{12}\mathrm{C}\mathrm{O}\leftrightarrow {}^{12}{\mathrm{C}}^{+}+{}^{13}\mathrm{C}\mathrm{O}+34.7\mathrm{K},$$(1)

which preferentially produces ¹³CO in the low temperature gas phase (Langer et al. 1984) before it freezes out on dust grains, where it can then be involved in reactions that produce more complex species such as CH₃OH. Conversely, molecules formed from C⁺, such as HCN or CN, would be depleted in ¹³C. However, low ¹²C/¹³C ratios have also been measured for these species (see Fig. 1). New reactions have subsequently been explored and found to be efficient in enriching these molecules in ¹³C (e.g. Roueff et al. 2015; Colzi et al. 2020; Loison et al. 2020), for example: $${}^{13}{\mathrm{C}}^{+}+{}^{12}\mathrm{C}\mathrm{N}\leftrightarrow {}^{12}{\mathrm{C}}^{+}+{}^{13}\mathrm{C}\mathrm{N}+31.1\mathrm{K},$$(2) $${}^{13}\mathrm{C}+{}^{12}\mathrm{C}\mathrm{N}\leftrightarrow {}^{12}\mathrm{C}+{}^{13}\mathrm{C}\mathrm{N}+31.1\mathrm{K},$$(3) $${}^{13}\mathrm{C}+{\mathrm{H}}^{12}\mathrm{C}\mathrm{N}\leftrightarrow {}^{12}\mathrm{C}+{\mathrm{H}}^{13}\mathrm{C}\mathrm{N}+48\mathrm{K}.$$(4)

In addition to isotopic fractionation by exchange reactions, there are other processes that can locally alter the isotopic ratio, for example, selective photodissociation (e.g. Visser et al. 2009).

These local variations in the fractionation of molecules are well suited to shedding light onto the formation of the species in general but also on the physical conditions that certain isotopic ratios may be associated with. In particular, isotopologue fractionation can be used to address the question of inheritance or in situ formation of molecules during the various stages of the star formation process, that is, from cold molecular clouds to protostellar environments and, eventually, planets (e.g. Caselli & Ceccarelli 2012; Wampfler et al. 2014; Cleeves et al. 2014). For example, the deuterium fractionation of molecules and the D/H ratio have been studied, and evidence for the inheritance of deuterated molecules from the prestellar to later stages has been found (e.g. Ceccarelli et al. 2014; Drozdovskaya et al. 2021).

Recently, ¹²C/¹³C ratios have been determined in the atmospheres of exoplanets using CO isotopologues: the super Jupiters VHS 1256b (Gandhi et al. 2023) and TYC 8998-760-1b (Zhang et al. 2021) and the hot Jupiter WASP-77Ab (Line et al. 2021). For VHS 1256b a ¹²C/¹³C ratio close to the expected local ISM value of 62 ± 2 was derived, while the ratios for the other two are significantly lower: ${31}_{-10}^{+17}$ and 10.2-42.6, respectively (see Fig. 1). Low ¹²C/¹³C ratios for CO isotopologues were also reported for the protoplanetary disks DM TAU (20 ± 3), LkCa 15 (11 ± 2), and MWC 480 (8±2) by Piétu et al. (2007) as well as for the disk around TW Hya (~20 and 40; Zhang et al. 2017; Yoshida et al. 2022). The values are summarised in Fig. 1. All together, these observations indicate that exoplanets likely inherit the ¹²C/¹³C of the disk material from which they form. In contrast to CO, ratios similar to or even higher than the local ISM value were derived for HCN isotopologues in PDS 70 (Rampinelli et al. 2025) and for HCN and CN in the TW Hya disk (Hily-Blant et al. 2019; Yoshida et al. 2024). The values in the TW Hya disk were derived within a radius of ~110 au. A high ¹²C/¹³C ratio (>84) for CO was found at larger radii (Yoshida et al. 2022). All ratios that have been derived for TW Hya are summarised in Bergin et al. (2024). In order to get to the bottom of the apparent presence of two carbon isotopic reservoirs, these authors used thermochemical models to trace changes in the ¹²C/¹³C chemistry as a function of time in the disk. The observations can be reproduced with a model ofa disk at an early evolutionary stage that has C/O > 1, has no carbon depletion, and is exposed to cosmic rays (i.e. cosmic rays are not attenuated due to high densities in the disk). This effect of [C/O]_elem > 1 on the ¹²C/¹³C ratios has been confirmed in simulations by Lee et al. (2024). These simulations assume that the two observed carbon isotope reservoirs are the result of in situ chemistry in the disk.

However, low ¹²C/¹³C ratios have also been found towards a few low-mass pre- and protostellar sources¹ (see Fig. 1). Reliable measurements at these early stages remain scarce, as emission from simple species such as CO, HCN, or CN is optically thick, sometimes even emission from the rarer isotopologues. Instead, complex organic molecules (COMs; C-bearing molecules with >6 atoms) can be used to shed light on the evolution of the ¹²C/¹³C ratio. This has only been done for a variety of COMs towards the IRAS 16293-2422 (IRAS 16293 hereafter) low-mass protostellar system. With this work, we aim to increase the number of measurements of ¹²C/¹³C ratios towards young protostars to see whether low ratios are already observed at this earlier evolutionary stage. In addition, we test whether there is a difference between N- and O-bearing molecules in the ¹²C/¹³C ratios that may hint at the presence of two carbon isotope reservoirs, as indicated in TW Hya.

In addition, it is crucial to know of any deviations from the local ISM value at these stages of star formation to accurately derive molecular abundances. For example, the simplest and most abundant of the COMs, methanol (CH₃OH), is often used to normalise column densities of other COMs to compare molecular inventories of different sources (e.g. Belloche et al. 2020; van Gelder et al. 2020; Yang et al. 2021). However, the main isotopologue is likely optically thick towards young protostars, so its optically thin(ner) ¹³C isotopologue can be used to infer the main isotopologue’s column density assuming a ¹²C/¹³C ratio (e.g. van Gelder et al. 2022). Therefore, the ¹²C/¹³C ratio can be a crucial parameter in determining and comparing chemical inventories.

We made use of the data taken as part of the Northern Extended Millimetre Array (NOEMA) large program PROtostars & DIsks: Global Evolution (PRODIGE, PIs: P. Caselli and Th. Henning) at 1 mm. The project targets 30 Class 0/I systems in the Perseus molecular cloud, which is at a distance of 294 ± 17pc (Zucker et al. 2018). The survey covers a plethora of molecular lines thanks to its broad bandwidth of 16 GHz. The first results using this data set were published in Valdivia-Mena et al. (2022), Hsieh et al. (2023, 2024), and Gieser et al. (2024). We used the isotopologues² of methanol (CH₃OH and ¹³CH₃OH) and methyl cyanide (CH₃CN, ¹³CH₃CN, and CH₃¹³CN) to probe their ¹²C/¹³C ratios in the hot gas towards the selected sources. In addition, we include CH¹₃⁸OH in the analysis to address optical depth effects. This article is structured as follows: Sect. 2 briefly describes the observations, data reduction, and imaging. In Sect. 3, we describe the methods of the spectral line analysis and present the results. The results are compared with studies towards other sources and with predictions of astrochemical models in Sect. 4. Conclusions are provided in Sect. 5.

Fig. 1 Literature 12C/13C ratios towards sources at different evolutionary stages using different isotopologues. Simple O- and N-bearing species (<6 atoms) share a light blue and light yellow colour, respectively, but have different markers. Black rectangles around the orange triangle (N) and coloured circles (O) indicate complex organics (>6 atoms). Arrows indicate lower limits. References to the values can be found in Table A.9.

2 The PRODIGE data set

The PRODIGE large program is part of the MPG/IRAM Observatory Program (MIOP, Project ID L19MB) observed with NOEMA. A total of 30 Class 0/I protostellar systems were targeted, of which we study seven in this work (see Table 1 and Sect. 3). The observations were carried out using the Band 3 receiver and the PolyFix correlator, tuned to a local-oscillator (LO) frequency of 226.5 GHz. This setup provides a total bandwidth of 16 GHz with a channel width of 2 MHz, corresponding to ~2.5-2.8 km s⁻¹. There are four sidebands: lower outer (214.7-218.8 GHz), lower inner (218.8-222.8 GHz), upper inner (230.2-234.2 GHz), and upper outer (234.2-238.3 GHz). In addition, 39 high spectral-resolution windows (62.5 kHz channel width), each covering a 64 MHz bandwidth, were placed within the 16 GHz bandwidth.

Observations were conducted in antenna configurations C and D and finished in 2024. In this work, we only used data taken in 2019 to 2022. The maximum recoverable scale (MRS) is 16.9″ at 220 GHz, which corresponds to approximately 5000 au at the distance of Perseus. The data calibration was done using the standard observatory pipeline GILDAS/CLIC³ package. We used 3C84 and 3C454.3 as bandpass calibrators and 0333+321 and 0322+222 as phase and amplitude calibrators. The observations for these calibrators were taken every 20 min. For the flux calibration we used LkHα 101 and MWC 349. The uncertainty in absolute flux density is 10%. Continuum subtraction and data imaging were done with the GILDAS/MAPPING package using the uv_baseline and clean tasks, respectively. The imaging of the spectral line cubes was done with natural weighting to minimise noise, while for the continuum maps we used robust = 1, to improve the angular resolution. More details on data reduction and imaging can be found in Gieser et al. (2024). The achieved angular resolution is ~1″, which corresponds to a spatial resolution of ~300 au. The accurate synthesised beams for each source as well as root-mean-square (rms) noise levels can be found in Table A.1.

3 Spectral line LTE analysis and results

We detected emission from the rarer ¹³C isotopologues of methanol and methyl cyanide towards seven and six protostellar sources of the PRODIGE sample (see Table 1), respectively. In this section we first present the morphology of the molecular emission and the position selection (Sect. 3.1), then discuss the line widths and peak velocities (Sect. 3.2). The following spectral line analysis is performed assuming a uniform medium that is in local thermodynamic equilibrium (LTE), which is applicable given the high volume densities (Yang et al. 2021). The spectral line identification and analysis were done using Weeds, which is an extension of the GILDAS/CLASS software (Maret et al. 2011). Weeds produces synthetic spectra by solving the radiative transfer equation. More details on the modelling with Weeds are in Sect. 3.3. The molecular column densities as well as the rotational temperatures that were obtained from this method are validated by results from a population-diagram analysis, which is described in Sect. 3.4. The results are presented in Sects. 3.5 and 3.6.

Fig. 2 Overview of the PRODIGE data towards IRAS 4B and presentation of the spectral line analyses. Panel a: integrated intensity map for th CH3CN 123−113 transition at 220.708 GHz. The beam (half power beam width, HPBW) is shown in the bottom left. The dashed yellow contou is at 5σ, where σ is the rms noise level and is written in the top-right corner. Panel b: observed spectra towards the peak position, which is a (Δα, Δβ = (0″, 0″) (black cross in panel a), of CH3CN at low (~2 MHz, black) and high (~62.5 kHz, teal) spectral resolution, where the latte] was shifted by 4K for visualisation. The low spectral-resolution data were analysed in this work (see Sect. 3). The modelled CH3CN spectra from Weeds are shown in red (low resolution) and orange (high resolution), where the same input parameters were used for both. Panel c: population diagram for CH3CN. Teal circles show the observed data points while the modelled data points from Weeds are shown with orange squares (se Sect. 4.6 for possible origins of differences between the two). No corrections are applied in the left panel, while in the right panel, corrections fo opacity and contamination by other molecules have been considered (negligible effect in this case), and only transitions with optical depths <0.5 (i.e. all transitions in this case) are shown. The results of the linear fit to the observed data points are shown in the right panels.

3.1 COM emission morphology and position selection

As an example, Fig. 2a shows an integrated intensity map of CH₃CN towards IRAS 4B. Integrated intensity maps of the ¹²C and ¹³C isotopologues of CH₃OH and CH₃ CN for all sources are shown in Figs. B.1–B.7. In most cases, the emission is compact and remains spatially unresolved in the PRODIGE data. Only for IRAS 4B, emission outside the compact hot-corino region is identified, which can be associated with its protostellar outflow (e.g. Sakai et al. 2012).

For the LTE analyses, we extracted the spectra at the continuum peak position at (0″ ,0″ ), which corresponds to the coordinates in Table 1. Although the line data may be affected by continuum optical depth, we did not select a position off peak (as done, for example, for IRAS 16293 by Jørgensen et al. 2016).

If we did so, we would risk to detect fewer transitions, which might, in turn, result in larger uncertainties in the LTE analyses as both methods require as many lines as possible.

3.2 Line widths and peak velocities

For the spectral line analyses (Sects. 3.3 and 3.4), we made use of the low spectral-resolution (LR) data that covers the total 16 GHz bandwidth and hence more spectral lines of our molecules of interest. However, because of the narrow line widths of the spectral lines (full width at half maximum, FWHM, of <4 km s⁻¹, see Tables A.2–A.8) they are only marginally spectrally resolved if at all. To determine the true FWHM and peak velocity, we used the higher resolution (HR) data. The HR data only cover transitions of both methanol isotopologues in their vibrational ground state, while the LR data also cover vibrationally excited states. Figure 3 shows peak velocities, ν_peak, and FWHM, which we derived from 1D Gaussian fits to all available transitions, as a function of upper-level energy, E_u. The median values for both species are shown with dashed lines. Outliers may appear as a consequence of weak transitions or unidentified contamination by other species, which can lead to shifts. In addition, Fig. 3 only shows transitions with τ < 10, though some FWHM may be larger due to optical depth effects. For some sources, the presence of outflow wing emission may broaden the lines. In some cases, there are significant differences between ¹²C and ¹³C methanol, which may indicate that the two species probe different material along the line of sight. Higher spectral resolution data that cover a larger number of spectral lines per species will help draw conclusions on this.

In the following, we ignore possible kinematic differences between species and use an effective FWHM, that is the lower median between the two species, and the respective ν_peak value. Because of their unresolved hyperfine structure, methyl cyanide isotopologues are not used to determine line properties. However, the hyperfine structure (HFS) is considered later during the determination of column densities, except for ¹³CH₃CN, as the Cologne database of molecular spectroscopy (CDMS; Endres et al. 2016), from which all spectroscopic data were taken, does not include the HFS for frequencies beyond 220 GHz.

Fig. 3 Peak velocity, νpeak, and line width, FWHM, derived from 1D Gaussian fits (including 1σ error bars) to the high spectral-resolution (HR) PRODIGE data as a function of upper-level energy, Eu, for the vibrational ground states of CH3 OH (orange) and 13CH3OH (teal). The coloured dashed lines indicate the respective average values for νpeak and FWHM.

3.3 Radiative transfer modelling with Weeds

We classify a molecule as detected when there are at least five transitions per isotopologue in the observed LR spectra above an intensity threshold of 3σ, where σ is the rms noise level measured in the 3D data cubes (see Table A.1), and the modelled spectra reasonably reproduce these transitions. In order to derive a synthetic spectrum for a molecule, Weeds requires five input parameters: total column density, rotational temperature, source size, FWHM, and velocity offset, ν_off = ν_peak,HR - ν_sys, of the spectral line from the source systemic velocity, where ν_sys is taken from Table 1. The FWHM and ν_peak,HR present observed values, for which we ignore any unresolved kinematic substructures and velocity gradients, and were determined from 1D Gaussian fits in Sect. 3.2. For FWHM that are spectrally unresolved in the observed LR spectra, Weeds smooths the line width for the given channel width. The source size, that is the size (in arcseconds) of the emitting area is needed to compute the beam filling factor.

We adopted source sizes from the literature (Table 1), which were estimated from COM emission (Belloche et al. 2020; Chen et al. 2024) or, if this was not available, from continuum emission (Yang et al. 2021). Because Weeds does not perform a fit to the observed data, the total column density, N_col, and the rotational temperature, T_rot, were adjusted in an iterative way: starting from a first educated guess, population diagrams (see Sect. 3.4) were derived and the resulting values for N_col and T_rot were fed back into Weeds. This was repeated until both Weeds and the population diagrams provide matching results for N_col and T_rot within the uncertainties. In the process, Weeds computes line optical depths and takes it into account when creating the model spectrum, based on the input parameters. We used this opacity correction when deriving population diagrams (see Sect. 3.4). The parameters of all final models are summarised in Tables A.2–A.8. In Fig. 2b we show, as an example, a spectrum of CH₃CN towards IRAS 4B from the low (black) and high (teal) spectral resolution PRODIGE data. The former is used for the spectral line analysis because a wider bandwidth and hence a larger number of transitions was covered. The derived Weeds model is shown in red. For comparison, we also show the model of the high resolution case in orange.

It happens that the final model over- or underestimates the peak intensity for a few transitions, even if it (and the corresponding population diagram) matches well the rest of the observed spectrum. If the model underestimates the intensity, this likely indicates contamination by another species that adds to the observed intensity. In severe cases, the line is excluded from further analysis. It may also be that the intensity of a transition is not well estimated by the model because the observed emission is not described by a single excitation temperature. This cannot be accounted for in the Weeds models (nor the population-diagram analysis) and can result in an over-or underestimation by the model. Especially, transitions with lower upper-level energies suffer from this effect as they can be excited more easily. Deviations from Gaussian profiles caused by such a non-uniform medium or optical depth effects may be worsened due to spectral dilution in the LR data. Such limitations of this analysis and their consequences are taken up on in Sect. 4.6.

3.4 Population diagrams

To verify total molecular column densities, N_tot, and rotational temperatures, T_rot, used to derive the synthetic spectra, we derived population diagrams (PDs), which are based on the following formalism (Goldsmith & Langer 1999; Mangum & Shirley 2015): $$\ln \left(\frac{N_u}{g_u}\right)=\ln \left(\frac{8\pi k_B{\nu }^2\int J(T_B)\mathrm{d}v}{c^{3}h{A}_{\mathrm{u}}{g}_{u}f}\right)=\ln \left(\frac{N_{\mathrm{t}\mathrm{o}\mathrm{t}}}{Q(T_{\mathrm{r}\mathrm{o}\mathrm{t}})}\right)-\frac{E_u}{k_B{T}_{\mathrm{r}\mathrm{o}\mathrm{t}}},$$(5)

where N_u is the upper-level column density, g_u the upper-level degeneracy, E_u the upper-level energy, A_u the Einstein A coefficient, k_B the Boltzmann constant, c the speed of light, h the Planck constant, $f=\frac{{\mathrm{s}\mathrm{o}\mathrm{u}\mathrm{r}\mathrm{c}\mathrm{e}\mathrm{s}\mathrm{i}\mathrm{z}\mathrm{e}}^2}{{\mathrm{s}\mathrm{o}\mathrm{u}\mathrm{r}\mathrm{c}\mathrm{e}\mathrm{s}\mathrm{i}\mathrm{z}\mathrm{e}}^2+{\mathrm{b}\mathrm{e}\mathrm{a}\mathrm{m}\mathrm{s}\mathrm{i}\mathrm{z}\mathrm{e}}^2}$ the beam filling factor, and Q the partition function. Intensities in brightness temperature scale, J(T_B), are integrated over a visually selected velocity range, dν, in the continuum-subtracted spectra. When computing the left-hand term for each transition of a species and plotting it against E_u, all data points should ideally follow a linear trend, in which case the level populations can be explained by a single rotational temperature.

In Fig. 2c we show one PD of CH₃CN obtained for IRAS 4B as an example (all other PDs can be found in Figs. C.1–C.8 in the appendix). Both the observed data points as well as those derived from the Weeds models are shown. If detected, rotational transitions from vibrationally excited states (v = 1, v = 2 for CH₃OH and v₈ = 1 for CH₃CN) are used in addition to transitions from the vibrational ground state (v = 0). The left panel shows the original data obtained from Eq. (5), while two correction factors (i.e. optical depth and contaminating emission, see below) are applied to the modelled and observed data in the right panel (negligible effect in Fig. 2c). One of the correction factors considers line opacities. Besides computing the synthetic spectra for a molecule, Weeds also calculates the optical depth for each transition for the given input parameters. Using the respective peak opacity from Weeds, τ_i, for a given transition i, we computed a correction factor, $\frac{{\tau }_i}{1-{\mathrm{e}}^{-{\tau }_i}}$ (e.g. Mangum & Shirley 2015). This factor was then multiplied by the observed and modelled integrated intensities. However, the right-hand panels of the PDs only show transitions for which Weeds computes optical depths of <0.5. Plus, the linear fit is only applied to those data points. Since it is possible that we still underestimated optical depth for some lines (see Sect. 4.6), we did not use transitions with higher opacities to ensure that the fit parameters apply to the sufficiently optically thin lines.

Secondly, contaminating emission from other species was subtracted from the integrated intensities, if identified. In particular, this applies to CH₃CN and CH₃¹³CN, which have overlapping transitions at ~220.636 GHz. Assuming that the emission spatially coincides, this addition of spectra from different species and a subsequent subtraction of contaminating emission only works for optically thin transitions. Even with the two correction factors applied, some scatter between the observed and modelled data points remains. In the example in Fig. 2c, the model clearly underestimates the observations. This is a result of the narrow line widths. As can be seen from the spectra in Fig. 2b, the spectral lines are rarely broader than one channel. However, for the integration we used at least 2-3 channels, adding contaminating emission that seems to be an additional CH₃CN component in the case of IRAS 4B, possibly associated with its outflow. Observational limitations such as these are discussed in detail in Sect. 4.6.

The linear fit is applied to the corrected observed data without taking into account their displayed uncertainties, since not all of them can be quantified, such as possible contamination from unidentified species. The displayed error bars only include the standard deviation from integrating intensities and the 10% uncertainty in absolute flux density. From the linear fit, we obtain T_rot from the slope and N_tot from the y-intercept. The latter requires knowledge of Q(T_rot), which is taken from the CDMS⁴, together with other spectroscopic information, such as A_u, E_u, and the rest frequencies.

Fig. 4 Rotational temperatures used for the Weeds modelling (crosses) and derived from the population diagrams (circles with 1σ error bars). The dashed lines indicate the median values considering all isotopologues of a species.

3.5 Rotational temperatures and ¹²C/¹³C isotopologue ratios

Figure 4 compares the rotational temperatures derived from both the Weeds modelling and the PDs for the various isotopologues in each source. In all shown cases, the temperatures derived from both methods agree within the error bars. Large error bars result from a smaller number of mostly weak transitions that were available to derive the respective PD. Some PD temperatures are not shown either because the molecule was not detected (both ¹³C CH₃CN isotopologues towards L1448-IRS2) or because the PD fit was unreliable and T_rot was fixed in the PDs to a value from the other detected species (see Tables A.2–A.8).

Overall, the results show that, for most sources, CH₃CN has higher temperatures than CH₃OH. The median rotational temperature derived from methanol is 180K (140-400 K), that from methyl cyanide is 290K (180-420 K). In some cases, the temperature values for the respective isotopologues are different. For example, for B1-bS and IRAS 2A, CH₃ OH shows higher temperature than ¹³CH₃OH. This similarly applies to the CH₃CN isotopologues for IRAS 4B. This may indicate that the isotopologues probe different excitation conditions.

Table 2 summarises the ¹²C/¹³C isotopologue ratios computed from the Weeds column densities of CH₃OH and CH₃CN for the PRODIGE sample sources. The ratios are shown in comparison to other sources in Fig. 5, where we used an average value for the two ¹³C isotopologues of methyl cyanide. The error bars are those from the linear PD fits. The ratios range from 4 to 30 and are thus all significantly below the local ISM value of 68. There is no significant difference between CH₃ OH and CH₃ CN. In some cases, uncertainties are likely larger than depicted (see Sect. 4.6).

Fig. 5 Same as Fig. 1 but including the 12C/13C ratios derived from the PRODIGE data in this work (light red shaded). Bolometric luminosities of the protostars increase from left to right. Black rectangles around markers indicate complex organic molecules. Arrows indicate lower limits.

3.6 CH₃¹⁸OH

The ¹⁸O methanol isotopologue is detected towards four of the seven sources of our sample: SVS13A, IRAS 2A, IRAS 4B, and B1-c. We applied the same analysis to this molecule: synthetic spectra were derived by fixing all parameters to those of ¹³CH₃OH, except for the rotational temperature and column density. Population diagrams can be found in Fig. C.8.

The derived rotational temperatures (Fig. 4) remain below 200 K and are hence either comparable to or lower by a few tens of Kelvin than the other CH₃OH isotopologues. Figure 6 shows the resulting column density ratios between ¹³CH₃OH and CH₃¹⁸OH together with lower limits on this ratio for the sources towards which the latter species remains undetected. Based on the local ISM value of ¹⁶O/¹⁸O = 557 (e.g. Nomura et al. 2023), we expect ${}^{13}{\mathrm{C}\mathrm{H}}_3^{16}\mathrm{O}\mathrm{H}{/}^{12}{\mathrm{C}\mathrm{H}}_3^{18}\mathrm{O}\mathrm{H}\sim 8$, which is consistent with our results within the error bars.

Fig. 6 Column density ratio between 13CH3OH and CH318OH using the Weeds column densities. The expected local ISM value of ~8 is indicated with a dashed line. Arrows show lower limits.

4 Discussion

4.1 Low ¹²C/¹³C COM isotopologue ratios in other protostellar systems

Low ¹²C/¹³C ratios for (COM) isotopologues towards protostellar systems have been reported in the past. Using the PRODIGE data, Hsieh et al. (2023) derived a similarly low ¹²C/¹³C ratio (~15—25) for CH₃CN towards the binary protostellar system SVS13A. These authors tested whether modelling one or two excitation components would lead to a different outcome, which it did not, as low ratios were seen in either scenario. In a followup study at much higher angular resolution with the Atacama Large Millimetre/submillimetre Array (ALMA) the molecular emission was resolved and the low ¹²C/¹³C ratios for CH₃CN were confirmed (Hsieh et al. 2025) towards one of the protostars.

Based on data obtained with ALMA at ~0.5″ angular resolution as part of the Protostellar Interferometric Line Survey (PILS) towards the young protostellar triple system IRAS 16293-2422 (IRAS 16293 hereafter), Jørgensen et al. (2016) reported ¹²C/¹³C ratios of 24-35 for the glycolaldehyde (CH₂OHCHO) isotopologues towards IRAS 16293B. Based on the same observational data, low ¹²C/¹³C isotopologue ratios were also (tentatively) found for C₂H₅OH, CH₃OCHO, and CH₃OCH₃ (see Fig. 5 and Jørgensen et al. 2018). These three O-bearing molecules, together with glycolaldehyde and methanol, likely form in the solid phase during the cold prestellar stage as predicted by astrochemical models (e.g. Garrod et al. 2022). They are formed from CO that was previously enriched with ¹³C via the exchange reaction (Eq. (1)) in the cold gas phase prior to freeze-out. This similarly applies to CH₃CHO, which can be formed in the solid phase from CH₂CO, which in turn is formed from CO (e.g. Ibrahim et al. 2022; Garrod et al. 2022). However, these two species and t-HCOOH, whose formation also includes CO, have ¹²C/¹³C ratios close to the local ISM value in IRAS 16293B, based on the PILS data. At temperatures >100 K, CH₃CHO has an efficient formation route in the gas phase (C₂H₅ + O, e.g. Garrod et al. 2022), which might explain why its ¹²C/¹³C ratio is different from the other O-bearing COMs, that is, if the ratio for C₂H₅ were as high. Again using the PILS data, Calcutt et al. (2018) derived ¹²C/¹³C for the CH₃CN isotopologues of ~70-80 in IRAS 16293A and B, which is in contrast to our results. Differences between starforming regions in the ¹²C/¹³C ratio of a given species may be linked to the duration of the prestellar phase as indicated by models (see Sect. 4.5).

The ¹²C/¹³C COM isotopologue ratio has also been investigated within the framework of the ALMA Survey of Orion Planck Galactic Cold Clumps (ALMASOP; Hsu et al. 2020, 2022) at an angular resolution of −0.4″. Of the 72 targeted clumps (some containing multiple pre- and protostellar sources) 56 Class 0/I protostars have been identified (Dutta et al. 2020), of which eleven show warm methanol emission (Hsu et al. 2022). For another six of these sources, ¹³CH₃OH has been detected. The resulting ¹²C/¹³C ratios are all below 5, likely because not many methanol transitions are covered in the survey resulting in a strong underestimation of the main isotopologue’s optical depth.

Recently, low ¹²C/¹³C ratios of 20-30 have also been measured for several O-bearing COMs towards the V883 Ori disk (Yamato et al. 2024; Jeong et al. 2025), based on ALMA data at 3 mm. Similar to the considerations done for IRAS 16293, the authors argue that one possible reason may be the COMs’ formation from CO in ices at early stages, where CO was enriched in ¹³C prior to freeze-out, but further investigation is needed.

4.2 No evidence for an O/N dichotomy based on the ¹²C/¹³C ratio

We selected CH₃OH and CH₃CN for this work because they trace the warm gas towards the protostellar sources and to study a possible segregation between the O- and N-bearing molecules based on their ¹²C/¹³C isotopologue ratios. For the disk around the T Tauri star TW Hya (see Sect. 1 and, e.g., Bergin et al. 2024), low ¹²C/¹³C ratios were derived for CO inside a 100 au radius and ratios close to or higher than the local ISM value for C₂H, HCN and CN (see also Fig. 5). The presence of two different carbon isotopic reservoirs was proposed as the origin for this dichotomy. Disk simulations (Lee et al. 2024) were able to reproduce the various ¹²C/¹³C isotopologue ratios when considering C/O > 1. In this work, we explored whether there is any difference between CH₃OH and CH₃CN in terms of the ¹²C/¹³C ratio at the Class 0 stage. First evidence for this can be seen from the ratios derived towards IRAS 16293B (see Sect. 4.1 and Fig. 5), where (some) O-bearing COMs show an enhancement in ¹³C, while CH₃CN shows ISM-like values. Our results do not show significant differences between CH₃OH and CH₃CN (Fig. 5). It is worth mentioning that for some sources the rotational temperatures that we derived for the CH₃OH and CH₃CN isotopologues were different. This may indicate that they probe different excitation conditions. This effect on the ¹²C/¹³C ratio and others that arise from observational limitations (see Sect. 4.6) can be tackled in observations at higher angular and spectral resolution.

4.3 Implications of the ¹³CH₃OH/CH₃¹⁸OH ratio

The local ISM ¹⁶O/¹⁸O isotopic ratio is 557 ± 30 (Wilson 1999). Dividing this value by the ¹²C/¹³C ratio of 68, yields a theoretical ¹³C¹⁶O/¹²C¹⁸O ratio of –8, which agrees with the values that we derived for CH₃OH towards the PRODIGE sources within the error bars (see Fig. 6). Therefore, based on this isotopologue ratio, we do not expect an enrichment in ¹³C, which may imply that we are underestimating the column density of the main isotopologue. As mentioned in Sect. 4.2, differences in rotational temperatures between CH₃OH isotopologues may affect the ratios.

Regarding the robustness of the ¹³C¹⁶O/¹²C¹⁸O ratios, lower values would mean that the abundance of CH₃¹⁸OH is enhanced, which has not been reported so far. On the contrary, oxygen fractionation is proposed to be a consequence of isotope-selective photodissociation that depletes CO-related molecules in ¹⁸O (e.g. Nomura et al. 2023). Loison et al. (2019) modelled the ¹⁶O/¹⁸O ratio in dark cloud conditions and compared it with observations. They found that in order for molecules to be enriched in ¹⁸O, the isotope has to be abundant in the gas phase. Consequently, only a small fraction would be available in ices. They focussed on S-bearing species, but also computed the ¹⁶O/¹⁸O for CO and CH₃OH in the gas phase as a function of time and find no significant deviations from the local ISM value. More simulations are needed to ascertain whether this remains to be the case during the protostellar phase.

The ¹³C¹⁶O/¹²C¹⁸O ratios have recently been determined for Class II sources, which were observed as part of the PRODIGE - Planet-forming disks in Taurus with NOEMA (Semenov et al. 2024). They were mostly found to be close to the expected value of 8, in agreement with the ratios that we found at earlier stages, although some of the Class II disk systems are rather described by isotopologue ratios deviant from the ISM values, according to simulations (Franceschi et al. 2024).

4.4 Tentative evidence for an O/N temperature segregation

Segregation between O- and N-bearing molecules is frequently observed (e.g. Widicus Weaver et al. 2017; Allen et al. 2017; Jørgensen et al. 2018; Csengeri et al. 2019; Pagani et al. 2019; Busch et al. 2022, 2024; Bouscasse et al. 2024), either in terms of spatial variations in their emission or in terms of differing rotational temperature and abundance. We derived on average higher rotational temperatures (Fig. 4) for the CH₃CN (290 K) than for the CH₃OH isotopologues (180 K), implying that the two species probe different conditions in the gas, which may trace back to their dominant formation and destruction mechanisms.

Methanol gas-phase abundances are predicted to peak at ~ 100150K (Garrod et al. 2022). At these temperatures the COM thermally desorbs from the icy dust grains, where it was formed prior to the birth of the protostar via successive hydrogenation of CO. Due to a lack of efficient gas-phase formation routes, the COM was found to be destroyed in the high-temperature gas phase right after its desorption (Garrod et al. 2022; Busch et al. 2022). On the other hand, methyl cyanide can efficiently be produced at high temperatures in the gas phase (Garrod et al. 2022). Another process that may explain this difference is carbon-grain sublimation (van’t Hoff et al. 2020, see also Nazari et al. 2023 and van’t Hoff et al. 2024). It was proposed to result in enhanced abundances of N-bearing COMs at temperatures >300 K through release of more refractory C- and N-rich species from the grain mantle.

Some of our analysed sources were also covered in the analysis of the spectral line data obtained as part of the IRAM Plateau de Bure Interferometre (PdBI) large program CALYPSO (Continuum And Lines in Young ProtoStellar Objects; Belloche et al. 2020). This survey resembles the PRODIGE survey in angular and spectral resolution but covers a smaller bandwidth. The authors used the same methods of analysis as in this work and do not find a significant difference in rotational temperature between O- and N-bearing molecules in general. Their values are likely more uncertain because they cover fewer transitions and a smaller range of upper-level energies. The authors do not report column densities for the rarer isotopologues of CH₃ OH and CH₃CN. Their column densities of the main isotopologues agree within a factor of 2 with ours, except for SVS13A and IRAS 2A, for which we derived methanol column densities that are a factor of 5.5 and 8 higher, respectively, than reported by Belloche et al. (2020). We assumed a source size of 0.2″ (Chen et al. 2024) for IRAS 2A, whereas Belloche et al. (2020) used 0.35″, resulting in a threefold difference in the column densities. The remaining difference is likely attributed to optical depth. According to our Weeds models for the two sources, all methanol transitions with E_u < 600 K show opacities above 0.5. Hence, we did not include them in the PD fits. In contrast, Belloche et al. (2020) use all transitions, likely resulting in an underestimation of the column density.

4.5 Comparison with astrochemical models

We ran simulations of carbon and nitrogen isotopic fractionation using the gas-grain chemical code pyRate, whose basic properties are described in Sipilä et al. (2015). We employed the gas-phase and grain-surface chemical networks presented recently in Sipilä et al. (2023), where an extensive set of isotopic exchange reactions for C- and N-bearing species - also taking into account the possibility of simultaneous isotopic substitution in both elements (e.g. H¹³C¹⁵N) - is considered, though in the present paper we do not provide predictions for ¹⁵N-bearing species. The networks, however, only include species consisting of up to five atoms; hence, we cannot obtain direct estimates for the fractionation of CH₃OH and CH₃CN. Instead, we simulated the abundances of precursor species and used them to estimate the isotopic ratios of CH₃ OH and CH₃CN, based on their most likely formation pathways, to constrain the possible origins of the low ratios derived for these COMs from the present observations.

We have not attempted to simulate each observed source individually. Instead, we have adopted as a template the physical model for the protostellar core IRAS 16293-2422 presented by Crimier et al. (2010), and carried out chemical simulations to predict the ¹²C/¹³C ratios of various species. Here, we give a brief account of the model setup; an expanded description including the initial conditions is given in Appendix D. Following Brünken et al. (2014) and Harju et al. (2017), we have employed a two-stage approach where the chemistry first evolves over a period of 10⁶ yr in an intermediate-density dark cloud, and then for another period of 10⁴ yr in the protostellar conditions using the abundances obtained in the first stage as initial abundances for the second stage. For simplicity, the primary cosmic-ray ionisation rate was set to the canonical value of ζ_p = 1.3 × 10⁻¹⁷ s⁻¹ across the core⁵. The initial elemental ¹²C/¹³C ratio is set to 68 in the first stage (Milam et al. 2005). The model of Crimier et al. (2010) is not appropriate for describing the very innermost regions of the protostellar core (inside a few hundred au) but is applicable in the present context where the observations have not been carried out using high angular resolution.

Figure 7 shows the ¹²C/¹³C ratios for various CH₃OH and CH₃CN precursors as a function of distance from the centre of the protostellar core, taken at the simulation time t₄ = 10⁴ yr in the second (i.e. protostellar) stage of the simulation. We only show this time step because the ratios do not change significantly at later stages, except for H₂CO (see below). The complex features evident in the various ratios as a function of radius, especially prominent at radii in the range 1000-2000 au, are a result of complex interactions involving different binding energies between the species, isotopic exchange reactions etc. Quantifying these is outside the scope of the present paper, where we concentrate only on the innermost regions of the core (left of the vertical dotted line in Fig. 7) and analyse the general trends. In what follows, we discuss the results for species associated with methanol and methyl cyanide separately.

Fig. 7 Model predictions for selected gas-phase 12C/13C ratios as a function of distance from the protostar taken at a simulation time of 104 yr. The constant ratio of 68 is indicated by the horizontal dashed black line. The vertical dotted black line separates the spatial scales discussed in this paper (R < 300 au) from the larger scales, which are not largely discussed but shown for completeness.

4.5.1 Methanol (precursors)

Methanol is formed during the prestellar phase in the ice mantles of dust grains through successive hydrogenation of CO.

Figure D.1 shows the evolution of the ¹²C/¹³C ratio during the prestellar phase, indicating that, as time advances, the ¹²CO/¹³ CO ratio approaches the ISM value of 68 in both gas and ice. We would expect the ratio for methanol to follow a similar trend, so that, when CH₃OH sublimates from the ices, the value should not be too different from what is seen at the end of the prestellar stage. This is higher than what we observe. It may be that the observed protostellar systems spent a shorter time in the prestellar phase than assumed in the model, thereby inheriting lower ¹²C/¹³C ratios (Fig. D.1). Such a time dependence may also explain the differences in the ¹²C/¹³C ratios seen between the Perseus sources and, for example, IRAS 16293 (see Sect. 4.1). However, the ¹²CO/¹³CO ratios are significantly lower only for t < 10⁵ yr. Predictions by Ichimura et al. (2024) show that, between 10⁵−10⁶ yr, the ¹²C/¹³C ratio for CH₃OH in ices during the prestellar phase can be somewhere between 40-60. Hence, although a shorter prestellar phase may not solely explain our observed ratios of <30, it can still affect the observed ¹²C/¹³C ratios observed during the protostellar stage. In Section 4.6 we discuss that the column densities of especially the main isotopologue of methanol derived from observations can be underestimated due to an underestimation of optical depth, and hence the ¹²C/¹³C ratios may be higher. A final conclusion needs further observational and computational efforts.

Moreover, the ¹²C/¹³C ratios close to the local ISM value for CO as predicted by the current simulations are in agreement with observations of ices in the ISM (Brunken et al. 2024, see their Fig. 10, and references therein). In addition, we note that our simulation predictions for the ¹²CO/¹³CO ratios are similar to those presented recently by Bergin et al. (2024) and Ichimura et al. (2024); the latter authors predicted a CH₃OH ¹²C/¹³Cratio close to that of CO (below 68).

One intermediate product of the methanol formation in ices is formaldehyde (H₂CO), which is included in the model. Its gasphase ¹²C/¹³C ratio is shown in Fig. 7 and is ~62 at the scales traced by our observations (i.e. a few 100 au). The slight difference between the CO and formaldehyde ¹²C/¹³C ratios likely results from competition between the formation of H₂CO in the solid phase from CO and in the gas phase; an in-depth examination of the simulation results shows that the H₂CO ¹²C/¹³Cratio rises in the central core at later times (not shown) due to a chain of chemical reactions starting with desorbed CH₄.

4.5.2 Methyl cyanide (precursors)

The pathway(s) that dominate the formation of CH₃CN are much less certain. Recently, the gas-phase chemical network of the COM has been revised by Giani et al. (2023). They used quantum-mechanical calculations to test the relative importance of the various proposed formation routes of CH₃CN. Their results suggest that the COM is dominantly formed from protonated methyl cyanide (CH₃CNH⁺) that recombines with an electron or gives the extra proton to NH₃ upon reaction. Protonated CH₃CN is in turn a product of the reactions CH₃⁺ + HCN and CH₃OH₂⁺ + HNC. Methyl cyanide can also form on dust grain surfaces via CH₃ + CN and hydrogenation of CH₂CN (e.g. Garrod et al. 2022).

In Fig. 7, we show ¹²C/¹³Cratios for CN, HCN, HNC, and CH₃⁺. The first three follow a similar trend: the ratio is significantly below the local ISM value at close distances to the protostar and only increases to above the ISM value beyond 500-1000 au. This is a result of exchange reactions such as Eqs. (2)–(4), whose relative importance depends on time and the prevalent physical conditions. On the other hand, the simulated ¹²C/¹³C ratio of CH₃⁺ is approximately 100 in the inner core, and therefore, if the CH₃⁺ + HCN formation pathway is dominant, we may expect ¹³CH₃CN and CH₃¹³CN to show different fractionation levels. Only for IRAS 4B, is there tentative evidence for a different ratio (factor of 2, see Table 2); however, due to the small number of transitions (<5) available for the rare CH₃CN isotopologues, we cannot be certain.

Therefore, the CH₃OH₂⁺ + HNC formation pathway thus seems more likely, but we cannot quantify this expectation with the present model. Similarly, the model predicts a high (~ 100) ¹²C/¹³C ratio for CH₃ ice and a low ratio (~30) for CN ice at the end of the cloud stage (Fig. D.1), indicating a similar dichotomy between ¹³CH₃CN and CH₃¹³CN if they were predominantly formed on grains via the association of CH₃ and CN. An expanded chemical model, taking into account isotopic fractionation of COMs, is required to draw definite conclusions on the ¹²C/¹³C ratios, also keeping in mind that the physical model used in the present work is extremely simple.

Our simulation predictions for the ¹²C/¹³C ratios for CN, HCN, and HNC are contrary to those presented recently by Ichimura et al. (2024). In their Fig. B5, the authors show that the ratios for these species and CH₃CN remain elevated (80-90) throughout the protostellar phase. During the prestellar phase, the ¹²CN/¹³CN ratio behaves similarly in our and their models. It is not clear whether this discrepancy in the protostellar phase is due to the warm-up scenario considered in their model versus our static approach or to differences in the chemical network. A detailed investigation is outside the scope of the present paper.

4.6 Observational limitations

In Sect. 3.4 we already mentioned that scatter between the observed and synthesised data points in the population diagrams (PDs) shown in Figs. C.1–C.8 arises from the integration of intensities over at least three channels. Hence, by using the LR data (with a channel width of ~2MHz) for the analysis, integrated intensities for lines that are spectrally unresolved (i.e. FWHM ≲ channel width) will inevitably contain some contamination or extra noise. Contamination may come from unidentified species or unresolved additional velocity components from the same species (e.g. outflow emission; see Fig. 2b). In addition, the optical depth of optically thick lines may be underestimated by the model due to the spectral dilution as the true line profiles remain unknown. This fact together with line optical depths that generally exceed the capability of the Weeds model to correct for them (τ & 2) are the main reasons for the discrepancy between transitions with low and high upper-level energies for sources with highest molecular column densities (e.g. for CH₃OH E_u ≲ 600 K for B1-c, IRAS 2A, and SVS13A).

An overall underestimation of a molecule’s column density may be a consequence of high continuum optical depth. If the continuum was optically thick, it would partially obscure the line emission (De Simone et al. 2020). Recent ALMA observations towards IRAS 2A and B1-c (Chen et al. 2024) at 0.9 mm showed that the continuum is indeed optically thick towards the central position of IRAS 2A, however, the effect is negligible for B1-c. The PRODIGE data cover lower frequencies, therefore, continuum opacities are generally lower. We cannot exclude a contribution of this effect to the uncertainties for some (i.e. the brightest) sources, but it does not apply to all. Observations at lower frequencies will help constrain the continuum optical depth at our observed frequencies.

Another source of uncertainty comes from the spatially unresolved emission. An overestimation of the size of the emitting area will lead to an underestimation of the column density. If the emission is optically thick or includes unresolved structures for one but not the other isotopologue, this can also affect the column density ratios. Moreover, we do not resolve any morphological structures, such as binaries (e.g. SVS13A, Hsieh et al. 2025, and IRAS2A, Jørgensen et al. 2022), streamers (e.g. Valdivia-Mena et al. 2024), or outflows (e.g. Podio et al. 2021, see also the HR spectrum in Fig. 2b, which shows outflow wing emission for IRAS 4B). These may come with different velocity components that we are not able to spectrally resolve, either, but rather smooth over. Accordingly, spatial variations in the ¹²C/¹³C ratios can likely be expected. Therefore, we may only see an average ratio in the current observations, which may not reflect the true conditions at all.

Lastly, although we are correcting for line opacities and focus the analysis on the optically thin(ner) lines, there is still a chance that we are underestimating the values, as they are model-dependent. Ultimately, observations at higher angular and spectral resolution that simultaneously cover a broad bandwidth are needed to address these issues and to put better constraints on the ¹²C/¹³C ratios.

5 Conclusions

We aimed to increase the number of measurements of ¹²C/¹³C ratios in the warm gas around young embedded protostellar systems using the data obtained as part of the PRODIGE large program. The ¹³C isotopologues of CH₃OH and CH₃CN were detected towards seven and six sources, respectively, out of a total of 30 Class 0/I protostellar systems targeted by the PRODIGE survey. By using isotopologues of these complex organic molecules (COMs), we ensured that the molecular emission stems from the innermost hot gas surrounding the proto-star(s). Additionally, we could test whether there is a segregation between O- and N-bearing molecules in terms of the ¹²C/¹³C ratio. Our main findings are the following:

1. We derived excitation conditions and column densities using optically thin(ner) spectral lines assuming LTE;

2. The resulting ¹²C/¹³C column density ratios for CH₃ OH and CH₃CN range from 4 to 30, which is lower than the expected ISM value of 68 by factors of a few up to an order of magnitude;

3. There is tentative evidence for an O/N segregation based on rotational temperatures, where CH₃CN isotopologues trace warmer gas on average (290 K) than CH₃OH (180 K). This may hint at the dominant formation pathways of the two COMs on dust grain surfaces during the prestellar phase (CH₃OH) or in the hot gas phase (CH₃CN);

4. There is no significant difference between the ¹²C/¹³C ratios from CH₃ OH and CH₃ CN. This is in contrast to the ratios derived in previous works from COMs towards the Class 0 protostellar system IRAS 16293-2422. In this source, some O-bearing COMs show ¹²C/¹³C ratios as low as ours, however, the values from CH₃CN isotopologues are around the ISM value. The origin of differences between the various molecules among different star-forming regions is not yet understood and requires further observational and modelling efforts;

5. Astrochemical models that we conducted for potential precursor molecules of CH₃OH and CH₃CN predict ¹²C/¹³C ratios for CO and H₂CO that are close to or slightly below the local ISM value of 68 (in ice and gas). The ratios do not reach as low as in the observations. One reason may be the assumed duration of the prestellar phase in the model (10⁶ yr), but this requires further investigation. The gasphase ¹²C/¹³C ratios for CH₃CN precursors, in particular CN, HCN, and HNC, are predicted to be in the range 20-50 at the spatial scales probed by the observations (~300 au), more similar to the observational results for CH₃CN.

Our observational and modelling results indicate that there exists ¹³C enrichment at the Class 0 stage traced by COMs, which likely originates from the prestellar stage. There are some observational limitations, which we evaluated in detail, that come with the spatially and (mostly) spectrally unresolved PRODIGE data set and that can be addressed in future wideband observations at higher angular and spectral resolution. Nevertheless, our work opens the avenue to using the ¹²C/¹³C ratio to study the question of chemical inheritance, as protostellar and protoplanetary disks are likely to inherit at least some fraction of the parent cloud material.

In addition, our simulations predict that the ¹²C/¹³C ratios for simpler species, such as H₂CO, HCN, or HNC, greatly vary as a function of distance from the protostar (i.e. density and/or temperature) and depend on the duration of the preceding prestellar phase. An in-depth study of these isotopologues will shed light onto the dependence of the ¹²C/¹³C ratios on time and environment and may help put constraints on the chemical network of more complex species such as CH₃CN.

Acknowledgements

The authors thank the IRAM staff at the NOEMA observatory for their support in the observations and data calibration. This work is based on observations carried out under project number L19MB with the IRAM NOEMA Interferometer. IRAM is supported by INSU/CNRS (France), MPG (Germany) and IGN (Spain). L.A.B., J.E.P., O.S., P.C., M.J.M., C.G., M.T.V.M., Th.H., D.S., and Y.C. acknowledge the support by the Max Planck Society. A.F. and J.J.M.P. acknowledge support from ERC grant SUL4LIFE, GA No. 101096293. Funded by the European Union. Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council Executive Agency. Neither the European Union nor the granting authority can be held responsible for them. A.F. also thanks project PID2022-137980NB-I00 funded by the Spanish Ministry of Science and Innovation/State Agency of Research MCIN/AEI/10.13039/501100011033 and by “ERDF A way of making Europe”. Th. H. and D. S. acknowledge support from the European Research Council under the Horizon 2020 Framework Program via the ERC Advanced Grant Origins 83 24 28 (PI: Th. Henning). L.C. acknowledges support from the grant No. PID2022-136814NB-I00 by the Spanish Ministry of Science, Innovation and Universities/State Agency of Research MICIU/AEI/10.13039/501100011033 and by ERDF, UE.

Appendix A Additional tables

Table A.1 lists beam sizes and rms noise values for the seven sources of the PRODIGE sample studied in this work. Tables A.2–A.8 contain the parameters to derive the LTE radiative transfer models for each analysed molecule as well as the results from the population-diagram analysis. The literature ¹²C/¹³C ratios for all additional sources shown in Figs. 1 and 5 are listed in Table A.9 together with their respective reference and the species that were used to derive them.

Appendix B Additional figures: Integrated intensity maps

Figures B.1–B.7 show integrated intensity maps for selected transitions of the main and rarer isotopologues towards the seven sources of the PRODIGE sample that were analysed in this work. The transitions were selected because they are bright and unblended.

Fig. B.1 Integrated intensity maps using two transitions for both isotopologues of CH3 OH, the main CH3CN isotopologue, and one transition for each of the 13C isotopologue towards IRAS 2A. Respective rest frequencies (in GHz) are shown in the bottom-right corner, the beam (half power beam width, HPBW) on the bottom left. Multiplication factors (e.g. ×3, next to the species name) were applied to some species and transitions to be able to use the same intensity scale for all panels. The yellow dashed contour is at 5σ, where σ is the rms noise level (in K km s−1) and is written in the top-right corner.

Fig. B.2 Same as Fig. B.1 but for IRAS 4B.

Fig. B.3 Same as Fig. B.1 but for L1448-IRS2.

Fig. B.4 Same as Fig. B.1 but for B1-c.

Fig. B.5 Same as Fig. B.1 but for L1448C.

Fig. B.6 Same as Fig. B.1 but for SVS13A.

Fig. B.7 Same as Fig. B.1 but for B1-bS.

Appendix C Additional figures: Population diagrams

In Figs. C.1–C.8 all derived population diagrams are displayed.

Fig. C.1 Population diagrams for CH3 OH, 13CH3OH, CH3CN, 13CH3CN, and CH313CN towards IRAS 4B. Observed data points are shown with teal circles as indicated in the top-right corner of the left panel while the modelled data points from Weeds are shown as orange squares. No corrections are applied in the left panels while in the right panels corrections for opacity and contamination by other molecules have been considered for the observed and modelled populations. The results of the linear fit to the observed data points are shown in the right panels.

Fig. C.2 Same as C.1 but for L1448C.

Fig. C.3 Same as C.1 but for L1448-IRS2.

Fig. C.4 Same as C.1 but for B1-c.

Fig. C.5 Same as C.1 but for B1-bS.

Fig. C.6 Same as C.1 but for IRAS 2A.

Fig. C.7 Same as C.1 but for SVS13A.

Fig. C.8 Same as C.1 but for CH318OH towards IRAS 2A, IRAS 4B, SVS13A, and B1-c.

Appendix D Details of the chemical and physical models

Fig. D.1 Model predictions for selected gas-phase and ice 12C/13C ratios as a function of time in the first simulation stage, representing evolution in dark cloud conditions.

As noted in Sect. 4.5, we have performed the chemical simulations in two stages. In the first (dark cloud) stage, we have run a simple point (i.e. zero-dimensional) model to represent chemical evolution in a dark cloud, assuming canonical physical conditions n(H₂) = 10⁴ cm⁻³, T_gas = T_dust = 10K, A_V = 10 mag. The chemistry was let to evolve over a period of 10⁶ yr in these conditions. Figure D.1 shows the time evolution of the ¹²C/¹³C ratios of selected CH₃OH and CH₃CN precursors in the gas phase and in the ice. It is evident that the ratios not only depend strongly on the time but also on the species in question; similar behaviour is shown in our previous works (Colzi et al. 2020; Sipilä et al. 2023). As a result, the initial conditions of the chemistry in the prestellar stage is very sensitive to how long the chemical evolution has proceeded in the dark cloud stage. For example, if we chose a duration of 2 × 10⁵ yr for the dark cloud stage, then, the ¹²C/¹³C ratio of CO in the ice would (slightly) exceed 68, and we would expect a similar ratio for CH₃OH as it is a product of CO hydrogenation; assuming that no other processes impact the methanol ¹²C/¹³C ratio in the protostellar stage, a very short dark cloud stage would be implied. On the other hand, this would also lead to the two ¹³C variants of CH₃ CN having very different ¹²C/¹³C ratios in the protostellar stage (assuming that CH₃CN forms through CH₃ + CN). Figure D.1 also shows a dichotomy in the results in that some species are enriched in ¹²C while others are enriched in ¹³C; at 10⁶ yr, the two most abundant carbon carriers are H₂CO and CH₄ ice (not shown). The former has a ¹²C/¹³C ratio close to that of CO (62), while the latter is anti-fractionated with a ¹²C/¹³C ratio of 134.

For the protostellar stage we employ the physical model by Crimier et al. (2010) which gives the density and temperature profiles as a function of distance from the centre of the core. These profiles are displayed in Fig. D.2. We assume T_gas = T_dust, which is reasonable given the high volume density and hence efficient gas-dust thermal coupling. We note that this physical model is not appropriate for describing the very innermost regions of the core (inside a few hundred au) where substructure is known to be present but is applicable in the present context where the observations were not carried out using high angular resolution.

Fig. D.2 H2 density (left-hand y-axis) and kinetic temperature (Tgas = Tdust; right-hand y-axis) as a function of distance from the centre in the IRAS16293-2422 model (from Crimier et al. 2010).

The chemical simulations in the protostellar stage are carried out by selecting a set of points along the physical model profiles and running the chemical simulation separately for each of the selected points, taking the initial abundances obtained from the preceding dark cloud simulation. This very simple model setup hence neglects any physical evolution between the dark cloud and protostellar core stages, and assumes that the initial chemical conditions are equal across the protostellar core. While such an approach is certainly not realistic, it does allow us to obtain rough estimates of the gas-phase and ice ¹²C/¹³C abundance ratios across the protostellar core. Naturally, the simulation results are also in this case sensitive to the choice of evolutionary time when the abundances are extracted from the model output. We chose here an early time of 10⁴ yr to reflect the fact that the physical model is not describing well the very inner regions of the core where substructure has already formed.

References

All Tables

All Figures

Fig. 1 Literature 12C/13C ratios towards sources at different evolutionary stages using different isotopologues. Simple O- and N-bearing species (<6 atoms) share a light blue and light yellow colour, respectively, but have different markers. Black rectangles around the orange triangle (N) and coloured circles (O) indicate complex organics (>6 atoms). Arrows indicate lower limits. References to the values can be found in Table A.9.

In the text

Fig. 2 Overview of the PRODIGE data towards IRAS 4B and presentation of the spectral line analyses. Panel a: integrated intensity map for th CH3CN 123−113 transition at 220.708 GHz. The beam (half power beam width, HPBW) is shown in the bottom left. The dashed yellow contou is at 5σ, where σ is the rms noise level and is written in the top-right corner. Panel b: observed spectra towards the peak position, which is a (Δα, Δβ = (0″, 0″) (black cross in panel a), of CH3CN at low (~2 MHz, black) and high (~62.5 kHz, teal) spectral resolution, where the latte] was shifted by 4K for visualisation. The low spectral-resolution data were analysed in this work (see Sect. 3). The modelled CH3CN spectra from Weeds are shown in red (low resolution) and orange (high resolution), where the same input parameters were used for both. Panel c: population diagram for CH3CN. Teal circles show the observed data points while the modelled data points from Weeds are shown with orange squares (se Sect. 4.6 for possible origins of differences between the two). No corrections are applied in the left panel, while in the right panel, corrections fo opacity and contamination by other molecules have been considered (negligible effect in this case), and only transitions with optical depths <0.5 (i.e. all transitions in this case) are shown. The results of the linear fit to the observed data points are shown in the right panels.

In the text

	Fig. 3 Peak velocity, νpeak, and line width, FWHM, derived from 1D Gaussian fits (including 1σ error bars) to the high spectral-resolution (HR) PRODIGE data as a function of upper-level energy, Eu, for the vibrational ground states of CH3 OH (orange) and 13CH3OH (teal). The coloured dashed lines indicate the respective average values for νpeak and FWHM.
In the text

	Fig. 4 Rotational temperatures used for the Weeds modelling (crosses) and derived from the population diagrams (circles with 1σ error bars). The dashed lines indicate the median values considering all isotopologues of a species.
In the text

	Fig. 5 Same as Fig. 1 but including the 12C/13C ratios derived from the PRODIGE data in this work (light red shaded). Bolometric luminosities of the protostars increase from left to right. Black rectangles around markers indicate complex organic molecules. Arrows indicate lower limits.
In the text

	Fig. 6 Column density ratio between 13CH3OH and CH318OH using the Weeds column densities. The expected local ISM value of ~8 is indicated with a dashed line. Arrows show lower limits.
In the text

Fig. 7 Model predictions for selected gas-phase 12C/13C ratios as a function of distance from the protostar taken at a simulation time of 104 yr. The constant ratio of 68 is indicated by the horizontal dashed black line. The vertical dotted black line separates the spatial scales discussed in this paper (R < 300 au) from the larger scales, which are not largely discussed but shown for completeness.

In the text

Fig. B.1 Integrated intensity maps using two transitions for both isotopologues of CH3 OH, the main CH3CN isotopologue, and one transition for each of the 13C isotopologue towards IRAS 2A. Respective rest frequencies (in GHz) are shown in the bottom-right corner, the beam (half power beam width, HPBW) on the bottom left. Multiplication factors (e.g. ×3, next to the species name) were applied to some species and transitions to be able to use the same intensity scale for all panels. The yellow dashed contour is at 5σ, where σ is the rms noise level (in K km s−1) and is written in the top-right corner.

In the text

	Fig. B.2 Same as Fig. B.1 but for IRAS 4B.
In the text

	Fig. B.3 Same as Fig. B.1 but for L1448-IRS2.
In the text

	Fig. B.4 Same as Fig. B.1 but for B1-c.
In the text

	Fig. B.5 Same as Fig. B.1 but for L1448C.
In the text

	Fig. B.6 Same as Fig. B.1 but for SVS13A.
In the text

	Fig. B.7 Same as Fig. B.1 but for B1-bS.
In the text

Fig. C.1 Population diagrams for CH3 OH, 13CH3OH, CH3CN, 13CH3CN, and CH313CN towards IRAS 4B. Observed data points are shown with teal circles as indicated in the top-right corner of the left panel while the modelled data points from Weeds are shown as orange squares. No corrections are applied in the left panels while in the right panels corrections for opacity and contamination by other molecules have been considered for the observed and modelled populations. The results of the linear fit to the observed data points are shown in the right panels.

In the text

	Fig. C.2 Same as C.1 but for L1448C.
In the text

	Fig. C.3 Same as C.1 but for L1448-IRS2.
In the text

	Fig. C.4 Same as C.1 but for B1-c.
In the text

	Fig. C.5 Same as C.1 but for B1-bS.
In the text

	Fig. C.6 Same as C.1 but for IRAS 2A.
In the text

	Fig. C.7 Same as C.1 but for SVS13A.
In the text

	Fig. C.8 Same as C.1 but for CH318OH towards IRAS 2A, IRAS 4B, SVS13A, and B1-c.
In the text

	Fig. D.1 Model predictions for selected gas-phase and ice 12C/13C ratios as a function of time in the first simulation stage, representing evolution in dark cloud conditions.
In the text

	Fig. D.2 H2 density (left-hand y-axis) and kinetic temperature (Tgas = Tdust; right-hand y-axis) as a function of distance from the centre in the IRAS16293-2422 model (from Crimier et al. 2010).
In the text