© The Authors 2023

1 Introduction

In the early stages of the low-mass star formation process the envelope, disk, and outflow components play a key role. Material infalls from the envelope to the disk, part of this material is accreted from the protoplanetary disk into the central protostar and part of the parental cloud material goes on to be ejected through energetic jets and outflows (Reipurth & Bally 2001; Pudritz et al. 2007). Tychoniec et al. (2020) and Cridland et al. (2022) proposed that the first planetary embryos start to form in the embedded Class 0/I stages, therefore, the chemical content of the early disk will potentially influence the natal composition of emerging planets. In this sense, it is essential to identify and study the chemical composition of embedded Class 0/I sources.

One of the current enigmas is the sulfur chemistry, which is poorly understood in the process of low-mass star and planet formation. The sulfur content in dense molecular gas was found to be depleted by several orders of magnitude in comparison to diffuse clouds (Ruffle et al. 1999; Laas & Caselli 2019; Rivière-Marichalar et al. 2019), suggesting that most of the sulfur is locked into dust grains (Ruffle et al. 1999; Laas & Caselli 2019; Kama et al. 2019; Shingledecker et al. 2020; Cazaux et al. 2022). More recently, Fuente et al. (2023) proposed that this sulfur-depletion in the gas-phase might depend on the environment, finding a larger depletion factor for Taurus and Perseus, in comparison with Orion A. In the solid phase, it is also unclear what is the main sulfur carrier. Laas & Caselli (2019) predicted that most of the missing sulfur is in the form of organo-sulfur species, while Kama et al. (2019) proposed that the main sulfur carrier is much more refractory than water ice, consistent with sulfide minerals such as FeS. In addition, Shingledecker et al. (2020) and Cazaux et al. (2022) concluded that sulfur-chain molecules, such as S₈, should be abundant in dust grains. These results are not reflected in the detection of sulfur-bearing molecules in interstellar ices, where OCS is the only S-bearing species unambiguously detected (Geballe et al. 1985; Palumbo et al. 1995; McClure et al. 2023). SO₂ has been tentatively detected by Boogert et al. (1997), Yang et al. (2022), and McClure et al. (2023), while the dominant S-bearing ice in comets (H₂S; Calmonte et al. 2016) remains undetected (Jiménez-Escobar & Muñoz Caro 2011; McClure et al. 2023). Therefore, the main carrier of sulfur is still unknown in both, the gas and the solid-phase.

CS, SO, and SO₂ are simple S-bearing species commonly detected in the gas phase toward embedded sources (Artur de la Villarmois et al. 2019; Tychoniec et al. 2021). CS emission typically traces the envelope material and outflow cavity walls in young Class 0 sources (Tychoniec et al. 2021), SO emission has been detected in outflow cavity walls and the disk component of Class 0/I sources (Sakai et al. 2014; Harsono et al. 2014; Tychoniec et al. 2021) and SO₂ is commonly seen in compact emission toward Class 0 sources and has been proposed to trace accretion shocks at the disk-envelope interface in more evolved Class I sources (Oya et al. 2019; Artur de la Villarmois et al. 2022). In addition, Artur de la Villarmois et al. (2019) studied a sample of 10 Class I sources in Ophiuchus and found that SO₂ was only detected in sources with high bolometric luminosities (L_bol ≥ 1.4 L_⊙). In contrast to embedded sources, a handful of more evolved Class II disks show emission of CS, H₂CS (Le Gal et al. 2019), SO (Rivière-Marichalar et al. 2020; Booth et al. 2023; Huang et al. 2023), and H₂S (Phuong et al. 2018; Rivière-Marichalar et al. 2021, 2022), while SO₂ and CCS emission has only been detected in IRS 48 (Booth et al. 2021) and GG Tau (Phuong et al. 2021), respectively. In addition, only one sulfur-bearing complex organic molecule (COM) has been detected so far (CH₃SH toward the Class 0 source IRAS 16293-2422; Drozdovskaya et al. 2018). The chemical evolution of S-bearing species from Class 0 to Class II disks is not very well understood and there is a lack of a statistically unbiased survey of sulfur-bearing molecules toward embedded Class 0/I sources, where they are more easily detected.

In this work, we present the analysis of the simple sulfur-bearing molecules, CS, SO, ³⁴SO, and SO₂ toward a sample of 50 Class 0/I sources in the Perseus star-forming region. A special focus in the analysis of the SO₂ emission is given, as this molecule has not been observed as much as the others and it usually shows centrally peaked emission (not commonly seen in outflow cavity walls or extended structures). The Perseus molecular cloud is one of the most active nearby star-forming regions, with distances between 234 and 331 pc (Zucker et al. 2020). The data analyzed in this paper are part of the PErseus ALMA CHEmical Survey (PEACHES; Zhang et al. 2018, 2023; Yang et al. 2021), with an average angular resolution of 0″ 6 in ALMA band 6. A more detailed study of COM detections toward the PEACHES sample is presented in Yang et al. (2021) and the connection between sulfur-bearing species and dust polarization of nine of the sources has been studied by Zhang et al. (2023).

Section 2 describes the observational procedure, calibration, and the parameters of the observed molecular transitions. Results are presented in Sect. 3, with moment 0 maps, spectra, and detection rates. Section 4 is dedicated to the analysis and discussion, where fluxes and column densities are analyzed as a function of the protostellar physical properties, SO₂ abundances are compared with other sources in the literature and the detectability of COMs is investigated in those sources where S-bearing molecules are also detected. Finally, we end with a summary in Sect. 5.

2 Observations

The data analyzed in this work are part of the PEACHES survey (Yang et al. 2021) observed with ALMA (project codes: 2016.1.01501.S and 2017.1.01462; PI: N. Sakai). J0237+2848 was used as bandpass calibrator, while J0238+1636 and J0336+3218 were employed as flux and phase calibrators, respectively. The integration time corresponds to 10 min per source in ALMA band 6. The molecular transitions analyzed in this paper are listed in Table 1 and they are the only molecular lines detected in the PEACHES sample that contain sulfur. Self-calibration to the phase and amplitude was performed in most of the sources using the continuum data after the normal calibration and applied to the molecular line data. The exceptions were Per-emb 15, Per-emb 20, Per-emb 54, Per-emb 60, L1448 IRS2E, and L1455 IRS2, given that the continuum emission toward these sources was not bright enough to perform self-calibration.

A robust parameter of 0.5 and an auto-mask option were employed in the cleaning process, and primary beam responses of 0.1 were used to maximize the field of view. The synthesized beam is about 0″ 6 × 0″.4 averaged across all spectral windows and, assuming a distance of 300 pc to the Perseus molecular cloud, this is equivalent to 180 × 120 au. The final spectral resolution is 0.14 km s⁻¹ (121 kHz) for both SO transitions, 0.15 km s⁻¹ (122 kHz) for the CS and SO₂ transitions, and 1.2 km s⁻¹ (987 kHz) for the ³⁴SO transition. The source sample is listed in Table A.1, with other common names, coordinates, and physical properties. More detailed information on the PEACHES survey and other observed molecular transitions, as well as the calibration strategy is presented in Yang et al. (2021).

3 Results

In this section, we present moment 0 maps of CS, SO, ³⁴SO, and SO₂, integrated spectra, detection rates, the physical components where the emission of the different S-bearing species are observed, continuum and molecular fluxes, and column density calculations. A 2D Gaussian was employed to fit the continuum emission, retrieving the coordinates and continuum fluxes listed in Table A.1. Continuum emission maps are not presented in this paper but can be found in Fig. 1 of Yang et al. (2021), as well as convolved sizes, deconvolved sizes, and positions angles (PA) of 2D Gaussian fits (Table 4). Comparing the beam size with deconvolved sizes, most of the sources are not resolved. The only resolved sources in the continuum are Peremb 2, Per-emb 11 B/C, Per-emb 12 A/B, Per-emb 13, Per-emb 20, and Per-emb 44.

Fig. 1 Class 0 sources. Moment 0 maps of CS, SO, 34SO, and SO2. The contours start at 3σ, follow steps of 1σ, and the emission has been integrated over the velocity ranges listed in Table B.1. Maximum values (in units of Jy beam−1 km s−1) are indicated in the bottom-right corner of each panel. Red and blue arrows represent the outflow direction and the yellow star indicates the position of the continuum peak. Magenta stars represent binary components and the synthesized beam for each molecular transition is shown in the panels of the first row. Sources with no detections are not shown in this figure.

3.1 Moment 0 maps

