© The Authors 2026

1 Introduction

Molecular clouds are known to be the birthplaces of stars. Their structure consists of a complex network of elongated filaments produced by the interplay of turbulence, magnetic fields, and self-gravity (e.g., André et al. 2014). The gravitational fragmentation of these filaments leads to the formation of dense substructures known as starless cores, which do not yet contain an embedded protostar (e.g., Myers 2009). As their central densities increase (n_H > 10⁵ cm⁻³) and their temperatures remain low (T < 10 K), a subset of these starless cores becomes gravitationally bound. Once the core becomes gravitationally unstable, it enters the pre-stellar phase before undergoing gravitational collapse (e.g., Crapsi et al. 2005).

At the pre-stellar stage, most of the material that will later form stars and planetary systems is already present. While this concept has been demonstrated for water through both observations and models, it has only been suggested by observations for complex organic molecules (Caselli & Ceccarelli 2012; Ceccarelli et al. 2014; Cleeves et al. 2014; Drozdovskaya et al. 2021; van Gelder et al. 2020; Jensen et al. 2021a; Jensen et al. 2021b; Scibelli et al. 2021; Scibelli et al. 2025). Starless and pre-stellar cores therefore provide a unique window into the initial conditions just before the collapse and are essential for constraining the chemical budget inherited by forming stars and planetary systems (Herbst & van Dishoeck 2009).

In dense cores, the chemistry is dominated by an interplay between gas-phase reactions and grain-surface processes. The gas-phase chemistry is mainly driven by ion-molecule reactions, which are efficient at the low temperatures of dense cores. With the evolution of the core, the central density increases, leading to a decrease in temperature and the freezing-out of gas-phase molecules onto dust grains, which thus form icy mantles. Freeze-out onto dust grains was first observed for CO and CS (e.g., Caselli et al. 1999; Tafalla et al. 2002; Tafalla et al. 2006). Despite this depletion, the remaining gas-phase molecules still trace the ongoing chemistry. The resulting changes in molecular abundances make species such as CO and deuterated molecules powerful tracers of core evolution. CO is the main destruction partner of ${\text{H}}_3^{+}$. The freeze-out of CO promotes the reaction of ${\text{H}}_3^{+}$ with HD to form H2D⁺, initiating efficient deuterium fractionation. A common tracer of deuterium fractionation is the N₂D⁺/N₂H⁺ ratio (e.g., Crapsi et al. 2005; Emprechtinger et al. 2009).

Sulfur is one of the most abundant elements in the Universe, with a cosmic S/H abundance of ∼1.6 × 10⁻⁵ (Asplund et al. 2021). Sulfur-bearing molecules are widely observed in star-forming regions and are used to probe their physical conditions, such as the presence of shocks (e.g., Dutrey et al. 2024; Esplugues et al. 2013; Hatchell et al. 1998; Bachiller & Pérez Gutiérrez 1997; Zhang et al. 2023). Despite all the efforts in observations and theory (e.g., Fuente et al. 2023; Wakelam et al. 2004), sulfur chemistry remains one of the least understood areas in astrochemistry. In particular, the total sulfur abundance observed in pre-stellar cores represents only a small fraction (<1%) of the cosmic sulfur abundance observed in diffuse clouds (e.g., Tieftrunk et al. 1994; Woods et al. 2015). While part of the sulfur reservoir could be locked up in ices or refractory grain material (e.g., Ruffle et al. 1999; Martín-Doménech et al. 2016; Laas & Caselli 2019), the solid-phase sulfur molecules identified so far account for only a few percent of the total sulfur budget, indicating that the main carrier of sulfur remains unknown; a part could be locked up in refractory sulfur allotropes (Shingledecker et al. 2020; Cazaux et al. 2022) and salts (Poch et al. 2020; Slavicinska et al. 2025) that are not detectable via rotational spectroscopy.

Observations of sulfur-bearing molecules hint that sulfur depletion is linked to the evolutionary stage of the core. Nagy et al. (2019) find that the young starless core L1521E exhibits higher sulfur abundance levels than the more dynamically evolved pre-stellar core L1544. The difference in abundance depends on the sulfur-bearing molecules, with C₃S being roughly ten times more abundant in the more dynamically evolved core compared to the starless core. This trend suggests that a significant fraction of sulfur is depleted as cores evolve, consistent with the predictions by Laas & Caselli (2019). Hily-Blant et al. (2022) confirmed this by using the observed NS/N₂H⁺ ratio to infer the abundance of the gas-phase atomic sulfur from a sample of four starless and pre-stellar cores, and support a scenario in which sulfur becomes progressively depleted as dense cores evolve.

Whether this behavior represents a general evolutionary trend or depends on the specific tracers used remains unclear. In this work, we observed various sulfur-bearing molecules in a sample of seven starless and two pre-stellar cores to investigate whether different evolutionary indicators trace sulfur depletion. Section 2 introduces the source sample, observational setup, instruments, and configurations used, and Sect. 3 outlines the derivation of column densities and the determination of the CO depletion factor. In Sect. 4 we compare the observed abundances and selected molecular ratios according to different evolutionary tracers. We compare the observational abundances with those from chemical models, spanning different initial sulfur abundances, in Sect. 5. Finally, in Sects. 6 and 7 we discuss the implications of our findings and conclude.

2 Observations

Using the IRAM 30 m telescope, situated on the Pico Veleta in Spain, we targeted the dust peaks of seven starless and two pre-stellar cores. All sources are located in the Taurus Molecular Cloud. Table 1 summarizes the observed objects, including their coordinates.

We collected the data using the Eight MIxer Receiver (EMIR) with the E090 band. The individual subbands were connected to the fast Fourier Transform Spectrometer (FTS), offering a total spectral coverage of 7.2 GHz and a resolution of 50 kHz (0.15 km/s at 100 GHz). We employed frequency switching with a frequency throw of ± 3.9 MHz and observed the core sample over a noncontinuous frequency range between approximately 82 GHz and 114.5 GHz.

During the observations (project IDs: 103-15, 120-22, 073-23, and 020-25), the weather conditions remained stable, with the optical depth (τ) consistently being between 0.1 and 0.2. The pointing of the source was checked at least every two hours and corrected on nearby planets or pulsars, with an accuracy of 1″ to 4″. The resulting spectra reached typical root mean square (rms) noise levels of 2–25 mK per 50 kHz channel, depending on the source and frequency. The spectral intensities are given in units of the antenna temperature, $T_{\text{A}}^{*}$. This translates to the main beam temperature T_MB via $T_{\text{MB}}\text{\hspace{0.17em}via\hspace{0.17em}}{T}_{\text{MB}}=F_{\text{eff}}/B_{\text{eff}}\cdot T_{\text{A}}^{*}$. Here, F_eff and B_eff represent the forward efficiency, which was assumed to be 95% for all transitions, and the beam efficiency, which changed with the frequency of the telescope. The beam efficiency was determined through interpolation using the data provided on the IRAM 30 m web page¹. For a few lines in Table 2, the updated list of efficiencies was used, namely: SO (2,2–1,1) in L1495A-S and L1495B, SO (2,3–1,2) in TMC2, SO₂ (3_1,3– 2_0,2) in TMC2, ¹³CS (2–1) in L1512 and L1517B, CCS (7,6–6,5) in L1495AS and L1495B, CCS (8,7–7,6) in CB23, H₂CS (3_1,3–2_1,2) in L1495A-S and L1495, and H₂CS (3_0,3–2_0,2) in TMC2.

Within the observed frequency range, we focused on the sulfur-bearing molecules. All detected species contain hydrogen, carbon, and/or oxygen, in addition to sulfur. Table 2 summarizes the detected molecules and their spectroscopic information: E_u represents the upper level energy, Q_rot the rotational partition function, g_u the upper level degeneracy, and A_ul the Einstein coefficient for spontaneous emission.

3 Derivation of the column density and CO depletion factor

3.1 Molecular column densities

The initial reduction of the observational data were performed using the GILDAS software (Pety 2005). To correct for the baseline caused by frequency switching, we fitted and subtracted a third-degree polynomial. The individual lines were then fitted with a Gaussian line profile using python. Figure 1 shows the spectra of L1495A-S with the Gaussian fits overlaid. The spectra of the remaining cores are shown in Appendix A. For the calculation of the total column density $N_{\text{tot}}^{\text{thin}}$, we employed the formula reported in Mangum & Shirley (2015) for the optically thin emission under the assumption of local thermal equilibrium (LTE) and a constant excitation temperature (T_ex) of 10 K: $$N_{\text{tot}}^{\text{thin}}=\frac{8\pi k_{\text{B}}{\nu }^2}{h{c}^3}\frac{Q_{\text{rot}}(T_{\text{ex}})}{g_{\text{u}}{A}_{\text{ul}}}\frac{J(T_{\text{ex}})}{J(T_{\text{ex}})-J(T_{\text{bg}})}{e}^{E_{\text{u}}/k_{\text{B}}{T}_{\text{ex}}}{\displaystyle \int T_{\text{mb}}\text{dv},}$$(1)

where ν represents the transition frequency, k_B is the Boltzmann constant, h is the Planck constant and c is the speed of light. $J(T)=\frac{h\nu }{k_{\text{B}}}{\left(e^{h\nu /k_{\text{B}}T}-1\right)}^{-1}$ is the Rayleigh-Jeans equivalent temperature and T_bg= 2.73 K is the temperature of the cosmic microwave background. When assuming a Gaussian line profile, the integrated main beam temperature was determined via $\frac{1}{2}\sqrt{\pi /ln(2)}\text{Δv}\cdot T_{\text{mb}}$, with ∆v describing the full width half maximum of the line. To account for any deviations of the optically thin regime, we determined the optical depth (τ) and scaled the column density with the optical depth correction factor (Goldsmith & Langer 1999) via $$N_{\text{tot}}=N_{\text{tot}}^{\text{thin}}\frac{\tau }{1-e^{-\tau }},$$(2)

with $$\tau =-\ln \left[1-\frac{T_{\text{mb}}}{f\left[J\left(T_{\text{ex}}\right)-J\left(T_{\text{bg}}\right)\right]}\right].$$(3)

Here, f describes the beam filling factor, which we set to 1 for all cores. For the excitation temperature, we assumed 10 K for all molecules. Although the innermost region of a core can be a few kelvin lower, the envelopes are well characterized by temperatures of around 10 K, where most of the sulfur-bearing species are located (e.g., Tafalla et al. 2006). Varying the excitation temperature by ± 1 K had a negligible effect on the derived column density for all molecules except SO, which showed a change of up to 10%. Given that the critical densities of the observed transitions (see Table 2) are lower or comparable to typical densities of the core, LTE provides a reasonable first-order approximation for describing the observed emission. While this may not hold for all transitions in all cores, the limited number of detected optically thin transitions per molecule does not provide enough constraints for a meaningful non-LTE analysis. We therefore used an approach with a uniform excitation temperature and LTE to ensure a homogeneous comparison across all species.

Table 3 reports the Gaussian fit parameters, the optical depth, and the corrected column densities for all the detected lines for L1495A-S. The information of the remaining cores can be found in Appendix B. For non-detections, we determined upper limits for the integrated intensity, which were further used to derive upper limits on the column density. The uncertainties for the integrated intensities and linewidths were taken from the Gaussian fit. For the column densities under the optically thin assumption, the uncertainties were determined by propagating the uncertainty of the integrated intensity. For the uncertainty of τ and the corrected total column density, we then used Gaussian error propagation.

CS, one of the most abundant sulfur-bearing species in starless and pre-stellar cores, shows clear signs of high optical depth (e.g., self-absorption) in L1495A-S, L1495A-N, and L1544 (Figs. 1, A.4, and A.8) and was therefore excluded from further analysis in these sources. In L1495B (Figure A.5), wing-like structures affect the CS (2–1) and SO (2,3–1,2) transitions, leading to their exclusion. The weak transitions owing to OCS (8–7), SO₂ (3_1,3–2_0,2) in TMC2 (Figure A.2), HDCS (3_0,3–2_0,2) in L1495 (Figure A.3), OCS (8–7) in L1495A-N (Figure A.4), and HCS⁺ (2–1) in L1517B (Figure A.7) show features resembling a second velocity component over one to two channels. Given their low S/R, these features were attributed to noise and do not influence the Gaussian fits. Therefore, only HCS⁺ in L1517B was excluded due to a linewidth, which is ∼40% broader than the other lines in that core. Additionally, four molecular transitions in each TMC2 and L1495A-S (Figs. A.2 and 1) show one-sided wings, which may indicate infall motions (e.g., Lee et al. 1999; Lee et al. 2001; Crapsi et al. 2005). The presence of both red- and blueshifted asymmetries, which vary between molecular tracers, suggests that the features could arise from different regions within the cores or reflect different evolutionary stages. To further constrain their origin, maps to understand the spatial distribution and spectra with higher frequency resolution would be necessary. As the features were not accounted for in the Gaussian fit, the line parameters from the fit are unaffected and were thus included in the further analysis. The transitions with an asterisk in Table B.1 mark the excluded lines.

Excluding the specified transitions, the average linewidths across the cores range from approximately 0.3 km/s to 0.45 km/s. The pre-stellar cores L1544 and TMC2 have the highest average linewidths, and L1512 and CB23 the narrowest. For most sources, the average linewidths agree within uncertainties; only the linewidths of L1544 and TMC2 do not agree with those of CB23, L1512, L1495A-S, and L1495B.