Moment 0 maps of CS, SO, ³⁴SO, and SO₂ were created by integrating the data cubes over the velocity ranges shown in Table B.1. Figures 1 and 2 show moment 0 maps for Class 0 and Class I sources, respectively. Moment 0 maps of those sources where the value of the bolometric temperature (T_bol) is unknown, are presented in Fig. 3, given that some of them (SVS 13 A2 and RAC1999 VLA 20) show emission of CS and/or SO.

From moment 0 maps, the CS line usually shows centrally peaked and extended emission. In some cases, the CS extended emission is clearly tracing the outflow cavities, such as in Peremb 1, Per-emb 11A, Per-emb 12, Per-emb 13, Per-emb 26, Per-emb 27, and Per-emb 29, while in other cases it seems to be tracing emission perpendicular to the outflow direction, possibly related to a wider outflow cavity wall or an infalling envelope: good examples are Per-emb 8, Per-emb 18, Per-emb 20, Per-emb 25, Per-emb 36, L1448 NW, and L1448 IRS 3A. Only compact emission (comparable to the beam size) is seen around Per-emb 28 and SVS 13 A2, while the emission peaks offset from the continuum peak in the cases of Per-emb 2, Per-emb 12 A and B, Per-emb 13, IRAS 4B2, Per-emb 53, Per-emb 54, Per-emb 55, and L1455 IRS 2. The absence of line emission around the central protostar could be due to optically thick dust, and De Simone et al. (2020) showed that this is the case for Per-emb 12A (also known as IRAS 4A) where COMs are detected in centimeter wavelengths but not in the millimeter regime. Yang et al. (2021) showed that the other sources with mean continuum brightness temperature (T_b,cont) (a proxy of gas column density that correlated with the dust mass of the disk-like structure) higher than 10 K, Per-emb 5, Per-emb 12B, Per-emb 13, Per-emb 11A, Per-emb 33A, and B1-b S, feature detections of several molecular lines, suggesting that the continuum opacity is not significant. For a reference, the T_b,cont value of Per-emb 12A is 21.9 K (Yang et al. 2021), and this is the only source in this sample where the dust emission will be considered as optically thick.

Both SO transitions agree with each other and they show a combination of compact and extended emission among the sample. Outflow structures are clearly traced by SO toward Per-emb 1, Per-emb 11A, Per-emb 12, Per-emb 13, Per-emb 26, and Per-emb 27. We note that those sources that have SO tracing the outflow also have CS tracing the outflow – and they are all Class 0. Additionally, some arc-like structures are seen toward Per-emb 36, Per-emb 40, and Per-emb 44.

The ³⁴SO and SO₂ lines generally show compact emission around the protostar. Only a few exceptions could be seen: the emission extends along the outflow direction toward Per-emb 27, it presents a banana shape perpendicular to the outflow direction in Per-emb 36, it peaks toward the south in Per-emb 50, it peaks toward the east in Per-emb 53, and it is seen perpendicular to the outflow direction in L1448 NW. In addition to the compact emission around the protostar, the SO₂ transition also shows an east component perpendicular to the outflow direction toward Per-emb 18. It is interesting to note that, when ³⁴SO and SO₂ are detected, they trace the same structures. There are a few sources for which SO₂ is detected but ³⁴SO is not, such as Per-emb 1, Per-emb 20, Per-emb 22 A and B, Per-emb 25, and Per-emb 42.

A peculiar case is seen toward Per-emb 50, where SO, ³⁴SO, and SO₂ lines are detected but CS is not. One possible reason for the missing CS is a resolving out effect by the interferometer. For two sources, EDJ2009-237 and RAC1999 VLA20, it is unknown whether those are protostars or not because their evolutionary stages have not been determined. No sulfur-bearing species are detected toward EDJ2009-237, but CS emission is seen in RAC1999 VLA20. Tobin et al. (2018) reported that RAC1999 VLA20 might be an extragalactic source, however, given that CS emission is detected toward this sources, we cannot rule out the possibility that RAC1999 VLA20 is a protostar or that the CS emission is associated with nearby protostars (SVS13 sources). Since it is difficult to compare chemical characteristics with other protostars without knowing their physical parameters, we present continuum fluxes, moment 0 maps, and spectra of EDJ2009-237 and RAC1999 VLA20, but they will no longer be discussed in this paper.

In addition to EDJ2009-237, there are other four sources with no detections of sulfur-bearing species: Per-emb 41, Per-emb 60, EDJ2009-172, and EDJ2009-235, all of them Class I sources. Per-emb 41 lies in the border of the field of view that corresponds to B1-b N and B1-b S, and this could prevent weak lines from being detected. Larger scale (≥3″) CS emission is detected toward L1448 IRS 2E and B1-b N, shown in Fig. C.1, and this emission seems to be related to the outflow component. For B1-b S, only CS absorption is seen toward the source position (see Fig. 4) and given that no emission is seen in the surroundings of the protostar, the absorption feature is likely to be caused by foreground gas and we consider this as a non-detection.

The detection, or no detection, of CS, SO, ³⁴SO, and SO₂ is shown in Table B.4, together with the amount of complex organic molecules (COMs) detected by Yang et al. (2021).

The relationship between detectability of COMs and S-bearing species will be discussed in more detail in Sect. 4.3. In addition, large scale emission of CS is presented in the Appendix C for some sources with complex morphologies: Per-emb 12, Per-emb 13, IRAS 4B2, Per-emb 26, Per-emb 42, Per-emb 27, Per-emb 33, and L1448 IRS ЗА (Fig. C.2).

Fig. 2 Class I sources. Moment 0 maps of CS, SO, 34SO, and SO2. The contours start at 3σ, follow steps of 1σ, and the emission has been integrated over the velocity ranges listed in Table B.1. Maximum values (in units of Jy beam−1 km s−1) are indicated in the bottom-right corner of each panel. Red and blue arrows represent the outflow direction and the yellow star indicates the position of the continuum peak. Magenta stars represent binary components and the synthesized beam for each molecular transition is shown in the panels of the first row. Sources with no detections are not shown in this figure.

Fig. 3 Sources with unknown TЬо1. Moment 0 maps of CS, SO, 34SO, and SO2. The contours start at 3σ, follow steps of 1σ, and the emission has been integrated over the velocity ranges listed in Table B.1. Maximum values (in units of Jy beairr−1 km s−1) are indicated in the bottom-right corner of each panel. Red and blue arrows represent the outflow direction and the yellow star indicates the position of the continuum peak. The synthesized beam for each molecular transition is shown in the panels of the first row. Sources with no detections are not shown in this figure.

3.2 Spectra

A circular area that corresponds to the major axis of the synthesized beam (see Table B.1) around the continuum peak position was used to integrate the spectra of the CS, SO 6₆−5₅, and SO₂ lines (Fig. 4), and SO 6₆−5₅, SO 7₆−6₅, and ³⁴SO transitions (Fig. D.1). In most cases, all the sulfur-related species show emission within the same velocity range and the CS line presents absorption features in most of the cases, therefore, integrated CS fluxes calculated in the next section correspond to lower limits. The spectra of Per-emb 28 and Per-emb 41 look very noisy because they lie close to the edge of the field of view.

Five of the sources, Per-emb 12B, Per-emb 13, Per-emb 27, Per-emb 29, and Per-emb 44, show an offset emission (about +8 km s⁻¹ from v_sys) in the SO 6₆−5₅ and SO₂ spectra (see vertical arrows in Fig. 4). This offset emission is also seen in ³⁴SO but not in SO 7₆−6₅ (see Fig. D.1). As listed in Table 1, both SO transitions have similar A_ij and E_up values, therefore, a similar behavior is expected between the two transitions, unless SO 7₆−6₅ is much more optically thick than SO 6₆−5₅. If this offset component is related with redshifted SO 6₆−5₅ and SO₂ emission, it is intriguing to note that there is no blueshifted counterpart seen around −8 km s⁻¹. Another possibility could be that the offset emission is related with another molecular transition, most likely a COM, given that those five sources are the ones where most of the COMs have been detected in Yang et al. (2021; see Table B.4). Possible species for this offset emission are methyl cyanate (CH₃OCN at 244.2471 GHz), formic acid (c-HCOOH at 244.2479 GHz), acetaldehyde (CH₃CDO at 244.2466 GHz), vinyl cyanide (CH₂CHC¹⁵N at 258.2465), or ethylene glycol ((CH₂OH)₂ at 258.2506 GHz), which are the strongest lines found close to the offset frequencies (244.2477 and 258.2493 GHz) in SPLATALOGUE¹ and using the molecular data from the Cologne Database of Molecular Spectroscopy (CDMS; Müller et al. 2001, 2005; Endres et al. 2016). The origin of this redshifted offset emission is out of the scope of this paper and we focus on the central emission for the rest of this work.