Fig. 1 Spectra of molecules observed toward the dust peak of the starless core L1495A-S (in black). The red line indicates the best fit to a Gaussian line shape model. For each panel, the molecule is indicated on the top right, with the transition stated below. CS shows signs of self-absorption; hence, no Gaussian fit was performed.

3.2 CO depletion factor

Throughout most of their evolution, dense cores maintain a temperature of around 10 K. When their central number densities exceed ∼10⁵ cm⁻³, molecular freeze-out becomes catastrophic, as gas-phase molecules collide with and stick to dust grains, drastically reducing their gas-phase abundances. This process, which reflects the evolutionary stage of the core, was first observed for CO (e.g., Caselli et al. 1999; Tafalla et al. 2002; Willacy et al. 1998). The degree of the CO freeze-out can be measured using the CO depletion factor f_d. This measure compares the observed CO abundance X(CO) with the CO abundance in the local ISM via (e.g., Emprechtinger et al. 2009) $$f_{\text{d}}(\text{CO})=\frac{X_{\text{ref}}(\text{CO})}{X(\text{CO})}.$$(4)

For the reference CO abundance X_ref(CO), we used 9.5 · 10⁻⁵ (Frerking et al. 1982). We derived the observed CO abundance from the C¹⁷O (J = 1–0) transition by comparing its total column density (see the derivation in Appendix C) with the H₂ column density of the core while assuming X(C¹⁶O)/ X(C¹⁷O) = 2044 (Wilson & Rood 1994; Penzias 1981). To determine the observed CO abundance, we then used $$X(C^{16}\text{O})=\frac{\text{X}({\text{C}}^{17}\text{O})\cdot 2044}{N({\text{H}}_2)}.$$(5)

For all the starless and pre-stellar cores in the sample, we adopted the H₂ column densities derived within a 40″ beam from Herschel/SPIRE observations as reported in Chantzos et al. (2018). The resulting derived CO depletion factors are shown in Table 1.

4 Results

To compare the abundances, we selected one transition per molecule across the core sample to ensure a consistent, nonredundant comparison. Whenever more than one transition was detected, for example CCS, the transition detected in most cores was selected. When multiple transitions had the same number of detections, for example, H₂CS, the transition with the higher observed intensity was taken. For SO, the N, J = 2,3–1,2 transition was used for all cores except L1495B due to optical depth effects. Because of high optical depth, CS was also excluded for four cores. Only five CS abundances remained, which limits the statistical significance of the analysis. Given that the isotopic species ¹³CS and C³⁴S, which we expect to trace similar regions as CS, are detected in all cores, we excluded the main isotopolog CS from the analysis. The transitions analyzed in this work are listed in Table 4.

The following section is divided into two parts: molecular abundances and selected abundance ratios used to probe chemical trends among sulfur-bearing species. In both parts, abundances and various molecular ratios are shown as a function of three evolutionary tracers: the H₂ column density (taken from Chantzos et al. 2018), the N₂D⁺/N₂H⁺ ratio (taken from Crapsi et al. 2005) and the CO depletion factor (derived in Sect. 3.2), in order to search for possible correlations. Potential trends were quantified with the Spearman rank correlation coefficient (including both detections and non-detections). For molecules or ratios showing a trend (|ρ| ≥ 0.6 and p ≤ 0.1), a leave-one-out analysis was performed to identify possible outliers. Note that the small sample size of cores limits the statistical significance of the results, and the Spearman rank coefficient does not account for uncertainties in the observed abundances. This analysis was also applied to the correlations between the evolutionary tracers themselves. No strong correlations are found within the sample, with the strongest trend observed between the H₂ column density and the N₂D⁺/N₂H⁺ ratio, with a Spearman rank coefficient of ρ = 0.57 and p = 0.11.

4.1 Abundances

The observational (fractional) abundances are defined as the ratio of the molecular column density and the column density of H₂. Table 4 shows the derived abundances for each molecule in the individual cores, with upper limits indicating non-detections. The highest abundances for all observed molecules are found in L1495A-S and L1495A-N, with L1495A-S showing the largest number of detected molecules together with L1544. L1517B exhibits the lowest abundances, and together with L1495, CB23, and L1512, displays the highest rate of non-detections. For these four cores, only 7 out of 12 molecules were detected.

4.1.1 H₂ column density

As the core evolves, the H₂ column density increases, making it a proxy for the evolutionary stage of starless and pre-stellar cores. Since molecular abundances are defined as the ratio between the molecular and the H₂ column density, these two quantities are not fully independent. Consequently, any observed trends in the H₂ column density reflect relative changes in molecular abundances with evolution rather than fully independent correlations.

Figure 2 shows the molecular abundances as a function of H₂ column density for all cores, with each subplot representing a molecule and hatched bars indicating upper limits. Overall, no clear trends are observed, as reflected in the Spearman rank correlation coefficient, which is generally around ρ ∼ 0.6. We find the strongest positive correlation for ³⁴SO with ρ = 0.65, which cannot be counted statistically significant. A leave-one-out analysis finds that excluding either L1495B or L1517B increases ρ to ∼0.8 (p ∼ 0.01). A similar, but weaker behavior is seen for SO₂. HDCS also shows a possible correlation of ρ = 0.79 (p = 0.02), when L1517B is left out. While none of them is statistically significant, we find a slight positive correlation among ³⁴SO, SO₂, and HDCS when L1517B is excluded.

Even though there are no clear monotonic trends, several molecules show similar abundance patterns across the cores. In particular, H₂CS and HDCS display a comparable relative behavior: the cores with a lower H₂ column density show more variation in abundance, whereas higher N(H₂) cores cluster more closely in their observed abundance. However, SO₂ and ³⁴SO also show a similar pattern in their abundances across the cores. Besides these pairwise similarities, we also observe significant variation among the cores, with L1517B consistently showing the lowest abundances and L1495A-S and L1495A-N displaying the highest abundances for all molecules.

We find that the cores with low N(H₂) values (L1495, CB23, L1512, and L1517B) display non-detections for most oxygen-bearing species, except SO. Only L1495B deviates from this trend. Even though it has a similar N(H₂) value, the core exhibits high abundances of oxygen-bearing species. Moreover, unlike the other cores embedded in the L1495 filament (L1495, L1495A-S, and L1495A-N), which exhibit similar abundances for all molecules, L1495B shows enhanced levels of oxygen-bearing sulfur molecules relative to carbon-bearing species. Given their proximity within the L1495 filament and their similar abundance levels, which agree within uncertainties for six out of ten molecules, it is likely that L1495A-N and L1495A-S have formed under similar physical and chemical conditions.

4.1.2 N₂D⁺/N₂H⁺ ratio