An interesting spectral shape is seen for Per-emb 50, where the CS line is not detected and the SO spectrum does not show a Gaussian profile, but a broader and increasing spectral shape instead. Valdivia-Mena et al. (2022) proposed that a streamer is feeding the envelope-disk system of Per-emb 50 and SO and SO₂ could be enhanced in shocked regions. The existence of shocks in this source is also supported by Zhang et al. (2023), who studied the relationship between sulfur-bearing species and dust polarization toward some of the PEACHES sources. Per-emb 50 will be discussed in more detail in Sect. 4.3.

Fig. 4 Integrated spectra of CS (red), SO 66−55 (blue), and SO2 (green). The integration area for each source is shown in Table B.1. Black arrows in the panels that correspond to Per-emb 12 B, Per-emb 13, Per-emb 27, Per-emb 29, and Per-emb 44 indicate an offset emission, likely related with a transition from an unidentified COM.

Fig. 5 Detection rate of CS, SO, 34SO, and SO2 toward Class 0 and Class I sources. Detected refers to emission above 3σ that is observed in the vicinity of the source (not necessarily toward the continuum peak).

3.3 Detection rates

The whole sample consists of 36 Class 0 and 14 Class I sources, along with detection rates for CS, SO, ³⁴SO, and SO₂ shown in Fig. 5. A given transition is considered to be detected if its integrated emission is above a Зσ value (see Table B.2) or if it shows extended emission in moment 0 maps (see Figs. 1–3 and C.1). Both SO transitions are detected toward the same sources, therefore, we do not distinguish between them when discussing detection rates. As seen in Fig. 5, the CS line is almost ubiquitous toward Class 0 sources and is detected in most of the Class I systems. In this case, the detection rate decreases in around 30% from Class 0 to Class I sources. SO is detected in most of the Class 0 sources and around half of the Class I ones. The SO detection rate also shows a considerable decrease from 86% to 57% from Class 0 to Class I sources, similar to the case of CS. A different behavior is seen for ³⁴SO and SO₂, where detection rates are similar for Class 0 and I sources: around 35% for ³⁴SO and around 45% for SO₂.

Given that the envelopes of Class 0 sources are more massive than those from Class I systems, the decrease of the CS and SO detections rates from Class 0 to Class I sources seems to be related to a reduction in the envelope mass. On the other hand, the detection rates of ³⁴SO and SO₂ do not decrease from Class 0 to Class I sources, suggesting that these ³⁴SO and SO₂ lines are not sensitive to the envelope mass, but they seem to be sensitive to L_bol (see Sect. 4.1), which can be linked to the mass accretion rate (Ṁ_ace).

The sulfur bearing species, when detected, are tracing multiple components of the system as shown in Fig. 6. The CS transition typically traces extended emission, consistent with the outflow direction and, in a few cases, the emission is seen perpendicular to the outflow direction. In addition, five of the sources (11%) are associated with CS extended structures but no central emission is detected and only one source (B1-b S) shows pure absorption features. Almost half of the SO detections show compact emission (when compared to the beam size) and the rest present extended emission, where two of them (5%) reveal a clear arc-like structure (Per-emb 40 and Per-emb 44). Most of the extended emission detected in SO seems to arise from the outflow component. Compact emission is mainly seen when ³⁴SO and SO₂ are detected and, when the emission is extended, most of the sources present a structure perpendicular to the outflow direction. This would be discussed in more detail in Sect. 4.3.

Fig. 6 Different components where the CS, SO, 34SO, and SO2 transitions are detected. Extended emission is divided into three groups, those sources where the emission is seen perpendicular or parallel to the outflow direction, and those where the structure is unclear (labeled as Extended).

3.4 Fluxes

From moment 0 maps, we chose an integrated region that coincides with the beam size (Table B.1) and employed the casa task imf it to calculate the molecular integrated fluxes listed in Table B.2. Given that the SO 6₆−5₅ transition is slightly optically thinner than the SO 7₆−6₅ one (see Sect. 4.2.1), integrated fluxes for SO correspond to the 6₆−5₅ transition. An upper limit of Зσ is set for non detections and lower limits represent optically thick emission.

3.5 Column densities

Column densities were calculated for both SO transitions, as well as ³⁴SO and SO₂. Given that the CS spectra show absorption features in most of the cases (an indication of optically thick emission and/or self absorption due to infall motion), column densities for CS are lower limits in most of the cases.

By employing the non-local thermodynamic equilibrium (LTE) radiative transfer code RADEX (van der Так et al. 2007), we found that the SO 6₆−5₅ transition is thermalized for H₂ number densities (n_H) above 10⁷ cm⁻³. Assuming that the SO emission is coming from dense regions (n_н ≥ 10⁷ cm⁻³), a LTE analysis was employed to calculate the molecular column densities for SO and ³⁴SO. On the other hand, we found that the SO₂ 14_0,14−13_1,13 transition is thermalized only for very high densities (n_h ≥ 10⁹ cm⁻³); therefore, the molecular column densities for SO₂ were estimated from a non-LTE analysis using RADEX.

3.5.1 LTE analysis

Assuming optically thin conditions and that the gas is under LTE, for a given molecular transition, the column density of the upper level (N_u) is given by: $$N_{\text{u}}=\frac{8\pi k{v}^{2}W}{h{c}^3{A}_{\text{ul}}},$$(1)

where k is the Boltzmann constant, ν is the frequency of the molecular transition, W is the integrated line intensity, h is the Planck constant, с is the speed of light, and A_ul is the Einstein coefficient for spontaneous emission. Equation (1) can be rewritten as: $$N_{\text{u}}=2375\times {10}^6\left(\frac{W}{1\text{\hspace{0.17em}\hspace{0.17em}Jy\hspace{0.17em}km\hspace{0.17em}}{\text{s}}^{-1}}\right)\left(\frac{1{\text{s}}^{-1}}{A_{\text{ul}}}\right)\left(\frac{{\text{arcsec}}^2}{{\text{Θ}}_{\text{area}}}\right),$$(2)

where Θ_area is the area of integration (in arcsec; see Table B.1) and N_u will be obtained in units of cm⁻². Employing a range of excitation temperatures (T_ех from 50 to 150 K), we calculated the total column density of the molecule (N) following $$N=\frac{N_{\text{u}}Q\left(T_{\text{ex}}\right)\text{\hspace{0.17em}exp}\left(E_{\text{up}}/T_{\text{ex}}\right)}{g_{\text{u}}},$$(3)

where Q(T_ex) is the partition function that depends on the excitation temperature, E_up is the upper level energy in K, and ɡ_u is the level degeneracy. Total column densities for SO and ³⁴SO are listed in Table B.3.

3.5.2 Non-LTE analysis

Given that the SO₂ 14_0,14−13_1,13 transition is not thermalized for typical inner envelope and disk densities (10⁷−10⁹ cm⁻³), RADEX was employed to estimate SO₂ molecular column densities (see the example of Per-emb 29 in Fig. E.1). Table B.3 lists a range of values for $N_{{\text{SO}}_2}$ for T_ex between 50 and 150 K and n_H between 10⁷ and 10⁹ cm⁻³. Collisional rates for SO₂ were taken from the Leiden atomic and molecular database (LAMDA; Balança et al. 2016) and a broadening parameter (b) of 2 km s⁻¹ was employed for the calculations. Four sources, Per-emb 12B, Per-emb 17, Per-emb 27, and Per-emb 44, are associated with lower limits of $N_{{\text{SO}}_2}$, given that their molecular column densities can be much higher than 10¹⁶ cm⁻² for n_H = 10⁷ cm⁻³ and the SO₂ emission toward these four sources was found to be optically thick. A more detailed discussion of the optical thickness of the SO₂ transition is presented in Sect. 4.2.2.

4 Analysis and discussion

This section is dedicated to the analysis and discussion of integrated fluxes as a function of T_bol and L_bol, optically thickness of the lines, molecular column densities, and detectability of COMs and sulfur-bearing molecules.

4.1 Fluxes as a function of T_bol and L_bol

Continuum and molecular line fluxes are compared with the physical properties of the sources, such as T_bol and L_bol, and are shown in Fig. 7. The continuum flux decreases as T_bol increases, which is expected, given that the envelope material deceases as the source evolves. A less clear positive trend is seen for the continuum flux as a function of L_bol. Multiple sources and unresolved binaries are highlighted in red and blue colors, given that the flux of the components cannot be totally separated in most of the cases.