Due to the freeze-out of CO, the fraction of deuterated N₂H⁺ increases in dense cores. Since N₂D⁺ and N₂H⁺ form through analogous reactions, the N₂D⁺/N₂H⁺ ratio (taken from Crapsi et al. 2005) serves as a good proxy for the core evolution. Figure 3 shows the molecular abundances as a function of this ratio.

The N₂D⁺/N₂H⁺ ratios for all but one core are below 0.12. The exception, L1544, exhibits a ratio of 0.23, which is consistent with it being in a more advanced evolutionary state. Cores with lower N₂D⁺/N₂H⁺ ratios tend to have lower abundances of most oxygen-bearing species (except SO), similar to the behavior seen in the abundances as a function of the H₂ column density.

This is true for L1495, L1512, and L1517B. L1495A-N does not follow the trend and exhibits high abundances for OCS and ³⁴SO, despite having the lowest deuteration ratio.

Overall, the scatter in abundances makes it difficult to identify any clear trends with N₂D⁺/N₂H⁺. A weak negative trend is visible only for HCS⁺, but the low Spearman rank coefficient (|ρ| < 0.6 for all molecules) makes this behavior statistically insignificant.

Fig. 2 Molecular abundances as a function of the H2 column density. Each panel shows one molecule, with bars representing the derived abundances for the different cores and widths corresponding to half the uncertainty in N(H2). Hatched bars with a downward-facing triangle indicate upper limits (non-detections). Each panel also shows the Spearman rank correlation coefficient (ρ) and the corresponding p-value.

4.1.3 CO depletion factor

The chemical evolution of starless and pre-stellar cores can also be traced using the CO depletion factor, which provides a direct measure of CO freeze-out onto dust grains. Figure 4 shows the molecular abundances as a function of the CO depletion factor. Similar to the other evolutionary tracers, no clear trend is visible. L1495A-S, one of the cores with high sulfur abundances, exhibits the lowest CO depletion factor, while L1517B has the highest CO depletion and the lowest abundance, despite not being classified as a highly evolved core (Crapsi et al. 2005). Within the L1495 cores, CO depletion and molecular abundances decrease gradually from L1495A-S to L1495A-N, then to L1495, with similar relative values. In contrast, L1495B deviates from this behavior and exhibits higher abundances of oxygen-bearing sulfur molecules compared to the other cores from the L1495 filament relative to the carbon- and hydrogen-bearing sulfur species.

The Spearman rank correlation coefficient shows no significant monotonic trend in the CO depletion factor. Nevertheless, several molecules show a decreasing trend with increasing CO depletion. This can be seen, for example, in the case of the CS isotopologs, ¹³CS and C³⁴S, or OCS (excluding non-detections).

In contrast, the pure oxygen-bearing sulfur species, such as SO₂, SO, and its two isotopologs, exhibit a more complex behavior. A few molecules, such as H₂CS and C³⁴S, show two nearly parallel abundance sequences rather than a single evolutionary trend. Overall, these observations suggest weak negative trends with increasing CO depletion factor, but none are statistically significant.

4.2 Ratios

To probe various aspects of the chemistry in starless and pre-stellar cores, abundance ratios can be used. While some determine the dominant sulfur chemistry by comparing carbon- and oxygen-bearing molecules, others, such as the ³⁴SO/SO ratio, trace isotopic fractionation. Here, we focused on six ratios, which we divided into three categories: ratios that trace isotopic fractionation (HDCS/H₂CS, 13CS/C³⁴S, and 34SO/SO), ratios that compare carbon and oxygen chemistry (CCS/³⁴SO and C³⁴S/³⁴SO), and a ratio that probes the degree of carbon-chain growth within sulfur-bearing molecules (C₃S/CCS). The results are presented in Table 5 and are further analyzed as a function of three evolutionary tracers introduced in Sect. 4.1. Note that for all cores, except L1544 and L1495A-S, upper limits were included in the ratio calculations, so the derived ratios are presented as upper or lower limits, depending on whether the non-detected molecule appears in the numerator or denominator.

4.2.1 Ratios that trace isotopic and deuterium fractionation

Figure 5 shows the isotopic ratios HDCS/H₂CS, ¹³CS/C³⁴S, and ³⁴SO/SO as a function of the different evolutionary tracers. Among these, the HDCS/H₂CS column density ratio serves as a proxy for deuterium enrichment in the individual cores, and thus complements the N₂D⁺/N₂H⁺ ratio by tracing the deuteration within sulfur-bearing molecules.

In the core sample, the deuteration ratios for each evolutionary tracer agree within uncertainties and span values between 0.12 and 0.21, without showing any monotonic trend. This is also reflected in the Spearman rank correlation coefficient, which shows a flat line for the H₂ column density and the CO depletion factor.

The ³⁴SO/SO ratio can indicate the optical depth of the SO line: if SO becomes optically thick, the ratio will be below the canonical isotopic value. In the optically thin limit, the expected ³⁴SO/SO ratio is 1/22 = 0.045 (Wilson & Rood 1994). Most sources are consistent with this value. However, L1495B, L1495A-S, and L1544 display particularly high ratios. For these cores, this may indicate either an elevated ³⁴SO abundance or an underestimated SO abundance. Comparing the ³⁴SO/SO ratio against the different evolutionary tracers, there is no notable trend. Excluding L1495B, the Spearman rank coefficient indicates a nonsignificant positive trend with the H₂ column density, while a slightly negative trend is found with the CO depletion factor.

The idea behind the ¹³CS/C³⁴S ratio is similar to that of the ³⁴SO/SO ratio. Taking the canonical values of ¹²C/¹³C = 68 (Milam et al. 2005) and ³⁴SO/SO = 1/22 = 0.045 (Wilson & Rood 1994), we expect a ratio of ∼0.35, assuming an optically thin line. All cores agree with the expected value within their uncertainties, except L1495A-S and L1495, which show higher ratios. Here, a possible explanation is an enhanced ¹³CS abundance or an underestimated C³⁴S abundance. When analyzed as a function of the different evolutionary tracers, the ¹³CS/C³⁴S ratio shows no notable trend.

Fig. 3 Same as Figure 2 but for the N2D+/N2H+ ratio. Widths correspond to half the uncertainty in the N2D+/N2H+ ratio.

Fig. 4 Same as Figure 2 but for the CO depletion factor. Widths correspond to half the uncertainty in the CO depletion factor.

Fig. 5 Molecular abundance ratios probing isotopic fractionation as a function of the three evolutionary tracers. For clarity, the width of each bar reflects half the uncertainty in the corresponding evolutionary tracer. Hatched bars indicate ratios where one of the two abundances represented an upper limit (non-detections). Each panel also shows the Spearman rank correlation coefficient (ρ) and the corresponding p-value.

4.2.2 Ratios that trace carbon- and oxygen-bearing sulfur species

Carbon-chain chemistry is known as “early-type chemistry” in dense cores (Suzuki et al. 1992). As the core evolves, a larger fraction of carbon-bearing species is locked up in CO and frozen onto dust grains, reducing the gas-phase carbon and gradually slowing the formation of carbon chains. Consequently, oxygen-bearing molecules become relatively more abundant in the gas phase. Thus, comparing the carbon- and oxygen-bearing sulfur molecules, we would expect a decrease in both the ratios of C³⁴S/³⁴SO and CCS/³⁴SO as a function of the evolutionary tracer. While CCS/³⁴SO focuses more on the carbon chain chemistry, C³⁴S/³⁴SO traces the more general balance between the carbon- and oxygen-bearing sulfur molecules.

Figure 6 shows both ratios as a function of the H₂ column density, the N₂D⁺/N₂H⁺ ratio, and the CO depletion factor. For the CCS/³⁴SO ratio, a moderate negative trend is found with the H₂ column density (ρ = −0.58, p = 0.10). L1495B strongly influences this trend, as it shows a comparatively high abundance across all oxygen-bearing sulfur molecules. When excluding L1495B, the trend becomes more pronounced and yields a statistically significant anticorrelation (ρ = −0.81, p = 0.01). The Spearman rank coefficient shows a similar value for the N₂D⁺/N₂H⁺ ratio with ρ = −0.81, p = 0.01, which does not change when removing individual cores from the sample. For the CO depletion factor, no correlation was found.

The C³⁴S/³⁴SO ratio exhibits a similar behavior as the CCS/³⁴SO ratio: for the H₂ column density, a moderate anti-correlation is found (ρ = −0.60, p = 0.09). This correlation becomes more pronounced if L1495B is excluded from the analysis. For the N₂D⁺/N₂H⁺ ratio, the anticorrelation is even stronger (ρ = −0.86, p < 0.01) and remains mostly unchanged when individual cores are removed from the analysis. Again, no clear trend is observed with increasing CO depletion factor.

Fig. 6 Same as Figure 5 but for the carbon- and oxygen-bearing sulfur species.

4.2.3 Carbon chains

The C₃S/CCS ratio provides information on the length of the carbon chains in the core. While carbon chemistry is generally more prominent in the earlier stages, C₃S forms, among other pathways via the surface reaction CCS + C (Laas & Caselli 2019). Based on the column densities reported by Vastel et al. (2018) for L1544, we would expect a ratio of 0.2–0.3. This compares well to the derived ratio of 0.27 for the same core. For L1495A-S, the value is 0.14, which agrees with the idea that it is less evolved. For the remaining cores, C₃S is not detected.

Figure 7 shows the ratio as a function of evolutionary tracers. We observe an increase as a function of the CO depletion factor and the N₂D⁺/N₂H⁺ ratio, as suggested by the Spearman rank correlation coefficient. In contrast, for the H₂ column density, no trend is visible. Given the large number of upper limits, this trend must be interpreted with caution.

5 Single-point chemical models

To investigate whether our current knowledge of sulfur chemistry can explain the observed abundances and the lack of clear evolutionary trends, we employed gas-grain chemical simulations with the code pyRate, varying the initial sulfur abundance (X_i(S)) and identifying the best-fitting model for each core. The specific details of the model, including the chemical network, are described in Sipilä et al. (2015) and Sipilä et al. (2019) and are briefly summarized here. Note that the network also incorporates the sulfur reactions reported in Laas & Caselli (2019). We applied a 0D chemical model assuming a constant temperature of T_gas = T_dust = 10 K, a density of n(H₂) = 10⁴ cm⁻³, and a cosmic ray ionization rate of ζ_X = 1.3 · 10⁻¹⁷ s⁻¹. These values were chosen to represent the typical conditions in the outer envelopes of starless and pre-stellar cores, which are known to host sulfur-bearing species. We adopted three different values for X_i(S)= 8 · 10⁻⁶, 8 · 10⁻⁷, and 8 · 10⁻⁸, representing an un-depleted sulfur abundance and increasing degrees of depletion. The three values of X_i(S) are often referred to as high, intermediate, and low X_i(S) in the following text.

Fig. 7 Same as Figure 5 but for the C3S/CCS ratio.

5.1 Comparison to models

We evaluated the detected sulfur-bearing molecules at a chemical age of 10⁶ years, an approximate timescale for dense cores (Evans et al. 2009). Isotopes of sulfur, oxygen, and carbon are not included in the chemical network. Thus, HDCS was the only isotopolog considered in the models. Since CS was previously excluded due to a lack of data caused by optically thick effects, we inferred its observational abundance from its optically thin isotopologs. For cores where CS and both ¹³CS and C³⁴S were detected, the CS/¹³CS and CS/C³⁴S ratios were calculated and compared to the canonical elemental ratios in the local ISM of ¹²C/¹³C = 68 (Milam et al. 2005) and ³²S/³⁴S = 22 (Wilson & Rood 1994). On average, both ratios are around half of the canonical values. Since the ¹³CS line in CB23 exhibits a broader linewidth than the average of the core, and as the CS/C³⁴S value is better represented, we adopted C³⁴S together with the canonical isotopic ratio to determine the CS abundance. We divided the nine remaining sulfur-bearing molecules into three groups: O-bearing sulfur molecules (SO, SO₂, and OCS), C-bearing sulfur molecules (CS, CCS, and C₃S), and HC-bearing sulfur molecules (HCS⁺, H₂CS, and HDCS). Note that for this analysis, no upper limits were considered.

Figure 8 compares the observational abundances with the models at the three different X_i(S), evaluated at 10⁶ years. Molecules of the same group are presented in similar colors for better distinction. For oxygen-bearing molecules, the models generally over-predict the observed abundances, with SO and SO₂ differing by up to an order of magnitude across the cores. OCS is better reproduced and best represented by the intermediate X_i(S). In contrast, the carbon-bearing molecules CCS and C₃S are mostly under-predicted and thus best match the high X_i(S) model. CS acts as an outlier within the carbon group, with the best-fitting model varying between intermediate and high X_i(S). For the hydrogen-bearing species, HCS⁺ is generally under-predicted, whereas H₂CS and HDCS are well reproduced. In all cores, both H₂CS and HDCS agree with the same model, indicating that the models successfully reproduce the observed HDCS/H₂CS ratios.

5.2 Mean distance of disagreement

To determine which model best agrees with the observed abundances, we employed the mean distance of disagreement (D(t)) approach introduced in Iqbal & Wakelam (2018) and incorporated a method that accounts for observational uncertainties (Wakelam et al. 2006). We implemented this via $$D\left(t\right)=\frac{1}{N_{\text{obs}}}{\displaystyle \sum }_i\left|d_i\right|,$$(6)