Molecular fluxes do not show any trend as a function of T_bol (left column of Fig. 7); however, a special case must be taken when analyzing CS and SO given that these molecular transitions could be optically thick: CS absorption is seen in the spectra of almost all the sources and (as we explain in Sect. 4.2) when ³⁴SO is also detected, both SO transitions are optically thick. A slightly positive trend is seen for the SO₂ flux as a function of Lbol, mainly given by the two most luminous sources, Per-emb 27 and Per-emb 44. Without them, the positive trend disappears. Despite the absence of a strong trend, it is clear that ³⁴SO and SO₂ are only detected in those sources with L_bol > 1 L_⊙. This result could be a selection bias, given that 70% of the sources in this sample are associated with L_bol > 1 L_⊙, or ³⁴SO and SO₂ are hard to detect in sources with low L_bol. This is consistent with the results from Artur de la Villarmois et al. (2019), where they studied a sample of 10 Class 0/I sources and detected SO₂ emission only in those sources with L_bol ≥ 1.4 L_⊙. More sensitive observations of sources with low L_bol (< 1 L_⊙) are needed to test these detection trends.

4.2 Column densities

4.2.1 Optically thick lines

The CS line shows absorption features in most of the spectra, therefore, the CS integrated fluxes and column densities represent lower limits. For both SO transitions, it is not obvious if they are optically thin or optically thick. Therefore, when the SO and ³⁴SO lines are detected, the molecular column density of both SO lines are compared to the ³⁴SO column density (for T_ex between 50 and 150 K). The SO/³⁴SO ratio should approach a value of 22 (³²S/³⁴S abundance ratio; Wilson 1999) if the SO line is optically thin. Figure 8 shows the molecular column density ratios between the two SO transitions and ³⁴SO, revealing that both SO transitions are considerably optically thick in those sources where ³⁴SO is also detected. For this reason, future calculations of the SO column density will be performed from the ³⁴SO integrated flux and later on multiplied by 22. It is important to notice that the ³⁴SO emission could also be optically thick and, in this case, SO column densities should be considered as lower limits and the SO/³⁴SO ratios as upper limits. Unfortunately, we cannot assess the optical thickness of the ³⁴SO 5₆−4₅ transition with the current data set.

Per-emb 13 is the source with the highest detection of COMs (Yang et al. 2021), followed by Per-emb 12 B, Per-emb 27, Per-emb 29, and Per-emb 44. These five sources have at least 14 COMs detected (see Table B.4) and they are associated with the lowest SO/³⁴SO ratios (see Fig. 8), consistent with the most optically thick SO emission. In this sense, the SO/³⁴SO ratio might be a good tracer of the inner high-density envelope, and a low value (≤1.5) seems to be related with the detection of multiple COMs.

4.2.2 Optical thickness of the SO₂ transition

Given that the PEACHES survey did not target any ³⁴SO₂ transition, the same analysis as the one carried out for SO cannot be replicated for SO₂. The SO₂ emission could, therefore, be optically thick. To assess the optical depth of the SO₂ 14_0,14−13_1,13 transition, we employed RADEX (van der Tak et al. 2007) and generated models for H₂ number densities (n_H) between 10³ and 10⁹ cm⁻³ and three different excitation temperatures: 50, 100, and 150 K. A broadening parameter (b) of 2 km s⁻¹ was employed for the calculations. RADEX results are presented in Fig. E.2 and compared with calculated SO₂ column densities (see Sect. 3.5.2 and Table B.2). For a value of $N_{{\text{SO}}_2}$ below 10¹⁵ cm⁻², the SO₂ emission could be interpreted as optically thin emission (assuming that the emission fills the beam uniformly) and most of the sources lie within this range. It should be noticed that the SO₂ emission might originate in a smaller region compared to the beam size and, in this case, the emission could be optically thick. Higher angular resolution observations would be needed to assess the optical depth of the SO₂ emission at smaller scales. On the other hand, there are four sources, Per-emb 12A, Per-emb 17, Per-emb 27, and Per-emb 44, associated with lower limits of $N_{{\text{SO}}_2}$ (≥2 × 10¹⁵ cm⁻²); therefore, they have high probabilities for the SO₂ emission to be optically thick (see Fig. E.2).

4.2.3 $N_{{\text{SO}}_2}/N_{\text{SO}}$

Assuming that the SO₂ transition is optically thin in most of the sources, the column density ratio between SO₂ and SO is shown in Fig. 9 as a function of T_bol and L_bol. No clear trend is seen for $N_{{\text{SO}}_2}/N_{\text{SO}}$ as a function of T_bol. On the other hand, $N_{{\text{SO}}_2}/N_{\text{SO}}$ shows a tentative increasing trend for L_bol between 1 and 10 L_⊙ (red dots). Those sources where SO emission is detected but ³⁴SO emission is not, are shown with blue dots and represent upper limits, given that the SO emission could be optically thick. For L_bol higher than 10 L_⊙, a tentative decreasing trend is seen for the $N_{{\text{SO}}_2}/N_{\text{SO}}$ ratio, however, caution must be taken given that the SO₂ emission is optically thick toward the most luminous sources. As seen in Fig. 7, both SO and ³⁴SO, as well as SO₂ fluxes increase as a function of L_bol; however, the UV and X-ray radiation field of the protostar could efficiently photo-dissociate SO₂ to form SO and decreases the SO₂/SO ratio for higher luminosities (Booth et al. 2021; van Gelder et al. 2021). Future observations of SO, ³⁴SO, SO₂, and ³⁴SO₂ toward other sources, within a broad range of bolometric luminosities, are needed to assess the $N_{{\text{SO}}_2}/N_{\text{SO}}$ ratio for L_bol higher than 10 L_⊙ and confirm (or otherwise) whether SO₂ is being efficiently destroyed in those luminous sources.

Fig. 7 Continuum and molecular fluxes (above 3σ) as a function of Tbol (left) and Lbol (right). The gray dashed line represents Tbol = 70, the separation between Class 0 and Class I sources, while the black solid line in the last two panels of the right column indicates Lbol = 1 L⊙. Red circles represent multiple systems (binaries or triples) where the continuum ALMA emission is resolved but they do not have individual Tbol and Lbol measurements from far IR. Blue circles indicate unresolved binaries in continuum ALMA emission.

Fig. 8 Molecular column density ratio of SO over 34SO, showing that both SO transitions are optically thick. A range of Tex, between 50 and 150 K has been employed for the calculations. Dark blue and light blue bars indicate ratios calculated from SO 76−65 and SO 66−55, respectively. The red vertical line separates Class 0 (left) from Class I sources (right).

4.2.4 $N_{{\text{H}}_2}$

To calculate molecular abundances, the column density of molecular hydrogen ($N_{{\text{H}}_2}$) is estimated from the continuum peak flux and following Eq. (A.27) from Kauffmann et al. (2008), as follows: $$\begin{array}{l}{N}_{{\text{H}}_2}=2.02\times {10}^{20}{\text{cm}}^{-2}\left({\text{e}}^{1.439{\left(\lambda /\text{mm}\right)}^{-1}{\left(T_{\text{dust}}/10\text{\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}K}\right)}^{-1}}-1\right){\left(\frac{\lambda }{\text{mm}}\right)}^3\hfill \\ \cdot {\left(\frac{{\kappa }_v}{0.01{\text{cm}}^2{\text{g}}^{-1}}\right)}^{-1}\left(\frac{F_v^{\text{beam}}}{\text{mJy\hspace{0.17em}}{\text{beam}}^{-1}}\right){\left(\frac{{\theta }_{\text{HPBW}}}{10\text{\hspace{0.17em}\hspace{0.17em}arcsec}}\right)}^{-2},\hfill \end{array}$$(4)

where λ is the wavelength, T_dust is the dust temperature, κ_ν is the dust opacity, $F_v^{\text{beam}}$ is the observed flux per beam, and θ_HPBW corresponds to the beam size. For λ = 1.2 mm, we employed κ_ν = 0.01056 cm² g⁻¹ (for a gas-to-dust radio of 100; Ossenkopf & Henning 1994), and a dust temperature of 30 K (found to be appropriate for Class 0 sources; Dunham et al. 2014). Values for $F_v^{\text{beam}}$ and θ_HPBW are listed in Tables A.1 and B.1, respectively. Column densities of H₂ are listed in Table B.3 and they will be used in the following sections to infer the abundances of SO₂ and CH₃OH.