where d_i describes the distance between the observational and model abundances and N_obs is the number of points taken into account. The individual distance element (d_i) is calculated differently depending on whether the observed abundance is higher or lower than the model value. We used d_i = log(X_{obs, min}) − log(X_model) for the case X_obs > X_model and d_i = log(X_model) − log(X_{obs, max}) for the case X_obs < X_model. If the observational and model abundances agree, d_i = 0. The smaller the resulting number, the greater the agreement between the model and the data.

Table 6 represents the mean distance of disagreement for the individual cores evaluated at 10⁶ years with different X_i(S). For six cores, the best model (indicated by the bolded value in Table 6) is the one with intermediate X_i(S). The remaining three fit best with the depleted (low) X_i(S). Regarding the individual groups, the results are similar to those shown in Figure 8: while the O-bearing group matches the lowest X_i(S), the carbon-bearing group generally agrees best with the highest, except for L1517B, where the intermediate X_i(S) reproduces the observations the best. Only the HC-bearing group is consistent with different X_i(S) across the cores. For this group, the model that best reproduces the observations is generally the one with the lowest mean distance of disagreement. The only exception is L1495B. While the HC-bearing group matches the lowest X_i(S), considering all molecules, the intermediate sulfur value fits best according to the mean distance of disagreement.

6 Discussion

We observe clear variations in the molecular abundances across the cores, yet no statistically significant trends with the three evolutionary tracers. L1495B clearly stands out for its high abundances of oxygen-bearing sulfur species. Conversely, L1517B exhibits low abundances and a high CO depletion factor, even though it was classified as a starless core (Crapsi et al. 2005). This might be due to the dense and quiescent conditions under which the core formed, rather than an advanced evolutionary stage. Additionally, we find that ratios of carbon- and oxygen-bearing molecules decrease with increasing H₂ column density and N₂D⁺/N₂H⁺, but show no notable trend when compared to the CO depletion factor. This suggests that the environment or the location of the molecular emission strongly influences the observed chemistry.

Indeed, most sulfur-bearing molecules have their emission peak not at the dust continuum maxima, but rather in the outer envelope of the cores (e.g., Spezzano et al. 2017; Tafalla et al. 2006). As the outer layers are less tightly coupled to the evolutionary state of the core’s center, it might explain why no clear trends between the sulfur abundances and the adopted evolutionary tracers were observed. Conversely, the outer envelopes are more directly affected by their environment, which can explain the variation in the observed abundances between the cores (Fuente et al. 2023).

Several observational surveys have cataloged young stellar objects and energetic sources in the Taurus Molecular Cloud and thus provide additional context on the environments of the different filaments. These include, for example, catalogs of young stars from the Spitzer Space Telescope (Rebull et al. 2010) and the XMM-Newton Extended Survey of the Taurus molecular cloud (XEST; Güdel et al. 2007), which identified X-ray emitting sources. Notably, the XEST survey locates the strongest X-ray source of the Taurus Molecular cloud, V410 Tau ABC (23-032), at a projected distance of 0.01 pc from L1495A-N (assuming a distance of 129.9 pc for the L1495 filament; Galli et al. 2019). At such a close projected distance, X-rays can locally increase the ionization rate and the cascade of ion-molecule chemistry. Under favorable conditions, it elevates the abundance levels of ${\text{H}}_3^{+}$, N₂H⁺, and HCS⁺ while also contributing to the destruction of CO (Stäuber et al. 2005). To evaluate the potential impact on our cores, we estimated the X-ray ionization rate produced by V410 Tau at the positions of the two closest cores, L1495A-N and L1495A-S, by integrating the attenuated stellar X-ray photon flux over energy using photo-absorption cross sections and accounting for secondary electron ionizations (Dalgarno et al. 1999). The resulting ionization rates, summarized in Table 7 for different assumed N(H₂) values between source and core, indicate that for a typical column density of N_H = 10²² cm⁻², X-rays clearly dominate the ionization in the vicinity of L1495A-N, while their effect becomes negligible at the larger distance of L1495A-S.

Interestingly, L1495A-N shows the highest abundances of several sulfur-bearing molecules in our sample, including HCS⁺, CCS, C³⁴S, H₂CS, and HDCS. Motivated by the strong X-ray influence of V410 Tau on L1495A-N, we increased the ionization rate in our chemical models by one and two orders of magnitude to test whether this improves the agreement with the observed abundances across the cores. Assuming the low X(S_i), increasing the ionization rate by one order of magnitude improves the overall agreement between modeled and observed abundances for all cores, which may reflect that the different molecules trace different layers of gas within the cores. The effect differs for the different sulfur species: oxygen-bearing sulfur molecules are generally better reproduced with the standard ionization rate, whereas carbon-bearing sulfur species show better agreement with the enhanced rate. The HC-group behaves differently with HCS⁺ better reproduced at the one order of magnitude higher ionization, while H₂CS and HDCS depend on the core. However, at higher ionization rates, the models no longer reproduce the observed D/H ratios, which were well reproduced by the standard rate. The differences between the standard and enhanced ionization models are particularly high for L1495A-N and L1495A-S, where the higher ionization rate provides a better fit. In contrast, increasing the standard ionization rate by two orders of magnitude leads to deterioration between the model and observations, as all molecules in all cores are systematically under-predicted, particularly SO₂ and OCS.

In addition to X-ray sources, nearby protostars can impact the sulfur chemistry through heating, shocks from accretion/outflows, or UV radiation (e.g., Codella et al. 2025). Using the catalog presented in Rebull et al. (2010), we find that the cores L1512, TMC2, and L1517B are relatively isolated, with no protostar within a projected distance of ∼0.25 pc. In contrast, the L1495 filament, especially the northern part, where the cores of this sample are located, contains many protostars. This part has been identified as one of the most evolved substructures within the filament, due to its high number of protostars (Hacar et al. 2013). For most of the starless cores in this filament, one or more protostars are found within a projected radius of 0.2 pc, except for L1495B, which is more isolated (the closest protostar is at ∼0.31 pc). The proximity of young stars may therefore contribute to the enhanced sulfur abundance levels in L1495A-N and L1495A-S.

Because L1495B is more shielded from external UV radiation, it exhibits high abundances of oxygen-bearing sulfur species. A higher degree of shielding favors the conversion of atomic carbon into CO, thus reducing the availability of free carbon. At the same time, enough atomic oxygen remains to sustain oxygen-driven reactions, thus allowing molecules such as SO or SO₂ to form more efficiently (Lattanzi et al. 2020). The L1495 cores thus demonstrate that local conditions can cause differences in sulfur abundance, even though they are separated by only 0.3 pc.

Within the core sample, L1517B stands out for its relatively low abundances of sulfur-bearing species and its high CO depletion factor. Hacar & Tafalla (2011) has studied the L1517 region in further detail and found it to be quite quiescent and dense, with no embedded protostars. The gas surrounding the starless cores exhibits subsonic velocity fields, indicating the absence of strong shocks. The generally low abundance of sulfur-bearing molecules, as well as the numerous non-detections of oxygen-bearing sulfur species in L1517B, support the idea of a quiet environment. The high CO depletion factor and the low N₂D⁺/N₂H⁺ ratio are likely a consequence of physical and chemical evolution not always progressing simultaneously and being influenced by local environmental effects. Therefore, the CO depletion factor alone does not indicate that this core is more evolved than the other starless cores in the sample.

Regarding the molecular ratios, the two ratios comparing carbon and oxygen chemistry show a dependence on two of the three evolutionary parameters: the H₂ column density and the N₂D⁺/N₂H⁺ ratio. As these evolutionary tracers increase, the ratios decrease. This agrees with the general trend that carbon chemistry dominates in less evolved cores, while oxygen chemistry becomes more important at later stages and is further confirmed by the non-detections in oxygen-bearing sulfur species for cores with low N(H₂) and N₂D⁺/N₂H⁺. A possible reason the ratios comparing carbon- and oxygen-bearing sulfur species do not correlate with the CO depletion factor could be that these ratios are more sensitive to specific chemical pathways that do not strictly trace the CO depletion factor.

The chemical models considered density and temperature to represent the core’s envelope and showed variations in the best-fitting X_i(S) across the molecules. While OCS is reproduced within an order of magnitude for the best fitting X_i(S), other species, such as SO, SO₂, and C₃S, are consistently over-or under-predicted regardless of the model. Additional tests, in which the time step for extracting the modeled abundances was reduced to 10⁵ years (before 10⁶ years), showed a change in behavior: the oxygen-bearing species better agree with the models, whereas the carbon-bearing molecules are over-predicted. This indicates that simply adjusting the timescale cannot reproduce both the carbon- and oxygen-bearing sulfur molecules simultaneously.

A similar behavior was found in Vastel et al. (2018), who found that no single value of the elemental sulfur abundance matched the observations of 21 sulfur-bearing molecules simultaneously in the pre-stellar core L1544. When varying the time in the models, we find that carbon- and oxygen-bearing sulfur species are reproduced differently, indicating that no single chemical age can simultaneously reproduce the abundances of all observed molecules. An alternative explanation is that the assumed C/O ratio in the models does not accurately describe the chemical conditions in the envelopes of dense cores. Byrne et al. (2026) show that even slight variations in the C/O ratio can strongly change the relative abundances of carbon- and oxygen-bearing molecules. Such sensitivities imply that simplified model assumptions might not be sufficient. Instead of using 0D chemical models, a 1D or higher-dimensional model would likely better capture the core’s physical structure and thus provide more accurate results, as seen for example in Jensen et al. (2023), Priestley et al. (2023), and Priestley et al. (2024). Given the limited observational constraints on the formation pathways of sulfur-bearing species, such as non-detections or low S/R, our results further indicate that our current understanding of sulfur chemistry is incomplete.

Fig. 8 Comparison between the observational data and the model abundance at three different Xi(S). The colored bars represent the observations, while the different symbols show the models evaluated at 106 years.

7 Conclusion

In this work, we have presented single-point observations of 13 different sulfur-bearing molecules toward seven starless and two pre-stellar cores in the Taurus Molecular Cloud. By comparing the abundances and selected molecular ratios with three evolutionary tracers, we tested how sulfur chemistry behaves across cores at different evolutionary stages. Additionally, the CO-depletion factors for CB23, L1495A-N, L1495A-S, L1495B, and L1512 were determined for the first time.

The observed abundances show notable variations between the cores and display weak or absent trends with respect to three evolutionary tracers, N(H₂), N₂D⁺/N₂H⁺, and the CO depletion factor. This agrees with the idea that sulfur-bearing molecules mainly trace the outer envelopes of the cores rather than the central regions, which are better represented by these tracers. Two molecular ratios that probe the difference between carbon-and oxygen-bearing sulfur species, C³⁴S/³⁴SO and CCS/³⁴SO, exhibit a more pronounced dependence on two of the three evolutionary indicators. This indicates that, while individual sulfur-bearing species are weak tracers of evolution, the relative abundances of carbon- and oxygen-bearing sulfur molecules are sensitive to both physical and chemical conditions. Additionally, we ran 0D chemical models with varying initial sulfur abundances. At 10⁶ years, the models reproduce the observed abundances for some species, such as OCS, H₂CS, and HDCS, within an order of magnitude, but cannot replicate the abundances of all sulfur-bearing molecules simultaneously. In particular, the oxygen-bearing molecules SO and SO₂ are consistently over-predicted, while carbon-bearing species, such as CCS and C₃S, are generally under-predicted. This difference between the carbon- and oxygen-bearing becomes even more prominent when varying the chemical age, reinforcing that, with the current 0D chemical models, no single chemical age can simultaneously describe all the observed sulfur-bearing species.

The results indicate that no single evolutionary parameter can describe the sulfur chemistry in starless and pre-stellar cores. Instead, the chemistry is influenced by a combination of the chemical age and different local environmental conditions. We find that nearby protostars or X-ray sources can measurably affect the detected sulfur abundance levels, leading to clear differences in abundance levels between sources. A single evolutionary tracer or model parameter cannot capture these variations.

The inability of 0D chemical models to simultaneously reproduce the observed sulfur abundances of all detected molecules demonstrates both their limitations and the complexity of sulfur chemistry. While the models can capture general chemical trends, the 0D structure limits the performance. It thus fails to account for variations in the physical conditions of individual cores or the influence of the local environment. Additionally, several formation and destruction pathways of sulfur-bearing molecules remain unknown under cold-core conditions. Together, this suggests that the sulfur chemistry in dense cores is a combination of the physical structure, environmental effects, and incomplete chemical knowledge. Future studies should thus focus on targeted observations of individual cores in well-characterized regions, such as L1517, combined with more complex modeling approaches to obtain a more comprehensive picture of sulfur chemistry in dense cores.

Data availability

Table B.1 is available at the CDS via https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/710/A251.

Acknowledgements

The authors wish to thank the anonymous referee for their comments, which helped to improve the clarity and quality of this work. L.S., S.S., and S.S.J. wish to thank the Max Planck Society for the Max Planck Research Group funding.

References

Appendix A Plots of observed spectra

Figures A.1 - A.8 show the observed spectra of all sulfur-bearing molecules including a Gaussian fit toward all sources in the sample.

Fig. A.1 Same as Figure 1 but for the starless core CB23.

Fig. A.2 Same as Figure 1 but for the pre-stellar core TMC2.

Fig. A.3 Same as Figure 1 but for the starless core L1495.