4.3 S-bearing species versus COMs detection

In the study by Yang et al. (2021), the authors present the detection (or non-detection) of COMs toward the PEACHES sample. From a total of 16 COMs, the most frequently detected ones were CH₃OH and CH₃CN. The rest are NH₂CHO, C₂H₅CN, CH₂DCN, HCOCH₂OH, CH₃CHO, C₂H₅OH, CH₃COCH₃, CH₃OCH₃, CH₃OCHO v =1, CH₃OCHO, t-HCOOH, ¹³CH₃OH, CH₃¹⁸OH, and CH₂DOH. We note that none of these COMs contain sulfur.

The molecule SO₂ has been associated with warm chemistry, given that its desorption temperature is above 60 K (van Gelder et al. 2021), and it is efficiently formed in the gas phase for temperatures between 100 and 200 K (Charnley 1997). Therefore, some correlation between the detection of COMs and SO₂ would be expected. To test this correlation, we separated those sources for which SO₂ emission is detected into three groups, following the detection rate of COMs from Yang et al. (2021): (i) sources very rich in COMs (at least five different COMs); (ii) sources with only a few COMs detected, between one and three; and (iii) sources with no COMs detected. These three groups are presented in Fig. 10, where the origin of the emission from the sulfur-bearing species is also indicated.

Most of the sources that show CS and SO emission that is consistent with the outflow component are rich in COMs, with the exception of Per-emb 44, which shows an arc-like structure, Per-emb 35 A, where the emission seems to connect the binary system, and Per-emb 17, whose CS emission is seen in both directions: parallel and perpendicular to the outflow direction. The only source rich in COMs, where SO₂ (in addition to SO and ³⁴SO) has not been detected, is B1-b S. Nine different COMs were detected toward this source (Yang et al. 2021) and only CS absorption is seen in this dataset. Nevertheless, the non detection of ³⁴SO and SO₂ is consistent with the low L_bol value of B1-b S (0.32 L_⊙). The detection of COMs in the gas-phase seems to be related with the presence of collimated outflows, where a high column density of warm gas is present. This is consistent with the results from Yang et al. (2021) and van Gelder et al. (2020), where they found that the detectability of COMs depends neither on the continuum brightness temperature, nor the bolometric luminosity and the bolometric temperature. The detection of SO₂ around the source position, together with CS and SO emission tracing high collimated outflows, could be an important indicator of the presence of multiple COMs in the gas-phase.

When only a few COMs have been detected (mainly CH₃OH and CH₃CN), the emission of S-bearing species is seen toward multiple components. The emission could be perpendicular to the outflow direction, tracing the material between resolved binaries, consistent with the outflow component and, for two sources L1448 IRS 3A and Per-emb 42, it is not clear if the emission is related with those sources or it is being affected by nearby outflows (see Fig. C.2).

For those sources with SO₂ detection and CS and SO emission perpendicular to the outflow direction (L1448 NW, Per-emb 25, Per-emb 36, and Per-emb 50), no COMs have been detected (see Fig. 10). These sources might be associated with wider outflow cavity walls and, consequently, with a lower column density of warm gas, or with infalling material from the envelope to the disk-protostar system. It is important to notice that, at disk scales, the outflow cavity walls could show a flattened structure in the disk direction, given that these regions very close to the protostar belong to the base of the outflow structure. The two scenarios, infalling material or wider cavity walls, require a more detailed dynamic study and it is beyond the scope of this paper. The SO₂ emission toward these sources with no detection of COMS could be originating in regions that are exposed to UV radiation from the protostar, in disk winds (Tabone et al. 2017), shocked regions in the inner regions of the outflow cavity walls (Podio et al. 2015), or accretion shocks at the disk-envelope interface (Artur de la Villarmois et al. 2022; Zhang et al. 2023). Observations with higher angular resolution and a more detailed study of these four scenarios are needed to assess the origin of the SO₂ emission.

Per-emb 50, one of the sources where COMs have not been detected, has been associated with an infalling streamer by Valdivia-Mena et al. (2022). In addition, Zhang et al. (2023) proposed that the ³⁴SO and SO₂ emission toward Per-emb 50 is related with shocks at the centrifugal barrier (~0″.4), which is located beyond the radius of the disk (~0″.11) derived from dust continuum maps (Agurto-Gangas et al. 2019). Apart from Per-emb 50, there are another three sources (Per-emb 2, Per-emb 27, and Per-emb 44) that have been associated with infalling streamers. Pineda et al. (2020) proposed that large-scale accretion flows funnel material down to disk scales in Per-emb 2 by observing molecular transitions from HC3N, CCS, and ¹³CS. The binary system Per-emb 27 has been associated with accretion streamers caused by the interaction of the binaries by Jørgensen et al. (2022), and Hsieh et al. (2023) proposed that the system Per-emb 44 might be fed by an infalling streamer from envelope scales (~700 au). Toward Per-emb 44, the arc-like structure that we observe in SO is consistent with the streamer suggested by Hsieh et al. (2023).

From this data set, the sources that show CS and SO emission perpendicular to the outflow direction (Per-emb 18, Per-emb 20, Per-emb 25, Per-emb 36, Per-emb 53, and L1448 NW) seem to be the perfect candidates for studying their large-scale structure in more detail and to search for shocks or disk winds at disk scales.

Figure 11 shows the CH₃OH abundance as a function of the SO₂ abundance, for those sources where SO₂ has been detected. The CH₃OH column densities were taken from Yang et al. (2021), and blue triangles show the upper limits, given that CH₃OH is not detected toward these sources. As seen in Fig. 11, the abundance of CH₃OH is not correlated with the detectability of multiple COMs: both groups (detected COMs between one and three and more than five COMs detected) present CH₃OH abundances between 10⁻⁹ and 10⁻⁷. Nevertheless, a positive trend is seen between the abundances of SO₂ and CH₃OH, and CH₃OH is between 10 and 100 times more abundant than SO₂. This positive trend may suggest that both species, SO₂ and CH₃OH, are coming from warm gas and, for those sources where COMs are not detected, the SO₂ emission might originate from a different physical process, such as shocked regions or disk winds. Values for SO₂ and CH₃OH abundances are listed in Table B.3.

Fig. 9 Column density ratio between SO2 and SO as a function of Tbol (left) and Lbol (right). Red dots are values calculated from 34SO fluxes and later on multiplied by 22 (Wilson 1999), while blue dots indicate those sources where34 SO is not detected and the ratios are calculated from the SO 66−55 transition. Black triangles represent lower limits (optically thick SO2 emission) and green triangles indicate upper limits, given by those sources where SO emission is detected but SO2 is not. The gray dashed line in the left panel represents Tbol = 70, the separation between Class 0 and Class I sources. Error bars in the y-axis represent ranges for Tех between 50 and 150 K.

Fig. 10 S-bearing species vs. COMs detection. Sources with detection of SO2 are clustered by the amount of detected COMs in Yang et al. (2021). The numbers on the left of the source names indicate the amount of different COMs detected. The different colors indicate where the emission of CS or SO is seen. Those sources where the CS or SO emission is parallel to the outflow direction are usually rich in COMs, while sources where the CS or SO emission is perpendicular to the outflow direction do not show detection of COMs.

Fig. 11 CH3OH abundance as a function of the SO2 abundance for the sources in Fig. 10. The different colors represent those sources where more than five different COMs were detected (green), between one and three different COMs were detected (red), and no COMs (blue) were detected in Yang et al. (2021). Gray dashed lines show a 10:1 and 100:1 linear correlation.

4.4 SO₂ detection in other sources

SO₂ emission has been detected in a handful of other sources, including the Class 0 sources IRAS 16293-2422 (Jørgensen et al. 2016; Drozdovskaya et al. 2018), a sample of five Class I sources in the Ophiuchus star forming region (Artur de la Villarmois et al. 2019), in different components of four Class I/II sources in Taurus (Garufi et al. 2022), and toward the Class II disk IRS 48 (Booth et al. 2021). The ratio between $N_{{\text{SO}}_2}$ and $N_{{\text{H}}_2}$ for these sources and for the PEACHES sample is shown in Fig. 12, highlighting the different star-forming regions and evolutionary classes. For the PEACHES sample, an excitation temperature (T_ex) between 50 and 150 K has been employed for the $N_{{\text{SO}}_2}$ calculations and the region of integration corresponds to the beam size (~0″.5 or 150 au). Assumptions and calculations made with values from the literature are presented in Appendix F and Table F.1.