Fig. A.4 Same as Figure 1 but for the starless core L1495A-N. CS shows signs of self-absorption; hence, no Gaussian fit was performed.

Fig. A.5 Same as Figure 1 but for the starless core L1495B.

Fig. A.6 Same as Figure 1 but for the starless core L1512.

Fig. A.7 Same as Figure 1 but for the starless core L1517B.

Fig. A.8 Same as Figure 1 but for the pre-stellar core L1544. CS shows signs of self-absorption; hence, no Gaussian fit was performed.

Appendix B Line list of the core sample

Table B.1, showing the fitting parameters, the optical depth, and column density for the remaining cores, is available in electronic form at the CDS.

Appendix C Column densities with hyperfine structure

Since the ¹⁷O isotope has a nonzero nuclear spin, the C¹⁷O (1–0) transition is split into three different hyperfine components. The magnetic and electric quadrupole moments of the nucleus interact with the rotational motion of the molecule, which causes the hyperfine splitting. The CLASS software contains a method for analyzing these hyperfine components, providing information about their positions and relative intensities. The fit provides the optical depth as the sum of the individual components, with a resulting value of τ = 0.1 indicating an optically thin line. Any line with a τ above 0.1 cannot automatically be assumed optically thin, which is why the line has to be corrected with the optical depth correction factor (see Equation 2).

Six out of the nine cores show a C¹⁷O (J = 1–0) transition that can be considered optically thin. For those, we used a three-component Gaussian fit to retrieve the line information. To calculate the column density, we used the sum of the individual areas from the three components as integrated intensity and the line parameters of the un-split C¹⁷O (J = 1–0) transition together with Equation 1 under the optically thin assumption. For the excitation temperature, we used 10 K, as in the previous analysis. The uncertainties were determined by propagating the integrated-intensity uncertainty.

Table C.1 shows the resulting line parameters for the individual components and the total column density for the transition. As all three components stem from the same molecule, each transition should have a similar linewidth. This was the case for all cores except L1495A-N and L1495B, for which the linewidth had to be fixed to the core’s average to ensure consistency. Additionally, L1495B shows a peculiar line feature that resembles a second velocity component. As none of the sulfur-bearing molecules in that core exhibit such a feature, we attempted to fit the three main components as accurately as possible. The resulting value should thus be seen as an approximation.

For CB23, L1495A-S, and L1544, the fitted optical depths exceed 0.1; therefore, we needed to correct them. When exhibiting non-optically thin hyperfine structure, we cannot simply assume a constant excitation temperature, as was done in the previous analysis. Instead, the hyperfine structure (hfs) fit allows us to calculate T_ex from the fit parameters themselves via $$T_{\text{ex}}=\frac{h\nu /k_{\text{B}}}{ln\left(\frac{h\nu /k_{\text{B}}}{(T_{\text{MB}}\cdot \tau )/({\tau }_{\text{main}}f)+J(T_{\text{bg}})}+1\right)},$$(C.1)

where T_MB · τ is the product of the intensity and the optical depth summed over all three components. In the hfs fitting, this corresponds to the first output parameter, while τ_MB corresponds to the last. The calculated excitation temperature can then be used to determine the total column density via $$N_{\text{tot}}={\tau }_{\text{main}}\sqrt{\frac{16{\pi }^3}{\ln (2)}}\frac{{\nu }^{3}Q(T_{\text{ex}})\text{Δv}}{c^3{A}_{\text{u}}{g}_{\text{u}}}\frac{e^{\frac{E_{\text{u}}}{k_{\text{B}}{T}_{\text{ex}}}}}{e^{\frac{h\nu }{k_{\text{B}}{T}_{\text{ex}}}}-1}.$$(C.2)

Table C.2 shows the results of the hfs fitting and the column densities for C¹⁷O in CB23, L1495A-S, and L1544. The uncertainties were determined using the Gaussian error propagation method. The resulting fits are presented in Figure C.1.

Fig. C.1 Spectra of the C17O line observed toward the different sources in black. The red line indicates the best fit line for the specified source in the upper-left corner of each panel.

All Tables

All Figures

	Fig. 1 Spectra of molecules observed toward the dust peak of the starless core L1495A-S (in black). The red line indicates the best fit to a Gaussian line shape model. For each panel, the molecule is indicated on the top right, with the transition stated below. CS shows signs of self-absorption; hence, no Gaussian fit was performed.
In the text

Fig. 2 Molecular abundances as a function of the H2 column density. Each panel shows one molecule, with bars representing the derived abundances for the different cores and widths corresponding to half the uncertainty in N(H2). Hatched bars with a downward-facing triangle indicate upper limits (non-detections). Each panel also shows the Spearman rank correlation coefficient (ρ) and the corresponding p-value.

In the text

	Fig. 3 Same as Figure 2 but for the N2D+/N2H+ ratio. Widths correspond to half the uncertainty in the N2D+/N2H+ ratio.
In the text

	Fig. 4 Same as Figure 2 but for the CO depletion factor. Widths correspond to half the uncertainty in the CO depletion factor.
In the text

Fig. 5 Molecular abundance ratios probing isotopic fractionation as a function of the three evolutionary tracers. For clarity, the width of each bar reflects half the uncertainty in the corresponding evolutionary tracer. Hatched bars indicate ratios where one of the two abundances represented an upper limit (non-detections). Each panel also shows the Spearman rank correlation coefficient (ρ) and the corresponding p-value.

In the text

	Fig. 6 Same as Figure 5 but for the carbon- and oxygen-bearing sulfur species.
In the text

	Fig. 7 Same as Figure 5 but for the C3S/CCS ratio.
In the text

	Fig. 8 Comparison between the observational data and the model abundance at three different Xi(S). The colored bars represent the observations, while the different symbols show the models evaluated at 106 years.
In the text

	Fig. A.1 Same as Figure 1 but for the starless core CB23.
In the text

	Fig. A.2 Same as Figure 1 but for the pre-stellar core TMC2.
In the text

	Fig. A.3 Same as Figure 1 but for the starless core L1495.
In the text

	Fig. A.4 Same as Figure 1 but for the starless core L1495A-N. CS shows signs of self-absorption; hence, no Gaussian fit was performed.
In the text

	Fig. A.5 Same as Figure 1 but for the starless core L1495B.
In the text

	Fig. A.6 Same as Figure 1 but for the starless core L1512.
In the text

	Fig. A.7 Same as Figure 1 but for the starless core L1517B.
In the text

	Fig. A.8 Same as Figure 1 but for the pre-stellar core L1544. CS shows signs of self-absorption; hence, no Gaussian fit was performed.
In the text

	Fig. C.1 Spectra of the C17O line observed toward the different sources in black. The red line indicates the best fit line for the specified source in the upper-left corner of each panel.
In the text