The SO₂ abundance of the PEACHES sample seems to increase from Class 0 to Class I sources, probably due to a decrease in the continuum flux (and therefore $N_{{\text{H}}_2}$) as the system evolves. Most of the Ophiuchus Class I sources show that the SO₂ abundance is between one and two orders of magnitude higher than the PEACHES Class 0/I sample. One of the Ophiuchus Class I sources, IRS 44, has been associated with the presence of accretion shocks, from the analysis of multiple SO₂ and ³⁴SO₂ transitions (Artur de la Villarmois et al. 2022). IRS 48, also located in Ophiuchus, shows higher values than less evolved Class I/II sources in Taurus. IRS 48 is the only Class II disk where prominent SO₂ has been detected so far and Booth et al. (2021) proposed that the SO₂ emission is originating from the sublimation of ices at the dust cavity wall, observed as an asymmetric dust trap (van der Marel et al. 2013). Therefore, the high SO₂ column densities toward the Ophiuchus Class I sources and toward IRS 48 seem to be related to the sublimation of SO₂ from the dust grains, either through accretion shocks or the exposure to the UV radiation from the protostar. Intermediate Class I/II sources in Taurus show a large spread in SO₂ abundances, mainly given by two lower limits and T Tau S, a source with bright SO₂ emission (~60 mJy beam⁻¹ km s⁻¹; Garufi et al. 2022) and weak continuum emission (~10 mJy beam⁻¹; Garufi et al. 2021). Finally, the hot-corino source IRAS 16293 presents an upper limit in the SO₂ abundance, given by the optically thick continuum, and the fluxes were taken in an offset position (Jørgensen et al. 2016; Drozdovskaya et al. 2018).

A possibility for the observed SO₂ abundances in Fig. 12 is that the Ophiuchus star-forming region could be enriched in sulfur-bearing species, in comparison with Perseus (the PEACHES sample) and Taurus (Class I/II sources). In this scenario, the sulfur enrichment will be present in the solid phase and, when some physical mechanism heats the dust grains, S-bearing species get desorbed from the dust to the gas-phase. The heating mechanism (outflows, accretion shocks, or exposure to the protostellar UV radiation) is essential, given that there are other Class I sources in Ophiuchus (associated with low L_bol) where SO₂ is not detected (Artur de la Villarmois et al. 2019). This interpretation is consistent with resent results from Fuente et al. (2023), where they studied the sulfur elemental abundance in Taurus, Perseus, and Orion A, finding that Taurus and Perseus are more depleted in sulfur than Orion A. They suggest that dust grains are negatively charged in regions where S⁺ is abundant toward Taurus and Perseus, promoting the adsorption of S⁺ to the grain surface. With this adsorption of ions, Laas & Caselli (2019) and Shingledecker et al. (2020) predicted that organosulfur compounds and S₈ could be an important sulfur reservoir in the solid phase.

The main difference between Taurus/Perseus and Orion/Ophiuchus is the external UV intensity, G₀. This value is ~30 toward Perseus and ≤ 100 in Taurus, while it could reach values from 100 to 1000 in Orion A and Ophiuchus (Xia et al. 2022). Given that Ophiuchus is more similar to Orion A in terms of UV intensity, it could be possible that Ophiuchus is also undepleted in sulfur, supporting its high SO₂ column density observed in Fig. 12. Future observations of SO₂ transitions in a representative number of sources located in different star-forming regions would clarify if the Ophiuchus molecular cloud is, indeed, particularly enriched in sulfur-bearing species. In addition, solid-phase observations could demonstrate if the sulfur reservoir varies between star-forming regions that are exposed to different UV radiation fields.

The linear trend between the SO₂ and CH₃OH abundances shown in Fig. 11 suggests that the SO₂ emission toward most of the Class 0 sources in Fig. 12 arises from the high column density of warm material. In the Ophiuchus Class I sources, CH₃OH emission has not been detected (Artur de la Villarmois et al. 2019), and SO₂ emission toward these sources and more evolved Class I/II might originate in accretion shocks. Finally, the more evolved Class II disk IRS 48 presents SO₂ emission that originates from a dust trap exposed to the UV radiation of the protostar. Therefore, the SO₂ emission detected in different evolutionary stages seems to arise from different physical mechanisms.

Fig. 12 Column density ratio between SO2 and H2 for sources in different evolutionary stages. Red, green, and blue symbols represent sources from Perseus, Ophiuchus, and Taurus, respectively. Dots indicate Class 0 sources, squares represent Class I sources, and diamonds show Class I/II and Class II sources. Empty triangles indicate upper or lower limits. Values for IRAS 16295 were taken from Drozdovskaya et al. (2018) and Jørgensen et al. (2016). Perseus sources are from this work, Ophiuchus Class I sources correspond to IRS 67, GSS30-IRS1, IRS 44, IRS 43, and Elias 29, from left to right, studied by Artur de la Villarmois et al. (2019), and Taurus sources are taken from Garufi et al. (2022) and represent DG Tau (outer disk), HL Tau (inner disk), HL Tau (outer disk), IRAS 04302, and T Tau S, from left to right. The upper limit for IRAS 16293 is given by the optically thick continuum emission and lower limits for Perseus and Taurus sources are given by optically thick SO2 emission. The PEACHES sample and Class I Ophiuchus sources are ordered by Tbol.

5 Summary

This work presents high-resolution (~0″.6, 180 au) ALMA observations of a sample of 50 Class 0/I sources located in the Perseus star-forming region. A detailed study of the sulfur-bearing species (CS, SO, ³⁴SO, and SO₂) has been carried out and the main results are summarized below:

• Moment 0 maps show that the CS emission traces extended structures and follows the outflow direction in most of the cases. The SO emission is also seen toward extended structures, but an important percentage of sources (~40%) show compact emission, while the majority of the ³⁴SO and SO₂ detections are compact;

• A comparison of the detection rates between Class 0 and Class I sources indicates that the CS and SO lines are sensitive to the evolutionary state of the source and, therefore, to the amount of mass present in the envelope. On the other hand, the ³⁴SO and SO₂ transitions seem to be more sensitive to the bolometric luminosity, which can be linked to the mass accretion rate;

• The ³⁴SO and SO₂ transitions are only detected in those sources with bolometric luminosities above 1 L_⊙, which could be a real finding or a selection bias, given that 70% of the sources have L_bol ≥ 1 L_⊙;

• When ³⁴SO is detected, the ratio SO/³⁴SO is lower than the canonical value of 22 in all the cases, suggesting that the SO emission is optically thick toward these sources. The lowest values (≤1.5) are associated with those sources rich in COMs, where at least 14 different COMs were detected. Therefore, the SO/³⁴SO ratio seems to be a good tracer of the inner high-density envelope and it could be used in the future to infer the presence of multiple COMs;

• The column density ratio between SO₂ and SO shows a positive trend for sources with L_bol ≤ 10 L_⊙ and a tentative decreasing trend for sources with L_bol ≥ 10 L_⊙, possibly related to the photo-dissociation of SO₂ and enhancement of SO in those very luminous sources;

• The comparison between the morphology of S-bearing species and the detectability of COMs suggests that the detection of COMs is linked to the presence of high col-limated outflows. The detection of SO₂ around the source position, together with CS and SO emission tracing high collimated outflows, could be an important indicator of the presence of multiple COMs in the gas-phase;

• When the emission of CS and SO is seen perpendicular to the outflow direction, no COMs (or only a few of them) are detected. In these cases, the CS and SO emission could be tracing wider outflow cavity walls or infalling material, and higher angular resolution observations are needed to assess the origin of the SO₂ emission: UV irradiated regions, disk winds, accretion shocks at the disk-envelope interface, or shocked regions in the outflow cavity walls;

• Per-emb 18, Per-emb 20, Per-emb 25, Per-emb 36, Per-emb 53, and L1448 NW seem to be ideal candidates to investigate in more detail their large scale structures and to search for the origin of the SO₂ emission close to the protostar;

• In those sources where SO₂ and CH₃OH are detected, a positive trend is seen between the CH₃OH and the SO₂ abundance, suggesting that both species are tracing warm material in those sources. The CH₃OH abundance is between 10 and 100 times higher than the SO₂ abundance and sources with multiple COMs detection do not show any trend with the CH₃OH abundance;

• When SO₂ is detected, the SO₂ abundances toward the sources in Perseus seem to be comparable to that in Taurus sources and lower than Ophiuchus sources. The depletion of sulfur seems to depend on the external UV radiation field of the molecular cloud and Perseus would be more depleted in sulfur than Ophiuchus;

• The SO₂ emission toward young Class 0 sources seem to be associated with a high column density of warm material, while accretion shocks might desorb SO₂ molecules to the gas-phase in Class I and Class I/II. Finally, the exposure to the protostellar UV radiation is the preferred mechanism to explain the SO₂ in more evolved Class II disks.

This paper shows the power of a large survey with multiple targeted molecular transition, making it possible to compare the emission of sulfur-bearing species between them and with the detection of COMs. Future observations of sulfur-bearing species toward Class 0/I sources in other star-forming regions are essential to confirm the observed values and trends in detection rates, SO₂ column densities and abundances, ratios between SO₂ and SO, and the possible connection with the detection of COMs.

Acknowledgements

We thank the anonymous referee for a number of good suggestions that helped us to improve this work. This paper makes use of the following ALMA data: ADS/JAO.ALMA#2016.0.00391.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. E.A.dlV. acknowledges financial support provided by FONDECYT grant 3200797. V. V.G. gratefully acknowledges support from FONDECYT Regular 1221352, ANID BASAL projects ACE210002 and FB210003, and ANID, – Millennium Science Initiative Program – NCN19_171. Y.-L.Y. acknowledges support from a Grant-in-Aid from the Ministry of Education, Culture, Sports, Science, and Technology of Japan (22K20389). Y.-L.Y. and N.S. acknowledge support from a Grant-in-Aid from the Ministry of Education, Culture, Sports, Science, and Technology of Japan (20H05845, 20H05844), and a pioneering project in RIKEN (Evolution of Matter in the Universe).

Appendix A Coordinates and continuum fluxes

Table A.1 lists the PEACHES sample, divided into Class 0 sources, Class I sources, and those sources where T_bol is unknown. A 2D Gaussian fit was employed to the continuum emission to estimate the coordinates and continuum fluxes.

Appendix B Integrated values, detection of S-bearing species, fluxes, and column densities

Table B.1 lists the systemic velocity (v_sys), estimated from the CS spectra (see Fig. 4), the outflow position angle (PA; measured from north to east), and the integration areas and velocity ranges used to create moment 0 maps and to extract molecular fluxes, which are listed in Table B.2. Later on, molecular fluxes were employed to calculate molecular column densities using T_ex between 50 and 150 K, and are presented in Table B.3. We note that the SO column density for those sources that have also detection of ³⁴SO need to be considered as lower limits, given that they are optically thick (see Sect. 4.2.1).

Table B.4 lists the number of COMs detected from Yang et al. (2021) and the detection (or not) of sulfur-bearing species. Detection could refer to emission toward the source position (above 3σ) or to extended emission that seems to be related with the source.

Appendix C Large-scale CS emission

Two Class 0 sources, L1448 IRS 2E and B1-b N, do not show molecular emission around the continuum peak position. However, extended CS emission is detected at large scales (Fig. C.1), possibly related with the outflow component.

Per-emb 12, Per-emb 13, and Per-emb 27 show strong and extended outflow structures that are shown in Fig. C.2. Per-emb 42 and L1448 IRAS3A lie very close to Per-emb 26 and Per-emb 33, respectively, and their close environments seem to be affected by the presence of these nearby sources (see Fig. C.2).

Fig. C.1 Same as Fig. 1 but for large-scale emission of CS toward L1448 IRS 2E and B1-b N.

Fig. C.2 Same as Fig. C.1 for Per-emb 12, Per-emb 13 and IRAS 4B2, Per-emb 26 and Per-emb 42, Per-emb 27, and Per-emb 33 and L1448 IRS3A.

Appendix D Spectra

Spectra of SO 7₆−6₅, SO 6₆−5₅ and ³⁴SO are presented in Fig. D.1, showing that both SO transitions exhibit a similar behavior. In addition, the bump in SO 6₆−5₅ (seen in a few sources) could be related to a different molecular transition (maybe to a COM) and not to the blueshifted or redshifted emission of SO.

Fig. D.1 Integrated spectra of SO 76−65 (red), SO 66−55 (blue), and 34SO 56−45 (green). The integration area for each source is listed in Table B.1 and the 34SO 56−45 spectra have been multiplied by three, for a better comparison. Black arrows in the panels that correspond to Per-emb 12 B, Per-emb 13, Per-emb 27, Per-emb 29, and Per-emb 44 indicate an offset emission, likely related with a transition from an unidentified COM.

Appendix E Radiative transfer

Given that the SO₂ 14_0,14−13_1,13 transition is not thermalized below a H₂ number density of 10⁹ cm⁻³, RADEX was employed to estimate a range of molecular column densities. Figure E.1 shows the example of Per-emb 29 for an excitation temperature of 50 and 150 K. A range of $N_{{\text{SO}}_2}$ is estimated for n_H between 10⁷ and 10⁹ cm⁻³, given the integrated flux (white contour), listed in Table B.3.

RADEX models were also employed to assess to optical depth of the SO₂ 14_0,14−13_1,13 transition. Figure E.2 shows the SO₂ molecular column density as a function of the H₂ number density, for three different values of T_ex: 50, 100, and 150 K. The range of calculated Ns¾ (see Table B.3) lies between 2 × 10¹⁴ and 5 × 10¹⁵ cm⁻²; thus, the SO₂ emission could be taken as optically thin emission in most of the cases. Four sources, Per-emb 12, Per-emb 17, Per-emb 27, and Per-emb 44, are associated with very bright SO₂ emission and their estimated $N_{{\text{SO}}_2}$ could reach values higher than 10¹⁶ cm⁻². Therefore, the SO₂ emission toward these four sources is considered to be optically thick and only lower limits are presented for $N_{{\text{SO}}_2}$ in Table B.3.

Fig. E.1 RADEX models for Per-emb 29. Fluxes of the SO2 140,14−131,13 transition for two different kinetic temperatures: 50 K (left) and 150 K (right). The white contour represents the flux observed for this source and the dashed cyan contour shows an H2 number density of 107 cm−3 (lower limit for nH used in the calculations). In this example, the $N_{{\text{SO}}_2}$ present a minimum value of 3 × 1014 cm−2 (for Tex = 150) and a maximum value of 1 × 1015 cm−2 (for Tex = 50), in a range between 107 and 109 cm−3.

Fig. E.2 RADEX models. Optical depth of the SO2 140,14−131,13 transition for three different excitation temperatures: 50, 100, and 150 K. The cyan contour represents τ = 1, the two dashed white vertical lines indicate the range of values obtained for $N_{{\text{SO}}_2}$ for most of the sources, and the dashed white horizontal line shows an H2 number density of 107 cm−3 (lower limit for nH used in the calculations).

Appendix F [SO₂/H₂] from the literature

There are only a handful of protostar where SO₂ emission has been detected at high angular resolution (≤ 1″.0) in nearby star-forming regions (≤ 300 pc): the Class 0 protostar IRAS 16293-2429, a sample of five Class I sources in Ophiuchus, different components of four Class I/II disks in Taurus, and the Class II disk IRS 48 (Drozdovskaya et al. 2018; Artur de la Vil-larmois et al. 2019; Garufi et al. 2022; Booth et al. 2021). All these sources are located at distances between 120 and 160 pc; therefore, special care has to be taken when comparing literature data with the PEACHES sample (Fig. 12), located at ~300 pc.

The region of integration for the PEACHES sample is equivalent to the size of the beam (~0″.5), proving regions of ~150 au in diameter around the central protostar. To trace similar regions in Taurus and Ophiuchus, a beam size (or integration area) of ~1″.0 would be needed. Table F.1 lists values taken from the literature to calculate the abundance of SO₂ (for an integration area of 1″.0 × 1″.0). Given that in most of the cases only the SO₂ column density is available, N(H₂) values were calculated using Eq. 4 and the conversion: $$\frac{F_v^{\text{beam}}}{{\theta }_{\text{HPBW}}^2}=\frac{\pi }{4\text{\hspace{0.17em}\hspace{0.17em}ln}\left(2\right)}\frac{F_v}{\text{Ω}},$$(F.1)

where F_v in the integrated flux and Ω the source area. Dust opacities were taken from (Ossenkopf & Henning 1994), where κ_ν = 0.0189 cm² g⁻¹ for λ = 0.87 mm (Ophiuchus Class I sources), κ_ν = 0.00899 cm² g⁻¹ for λ = 1.3 mm (Taurus Class I/II sources), and κ_ν = 0.0181 cm² g⁻¹ for λ = 0.89 mm (IRS 48). Jørgensen et al. (2016) employed κ_ν = 0.0182 cm² g⁻¹ for λ = 0.85 mm in IRAS 16293.

The [SO₂/H₂] value for IRAS 16293 results in an upper limit, given that the continuum emission is optically thick (Jørgensen et al. 2016). Garufi et al. (2022) calculated SO₂ column densities by taken an integrated region equivalent to the beam size for the Taurus sources, which varies from 0″.27 to 0″.64. For these sources, we applied a factor (beam_min × beam_maj/1″.0²), where beam_min and beam_maj are the beam size in Garufi et al. (2022), to be consistent with the integration area of 1″.0. This results in large error bars where the lower values correspond to abundances for an integration area of 1″.0 and the upper values were calculated using an integration area equivalent to the beam size.

This section highlights the difficulties in comparing data from different observations with diverse beam sizes, strategies, results, and different kind of information available in the literature. An unbiased survey is, therefore, crucial for a more precise comparison.

References

All Tables

All Figures

Fig. 1 Class 0 sources. Moment 0 maps of CS, SO, 34SO, and SO2. The contours start at 3σ, follow steps of 1σ, and the emission has been integrated over the velocity ranges listed in Table B.1. Maximum values (in units of Jy beam−1 km s−1) are indicated in the bottom-right corner of each panel. Red and blue arrows represent the outflow direction and the yellow star indicates the position of the continuum peak. Magenta stars represent binary components and the synthesized beam for each molecular transition is shown in the panels of the first row. Sources with no detections are not shown in this figure.

In the text

Fig. 2 Class I sources. Moment 0 maps of CS, SO, 34SO, and SO2. The contours start at 3σ, follow steps of 1σ, and the emission has been integrated over the velocity ranges listed in Table B.1. Maximum values (in units of Jy beam−1 km s−1) are indicated in the bottom-right corner of each panel. Red and blue arrows represent the outflow direction and the yellow star indicates the position of the continuum peak. Magenta stars represent binary components and the synthesized beam for each molecular transition is shown in the panels of the first row. Sources with no detections are not shown in this figure.

In the text

Fig. 3 Sources with unknown TЬо1. Moment 0 maps of CS, SO, 34SO, and SO2. The contours start at 3σ, follow steps of 1σ, and the emission has been integrated over the velocity ranges listed in Table B.1. Maximum values (in units of Jy beairr−1 km s−1) are indicated in the bottom-right corner of each panel. Red and blue arrows represent the outflow direction and the yellow star indicates the position of the continuum peak. The synthesized beam for each molecular transition is shown in the panels of the first row. Sources with no detections are not shown in this figure.

In the text

	Fig. 4 Integrated spectra of CS (red), SO 66−55 (blue), and SO2 (green). The integration area for each source is shown in Table B.1. Black arrows in the panels that correspond to Per-emb 12 B, Per-emb 13, Per-emb 27, Per-emb 29, and Per-emb 44 indicate an offset emission, likely related with a transition from an unidentified COM.
In the text

	Fig. 5 Detection rate of CS, SO, 34SO, and SO2 toward Class 0 and Class I sources. Detected refers to emission above 3σ that is observed in the vicinity of the source (not necessarily toward the continuum peak).
In the text

	Fig. 6 Different components where the CS, SO, 34SO, and SO2 transitions are detected. Extended emission is divided into three groups, those sources where the emission is seen perpendicular or parallel to the outflow direction, and those where the structure is unclear (labeled as Extended).
In the text

Fig. 7 Continuum and molecular fluxes (above 3σ) as a function of Tbol (left) and Lbol (right). The gray dashed line represents Tbol = 70, the separation between Class 0 and Class I sources, while the black solid line in the last two panels of the right column indicates Lbol = 1 L⊙. Red circles represent multiple systems (binaries or triples) where the continuum ALMA emission is resolved but they do not have individual Tbol and Lbol measurements from far IR. Blue circles indicate unresolved binaries in continuum ALMA emission.

In the text

	Fig. 8 Molecular column density ratio of SO over 34SO, showing that both SO transitions are optically thick. A range of Tex, between 50 and 150 K has been employed for the calculations. Dark blue and light blue bars indicate ratios calculated from SO 76−65 and SO 66−55, respectively. The red vertical line separates Class 0 (left) from Class I sources (right).
In the text

Fig. 9 Column density ratio between SO2 and SO as a function of Tbol (left) and Lbol (right). Red dots are values calculated from 34SO fluxes and later on multiplied by 22 (Wilson 1999), while blue dots indicate those sources where34 SO is not detected and the ratios are calculated from the SO 66−55 transition. Black triangles represent lower limits (optically thick SO2 emission) and green triangles indicate upper limits, given by those sources where SO emission is detected but SO2 is not. The gray dashed line in the left panel represents Tbol = 70, the separation between Class 0 and Class I sources. Error bars in the y-axis represent ranges for Tех between 50 and 150 K.

In the text

Fig. 10 S-bearing species vs. COMs detection. Sources with detection of SO2 are clustered by the amount of detected COMs in Yang et al. (2021). The numbers on the left of the source names indicate the amount of different COMs detected. The different colors indicate where the emission of CS or SO is seen. Those sources where the CS or SO emission is parallel to the outflow direction are usually rich in COMs, while sources where the CS or SO emission is perpendicular to the outflow direction do not show detection of COMs.

In the text

	Fig. 11 CH3OH abundance as a function of the SO2 abundance for the sources in Fig. 10. The different colors represent those sources where more than five different COMs were detected (green), between one and three different COMs were detected (red), and no COMs (blue) were detected in Yang et al. (2021). Gray dashed lines show a 10:1 and 100:1 linear correlation.
In the text

Fig. 12 Column density ratio between SO2 and H2 for sources in different evolutionary stages. Red, green, and blue symbols represent sources from Perseus, Ophiuchus, and Taurus, respectively. Dots indicate Class 0 sources, squares represent Class I sources, and diamonds show Class I/II and Class II sources. Empty triangles indicate upper or lower limits. Values for IRAS 16295 were taken from Drozdovskaya et al. (2018) and Jørgensen et al. (2016). Perseus sources are from this work, Ophiuchus Class I sources correspond to IRS 67, GSS30-IRS1, IRS 44, IRS 43, and Elias 29, from left to right, studied by Artur de la Villarmois et al. (2019), and Taurus sources are taken from Garufi et al. (2022) and represent DG Tau (outer disk), HL Tau (inner disk), HL Tau (outer disk), IRAS 04302, and T Tau S, from left to right. The upper limit for IRAS 16293 is given by the optically thick continuum emission and lower limits for Perseus and Taurus sources are given by optically thick SO2 emission. The PEACHES sample and Class I Ophiuchus sources are ordered by Tbol.

In the text

	Fig. C.1 Same as Fig. 1 but for large-scale emission of CS toward L1448 IRS 2E and B1-b N.
In the text

	Fig. C.2 Same as Fig. C.1 for Per-emb 12, Per-emb 13 and IRAS 4B2, Per-emb 26 and Per-emb 42, Per-emb 27, and Per-emb 33 and L1448 IRS3A.
In the text

Fig. D.1 Integrated spectra of SO 76−65 (red), SO 66−55 (blue), and 34SO 56−45 (green). The integration area for each source is listed in Table B.1 and the 34SO 56−45 spectra have been multiplied by three, for a better comparison. Black arrows in the panels that correspond to Per-emb 12 B, Per-emb 13, Per-emb 27, Per-emb 29, and Per-emb 44 indicate an offset emission, likely related with a transition from an unidentified COM.

In the text

Fig. E.1 RADEX models for Per-emb 29. Fluxes of the SO2 140,14−131,13 transition for two different kinetic temperatures: 50 K (left) and 150 K (right). The white contour represents the flux observed for this source and the dashed cyan contour shows an H2 number density of 107 cm−3 (lower limit for nH used in the calculations). In this example, the $N_{{\text{SO}}_2}$ present a minimum value of 3 × 1014 cm−2 (for Tex = 150) and a maximum value of 1 × 1015 cm−2 (for Tex = 50), in a range between 107 and 109 cm−3.

In the text

Fig. E.2 RADEX models. Optical depth of the SO2 140,14−131,13 transition for three different excitation temperatures: 50, 100, and 150 K. The cyan contour represents τ = 1, the two dashed white vertical lines indicate the range of values obtained for $N_{{\text{SO}}_2}$ for most of the sources, and the dashed white horizontal line shows an H2 number density of 107 cm−3 (lower limit for nH used in the calculations).

In the text