PRODIGE – envelope to disk with NOEMA

J. Miranzo-Pastor, A. Fuente, D. Navarro-Almaid a, J. Pineda, D. Segura-Cox,
P. Caselli, R. Martin-Domenech, M. Valdivia-Mena, T. Henning, T.-H. Hsieh,
et al.

► To cite this version:

J. Miranzo-Pastor, A. Fuente, D. Navarro-Almaid a, J. Pineda, D. Segura-Cox, et al. PRODIGE – envelope to disk with NOEMA. Astronomy & Astrophysics – A&A, 2025, 700, pp.A251. 10.1051/0004-6361/202554018. [hal-05367634](https://hal.science/hal-05367634v1)

HAL Id: hal-05367634

https://hal.science/hal-05367634v1

Submitted on 19 Nov 2025

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire **HAL**, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# PRODIGE – envelope to disk with NOEMA

## VI. The missing sulfur problem*

J. J. Miranzo-Pastor<fnref n="1" />,<fnref n="2" />,*, A. Fuente<fnref n="1" />, D. Navarro-Almaid<fnref n="1" />, J. E. Pineda<fnref n="3" />, D. M. Segura-Cox<fnref n="4" />,<fnref n="3" />, P. Caselli<fnref n="3" />, R. Martin-Domenech<fnref n="1" />, M. T. Valdivia-Mena<fnref n="3" />,<fnref n="5" />, T. Henning<fnref n="6" />, T.-H. Hsieh<fnref n="3" />,<fnref n="7" />,<fnref n="8" />, L. A. Busch<fnref n="3" />, C. Gieser<fnref n="6" />, Y.-R. Chou<fnref n="3" />, B. Commerçon<fnref n="9" />, R. Neri<fnref n="10" />, D. Semenov<fnref n="6" />, A. Lopez-Sepulcre<fnref n="11" />,<fnref n="10" />, N. Cunningham<fnref n="12" />, L. Bouscasse<fnref n="10" />, and M. Maureira<fnref n="3" />

Received 4 February 2025 / Accepted 2 July 2025

## ABSTRACT

*Context.* Determining the amount of sulfur in volatiles and refractories in the interstellar medium remains one of the main problems in astrochemistry. The detection of H₂S ices, which are thought to be one of the main sulfur reservoirs, is still a great challenge and has not been achieved yet, and the only sulfur-bearing species detected in the ices to date is OCS. The PROtostars and DIsks: Global Evolution (PRODIGE) large survey observations with the NOrthern Extended Millimeter Array (NOEMA) of several Class 0/I protostars in the Perseus Molecular Cloud provide a perfect opportunity to study the H₂S and OCS composition of the ices through the volatiles sublimated in the warm inner core (T>100K, n~10⁶ cm⁻<fnref n="3" />) of these protostars.

*Aims.* Our aim is to determine the H₂S/OCS ratio in the warm inner core of the protostars of our sample in order to study how it is affected by different factors during its evolution.

*Methods.* We used the NOEMA millimeter observations from the PRODIGE program of H₂S, H₂³S, OCS, OC³³S, and OC³⁴S to estimate the H₂S and OCS column densities in the warm inner core of 24 protostars of Perseus. In addition, we used SO and SO₂ data from the Atacama Large Millimeter/submillimeter Array (ALMA) archive to complete the sulfur budget and give a rough estimate of the total sulfur abundance in each of the sources. We explored the chemistry of H₂S and OCS in the warm cores using chemical and dynamical simulations of the collapse of a dense core to form a protostar.

*Results.* The compound H₂S is detected in 21 protostars and OCS in 17 protostars of our sample. The estimated H₂S/OCS ratio reveals a segregation of the sources into “OCS-poor” and “OCS-rich” protostars, where the OCS-poor protostars present higher H₂S/OCS ratios than the OCS-rich ones. The total sulfur abundance, which is always dominated by either H₂S or OCS, grows with evolution during the Class 0 phase, reaching a minimum depletion of a factor less than eight in the Class 0/I objects and decreasing again in the Class I. Chemo-MHD simulations show that temperature changes in the pre-stellar phase and during the collapse can produce substantial differences in the H₂S and OCS (ice and gas-phase) abundances and in the H₂S/OCS ratio.

*Conclusions.* Our analysis shows that the H₂S/OCS ratio is strongly influenced by the environment and the initial conditions of the cloud.

**Key words.** stars: formation – stars: protostars – ISM: abundances – evolution – ISM: molecules

## 1. Introduction

Sulfur is the tenth most abundant element in the Universe, and it is known to play a significant role in biological systems

(Leman et al. 2004; Chen & Yu 2019). Sulfur is found in a wide variety of biomolecules, such as amino acids, nucleic acids, sugars, and vitamins. In fact, along with hydrogen, carbon, oxygen, nitrogen, and phosphorus, it is considered to be one of the six elements fundamental to life. Moreover, some sulfur compounds, such as hydrogen sulfide (H₂S), have been proposed as a necessary catalyst to form amino acids in the interstellar medium (ISM) (Olson & Straub 2016). Nowadays, we are

<footnote n="1">¹ Centro de Astrobiología (CAB), CSIC-INTA, Ctra. de Torrejón a Ajalvir km 4, 28850 Torrejón de Ardoz, Spain</footnote>

<footnote n="2">² Departamento de Física de la Tierra y Astrofísica, Facultad de Ciencias Físicas, Univ. Complutense de Madrid, 28040 Madrid, Spain</footnote>

<footnote n="3">³ Max-Planck-Institut für extraterrestrische Physik, Giessenbachstrasse 1, 85748 Garching, Germany</footnote>

<footnote n="4">⁴ Department of Physics and Astronomy, University of Rochester, Rochester, NY 14627-0171, USA</footnote>

<footnote n="5">⁵ European Southern Observatory, Karl-Schwarzschild-Strasse 2, 85748 Garching, Germany</footnote>

<footnote n="6">⁶ Max-Planck-Institut für Astronomie, Königstuhl 17, 69117 Heidelberg, Germany</footnote>

<footnote n="7">⁷ Taiwan Astronomical Research Alliance (TARA), Taiwan</footnote>

<footnote n="8">⁸ Academia Sinica Institute of Astronomy and Astrophysics, No. 1, Section 4, Roosevelt Road, Taipei 10617, Taiwan</footnote>

<footnote n="9">⁹ Centre de Recherche Astrophysique de Lyon/ENS, Lyon, France</footnote>

<footnote n="10">¹⁰ Institut de Radioastronomie Millimétrique (IRAM), 300 rue de la Piscine, 38406 Saint-Martin d'Hères, France</footnote>

<footnote n="11">¹¹ Univ. Grenoble Alpes, CNRS, IPAG, 38000 Grenoble, France</footnote>

<footnote n="12">¹² SKA Observatory, Jodrell Bank, Lower Withington, Macclesfield SK11 9FT, UK</footnote>

<footnote>* Based on observations carried out under project number L19MB with the IRAM NOEMA Interferometer. IRAM is supported by INSU/CNRS (France), MPG (Germany) and IGN (Spain).</footnote>

<footnote>** Corresponding author.</footnote>
aware of the existence of more than 300 molecules in interstellar regions (Müller et al. 2005; McGuire 2022, CDMS<fnref n="1" />). Out of the currently detected interstellar molecules, only 33 contain sulfur atoms. The apparent paucity of sulfur-bearing molecules detected in the ISM is somewhat reflective of a great problem in astrochemistry: While the observed gaseous sulfur seems to account for its total cosmic abundance (S/H~1.5 × 10⁻⁵, Asplund et al. 2009; Daflon et al. 2009) in diffuse clouds and photodissociation regions (PDRs) (Neufeld et al. 2015; Goicoechea & Cuadrado 2021; Fuente et al. 2024), the sum of the abundances of the gas-phase sulfur-bearing molecules detected constitute less than 1% of the expected amount (Ruffle et al. 1999; Vastel et al. 2018; Rodríguez-Baras et al. 2021). In protoplanetary disks, observations also suggest that only 1% of sulfur is in volatiles (Semenov et al. 2018; Le Gal et al. 2019; Rivière-Marichalar et al. 2020). One could think that in these cold and dense regions, most of the sulfur is locked in the icy mantles that cover the dust grain surfaces. However, the detection of sulfur-bearing species in interstellar ices remains elusive. Nowadays, s-OCS (“s-” means that the molecule is in the icy grain mantles) is the only compound firmly detected in interstellar ices (Palumbo et al. 1995; Boogert et al. 2022). Different authors have published tentative detections of the 7.5 µm band of s-SO₂ (Rocha et al. 2024), but this band is overlapped by intense bands of abundant complex organic molecules (COMs; C-bearing molecules containing six atoms or more) and CH₄, which hinders its conformation (Taillard et al., in prep). Although s-H₂S is the most abundant sulfur species in comets (Calmonte et al. 2016) and is predicted to be the main sulfur reservoir in interstellar ices, it has not been detected in the ISM, yet. One main problem is that its band at 3.9 µm is overlapped by an intense methanol band (Taillard et al., in prep.). By now, only upper limits to its abundance relative to water, which are within the range of ~0.1–1%, have been determined (Jiménez-Escobar & Muñoz Caro 2011; Rocha et al. 2024). These upper limits are roughly consistent with some chemical models (see Taillard et al. 2025) and suggest that H₂S might not be as abundant as predicted, at least in the observed environments. It has been proposed that H₂S is destroyed within the ice to eventually form other sulfur compounds such as OCS (el Akel et al. 2022). Laas & Caselli (2019) modeled sulfur chemistry assuming enhanced accretion of cations (in particular, S⁺) on negatively charged grains and found that in timescales similar to the free-fall time, some species such as s-SO, s-OCS, or s-HSO could be more abundant than s-H₂S.

All in all, the abundances of the species detected in gas and ice contribute to less than 5% of the total sulfur, leaving around 95% unaccounted for. Several theories have been proposed to explain this “missing” sulfur. One possibility is that sulfur exists as neutral atomic sulfur, which is undetectable under the conditions typical of molecular clouds. Indeed, high abundances of atomic sulfur have been detected toward the Orion Bar by Goicoechea & Cuadrado (2021); Fuente et al. (2024). Also, Hily-Blant et al. (2022) showed, indirectly with the NS/N₂H⁺ ratio, that atomic sulfur could be the main carrier in dense cores. Another compelling hypothesis is that sulfur is sequestered in sulfur allotropes and hydrogen sulfides (Sₓ, H₂Sₓ), which could serve as semi-refractory sulfur reservoirs (Jiménez-Escobar et al. 2012; Shingledecker et al. 2020; Fuente et al. 2023; Cazaux et al. 2022; Carrascosa et al. 2024). This idea is supported by the detection of these compounds in comets (Calmonte et al. 2016)

and meteorites (Aponte et al. 2023) as well as by the detection of S₂H toward the Horsehead nebula (Fuente et al. 2017a). Other proposed reservoirs include ammonium hydrosulfide (Altwegg et al. 2022) and iron sulfides that have been detected in protoplanetary disks and meteorites (Keller et al. 2002). By now, the primary sulfur reservoir in cores – where the depletion timescale is shorter than the dynamical timescale – remains unidentified, which highlights a critical gap in our understanding of sulfur chemistry.

Millimeter observations of gas-phase molecules can help disentangle which is the most abundant sulfur compound in ices. We know that ices are formed in the cold pre-stellar phase where grain temperature decreases to less than 10 K and most molecules become frozen on their surfaces (Walmsley et al. 2004; Pineda et al. 2022; Caselli et al. 2022). As the collapse proceeds, a young protostellar object is formed in the interior of the dense core and heats the surrounding material. Dust and gas temperatures of hundreds of K can be reached in the inner ~100 au around the central protostar (Martín-Doménech et al. 2021; Bianchi et al. 2022). As the temperatures increase above 100 K, the water-rich ice mantles sublimate, injecting molecules into the gas phase (see, e.g., Collings et al. 2003). In some cases these warm regions present a rich chemistry with the detection of several COMs. In these cases, we refer to this compact warm region as the hot corino (Ceccarelli 2007; Jørgensen et al. 2020). The chemical composition of these warm regions should reflect to some extent the composition of the ices formed during the collapse, modulated by the effect of the different desorption mechanisms at play in the region. Binding energies determine the sublimation temperature of the different species, and one would expect a layered structure where the most volatile compounds appear at larger radii, while the less volatile ones are only detected in the hot interior (Ruaud et al. 2016). Grain surface molecules can also be released to the gas-phase due to the accretion shocks formed when the envelope material falls onto the disk (Artur de la Villarmois et al. 2022) and/or due to the energetic shocks produced by the bipolar outflow Caselli et al. (1997); Schilke et al. (1997); Jiménez-Serra et al. (2008); Holdship et al. (2016). Interferometric observations are needed to spatially resolve the interior of protostars and to determine the main desorption mechanism at work, a mandatory step to infer the connection between the gas and ice chemistry.

In Sections 2 and 3, we present high spatial resolution observations (~, corresponding to 1" ~ 300 au) of H₂S and OCS of 24 protostars located in Perseus. These compounds are expected to be the most abundant sulfur species in ices and present similar binding energies ($E_{bin}$ = 2400 K for OCS and 2700 K for H₂S, Wakelam et al. 2017), which ensures a fair comparison. In Section 4, we estimate the column density and abundances of H₂S, OCS and their isotopologues in all the sources of the sample, and in Section 5, we determine the H₂S/OCS to explore possible variation of the ice composition along the sample. In Section 6, we combine our observations with previous SO and SO₂ observations to determine the sulfur budget. In Section 7, we run a set of theoretical simulations to explore the parameters that could affect the different H₂S/OCS ratios. Finally, in Section 8, we discuss the astrochemical implications of our results.

## 2. Observations

This paper is based on observations of the MPG-IRAM Observing Program PRODIGE (PROtostars and DIsks: Global

<footnote n="1">¹ https://cdms.astro.uni-koeln.de/</footnote>
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Map of the Perseus Molecular Cloud and its several subregions. Represented with stars are the positions of the 24 protostars observed in this work. Their color shows the evolutionary stage of the source. Gray-scale background shows far-infrared and submillimeter emission from the Herschel Gould Belt Survey (André et al. 2010; HGBS team 2020).</figcaption></figure>

Evolution, Project ID: L19MB)<fnref n="2" /> PIs: Paola Caselli, Thomas Henning. Within this program, a series of 32 Class 0/I protostars have been observed with the NOrthern Extended Millimeter Array (NOEMA). In this program, the Band 3 receiver and the PolyFix correlator were used. PolyFix provides ~15.5 GHz of bandwidth (divided into two 7.744 GHz wide sidebands separated by 15.488 GHz). The whole band was observed with a channel width of 2 MHz. Moreover, several windows with high spectral resolution, 62.5 kHz channel width, were placed to observe selected lines.

We processed the NOEMA observations from uncalibrated data using the standard observatory pipeline in the Grenoble Image and Line Data Analysis Software (GILDAS) package Continuum and Line Interferometer Interferometer Calibration (CLIC). We used 3C84 and 3C454.3 as bandpass calibrators. 0333+321 and 0333+322 were used for phase and amplitude calibration, with observations for these calibrators taken every

20 min. Flux calibration sources were LKHA101 and MWC349. The uncertainty in flux density was 10%. The continuum was bright enough to allow for phase self-calibration. More details about calibration and imaging can be found in Section 2.1 of Gieser et al. (2024).

This work relies in a subset of 24 Class 0/I protostars from the previously mentioned observations (see Table D.1, Fig. 1), which were the ones observed and reduced between 2019 and 2023. This includes all Class 0 (the younger) sources and also some of the Class 0/I and I protostars.

## 3. Results

In this section, we show the moment 0 maps of the 24 sources in our sample, the spectra of H₂S, OCS(18–17) and OC³S (see Table 1) and their Gaussian fits, and we list the detection rates of each species, which are also shown in Fig. 2.

In Fig. 3, we present the moment-0 maps of the 24 protostars, showing the emission of the main species H₂S, OCS(18–17) and

<footnote n="2">² https://iram-institute.org/science-portal/proposals/lp/miop/miop-002/</footnote>
<table><caption>**Table 1.** List of studied NOEMA lines in order of increasing wavelength with atomic data retrieved from Splatalogue<fnref n="3" /> and CDMS<fnref n="4" />.</caption><thead><tr><th>ν (GHz)</th><th>Species</th><th>Transition Upper-lower</th><th>A<sub>ij</sub></th><th>E<sub>u</sub> (K)</th><th>Δν (kHz)</th><th>RMS (mJy/beam)</th></tr></thead><tbody><tr><td>215.5029</td><td>H<sub>2</sub><sup>33</sup>S</td><td>2<sub>2,0,4</sub>–2<sub>1,1,4</sub></td><td>–4.374</td><td>81.6</td><td>2000</td><td>[2.0–2.9]</td></tr><tr><td>216.7104</td><td>H<sub>2</sub>S</td><td>2<sub>2,0</sub>–2<sub>1,1</sub></td><td>–4.312</td><td>84.0</td><td>2000</td><td>[2.0–2.9]</td></tr><tr><td>216.1474</td><td>OC<sup>33</sup>S</td><td>J = 18–17</td><td>–4.420</td><td>98.6</td><td>62.5</td><td>[15–60]</td></tr><tr><td>218.9034</td><td>OCS</td><td>J = 18–17</td><td>–4.517</td><td>99.8</td><td>62.5</td><td>[10–50]</td></tr><tr><td>231.0615</td><td>OCS</td><td>J = 19–18</td><td>–4.446</td><td>111</td><td>62.5</td><td>[20–80]</td></tr><tr><td>237.2736</td><td>OC<sup>34</sup>S</td><td>J = 20–19</td><td>–4.411</td><td>120</td><td>2000</td><td>[2.5–3.5]</td></tr></tbody></table>

**Notes.** For H<sub>2</sub><sup>33</sup>S and OC<sup>33</sup>S, the A<sub>ij</sub> were calculated from the available data in CDMS. Δν is the channel width with which the transition was observed. In the last column (RMS) we give a range for the typical noise in the datacubes of each transition.

<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Detection rate of the different molecules in the warm inner core of the 24 sources of our sample. Of the total 24 protostars, 17 of them are Class 0 sources, 3 of them are classified as Class 0/I, and the last 4 sources are Class I.</figcaption></figure>

OCS(19–18) lines. We created the moment-0 maps by integrating the data cubes over specific velocity ranges for each map. Integration was done over a threshold of 3σ, for each channel individually. For a more accurate description of the maps, the different velocity ranges and σ values are detailed in the corresponding table of the repository referenced in Appendix D. The emission is presented in the form of a contour map, within a window of 4''×4'', centered in the peak of the continuum (see Table D.1). The contour lines delineate regions were emission is within 10%, 30%, 50%, 70% and 90% of the peak emission. In addition, to complement this representation, we also give the maximum value of the emission and the size of the beam, represented by a purple ellipse in the bottom-left corner. We have also marked the position of the peak of the continuum, which is represented with a purple star in the center of each image.

In the maps, we see that the H<sub>2</sub>S line usually appears centered in the peak of the continuum, showing an approximately symmetric distribution. In the sources where OCS was not detected, H<sub>2</sub>S emission may have a more irregular shape. We also find that OCS is only observed in sources where H<sub>2</sub>S was detected, and that OCS(18–17) and OCS(19–18) are always detected together, although sometimes they may have rather different distributions.

Regarding the isotopologues, we find that many of the OCS detections are also detected in OC<sup>34</sup>S. This is a clear indicator

of optically thick OCS lines, as we show later in Section 4. We assumed an ³²S/³⁴S isotopic ratio of 22.5 (Anders & Grevesse 1989). ³³S is even scarcer than its 34 isotope counterpart, but the value of the ³⁴S/³³S ratio has not been uniformly established, as it may change depending on the region. Standard values are 5.61 (solar, Anders & Grevesse 1989) and 6.27±1.01 (Chin et al. 1996).

In this work, we are interested in the H<sub>2</sub>S/OCS ratios in the warm inner core (T>100K, n~10⁶cm⁻<fnref n="3" />) of 24 Class 0/I protostars in Perseus. To carry out this analysis, we needed to integrate the spectra of only the most interior region of the protostar, where the core is found. By comparing with the typical extension of the H<sub>2</sub>S (2<sub>2,0</sub>–2<sub>1,1</sub>) emission, we decided to extract the spectra of a circular region of 1.5'' in diameter around the peak of the continuum for every source. After extracting all the spectra, we fitted a Gaussian profile to the detected lines using the CLASS software from GILDAS<fnref n="5" />.

In Fig. 4, we present the integrated spectra of the H<sub>2</sub>S 2<sub>2,0</sub>–2<sub>1,1</sub>, the OCS (18–17) and the OC<sup>34</sup>S (20–19) lines in the warm inner core (1.5''×1.5'') of the 24 protostars from our sample. We find clear detections (>5σ) of both H<sub>2</sub>S and OCS in most of the sources while OC<sup>34</sup>S detections are usually more faint (>3σ) or tentative (<3σ). We note that, after integrating the spectra over the 1.5'' region, measured emission can be over 3–5σ due to the increased signal to noise ratio after the integration, even if emission inside one beam was below 3σ for some channels.

We have also extracted the spectra of the OCS (19–18), the OC<sup>33</sup>S (18–17) and the H<sub>2</sub><sup>33</sup>S 2<sub>2,0</sub>–2<sub>1,1</sub> lines. H<sub>2</sub><sup>33</sup>S has hyperfine structure and its emission is divided in 10 different lines. We show an example of the H<sub>2</sub><sup>33</sup>S 2<sub>2,0</sub>–2<sub>1,1</sub> spectral distribution in Fig. A.1. OC<sup>33</sup>S also has hyperfine structure; however, for high J transitions, the line splitting cannot be resolved, and we do not account for it in our fits and calculations.

Table D.2 summarizes the detected and undetected lines in the 24 protostars of Perseus observed in this work. We can see that in all the sources where OCS was detected, H<sub>2</sub>S emission could also be observed, but not vice versa (see B5-IRS1, L1448NW, Per-emb-08). Additionally, all OCS detections were observed in both J=18–17 and J=19–18 transitions.

Regarding the 17 Class 0 sources, we observed the H<sub>2</sub>S 2<sub>2,0</sub>–2<sub>1,1</sub> line in 16 sources – although B1bN, the only Class 0 source not observed in H<sub>2</sub>S, was actually detected in its isotopologue, H<sub>2</sub><sup>33</sup>S – and we detected both OCS lines in 14 of them. That is a ~94% detection rate of H<sub>2</sub>S and an ~82% detection rate of OCS in Class 0 protostars (or ~88% detection rate of OCS in Class 0 protostars with H<sub>2</sub>S). On the other hand, for the seven Class 0/I and Class I protostars, we detected H<sub>2</sub>S in four of the sources and we observed OCS emission in only three of them. That corresponds to a ~57% detection rate of H<sub>2</sub>S and a ~43% for Class 0/I and Class I protostars (or ~75% detection rate of OCS in sources with H<sub>2</sub>S emission). Globally, we have a ~83% detection rate of H<sub>2</sub>S and a ~71% detection rate of OCS (or, equivalently, a ~85% detection rate of OCS in protostars with H<sub>2</sub>S).

When it comes to the isotopologues, in the Class 0 protostars, we found OC<sup>34</sup>S emission toward 10 of the 14 sources with OCS emission, and we detected OC<sup>33</sup>S in seven of them. This is a ~71% detection rate of OC<sup>34</sup>S and a ~50% detection rate of OC<sup>33</sup>S in Class 0 protostars with OCS (or ~70% detection rate of OC<sup>33</sup>S where detected in OC<sup>34</sup>S). For the Class 0/I and Class I objects, we found OC<sup>34</sup>S in only three of the seven

<footnote n="3">³ https://splatalogue.online</footnote>

<footnote n="4">⁴ https://cdms.astro.uni-koeln.de/</footnote>

<footnote n="5">⁵ https://www.iram.fr/IRAMFR/GILDAS</footnote>
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Emission of the main species ($H_2S$ and both OCS lines) in the inner core of protostars B1bS, HH211MMS, IRAS2A, and IRAS4B. The color map represents the $>3\sigma$ emission integrated images in a 4''x4'' square region. The color scale, shown at the right of each map, is the brightness temperature in Kelvin. The white contours represent 10%, 30%, 50%, 70%, and 90% of the peak temperature. The pink star shows the position of the protostar, determined by the position of the maximum emission in the continuum (Tobin et al. 2016b). The green cross marks the point with maximum emission of the line.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Spectra of the $H_2S$ (2,0–2,1), OCS (18–17), and OC$^{34}$S (20–19) lines in protostars B1bS, HH211MMS, IRAS2A, and IRAS4B. The red dotted lines represent the values of a Gaussian fit to each of the detected lines. In some cases, the combination of two Gaussian fits was necessary for the fit. The OCS (18–17) line was observed with a spectral resolution of 62.5 kHz, while the $H_2S$ and OC$^{34}$S lines were observed with a spectral resolution of 2 MHz. Gaussian fits have been plotted for the $>3\sigma$ detections except for the OC$^{34}$S, for which we also show the $<3\sigma$ fits when the OC$^{33}$S counterpart was detected with $>3\sigma$.</figcaption></figure>

protostars, and OC$^{33}$S in just two of them. This makes a ~38% detection rate of OC$^{34}$S and a ~25% detection rate of OC$^{33}$S in Class 0/I and Class I sources where OCS was observed (or ~67% detection rate of OC$^{33}$S in the sources with OC$^{34}$S emission). In general, we find that OC$^{34}$S has a detection rate of ~76% and that OC$^{33}$S has a ~53% detection rate in protostars with detected OCS (which is a ~69% detection rate of OC$^{33}$S in sources with OC$^{34}$S). Additionally, we found that six sources presented $H_2^{33}S$ emission in the warm inner core, which is a 30% detection rate of $H_2^{33}S$ in the 20 sources where $H_2S$ was detected. Also, one more $H_2^{33}S$ detection corresponding to the $<3\sigma H_2S$ emission of B1bN was found. This makes a total of seven detections and a detection rate of ~29% in the full sample. Of these seven detections, five of them were in Class 0 protostars and two in Class 0/I sources, making it a ~29% detection rate in both Class 0 and Class 0/I + I sources.

Here, we make a brief analysis of the sources in our sample. All 24 protostars are located in the Perseus Molecular Cloud (Fig. 1), and segregated in several smaller denser regions: Barnard 1, IC348, NGC1333, L1448 and Barnard 5.

### 3.1. B1bN

B1bN is the northern member of the multiple system embedded in B1b (Hirano et al. 1999), inside the Barnard 1 region (comprised by objects B1bN, B1bS and Per-emb-41). The low velocity of the outflow and the high degree of turbulence categorize the object as a first hydrostatic core (FHSC) (Gerin et al. 2015), even younger than its companion B1bS, as B1bN almost free-of-emission spectra suggests (Marcelino et al. 2018). There is a tentative detection of $H_2S$ below $3\sigma$, but the detection of the isotopologue $H_2^{33}S$ (with ~$3\sigma$) suggests a high opacity for the main species (see Fig. E.2a). OCS or its isotopologues are not detected in the protostar.

### 3.2. B1bS

B1bS is the southern member of B1b (Hirano et al. 1999), and is separated from B1bN by 17.4'' (~5218 au; Tobin et al. 2018). Interferometric NOEMA and ALMA observations identified B1bS as a very early Class 0 object due its young and slow outflow (Hirano & Liu 2014; Gerin et al. 2015) and small
and compact structures, including an incipient disk (Gerin et al. 2017). Additionally, Fuente et al. (2017a) provided evidence of pseudo-disk rotation by imaging this source in NH<sub>2</sub>D emission using the NOEMA interferometer. Unlike its northern counterpart, B1bS is considered a hot corino, showing a high chemical richness in COMs (Marcelino et al. 2018). This richness is also present in our analysis. In Fig. 3a, we find the observed S-bearing molecules emission centered in the peak of the continuum. All species, both H<sub>2</sub>S and OCS and all three isotopologues, are detected. Similarly to B1bN, the lines of the main species are optically thick (Table D.3). The OCS line shown in Fig. 4a seems to present very faint wings, which can be evidence of the molecular outflow.

### 3.3. B5-IRS1

B5-IRS1 is a protostellar object in the Barnard 5 star-forming core. It is identified as a Class I protostar from its spectral energy distribution (Enoch et al. 2009). The protostar presents a wide-angle bipolar outflow (Zapata et al. 2014), and a streamer injecting new material from the filaments into B5-IRS1, characterized by Valdivia-Mena et al. (2023). They also found processed and unprocessed material from the core infalling toward the filaments. Fig. E.1b shows an asymmetric distribution of the H<sub>2</sub>S emission, perpendicular to the outflow direction, which could be explained by the material flow from the streamer (see Fig. 14 from Valdivia-Mena et al. 2023). H<sub>2</sub>S spectrum presents a rather irregular profile (see Fig. E.2b). This, again, may be produced by the streamer or by the outflow present in B5-IRS1 – or both phenomena combined. We did not find emission of OCS or from the H<sub>2</sub>³S isotopologue coming from the warm inner core.

### 3.4. HH211MMS

HH211MMS, located in the IC348 complex, is one of the most studied and well-characterized outflows, featuring a highly collimated jet and a molecular outflow (Gueth & Guilloteau 1999; Lee et al. 2009; Ray et al. 2023; Caratti o Garatti et al. 2024). It shows a wide variety of atomic and ionic line emission, as well as molecular and PAH features (Dionatos et al. 2010, 2018). The central source has been classified as a low-mass, low-luminosity Class 0 protostar (Froebrich 2005). In Fig. 3b, we see the warm inner core of HH211MMS, where we detect irregular distributions of H<sub>2</sub>S and OCS. H<sub>2</sub>S emission is extended in the outflow direction, while OCS spatial distributions are much more compact. We also find a tentative detection (<3σ) of the OC³S isotopologue in the warm inner core (Fig. 4b).

### 3.5. IC348MMS

IC348MMS is another Class 0 object in the IC348 region. It has a close companion, JVLA3a (Rodríguez et al. 2014), 2.95" south-west, but we have not been able to observe it in any of the lines. Palau et al. (2014) suggests that JVL3a may be coincident with the origin of the outflow. Previous observations have not found any substructure in IC348MMS (Tobin et al. 2015b). We detect emission of all main species and isotopologues. The spatial distribution of all the species is centered in the peak of emission, as shown in Fig. E.1c. The detection of high column densities of OCS and H<sub>2</sub>S isotopologues suggests a high opacity in the lines of the main species (Table D.3).

### 3.6. IRAS2A

IRAS2A (see, e.g., Jennings et al. 1987) is a Class 0/I protostellar system (Brinch et al. 2009) in the NGC1333 region that hosts a protobinary as revealed by the high-resolution observations of the VANDAM survey (Tobin et al. 2015a), separated by only 0.62" (~186 au; Tobin et al. 2018). These two sources produce the two outflows, orthogonal to one another, known to originate from this system (see, e.g., Jørgensen et al. 2004; Tobin et al. 2015a). Molecular emission in COMS has been observed toward this system (Jørgensen et al. 2005; Maury et al. 2014). IRAS2A also forms a wide binary with IRAS2B. The moment 0 maps in Fig. 3c show largely extended emission for all the lines, especially in the south-north outflow direction (see also Fig. B.1), and centered at the emission peak. The spectra of the OCS lines, Fig. 4c, clearly reveal the presence of a bipolar outflow. The H<sub>2</sub>S line may present this same characteristic, but the lower spectral resolution complicates the interpretation. We also detect emission in the inner core from all three isotopologues.

### 3.7. IRAS2B

IRAS2B is a Class I protostar in the NGC1333 region, with the wide companion IRAS2A separated by 31.4" (~9426 au; Tobin et al. 2018). IRAS2B also has a close companion, separated by 0.31" (~93 au; Tobin et al. 2018), that we cannot resolve in our observations. In the observations, we can only detect the main species, H<sub>2</sub>S and the two OCS lines. The H<sub>2</sub>S moment 0 map in Fig. E.1d shows a barely resolved emission along the beam's major axis direction and an asymmetric spatial distribution along the beam's minor axis direction. The OCS emission in this source is faint, but above the 5σ threshold (Fig. E.2d). The spatial distributions of these species are also very different which suggests that the OCS and the H<sub>2</sub>S may trace different regions of the protostar (e.g., disk and core).

### 3.8. IRAS4B

IRAS4B is a Class 0 young stellar object (Sandell et al. 1991; Hirota et al. 2008; Enoch et al. 2009) characterized by a deeply embedded 0.24 M<sub>☉</sub> disk (Jørgensen et al. 2009) inside a 2.9 M<sub>☉</sub> envelope (Jørgensen et al. 2002). It is located in the NGC1333 region and there are several protostars in its near vicinity: IRAS4A (not studied in this work), IRAS4B and IRAS4C, and their close companions (IRAS4A' and IRAS4B'). IRAS4B is located 29.74" away from IRAS4A (~8922 au; Tobin et al. 2018) and ~54" away from IRAS4C. IRAS4B close companion, IRAS4B', is separated by 10.65" (3196 au; Tobin et al. 2018; Looney et al. 2000) east of the main source. This embedded system exhibits hot corino chemistry, with emission of COMS as shown in, for instance, the Perseus ALMA Chemistry Survey (PEACHES) (Yang et al. 2021). Water emission has also been identified in its outflow (Herczeg et al. 2012). IRAS4B has very spread emission of the two main species. The moment 0 maps in Fig. 3d show extended emission of H<sub>2</sub>S, OCS(18–17) and OCS(19–18) in the core but also in two lobes in the outflow direction (see also Fig. B.2): both OCS maps show a main emission centered in the peak of the continuum and second emission clearly visible in the northern part of the map; in the case of the OCS lines, also a third emission region appears in the south. These additional outer detections represent material that belongs to the outflow. On the other hand, the isotopologues are only detected in the warm inner core. High column densities have been estimated for all isotopologues, but specially H<sub>2</sub>³S,
revealing high opacity particularly in the H₂S (Table D.3). The irregular but slim shape of the spectra in Fig. 4d suggests that the velocity of the outflow is low. These narrow lines could also appear if the outflow is almost on the plane of the sky or if the emission is coming from the inner cavity walls.

### 3.9. IRAS4C

IRAS4C is a Class 0 protostar in the NGC1333 nebula (Looney et al. 2000). Although it is in the IRAS4 subregion, it is quite far from the nearest protostar, IRAS4B, which is separated by 54'' (Fig. E.1e), and the OCS emission is faint, but above the 5σ threshold too (Fig. E.2e). None of the isotopologues are observed in IRAS4C.

### 3.10. L1448-IRS3A

L1448-IRS3A is a Class I protostar (see, e.g., Reynolds et al. 2021) in the L1448 star-forming region. It forms a multiple system with other 5 sources: L1448-IRS3B's three components and L1448NW binary; it is located 7.32'' (~2195 au; Tobin et al. 2018) away from L1448-IRS3B. Recently, Gieser et al. (2024) found an elongated gas bridge with bright emission that connects the IRS3A and IRS3B protostellar systems. They also revealed infalling material toward IRS3A from both the red- and blueshifted sides. In Fig. E.1f we present the moment 0 maps of L1448-IRS3A, where we observe the object in H₂S and OCS, with the OCS emission localized in a compact region around the peak of the continuum. The H₂S emission is also compact but more stretched, in the direction of the outflow, and displaced from the peak of the continuum. In Fig. E.2f, we see a tentative detection of OC³⁴S but we do not observe either OC³³S nor H₂³S. The four lines observed in this source present a double peak (tentatively for the OC³⁴S line), which could be the product of a bipolar outflow or could be tracing the infalling material revealed by Gieser et al. (2024). The double peaked emission could also appear due to the disk/envelope rotation.

### 3.11. L1448-IRS3B

L1448-IRS3B is a Class 0 triple protostar in L1448 region (Tobin et al. 2016b), with the separation of the two companions from the main component being ~79 and ~238 au (Tobin et al. 2018). The formation of the triple system has been associated with the fragmentation of a gravitationally unstable disk, supported by the detection of prominent spiral arms surrounding the three protostars (Tobin et al. 2016a). The protostellar system is embedded, together with L1448-IRS3A, in extended molecular gas that extends well beyond 6600 au (Gieser et al. 2024). The moment-0 maps of L1448-IRS3B, Fig. E.1g, show extended emission of H₂S in the outflow direction, but a compact spatial distribution of OCS. In Fig. E.2g, we present the spectra of the six observed lines, with the detection of the two main species, H₂S and OCS, and also a tentative detection (<3σ) of the OC³⁴S isotopologue.

### 3.12. L1448C

L1448C is a Class 0 object in the L1448 region (Tobin et al. 2007) and a well known example of an outflow-driving source with extremely high velocity (EHV) jets. These jets have been detected in CO (Bachiller et al. 1990) and in several transitions of SiO (Bachiller et al. 1991; Dutrey et al. 1997). High angular resolution observations (Hirano et al. 2010) revealed that L1448C consists of two sources, L1448C(N) and L1448C(S),

with L1448C(N) being responsible for the EHV jets. In Fig. E.1h we show the 0-moment maps of the H₂S and the two OCS lines observed toward this source. The peak of emission of the three lines is centered in the peak of the continuum. In the spectra, shown in Fig. E.2h, wings are present in the OCS 18–17 line and linewidths are much larger than most objects in the sample. This suggests that the emission of these molecules may be perturbed by the extremely high-velocity jets coming from the source and the close binary.

### 3.13. L1448NW

L1448NW is a binary Class 0 protostar in the L1448 region (Barsony et al. 1998). Its two components are separated by 0.25'' (~58 au, Tobin et al. 2016b) and the binary is 21.5'' (~6450 au) away from L1448-IRS3B (Tobin et al. 2018). The source presents a bipolar outflow, were the blueshifted emission is much less visible possibly due to an asymmetric distribution of the surrounding gas (Lee et al. 2015). In the observations of the source (Fig. E.1i, we have only detected H₂S emission in a stretched and compact distribution perpendicular to the outflow direction. We also present the spectra of the lines in Fig. E.2i.

### 3.14. Per-emb-02

Per-emb-02 is a binary protostar in the Barnard 1 region and has its two components separated by only 0.08'' (~24.0 au; Tobin et al. 2018). It is a chemically poor Class 0 YSO (Yang et al. 2021), as no COMs have been detected besides CH₃OH. A streamer feeding this source with chemically fresh gas has been detected in carbon chain emission (Taniguchi et al. 2024; Pineda et al. 2020), and outflows in the northwest-southeast direction were detected (Stephens et al. 2019). This protostar is located in the Barnard 1 region. In Fig. E.1j, we present the moment-0 maps of Per-emb-02, where we find extended emission of H₂S. In Fig. E.2j, we also show the detection of the two OCS lines in the warm inner core of the source. There are no signs of emission of any of the isotopologues.

### 3.15. Per-emb-05

Per-emb-05 is also a Class 0 binary protostar, with both components separated by 0.10'' (~29.1 au; Tobin et al. 2018). The evidence of a companion embedded within the extended dust continuum structure is interpreted as evidence of ongoing disk fragmentation by Tobin et al. (2016b), though the surrounding structure is not confirmed to be rotationally supported. We present the H₂S and OCS moment-0 maps in Fig. E.1k. OCS emission is localized in a compact region around the peak of the continuum while H₂S emission is slightly more extended. Part of this OCS emission is stretched perpendicular to the south-eastern outflow, and is specially visible for the J=19–18 transition. OC³⁴S is clearly detected in the spectrum, displayed in Fig. E.2k. We think the emission corresponds to OC³⁴S because it is centered at the same velocity as the other lines; however, it is not clear if the larger width of the line is only produced by the isotopologue or also by other nearby transitions that may have contaminated the spectrum. OC³³S is tentatively detected within ~3σ.

### 3.16. Per-emb-08

Per-emb-08 is a Class 0 protostar in the IC348 region that forms a multiple system with binary Class I object Per-emb-55, separated by 9.3'' (~2870 au; Tobin et al. 2018), which is not
studied in this work. It is associated with a north-south bipo-
lar outflow and is surrounded by a warped protostellar envelope,
suggesting a highly turbulent environment or strong interactions
between this source and its companions in Per-emb-55 (Lin et al.
2024). It is chemically poor, as COMs have not been detected
in this source (van Gelder et al. 2022). The moment-0 maps,
Fig. E.11, reveal H₂S emission in the central region of the proto-
star. The spatial distribution is very asymmetric, with two peaks
of emission, west and east of the source's position, and a third
weaker peak north-east of the source. The peculiar distribution
might be related with the important disk present in this source,
which is very differentiated from the warm inner core (Tobin
et al. 2018). This could also explain the potentially double-
peaked emission in the H₂S spectrum, displayed in Fig. E.2l –
note that it is very uncommon to find a double peaked emis-
sion for H₂S in this work due to the lower spectral resolution for
this line. A large velocity gradient in Per-emb-08 (>10 km/s) is
expected.

3.17. Per-emb-18

Per-emb-18 is a Class 0 binary protostar in the NGC1333
region, with both components separated by only 0.085" (~25.6
au; Tobin et al. 2018). It forms a wider multiple system with
Per-emb-21, separated by 13.2" (~3976 au), and Per-emb-49,
which is 27.5" (~8240 au) away from Per-emb-18 (Tobin et al.
2018) – the two additional sources are not studied in this work.
It was previously observed in CS, SO and SO₂ in Zhang et al.
(2023), were SO was found to trace asymmetric accretion shocks
near the edge of the circumbinary disk. Fig. E.1m displays the
moment-0 maps of the source, where we can find similar emis-
sion of OCS and H₂S, with the peak of emission somewhat
displaced to the south-west from the peak of the continuum.
In the spectra, Fig. E.2m, we find a double-peaked emission
in both OCS lines, supported by the broader width of the H₂S and
OC³⁴S spectra. This line profile could be produced by the close
binary, by a strong outflow driven by the protostar, or enhanced
by the accretion shocks and an alignment with magnetic fields,
as suggested in Zhang et al. (2023).

3.18. Per-emb-22

Per-emb-22 is a Class 0 binary protostar in the L1448 region, with both components separated by 0.75" (~225 au; Tobin et al. 2018), which we cannot spatially resolve. It is a source with observed outflow activity in continuum emission whose shape resembles that of outflow cavities (Yang et al. 2021). Emission of H2S is detected in the inner core of the source, but rather displaced to the west of the peak of the continuum, as shown in Fig. E.1n. OCS emissions show a smaller but similar distribution, with the addition of some more extended emission in the NE-SW direction, perpendicular to the outflow. The spectra in Fig. E.2n might suggest the presence of an outflow, but the low level of OCS emission, and the low spectral resolution of the H2S observations make it difficult to determine certainly. None of the isotopologues are observed in this object.

3.19. Per-emb-29

Per-emb-29 is a Class 0 object in the Barnard 1 region. It is a hot
corino that has been observed to have rich emission lines of gas-
phase COMs (van Gelder et al. 2020), as well as a high velocity
outflow (Jørgensen et al. 2006). Chen et al. (2024) recently
detected several COMs also in the form of ices. We present the

moment-0 maps of Per-emb-29 in Fig. E.1o, were we find a reg-
ular and round distribution around the peak of the continuum
for the three main lines. The spectra (Fig. E.2o) show a self-
absorption in OCS, which is produced by a high opacity. The
detection of the OC³⁴S isotopologue, which peaks inbetween the
two peaks of the OCS line, supports this fact. In this source, we
also detected emission of the H₂³³S isotopologue.

3.20. Per-emb-30

Per-emb-30 is a very faint Class 0/I protostar in the Barnard 1 region. We have not detected any of the transitions in this object.

3.21. Per-emb-50

Per-emb-50 is a Class I protostar in the NGC1333 region. It is a very bright and massive object, compared with other Class I protostars (around ×10, Enoch et al. 2009; Fiorellino et al. 2021), and emission of some sulfurated species has been observed in the warm inner core (e.g., Valdivia-Mena et al. 2022). However, emission of H₂S and OCS has not been detected in our work (see Figs. E.1q, E.2q). A streamer depositing material close to the edge of the gas disk from roughly 1500 to 3000 au from the protostar was also identified in Valdivia-Mena et al. (2022).

3.22. Per-emb-62

Per-emb-62 is a Class I protostar in the IC348 region. We have
not detected any of the species in this source either.

3.23. SVS13A

SVS13A is a Class 0/I protostar in the NGC1333 region. It belongs to the SVS13 multiple system, whose components are SVS13A, SVS13B and SVS13C – this last one not being studied in this work. SVS13A is a triple system where the two main components are separated by 0.3" (~90 au, Tobin et al. 2018). The third component, SVS13A2, is located 5.3" south-west of the main source, SVS13A1; however, we do not detect it in our observations. The three components are surrounded by a molecular envelope (Lefloch et al. 1998) and drive a molecular outflow associated with the Herbig-Haro object HH711 (Reipurth et al. 1993) as well as other younger flows (Lefèvre et al. 2017). The system has been associated with a hot corino using deuterated water observations (Codella et al. 2016), and, more recently, many sulfurated species were observed in the warm inner core (Codella et al. 2021). Several COMs were also detected in this source, likely not only arising from the hot corino, but also from shocked gas at disk scales (Hsieh et al. 2024). A streamer, possibly infalling, with a length of ~700 au was identified by Hsieh et al. (2023) and can be observed with multiple tracers. SVS13A moment-0 maps, in Fig. E.1s, present extended emission in both H₂S and OCS, centered in the peak of the continuum. OCS(18–17) spectrum (Fig. E.2s) shows a self-absorption feature which, supported by the detection of the OC³⁴S isotopologue that peaks inbetween the two peaks of the OCS line. H₂³³S is also detected in this source. Apart from this, OCS(18–17) slightly asymmetric spectrum could be a sign of the present outflows and flows in the system.

3.24. SVS13B

SVS13B is a Class 0 protostar in the NGC1333 region (Grossman et al. 1987) (Sadavoy et al. 2014). It is located 14.9" south-west
of SVS13A (~4479 au, Tobin et al. 2018) and we can find both sources in the same field of view. The object presents a collimated outflow that could contribute to other flows around the HH711 (Bachiller et al. 1998). We present the moment-0 maps of the three main transitions for SVS13B in Fig. E.1t, and we can see a compact spatial distribution around the continuum peak for all of them. We only find emission of the two main species and the OC³⁴S isotopologue (Fig. E.2t). Neither OC³⁴S nor H₂³S were detected in the warm inner core.

## 4. Column densities and abundances

We calculated column densities for both OCS transitions, as well as for the H₂S transition. In the cases where H₂S was detected but OCS emission was not detected, we estimated an upper limit to the column density of OCS. We proceeded similarly with the 3 isotopologues. For the typical conditions of the warm inner cores (T>100 K, n~10⁶ cm⁻³) we assumed that all of the transitions are thermalised, and so, we did a local thermodynamic equilibrium (LTE) analysis to estimate the column densities.

Assuming that emission is optically thin and that the emission fills the beam, the column density of the upper level ($N_u$) of a given transition is (see, e.g., Goldsmith & Langer 1999)

$$N_u = \frac{8\pi k v^2 W}{hc^3 A_{ul}}, \qquad (1)$$

where k is Boltzmann's constant, v is the frequency of the molecular transition, W = $\int T_a$ dv is the integrated line intensity after the spectral fitting (with $T_a$ the brightness temperature and v the velocity), h is Planck's constant, c is the speed of light, and $A_{ul}$ is Einstein's coefficient for spontaneous emission. We note that, since $N_u$ depends on the integrated intensity, the spectral dilution in the lower spectral resolution lines has no effect in the result, as the area of a gaussian remains constant under dilution effects. Then, the total column density using the partition function, $Q(T_{ex})$, is

$$N = \frac{N_u Q(T_{ex}) e^{E_u / kT_{ex}}}{g_u}. \qquad (2)$$

The total column densities were calculated assuming a $T_{ex}$ of 100 K or $T_{ex} = T_{bol}$ if $T_{bol} > 100$ K in the warm inner core for all species and sources. We calculated two different N(OCS) values, one from the J=18–17 transition and another from the J=19–18 one (which should actually be the same). We decided the total column density of OCS to be the mean of both values, given that the two quantities are similar in all cases (within less than a factor of two). We estimated an additional value for the column densities of OCS and H₂S, when possible, using their isotopologues OC³⁴S and H₂³S (see Section 5). The column densities of H₂S and H₂³S have been calculated assuming the statistical value of the orto/para ratio of 3. Regarding uncertainties in the temperature, a factor of 2 in the temperature translates in a factor of ~1.5 in OCS, OC³⁴S and H₂S column densities, and slightly major changes in the H₂³S column density. The effect of the temperature is smaller in the H₂S/OCS ratio, as long as the same temperature is assumed for both species. In fact, the maximum change of the H₂S/OCS ratio in all the sample is of a factor ~1.5, usually within the uncertainties.

We then proceeded to the calculation of the abundances of each species, with respect to molecular hydrogen (H₂) gas-phase abundance. We calculated the column density of H₂ using the

data from our observations of the continuum, following the process in Kauffmann et al. (2008). We start from the expression of the intensity emitted by a source of temperature T, optical depth $\tau_v$ and frequency v, given by the equation of radiative transfer:

$$I_v = B_v(T)(1 - e^{\tau_v}), \qquad (3)$$

where $B_v$ is Planck's blackbody function, and the optical depth can be expressed as

$$\tau_v = \int \kappa_v \rho \, ds, \qquad (4)$$

with $\kappa_v$ as the absorption coefficient and $\rho$ as the gas density.

If most of the hydrogen is in its molecular form, H₂, the column density can be related with the optical depth using the relation

$$N_{H_2} = \int n_{H_2} \, ds = \int \frac{\rho}{\mu_{H_2} m_H} \, ds = \frac{\tau_v}{\mu_{H_2} m_H \kappa_v}, \qquad (5)$$

where $\mu_{H_2}$=2.8 is the molecular weight per hydrogen molecule and $m_H$ is the H-atom mass.

On the other hand, $B_v$ cannot be simplified using the Rayleigh-Jeans limit as the necessary $\lambda$ regime is not achieved in these conditions. For ~1.2 mm observations, in most cases (Kauffmann et al. 2008),

$$\lambda \gg 1.44 \, \text{mm} \left(\frac{T}{10 \, \text{K}}\right)^{-1}. \qquad (6)$$

Apart from this, we also look for an easy-to-use relation between flux per beam and intensity. In this case, we started with

$$F_{\text{beam},v} = \int I_v P \, d\Omega, \qquad (7)$$

where P is the normalised power pattern of the telescope and the solid angle of the antenna is defined as $\Omega_A = \int P d\Omega$. If we assume that the beam solid angle of the telescope is similar to a Gaussian function, we have that

$$P = P(\theta) = e^{-\theta^2 / 2\theta_0^2}, \qquad (8)$$

with $\theta_0 = \theta_{\text{HPBW}} / \sqrt{8 \ln 2}$. With this condition, we can then integrate the antenna's solid angle:

$$\Omega_A = \frac{\pi}{4 \ln 2} \theta_{\text{HPBW}}^2. \qquad (9)$$

Eventually, we can rewrite Eq. (5) using the optically thin approximation in Eq. (3) and using the average intensity from Eq. (7) as $\langle I_v \rangle = F_{\text{beam},v} / \Omega_A$ to obtain

$$N_{H_2} = \frac{\tau_v}{\mu_{H_2} m_H \kappa_v} = \frac{I_v}{\mu_{H_2} m_H \kappa_v B_v(T)} = \frac{F_{\text{beam},v}}{\mu_{H_2} m_H \kappa_v B_v(T) \Omega_A} \qquad (10)$$

Using useful units, we can recover Eq. (A.27) from Kauffmann et al. (2008)

$$N_{H_2} = 2.02 \times 10^{20} \text{cm}^{-2} \times (e^{1.439(\lambda/\text{mm})^{-1}(T/10\text{K})^{-1}} - 1) \times \left(\frac{\lambda}{\text{mm}}\right)^3 \times \left(\frac{\kappa_v}{0.01\text{cm}^2/\text{g}}\right)^{-1} \times \left(\frac{F_{\text{beam},v}}{\text{mJy/beam}}\right) \times \left(\frac{\theta_{\text{HPBW}}}{10 \text{ arcsec}}\right)^{-2}, \qquad (11)$$
<table><caption>**Table 2.** Column densities and H₂S/OCS ratios in Perseus.</caption><thead><tr><th>Protostar</th><th>Class</th><th>T<sub>bol</sub><br>(K)</th><th>T<sub>kin</sub><br>(K)</th><th>N(OCS)<br>(cm<sup>-2</sup>)</th><th>N(H<sub>2</sub>S)<br>(cm<sup>-2</sup>)</th><th>H<sub>2</sub>S / OCS</th></tr></thead><tbody><tr><td>B1bN</td><td>0</td><td>14.7</td><td>100.0</td><td>0.0</td><td>(1.54±1.31) ×10<sup>15</sup> <sup>(b)</sup></td><td>—</td></tr><tr><td>B1bS</td><td>0</td><td>17.7</td><td>100.0</td><td>(2.88±1.12) ×10<sup>15</sup> <sup>(a)</sup></td><td>(2.65±0.08) ×10<sup>16</sup> <sup>(b)</sup></td><td>9.21<sup>+4.33</sup><sub>-2.79</sub></td></tr><tr><td>B5-IRS1</td><td>I</td><td>287.0</td><td>287.0</td><td>&lt;7.15 ×10<sup>13</sup></td><td>(5.04±1.26) ×10<sup>14</sup></td><td>&gt;7.04</td></tr><tr><td>HH211MMS</td><td>0</td><td>27.0</td><td>100.0</td><td>(1.60±0.87) ×10<sup>15</sup> <sup>(a)</sup></td><td>(2.10±0.17) ×10<sup>14</sup></td><td>0.132<sup>+0.125</sup><sub>-0.053</sub></td></tr><tr><td>IC348MMS</td><td>0</td><td>30.0</td><td>100.0</td><td>(1.97±0.30) ×10<sup>15</sup> <sup>(a)</sup></td><td>(2.94±0.01) ×10<sup>16</sup> <sup>(b)</sup></td><td>14.9<sup>+1.9</sup><sub>-2.0</sub></td></tr><tr><td>IRAS2A</td><td>0/I</td><td>69.0</td><td>100.0</td><td>(1.85±1.16) ×10<sup>16</sup> <sup>(a)</sup></td><td>&gt;5.69 ×10<sup>16</sup> <sup>(b)(c)</sup></td><td>&gt;3.08</td></tr><tr><td>IRAS2B</td><td>I</td><td>106.0</td><td>106.0</td><td>(1.06±0.21) ×10<sup>14</sup></td><td>(1.03±0.15) ×10<sup>14</sup></td><td>0.973<sup>+0.333</sup><sub>-0.283</sub></td></tr><tr><td>IRAS4B</td><td>0</td><td>28.0</td><td>100.0</td><td>(1.44±0.39) ×10<sup>16</sup> <sup>(a)</sup></td><td>(2.06±0.01) ×10<sup>17</sup> <sup>(b)</sup></td><td>14.3<sup>+3.6</sup><sub>-3.1</sub></td></tr><tr><td>IRAS4C</td><td>0</td><td>31.0</td><td>100.0</td><td>(3.89±1.38) ×10<sup>13</sup></td><td>(1.15±0.13) ×10<sup>14</sup></td><td>2.97<sup>+1.53</sup><sub>-1.02</sub></td></tr><tr><td>L1448-IRS3A</td><td>0</td><td>47.0</td><td>100.0</td><td>(2.17±1.16) ×10<sup>15</sup> <sup>(a)</sup></td><td>(2.18±0.46) ×10<sup>14</sup></td><td>0.101<sup>+0.115</sup><sub>-0.049</sub></td></tr><tr><td>L1448-IRS3B</td><td>0</td><td>57.0</td><td>100.0</td><td>(1.61±0.76) ×10<sup>15</sup> <sup>(a)</sup></td><td>(2.71±0.12) ×10<sup>14</sup></td><td>0.168<sup>+0.112</sup><sub>-0.059</sub></td></tr><tr><td>L1448C</td><td>0</td><td>47.0</td><td>100.0</td><td>(3.82±3.69) ×10<sup>15</sup> <sup>(a)</sup></td><td>(1.31±0.03) ×10<sup>15</sup></td><td>0.34<sup>+6.44</sup><sub>-0.17</sub></td></tr><tr><td>L1448NW</td><td>0</td><td>22.0</td><td>100.0</td><td>&lt;3.43 ×10<sup>13</sup></td><td>(9.61±2.63) ×10<sup>13</sup></td><td>&gt;2.81</td></tr><tr><td>Per-emb-2</td><td>0</td><td>27.0</td><td>100.0</td><td>(2.79±1.11) ×10<sup>13</sup></td><td>(3.00±0.15) ×10<sup>14</sup></td><td>10.8<sup>+5.5</sup><sub>-3.4</sub></td></tr><tr><td>Per-emb-5</td><td>0</td><td>32.0</td><td>100.0</td><td>(1.91±0.58) ×10<sup>15</sup> <sup>(a)</sup></td><td>(1.48±0.11) ×10<sup>14</sup></td><td>0.078<sup>+0.030</sup><sub>-0.022</sub></td></tr><tr><td>Per-emb-8</td><td>0</td><td>43.0</td><td>100.0</td><td>&lt;2.92 ×10<sup>13</sup></td><td>(1.16±0.30) ×10<sup>14</sup></td><td>&gt;3.98</td></tr><tr><td>Per-emb-18</td><td>0</td><td>59.0</td><td>100.0</td><td>(3.65±0.55) ×10<sup>15</sup> <sup>(a)</sup></td><td>(2.77±0.13) ×10<sup>14</sup></td><td>0.076<sup>+0.013</sup><sub>-0.013</sub></td></tr><tr><td>Per-emb-22</td><td>0</td><td>43.0</td><td>100.0</td><td>(2.37±0.15) ×10<sup>14</sup></td><td>(3.79±0.13) ×10<sup>14</sup></td><td>1.60<sup>+0.13</sup><sub>-0.15</sub></td></tr><tr><td>Per-emb-29</td><td>0</td><td>48.0</td><td>100.0</td><td>(4.51±2.03) ×10<sup>15</sup> <sup>(a)</sup></td><td>(5.36±0.09) ×10<sup>16</sup> <sup>(b)</sup></td><td>11.9<sup>+6.8</sup><sub>-3.8</sub></td></tr><tr><td>Per-emb-30</td><td>0/I</td><td>78.0</td><td>100.0</td><td>0.0</td><td>0.0</td><td>—</td></tr><tr><td>Per-emb-50</td><td>I</td><td>128.0</td><td>128.0</td><td>0.0</td><td>0.0</td><td>—</td></tr><tr><td>Per-emb-62</td><td>I</td><td>378.0</td><td>378.0</td><td>0.0</td><td>0.0</td><td>—</td></tr><tr><td>SVS13A</td><td>0/I</td><td>188.0</td><td>188.0</td><td>(4.33±2.70) ×10<sup>16</sup> <sup>(a)</sup></td><td>&gt;1.48 ×10<sup>17</sup> <sup>(b)(c)</sup></td><td>&gt;3.42</td></tr><tr><td>SVS13B</td><td>0</td><td>20.0</td><td>100.0</td><td>(1.39±0.64) ×10<sup>15</sup> <sup>(a)</sup></td><td>(2.68±0.19) ×10<sup>14</sup></td><td>0.19<sup>+0.13</sup><sub>-0.07</sub></td></tr></tbody></table>

Notes. (a) Estimated using the OC³⁴S isotopologue column density. (b) Estimated using the H₂³S isotopologue column density. (c) The broad CH₃CHO line (see Appendix A) complicates the fit of the H₂³S lines, but a lower limit could be estimated.

which we used to calculate the H₂ column density from the flux of our continuum observations.

The magnitudes represented in Eq. 11 are the following: λ is the wavelength of the continuum observations, 1.29 mm; the temperature is the respective T<sub>dust</sub>, which we have estimated to be the T<sub>bol</sub> of each source; κ<sub>v</sub> is the dust opacity, for which we used the value of 0.01056 cm g⁻¹, following Artur de la Villarmois et al. (2023) and Ossenkopf & Henning (1994), for a gas-to-dust ratio of 100; θ<sub>HPBW</sub> is the width of the region where we integrate the flux, which in this case is a circular region of 1.5" of diameter; F<sub>beam,v</sub> is the integrated flux per beam of the dust continuum map inside the 1.5" region around the peak in millijanskys per beam.

## 5. H₂S/OCS ratio

We calculated the averaged column densities of H₂S, H₂³S, OCS, OC³S, and OC³S in a 1.5" diameter circle around the continuum peak using the averaged interferometric spectra and following the procedure described in Sect. 4. One important problem is the impact that the possible high opacities of the main isotopologue lines can have in our column density estimates, leading to a severe underestimation of the real ones. To avoid this problem, we used the observations of the less abundant isotopologues to estimate the column density of the most abundant one. We have detected the OC³S line in 13 of the total 17 targets detected in OCS. Significant deviations of the ³²S/³⁴S ratios from the solar value are not detected in the

interstellar medium and/or comets (Gratier et al. 2016; Calmonte et al. 2017). Therefore, we used OC³S as a proxy for OCS by calculating N(OCS)=22.5 ×N(OC³S). These are the values shown in Table 2 and used for the following discussion. We have also detected the OC³S line in nine protostars. The values of N(OC³S)/N(OC³S) in these targets are >2.5, proving that the OC³S line is not optically thick. In the case of H₂S, we have detected the H₂³S S<sub>2,0</sub>-S<sub>2,1</sub> in seven sources. In these sources, we used the H₂³S column density to compute that of H₂S assuming ³²S/³³S = 126 (Anders & Grevesse 1989). The column densities thus calculated are significantly higher than those estimated from the main isotopologue line. The most extreme case is B1bS, in which the column density estimated from H₂³S is ~144 times higher than that obtained from the main isotopologue observations. However, the difference is lower, ~15 to ~85, for the other protostars. We think that the extreme case of B1bS is caused by a very high opacity in the source, due to an extense and thick envelope surrounding the very young protostar, similarly to B1bN. This is supported by the very low ³⁴S/³³S ratio (see Table D.3) which suggests that even the OC³S isotopologue is optically thick in this source.

In Fig. 5, we plot N(OCS) versus N(H₂S) for the targets where H₂S has been detected. There is a large scatter between the two magnitudes, specially at low values of N(H₂S), reflecting significant changes in the N(H₂S)/N(OCS) ratio (hereafter, H₂S/OCS). The values of H₂S/OCS range from ~0.1 to ~10, thus expanding over two orders of magnitude (see Table 2 and Fig. 11). In fact, we can differentiate two groups. The first
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** OCS versus H₂S comparison in the sources where H₂S was detected. The color of the markers represents bolometric temperature (side colorbar). Left: column densities. The OCS column density increases with increasing H₂S column density. We see two different trends in this graph: (1) Several sources have H₂S/OCS ratios around ~7–10, which we classify as OCS-poor sources. (2) The rest of the sources present H₂S/OCS ratios in the ~0.1–1 range, and most of them are arranged in a cluster of several points in the {N(OCS)~2×10¹⁵ cm⁻², N(H₂S)~2×10¹⁴ cm⁻²} range. They are classified as OCS-rich sources. Right: abundances. We find similarities with the first graph: The OCS abundances grow together with H₂S abundances, the lower H₂S/OCS ratio sources follow a linear (log) trend, and the rest are very differentiated from the fitted line. In this case, the fit corresponds to H₂S/OCS ratio of 0.147±0.022.</figcaption></figure>

one, characterized by H₂S/OCS ~7–10, is formed by the seven sources detected in H₂³S and one additional source with low H₂S column density, and we refer to them as “OCS-poor” targets. The second one is composed of 8 sources with N(H₂S) < 2×10¹⁵ cm⁻² and presents H₂S/OCS ~0.1–1. We refer to these targets as “OCS-rich” protostars. In four sources, we have not detected OCS and only upper limits to the H₂S/OCS ratio can be derived. These limits imply H₂S/OCS >5 and consequently, they can be considered as OCS-poor protostars. One could think that this difference can be related to a different evolutionary stage but, as shown in Fig. 5, there is no correlations between H₂S/OCS and the bolometric temperature. We consider that this differentiated chemistry is more likely related with a different composition of the sublimated ice. In Sect. 7, we discuss the origin of this chemical diversity.

In the following, we discuss some observational uncertainties that can influence the derived values of H₂S/OCS and therefore, our classification in OCS-poor and OCS-rich protostars. One could think that the limited angular resolution of our observations can bias the estimated values. The unknown size of the emitting area would definitely have a significant impact on the estimated column densities but we do not expect a big impact on the calculated H₂S/OCS since the transitions observed present similar excitations conditions (E<sub>u</sub> ~ 80–100 K) and more likely come from the same region. A second problem would be the possible impact of the opacity in our column density estimates. As explained above, we avoided this problem through the observation of the less abundant isotopologues. In particular, we observed transitions of OCS, OC³⁴S, and OC³³S, which allowed us to derive accurate OCS column density estimates. In the case of H₂S, there is a significant fraction of sources in which we have not detected the H₂³S S₂,₀–S₂,₁ line and our estimates rely on the observations of the main isotopologue. Looking at the ratios between the H₂S column densities derived from the ³³S isotopologue over those using ³²S (see Table D.3), we expect that the opacity of the main isotopologue line toward the targets without H₂³S detection is less than ~15. This would mean that

our estimates are reliable within an order of magnitude in these sources. This is also corroborated by the lower column densities and opacities of the OCS lines, which suggest that these protostars host less massive warm regions. Even considering the worst case scenario where all OCS-rich protostars with OC³⁴S detection would scale up their H₂S by a factor 37 – which is the geometrical mean of the scalling factors for the sources detected in H₂³S (without B1bS; we note that the opacity scales exponentially, so the appropriate mean is the geometrical and not the arithmetic one) –, we would still find the two-family differentiation between OCS-rich and OCS-poor sources. In this scenario, which we show in Figure C.1, L1448C would be the only case where we have classified a source as OCS-rich that would prove to be an OCS-poor one, but the high uncertainty in the H₂S/OCS ratio makes this comparison a rather bad example.

## 6. Completing the sulfur budget

To determine sulfur depletion is an important open question in astrochemistry. In the dense and cold regions of molecular clouds, most of the sulfur is thought to be locked in the icy mantles that cover the dust grain surfaces. Chemical models (Navarro-Almaida et al. 2020) and comet observations (Calmonte et al. 2016) suggest that H₂S is the most abundant sulfur compound in the ice. However, its detection in interstellar ices remains elusive even with the high sensitivity provided by the JWST. One factor contributing to this is its relatively low binding energy as calculated by Bariosco et al. (2024). Currently, OCS is the only compound firmly detected in interstellar ices (Palumbo et al. 1995; Boogert et al. 2022). Different authors have published tentative detections of the 7.5 µm band of SO₂ (Rocha et al. 2024), but this band is overlapped with intense bands of abundant complex organic molecules (COMs) and CH₄, which hinder its confirmation (Taillard et al., in prep.). Ices are sublimated in the warm innermost regions of protostars where the dust temperature is greater than 100 K. The chemical composition of this warm region can therefore inform on the ice
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Abundances of the main sulfuretted species (H₂S, OCS, SO, SO₂) in each of the protostars in our sample. The purple marks represent the total sulfur abundance produced by these four sulfuretted molecules. Sources are sorted first by their Class and then by T<sub>bol</sub>. Total abundances seem to present a slight growing trend with T<sub>bol</sub> and Class. We find that with the data available, SO and SO₂ are tipically one to two orders of magnitude below H₂S and OCS, but optically thick lines (specially the H₂S lines not detected in H<sup>33</sup>S) may play a role in refuting or supporting these findings. We note that H₂S and OCS abundances are very similar in IRAS2B and hardly distinguishable in the figure.</figcaption></figure>

composition as long as the gas-phase chemistry has not signifi-
cantly altered it. Even if it were the case, these warm regions,
where all the volatiles are in gas phase, provide an excellent
opportunity to determine the total sulfur budget in volatiles, i.e.,
sulfur depletion. In Sect. 5, we have calculated the abundances
of H₂S and OCS, that are expected to be the most abundant
sulfur-bearing species in ices. Here, we use the SO and SO₂
column densities reported from the data of the PEACHES sur-
vey (Yang et al. 2021) (ALMA project codes: 2016.1.01501.S
and 2017.1.01462; PI: N. Sakai) by Artur de la Villarmois et al.
(2023), obtained to complete the sulfur budget in most of our tar-
gets. Unfortunately, observations of Per-emb-30 and Per-emb-62
in SO and SO₂ are not included in the sample studied by Artur
de la Villarmois et al. (2023).

The synthesized beam of PEACHES observations, $\theta_{\text{HPBW}} \sim 0.5''$, is three times smaller than the region size ($\theta_{\text{HPBW}} = 1.5''$) used for our calculations, and it is also beyond the maximum resolution of our continuum images ($\theta_{\text{HPBW}} \sim 1.3''$). Therefore, the comparison of the column densities derived from PEACHES observations with those listed in Table 2 is not straightforward, since it depends on the unresolved chemical and physical structure toward each target. To account for that, we have considered 2 limiting cases to estimate the SO and SO₂ abundances in our sample. In the first case, we assume that the emission of SO and SO₂ is extended and the continuum emission uniformly fills the 1.5" integration region. In this case, we calculate the SO and SO₂ abundances by adopting the column densities reported by Artur de la Villarmois et al. (2023) and the values of N(H₂) calculated as described in Sect. 4. In the second case, we assume that the continuum emission is unresolved, and that all of it comes from the innermost $\sim 0.5''$ region, similar to the beam used to calculate the SO and SO₂ column densities; then, we calculate the H₂ column density using the peak flux and using equation (11) with $\theta=0.5''$. By considering these two scenarios, we define an upper and lower limit to the SO and SO₂ abundances. Every other possibility should fall within this range of values. This uncertainty in the spatial distribution of the different emission lines and dust emission would introduce a large uncertainty, up to a factor of

~9, in the derived abundances. However, as commented below,
it does not have a big impact on the total sulfur budget.

In Fig. 6, we present the abundances of the main sulfuretted
species – H₂S, OCS, SO, SO₂ – in the 24 protostars of our sam-
ple, as well as the total sulfur abundance calculated as the sum
of the abundances of these species. There are large differences
between sources in the abundances of all S-bearing molecules:
SO and SO₂ abundances spread over almost three orders of mag-
nitude, while OCS and H₂S cover ranges of more than four orders
of magnitude. H₂S appears to be the most abundant species in
most of the sources, with a maximum >1.3×10⁻⁶ in SVS13A,
which is in accordance with the value measured in Orion KL
by Crockett et al. (2014) of (3.1±1.1)×10⁻⁶. The second most
abundant sulfuretted molecule is OCS, which in some cases is
greater than the H₂S abundance (OCS-rich sources). However,
these OCS-rich sources have less H₂S in the gas phase than the
OCS-poor ones, and the sum of H₂S and OCS abundances is
less than ~10⁻⁷ (see Fig. 6). When it comes to SO and SO₂, we
have given a range of values for their abundances, taking into
account the range of column densities calculated as explained
above. In all the sources, the abundances of these two species
are far less abundant than H₂S and OCS, <10⁻⁸. If the gas phase
abundances we are measuring reflects the ice composition in the
protostellar phase, the main sulfur reservoir in the Perseus Cloud
would be H₂S ice in ~ 60%, of the sources and OCS in the rest.
In addition, we find that the Class I sources gas-phase sulfur total
abundance declines. The SO and SO₂ remain with higher abun-
dances in the Class I, starting to gain more weight in the total
sulfur abundance.

According with the data shown in Fig. 6, the amount of sulfur
in gas phase increases during the Class 0 stage with the max-
imum in the Class 0/I sources (IRAS2A, SVS13A). One could
think that this is the consequence of missing some important
sulfur species, mainly CS, CCS, C₃S, and H₂CS. A more likely
explanation is that the warm region, i.e., the region in which ices
mantles are evaporated, is smaller in these young objects. Higher
angular resolution observations and the observation of a wide
variety of sulfur species are needed to probe these hot environ-
ments. Also interesting, the amount of sulfur seems to decrease
<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Evolution of the models with different warming curves. Left: initial ice and gas-phase components of the H₂S, OCS, SO, and SO₂ species in the inner core before the collapse. Comparison with the final abundance of these species in the gas-phase after the fiducial MHD collapse. Center: comparison of the fiducial MHD warming curve with the rest of warming curves explored in the section. Right: comparison of the fiducial MHD density curve with the rest of warming curves explored in the section. The different density curves are obtained using the barotropic equation of state from Machida et al. (2006).</figcaption></figure>

again in the Class I stage. It is true, however, that our sample of Class I objects is too small to draw firm conclusions.

1 $M_{\odot}$ of gas reaching up to the first hydrostatic core (FHSC) phase ($t \sim 10^4$ yr).

# 7. Chemical modeling

In this section, we describe a chemical model for evolved Class 0 and Class I warm inner cores aimed at understanding the chemical variety found in our sample and, in particular, the factors that could decide whether a source is H₂S or OCS-rich. Class 0 and Class I stages are the bridge between the pre-stellar phase and protoplanetary disks. While the degree of chemical processing along these stages is still unclear, chemical surveys (see, e.g., Bianchi et al. 2019; Mercimek et al. 2022), cometary abundances (Le Roy et al. 2015; Drozdovskaya et al. 2018; Altwegg et al. 2022), and chemo-dynamical models (Navarro-Almaida et al. 2024) suggest that there is a non-negligible fraction of the chemical composition that survives throughout this evolutionary process. Following this nature versus nurture idea, in this section we investigate whether the chemical diversity of our sample is a product of properties such as the observed bolometric temperature, whether it is inherited from the pre-stellar core from which they were born, or both.

Our chemical model for a Class 0 warm core consists of two steps. First, we set a pre-stellar phase model based on the Barnard-1b model presented in Navarro-Almaida et al. (2020) by taking its physical structure. With it, we then ran the 3-phase gas-grain chemical model NAUTILUS (Ruaud et al. 2016) in its latest version (Wakelam et al. 2024) for $10^6$ yr. We note that, in this chemical model, the gas-grain exchange processes considered are: thermal desorption, UV photodesorption, cosmic-ray induced photodesorption, cosmic-ray induced heating, chemical (reactive) desorption, cosmic-ray sputtering and diffusive chemistry. We then extracted the chemical abundances at the innermost point in the model as the representative chemical composition of the pre-stellar core and as the initial chemical abundances for the protostellar collapse. In a second step, the gravitational collapse of the pre-stellar core into the young Class 0 object was described by the density and temperature evolution of Lagrangian tracer particles present in the MHD simulation reported in Hennebelle et al. (2016) and Gerin et al. (2017). This simulation has also been used to investigate the impact of grain-growth on the chemical composition of the hot corino in Navarro-Almaida et al. (2024) and describes the collapse of

## 7.1. The role of collapse in the H₂S/OCS ratio

After obtaining the pre-stellar chemical composition of both gas and ice phases, we ran a core collapse model taking the density and temperature evolution of a tracer particle that belong to the warm core ($T_{gas} \sim 500$K) at the end of the MHD simulation. As discussed in Navarro-Almaida et al. (2024), the resulting chemical abundances do not change significantly among tracer particles that belong to the warm core by the end of the simulation and so we can use a random representative one to investigate the effect of the dynamical history on the H₂S/OCS ratio.

The gas and ice-phase abundances of H₂S, OCS, SO, and SO₂ in the pre-stellar phase (first step) are shown in the left panel of Fig. 7. As expected from cold core chemistry, most of molecular abundances are found in ices, with H₂S abundance being more than two orders of magnitude higher than that of OCS. SO₂ ice is also more abundant than that of SO. We then ran the MHD model of core collapse using these abundances as initial input. The evolution of density and temperature in this collapse is shown in the central and right panels of Fig. 7 as the fiducial MHD model. After the collapse, the resulting gas-phase abundances are shown in the left panel of Fig. 7. Comparing the ice-phase abundances of the pre-stellar phase and the gas-phase abundances at the end of the collapse, it is apparent that the H₂S in the warm core is mainly coming from the thermal desorption of the ice, with no significant chemical processing. This is not the case for OCS, with an enhancement of two orders of magnitude with respect to the total ice abundance. OCS ice abundance is increased when temperatures reach ~20 K, before being too high and thus thermally desorbing it. This is due to the enhanced diffusion of the reactants to form icy OCS. When temperatures reach several hundreds K, OCS ice is thermally desorbed and it is mainly found in the gas-phase, with an abundance that matches that of the icy OCS right before its desorption. The degree of chemical processing of OCS is therefore higher, in line with the results of Navarro-Almaida et al. (2024). Chemical processing is also found in SO and SO₂, molecules whose abundances are highly dependent on hot gas chemistry (see, e.g., Wakelam et al. 2011; Esplugues et al. 2013). In this setting, the H₂S/OCS ratio yields H₂S/OCS ~ 5.32. Compared to the H₂S/OCS ratio
<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Comparison of the H₂S and OCS chemical evolution during the final stages of the collapse with different warming curves. Left: H₂S ice and gas-phase abundances during the collapse with the different warming curves. The starting time of the figure has been set to 4×10⁴ years to show a closer detail of the major changes. The products after 5.5×10⁴ are very similar in all the models but, in the later stages, the model with a lower final temperature, “Lower T”, shows a drop of the gas-phase H₂S abundance. Center: OCS ice and gas-phase abundances during the collapse with the different warming curves. The starting time of the figure has been set to 4×10⁴, as in the left panel. The products after 5.5×10⁴ are very similar in all the models. Right: H₂S/OCS ratio after the collapse with the different warming curves. The models show a ratio close to five except for the model with a lower final temperature.</figcaption></figure>

measured toward the members of our sample, this ratio is high, only compatible with the OCS-poor sources, and close to what it is observed toward IRAS4C.

Column densities of OCS and SO₂ have been derived toward a sample of 26 massive hot cores by Santos et al. (2024). Moreover, they compared their abundances with that of methanol and concluded that the OCS/CH₃OH ratio is quite uniform in the sample, and similar to the values measured in the ices toward low-mass and massive star forming regions. They interpreted this similarity as an evidence that the abundances of these two species are inherited from the pre-stellar phase, which is in disagreement with our model results (see Navarro-Almaida et al. 2020). As they commented, this conclusion can be biased by the fact that OCS ice has only been detected toward two low-mass starless regions. Further observations are needed to understand the variety of ice composition in pre-stellar cores and test our model predictions.

Aiming at investigating why the H₂S/OCS ratio varies among the members of our sample and the physical or chemical factors behind these changes, we ran the chemical code Nautilus in several alternative collapse scenarios. These are shown in the middle and right panels of Fig. 7. Modifications were introduced in the final temperature of the warm core and the warming rate to make it slower. The corresponding changes in density were computed using the piece-wise barotropic equation of state presented in Machida et al. (2006):

$$T = T_0 \sqrt{1 + \left(\frac{n}{n_1}\right)^{2g_1} \left(1 + \left(\frac{n}{n_2}\right)^{g_2} \left(1 + \left(\frac{n}{n_3}\right)^{g_3}\right)\right)}, \quad (12)$$

with n as the total density, T₀=10 K, and

$$n_1 = 10^{11} \text{ cm}^{-3}, \quad n_2 = 10^{16} \text{ cm}^{-3}, \quad n_3 = 10^{21} \text{ cm}^{-3}, \quad (13)$$

$$g_1 = 0.4, \quad g_2 = -0.3, \quad g_3 = 0.56667. \quad (14)$$

The evolution of gas and ice-phase abundances of H₂S and OCS for the different collapse scenarios is shown in the left and central panels of Fig. 8, and the H₂S/OCS ratio is in the right panel of Fig. 8. While higher temperatures or a slower warming rate of the warm core do not produce significant changes in the H₂S/OCS abundance ratio, it decreases when the final temperature is lower than the ~500 K of the fiducial model.

Therefore, gas temperature seems to be an important factor governing the H₂S/OCS abundance ratio. The gas-phase abundance of H₂S drops from the value set in the pre-stellar phase only in the colder scenario, while OCS is mainly unaffected (see Figs. 8 and 9). We recall that the angular resolution of our observations does not allow us to resolve the protoplanetary disk. One would expect to have regions with different temperature within the beam, and the exact distribution of these regions could have a strong impact on the average value of the H₂S/OCS ratio.

To investigate why the gas-phase H₂S abundance presents different behaviors, we show the evolution of the main sulfur carriers in Fig. 9. We noted the drastically different evolution of molecular sulfur S₂ between the fiducial MHD model and the model with a lower final temperature. It is the main culprit in the decline of gas-phase H₂S abundance, as it is the only sulfur bearing molecule whose abundance is highly enough enhanced by the time H₂S abundance starts to drop. Comparing the fiducial model with the low temperature model (Fig. 9), the abundance of S₂ in the latter is up to four orders of magnitude higher than in the former scenario. On the one hand, according to the chemical network in the lower temperature case, H₂S is mainly destroyed by its neutral-neutral reaction with H such that H + H₂S → H₂ + HS, with a reaction rate of k = 3 × 10⁻⁴ cm³ s⁻¹. H₂S formation is less efficient, produced on the surface of grains by the radical-radical reaction between sulfanyleths HS + HS → S + H₂S, with a reaction rate of k = 1 × 10⁻⁴ cm³ s⁻¹. On the other hand, S₂ is being formed by the reaction of HS with atomic sulfur, such that S + HS → H + S₂ (k = 1 × 10⁻⁴ cm³ s⁻¹). Therefore, the sequestration of HS in the formation of S₂ at a similar rate as the formation of H₂S and the more efficient destruction of H₂S are the factors responsible for the decline of gas-phase abundance of H₂S and the enhancement of S₂ abundance. In the fiducial MHD scenario, H₂S production and destruction are again mediated by HS in the neutral-radical reaction of HS with molecular hydrogen H₂ + HS ⇌ H₂S, with its creation (forward) being more efficient than its destruction (backward). Furthermore, in this scenario, S₂ is no longer formed by reactions with HS, but in the neutral-ion reaction NO + S₂⁺ → S₂ + NO⁺ with the rate k = 3.73 × 10⁻⁸ cm³ s⁻¹. This leads to the less efficient formation of S₂ and the stable gas-phase abundance of H₂S we see in Fig. 9.
<figure><img src="image_10.png" /><figcaption>**Fig. 9.** Species abundance comparison between the fiducial MHD model and the lower temperature model. The starting time of the figure has been set to $5 \times 10^4$ years to show a closer detail of the major changes. Left: fiducial MHD model. In the fiducial case, we observed that the chemistry is mostly inherited from the ices. Once the molecules sublimate (at $\sim 5.4 \times 10^4$ years) their abundances remain almost constant, except for $CS_2$ and $S_2$. Right: lower temperature model. In this model, different molecules sublimate at different temperatures, and some need a few thousand years to stabilize their gas-phase abundance (e.g., $CS_2$, $SO$, $SO_2$ and $H_2CS$). After $\sim 5.8 \times 10^4$ yr, $H_2S$ abundance starts to fall, while $H_2S_2$ and $S_2$ start to form. In particular, $S_2$ increases its abundance abruptly by several orders of magnitude, reaching an abundance of $\sim 10^{-6}$, which is $10^4$ times more than in the fiducial case.</figcaption></figure>

Polysulfides $H_2S_x$ have recently been proposed as compounds that could be involved in the synthesis of sulfur chains. Large sulfur chains are refractory, remaining locked on grain surfaces even at high temperatures. These compounds are candidates for sulfur reservoirs and might be where the missing sulfur in molecular clouds is (Cazaux et al. 2022; Carrascosa et al. 2024). We noted the high abundance of the polysulfides $H_2S_2$ and $H_2S_3$ in the different models of collapse. Of particular importance is the lower temperature scenario, in which $H_2S_2$ abundance is enhanced one order of magnitude with respect to the fiducial case. In both cases, however, the $H_2S_3$ abundance is similar.

## 7.2. $H_2S/OCS$ ratio and the initial conditions

After analyzing the effect of different warming rates and final temperatures in the final chemical composition surrounding the protostar, we explored the role that initial conditions play in the chemical abundances of the molecules studied thus far. This is motivated by how $H_2S$ is inherited from the pre-stellar phase since the main contribution to the total $H_2S$ in gas phase after the collapse comes from the sublimation of $H_2S$ ices, as shown in the left panel of Fig. 7.

To do so, we ran a set of zero-dimensional (0-D) pre-stellar models with Nautilus for $10^6$ years. These models are warmer alternatives to the pre-phase considered in Sect. 7.1 (see Table 3), that is, the innermost point of the Barnard 1b model presented in Navarro-Almaida et al. (2020), but covering a range of dust and gas temperatures between 10 K and 25 K, spaced by 5 K. These models would be appropriate to describe, for instance, intermediate-mass star forming regions illuminated by a stronger interstellar FUV field or clusters subjected to intense stellar feedback. After $10^6$ years, the chemical abundances become the input of the corresponding core collapse. We modified the density and temperature evolution of the collapsing core in such a way that they take into account the rising temperature of the pre-phase: if

<table><caption>**Table 3.** Comparison of physical properties between models.</caption><thead><tr><th>Model</th><th>n<sub>H</sub><br>(cm<sup>-3</sup>)</th><th>T<sub>0</sub><br>(K)</th><th>A<sub>v</sub><br>(mag)</th><th>&chi;<sub>UV</sub></th><th>&zeta;<sub>H<sub>2</sub></sub><br>(s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>Model 7.1</td><td>2.0&times;10<sup>4</sup></td><td>10</td><td>15</td><td>25</td><td>4.0&times;10<sup>-17</sup></td></tr><tr><td>Model 7.2</td><td>2.0&times;10<sup>4</sup></td><td>15–25</td><td>15</td><td>100</td><td>5.0&times;10<sup>-17</sup></td></tr></tbody></table>

**Notes.** The five magnitudes compared in the table are, in order from left to right: Initial Gas Density, Initial Gas Temperature, Initial Visual Extinction, UV Flux and Cosmic-rays Ionization Rate.

at any given time step the temperature is lower than the temperature set at the pre-phase, it gets updated to the pre-phase value. The density is then updated accordingly following the barotropic equation of state. The results of the models with $T = 15$ K and $T = 25$ K are shown in Fig. 10. On one hand, we found that, even at 15 K, a temperature slightly higher than in the fiducial MHD model, icy $H_2S$ at the end of the pre-phase is less abundant than in the fiducial case. The final gas-phase abundance of $H_2S$ is not directly inherited as the thermal desorption of the ice content, but it is increased in the warm up of the collapsing core, reaching a comparable value with respect to the fiducial model. On the other hand, OCS ice-abundance is higher than its predicted value by the fiducial model and comparable to that of $H_2S$. This also results in a higher gas-phase abundance of this molecule at the end of the collapse, yielding a $H_2S/OCS$ ratio of $\sim 2.4$. This is no longer the case as we increase the temperature of the pre-phase. In the right panel of Fig. 10 we show the results of the warmest model among the set of 0-D models. The ice-phase abundance of $H_2S$ drops by one order of magnitude with respect to the 15 K case, while icy OCS abundance declines significantly, by three orders of magnitude. The corresponding gas-phase abundances at the end of the collapse are lower compared to the fiducial or the 15 K case, resulting in a $H_2S/OCS$ ratio of $\sim 53$.

The general decreasing $H_2S$ abundances in warmer scenarios is accompanied by the abundant formation of $H_2S_3$ ices during
<figure><img src="image_11.png" /><figcaption>Fig. 10. Comparison between the 15 K and 25 K models. In the panels, we show the initial ice and gas-phase components of the H₂S, OCS, SO, SO₂, and H₂S₃ species in the inner core before the collapse, and we compare them with the final abundance of these species in the gas-phase after the collapse. There is a significant fall of H₂S, OCS, SO and SO₂ in the ices after the 25 K pre-phase, compared with the 15 K pre-stellar phase. On the contrary, H₂S₃ is formed in the ices during the warmer pre-phase, something that did not happen in the cooler one. After the collapse of the 25 K model, the final H₂S, SO and SO₂ drops in ~1 order of magnitude with respect to the 15 K simulation. This is not true for OCS, whose abundance falls more than 2 orders of magnitude. H₂S₃ final gas-phase abundance is greater in the warmer model, reaching an abundance of ~5×10⁻⁶.</figcaption></figure>

the pre-stellar phase (right panel of Fig. 10) that were not present at lower temperatures (left panel of Fig. 10). These high abundances then translate into higher abundances in the gas-phase after the MHD collapse. H₂S₃ abundances can be even higher than those of H₂S (right panel of Fig. 10). With a final gas-phase abundance of 4.5 × 10⁻⁶ and given that this molecule contains three sulfur atoms, it accounts for a great fraction of the cosmic sulfur abundance. According to the chemical network, icy H₂S₃ is formed by the diffusion and reaction of HS with H₂S₂ in the ice matrix. This diffusion process is only efficient when dust temperature is higher than ~22 K and, consequently, only appears in the ice of the warmest models considered here. The warm-up during the collapse present in each model also triggers this diffusion of the reactants to form H₂S₃. This is why the gas-phase H₂S₃ abundance is high at the end of the collapse in colder models even when its was not present in the pre-phase. Nevertheless, this result should be considered with caution since the production and destruction of sulfur allotropes are not yet well understood and the family of sulfides H₂Sₓ is included in the chemical network only for x ≤ 3. However, Shingledecker et al. (2020) included the formation of allotropes in their code, finding that S allotropes can become a major sink of sulfur, which agrees with our results for the warmer models.

In conclusion, our simulations show that both the absolute abundances of H₂S and OCS as well as the H₂S/OCS are very sensitive to the gas and dust temperature in the pre-stellar phase and also the final temperature in the inner core. High temperatures in the pre-stellar phase drive to higher H₂S/OCS ratio but lower absolute abundances of these molecules, since S is sequestered in hydrogen sulfides. The final temperature of the warm inner core has also a strong impact on the chemical composition of the warm gas. Low values of the H₂S/OCS ratio are consistent with regions with T ~ 100 K. Summarizing, the H₂S and OCS are very sensitive to the thermal history of the gas and dust, producing a high range of these molecules abundances and changing the H₂S/OCS ratio.

## 8. Discussion: The missing sulfur

This section is dedicated to the discussion of the gas-phase depletion of sulfur-bearing species in the Class 0/I protostars

of Perseus Molecular Cloud and how it compares to other environments. Sulfur depletion is defined as the ratio of the cosmic sulfur abundance over the amount of sulfur in volatiles:

$$D_S = \frac{[S/H]_{\text{cosmic}}}{[S/H]_{\text{volatiles}}}. \qquad (15)$$

Since some of the most abundant sulfur bearing species cannot be easily observed, the value of S/H<sub>volatiles</sub> needs to be indirectly estimated through chemical modeling (see, e.g., Laas & Caselli 2019; Navarro-Almaida et al. 2020; Fuente et al. 2023). Warm cores in the nuclei of young protostars are specially favorable environments to estimate D<sub>S</sub> because we can assume that the ices have sublimated and all the volatile sulfur is in gas phase. Thus, we can estimate the amount of sulfur in volatiles by adding the abundances of the most abundant sulfur-bearing molecules in this kind of environment. We have used the H₂S and OCS abundances derived with PRODIGE data and the SO and SO₂ abundances estimated from PEACHES observations, to have an estimate of sulfur depletion in our sample (see Table D.4). These species are expected to be the most abundant sulfur compounds in warm inner cores and lock essentially all the sulfur atoms (see, e.g., Vidal & Wakelam 2018), their sum providing a good estimate of D<sub>S</sub>. We note, however, that this is true only if a negligible fraction of S is in allotropes (Shingledecker et al. 2020). Thus calculated, we found that the highest sulfur abundance in our sample is measured toward SVS13A, of around ~1.7×10⁻⁶, which corresponds to a sulfur depletion of D<sub>S</sub> ~8. This value is comparable to the depletion observed in sources in Orion A (Fuente et al. 2023, 2024). However, the rest of the sources present higher depletion of sulfuretted species, with most of the values between D<sub>S</sub> ~500–5000, and the least abundant ones reaching values up to D<sub>S</sub> ~5×10⁵ (e.g., Per-emb-02, L1448NW). These values are higher than those found by Fuente et al. (2023) in Orion, Taurus and Perseus molecular clouds, suggesting that sulfur depletion increases during the collapse of a starless core. Yet, there are some observational and theoretical factors that can contribute to these large values of D<sub>S</sub>.

The first one implies that the sulfur-bearing species may be still in the ices in some regions within our beam. We remind that the beam of our observations, ~1.3" (~390 au) may not be enough to resolve the smaller disks (Tobin & Sheehan 2024).
Even in the case of Class I larger protoplanetary disks, it may
not be enough to resolve the warmest areas around the proto-
star. The majority of the sources in our sample are very young
(Class 0), and the temperature in the warm inner core might
not be as high as in the more evolved ones. This would have a
strong effect on the position of the snow line and would change
the amount of ices – and, with them, the amount of sulfuretted
molecules – sublimated into the gas phase, which is the only
component we observed in this work. This hypothesis is sup-
ported by the results shown in Fig. 6, where the protostars’ total
abundance of sulfur is growing with evolution (given by T<sub>bol</sub>)
until they reach Class 0/I. A lower bolometric temperature in the
core entails a closer snow line, as we explained previously, which
would produce a concentration of gas-phase species in the most
inner part that would remain unresolved by our observations.
Possible examples could be L1448-IRS3A or L1448NW, among
others. As the protostar evolves, the snow lines would move far-
ther from the protostar. In the opposite side of the evolutionary
track, in a evolved disk (Class I), where the protoplanetary disk
is thinner than in the Class 0 sources and the hot material is
concentrated in the innermost region (R < 50 au) of the disk
and its surfaces, while the outer disk is mostly cold. Again, our
observations will not resolve the inner region and we would be
essentially detecteting the colder outer disk. This phenomena
could give an explanation to the later fall of sulfur abundance
in the Class I objects (e.g., B5-IRS1) and the lack of detec-
tions of the sulfuretted species (e.g., Per-emb-50, Per-emb-62).
Our limited spatial resolution would affect the absolute value of
the H<sub>2</sub>S and OCS estimated abundances but not the H<sub>2</sub>S/OCS
ratio since all the observed transitions are expected to come from
the warm region. The large scattering in the values of H<sub>2</sub>S/OCS
demonstrates important chemical differences within our sample.

It has been suggested that sulfur depletion depends on the
environment (Fuente et al. 2016; Navarro-Almaida et al. 2020;
Fuente et al. 2023) in molecular clouds, and therefore, the
estimated values of D<sub>S</sub> could depend on the location of the pro-
tostars of our sample. In a more crowded region, protostars are
exposed to greater fluxes of UV rays, which can lead to higher
levels of photo-desorption, and shocks associated with bipolar
outflows, which can release sulfur-bearing species to the gas-
phase. Although all the protostars in our sample are located in
Perseus, there is a patent diversity of conditions in which these
objects are present. The most clustered zone of Perseus is SVS13,
where, apart from the cluster of protostars, there is a Herbig
Haro object (HH 711) in its vicinity (Hatchell & Fuller 2008;
Hatchell & Dunham 2009). IRAS2A, IRAS4B, Per-emb-18, and
SVS13A are located in this active star forming region (Hatchell
et al. 2013; Codella et al. 2021). Another interesting region is
IC 348, where the UV flux is expected to be higher because its
proximity to a cluster of massive stars and a second Herbig Haro
object is located, HH 211. In these regions, we find values of D<sub>S</sub>
~10 to 300, with variations of more than one order of magni-
tude. Although environment is surely playing a role, we do not
find a clear correlation between the protostar's location and D<sub>S</sub>
because other factors, such as the evolutionary stage, have also
an important impact.

The question here is whether the total sulfur abundance – in
gas phase and in the ices – is the same, or similar, in all these
sources in spite of the large dispersion of D<sub>S</sub> values. Chemi-
cal models predict the presence of large sulfur compounds that
possess sublimation temperatures higher than that of water and
thus are not expected to be sublimated in warm inner enve-
lope of young protostars, where the dust temperature is, at
most, ~500 K. In particular, both laboratory experiments and

theoretical work show that sulfur allotropes, such as S<sub>8</sub> and
large hydrogen sulfides (H<sub>2</sub>S<sub>x</sub>), could be important sulfur reser-
voirs (Jiménez-Escobar et al. 2012; Shingledecker et al. 2020;
Cazaux et al. 2022; Fuente et al. 2023; Carrascosa et al. 2024).
The sublimation temperature of these compounds when x > 3
is ≫100 K and can be considered as semi-refractory material
(Perrero et al. 2024; Carrascosa et al. 2024). Our simulations
show a very efficient formation of H<sub>2</sub>S<sub>3</sub> under some physical
conditions, locking most of the sulfur atoms (see Section 7)
while the abundances of H<sub>2</sub>S and OCS remain lower than 10<sup>-8</sup>
The formation and destruction processes of allotropes are not
well understood yet but our predictions demonstrate the forma-
tion of these compounds in molecular clouds. Several laboratory
experiments confirm that large hydrogen sulfides are formed dur-
ing the irradiation of H<sub>2</sub>S and H<sub>2</sub>S:H<sub>2</sub>O ices (Cazaux et al. 2022;
Carrascosa et al. 2024). More recently, experiments reported by
Martín-Doménech et al. (2024) suggest that sulfur chains (S<sub>x</sub>)
can be formed in CO:CS<sub>2</sub> and CO<sub>2</sub>:CS<sub>2</sub>. Morevoer, sulfur chains
have been detected in comets (Calmonte et al. 2016) and the
compound S<sub>2</sub>H was detected in the Horsehead nebula by Fuente
et al. (2017b). The existence of these molecules in the interstellar
medium is proved and the relative importance of these com-
pounds as sulfur reservoirs in different environments deserves
further study.

Other compounds such as NH<sub>4</sub>SH have also been suggested as important sulfur reservoirs in dense regions (Vitorino et al. 2024). This compound sublimates together with water, and is released to the gas phase increasing the abundances of H<sub>2</sub>S and NH<sub>3</sub>. Therefore, the formation of this molecule cannot be responsible for the low abundance of H<sub>2</sub>S in our sample of warm cores. Finally, some molecules not considered in this work such as CS, CCS, and H<sub>2</sub>CS could lock a significant fraction of sulfur atoms. Recent observations of these compounds toward SVS 13A unveiled that CCS, and H<sub>2</sub>CS column densities are one order of magnitude lower than those of the molecules in this work. Regarding CS, the column density is very uncertain and could be of the same order, even higher, than that of H<sub>2</sub>S and OCS (Codella et al. 2021). High angular observations of these molecules could help disentangle the sulfur chemistry in these cores.

Our observations and simulations show that the composition
of sulfur species in the ice formed during the pre-stellar and pro-
tostellar phase, and also later, in the warm gas enriched in sulfur
compounds because of ice evaporation, are extremely sensitive
to the thermal history of the gas and dust. We note the variations
of more than one order of magnitude in the abundances of rou-
tinely observed species such as H<sub>2</sub>S and OCS. This challenges
the estimation of sulfur depletion based on the observations of
only these molecules. The knowledge of the ice composition at
the end of the pre-stellar phase would be essential to obtain accu-
rate predictions of the abundances of all sulfur compounds in the
warm phase. Observation of the chemical composition of ices
would be a valuable input to constrain chemical predictions.

The detection of sulfur species in solid phase remains chal-
lenging even in the James Webb Space Telescope (JWST) era.
Thus far, only OCS and SO₂ have been detected in interstellar
ices (Geballe et al. 1985; Palumbo et al. 1997; Boogert et al.
1997, 2022; McClure et al. 2023), while solid H₂S remains
undetected. Recent measurements carried out with the JWST
have improved previous upper limits of the H₂S/OCS in the ice.
Yet, the limits are far from constraining our knowledge of sul-
fur chemistry. In Fig. 11, we compare the values of H₂S/OCS
obtained toward the warm protostellar cores using millimeter
interferometers with the upper limits obtained for the solid phase
<figure><img src="image_12.png" /><figcaption>**Fig. 11.** H₂S/OCS ratio in the 24 protostars of our sample. The class of each protostar has been represented using the color code in the legend. Sources are sorted by class, and then ordered by increasing bolometric temperature within each class. The discontinuous horizontal lines represent different upper limits to the H₂S/OCS ratio in ices from different regions: in black, W33A; in green, Mon R2 IRS2; in red, ChaMMS1. Only the OCS-rich sources are compatible with the more restrictive limit, ChaMMS1, which could give us some hints about the environment where the ice was formed.</figcaption></figure>

using infrared telescopes. Some of the gas-phase measurements are higher than the H₂S/OCS ratio obtained by McClure et al. (2023) for the ice toward ChaMMS1, but most of them remain compatible with it. Taking into account the great sensitivity if the icy H₂S/OCS ratio to the local physical conditions and past thermal history of the interstellar grains, current observations are far to be conclusive and new searches of solid H₂S are necessary to increase our knowledge of sulfur chemistry and further contrast our models.

## 9. Summary and conclusions

In this work, we have studied several sulfur-bearing species – H₂S, H³³S, OCS, OC³³S, OC³⁴S, SO, and SO₂ – in the warm inner core of 24 Class 0/I protostars in the Perseus Molecular Cloud. We derived the column densities of the first five by fitting the gas emission lines in the NOEMA Band 3 receiver and the PolyFix correlator spectra from the PRODIGE large program observations, and we used the SO and SO₂ data from Artur de la Villarmois et al. (2023). We also estimated the abundances of all seven molecules using our observations of the continuum, and we analyzed the variations of the H₂S/OCS ratio over the 24 sources. Our main conclusions are the following:

- We calculated the gas-phase column densities of H₂S, H³³S, OCS, OC³³S, and OC³⁴S NOEMA Band 3 data and the N(H₂S) ratio over N(OCS). We detected H₂S in 20 of the 24 sources in our sample and OCS in 17 of them;

- The protostars can be characterized by their H₂S and OCS composition. We differentiated two kinds of objects: OCS-poor protostars, with H₂S/OCS~7–10, and OCS-rich protostars, with H₂S/OCS~0.1–1;

- The total sulfur abundance grows with evolution until the Class 0/I stage, where the minumum depletion is achieved: Dₛ ~8, with a total S abundance of (1.7±0.1)×10⁻⁶. This value is in accordance with the H₂S abundance measured in Orion KL: (3.1±1.9)×10⁻⁶ (Crockett et al. 2014). In more evolved stages (Class I), the total sulfur abundance decreases, though more observations of Class I/II sources

would be needed to confirm this. The total sulfur abundance is usually dominated by the H₂S molecule, followed by OCS. The SO and SO₂ abundances are typically two orders of magnitude lower than H₂S and OCS; however, SO and SO₂ abundances grow slowly toward the Class I stage, gaining relevance in the sulfur budget in the more evolved sources;

- Theoretical numerical simulations showed that H₂S and OCS gas-phase and ice abundances are very responsive to temperature changes both in the pre-stellar phase and during the collapse. H₂S is less abundant when the initial temperature of the cloud is warmer, resulting in lower amounts of OCS and eventually greater H₂S/OCS ratios after the collapse. Also, a lower temperature of the inner core during the collapse showed a fall in the H₂S abundance, leading to lower H₂S/OCS ratios. Furthermore, the 7.1 and 7.2 (15 K) models (Table 3) return significant variations of the H₂S and OCS in ices after the pre-phase and a factor of approximately two difference of the H₂S/OCS ratio after the collapse, with only slightly different conditions of the initial cloud. Results suggest that in some scenarios, H₂Sₓ molecules, for x ≥ 3, might be a significant sulfur reservoir;

- We found a wide range of values for the depletion of sulfur (8 to >100) in our sample of protostars. This behavior only seems to be explainable with the differences in the initial conditions and in the first evolutionary stages. We consider that the details of the ice composition after the pre-stellar phase (before the collapse) are essential to understanding the chemical composition of a protostar after the collapse. Also, higher angular resolution observations are required to resolve the circumstellar disk to probe the innermost and warmest region.

The NOEMA interferometer observations have allowed us to study the sulfur chemistry of the warm inner core of 24 protostars in Perseus. Our results show that the evolution history of a protostar plays a decisive role in the chemistry of its later stages. However, the low spatial resolution of the images and the low spectral resolution of the H₂S, H³³S, and OC³⁴S bands invites the use of higher resolution observations in order to better resolve the inner core and better estimate the column densities
of some species. Nevertheless, the PRODIGE data reveal that the sulfur chemistry of young protostars is still not well understood and that further study is needed to comprehend how the different environmental and evolutionary factors affect the composition of these young stellar objects.

**Acknowledgements.** This work is supported by ERC grant SUL4LIFE, GA No. 101096293. Funded by the European Union. Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council Executive Agency. Neither the European Union nor the granting authority can be held responsible for them. J.E.P., D.M.S.-C., P.C., M.T.V.-M., T.-H.H., L.B., C.G. Y.-R.C. and M.M. are grateful for support from the Max Planck Society. R.M.-D. was supported by a La Caixa Junior Leader grant under agreement LCF/BQ/PI22/11910030. This work is based on observations carried out under project number L19MB with the IRAM NOEMA Interferometer. IRAM is supported by INSU/CNRS (France), MPG (Germany) and IGN (Spain). 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), NSTC 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.

## References

Altwegg, K., Combi, M., Fuselier, S. A., et al. 2022, MNRAS, 516, 3900

Anders, E., & Grevesse, N. 1989, Geochim. Cosmochim. Acta, 53, 197

André, P., Men'shchikov, A., Bontemps, S., et al. 2010, A&A, 518, L102

Aponte, J. C., Dworkin, J. P., Glavin, D. P., et al. 2023, Earth Planets Space, 75, 148

Artur de la Villarmois, E., Guzmán, V. V., Jørgensen, J. K., et al. 2022, A&A, 667, A20

Artur de la Villarmois, E., Guzmán, V. V., Yang, Y. L., Zhang, Y., & Sakai, N. 2023, A&A, 678, A124

Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481

Bachiller, R., Cernicharo, J., Martin-Pintado, J., Tafalla, M., & Lazareff, B. 1990, A&A, 231, 174

Bachiller, R., Martin-Pintado, J., & Fuente, A. 1991, A&A, 243, L21

Bachiller, R., Guilloteau, S., Gueth, F., et al. 1998, A&A, 339, L49

Bariosco, V., Pantaleone, S., Ceccarelli, C., et al. 2024, MNRAS, 531, 1371

Barsony, M., Ward-Thompson, D., André, P., & O'Linger, J. 1998, ApJ, 509, 733

Bianchi, E., Codella, C., Ceccarelli, C., et al. 2019, MNRAS, 483, 1850

Bianchi, E., Ceccarelli, C., Codella, C., et al. 2022, A&A, 662, A103

Boogert, A. C. A., Schutte, W. A., Helmich, F. P., Tielens, A. G. G. M., & Wooden, D. H. 1997, A&A, 317, 929

Boogert, A. C. A., Brewer, K., Brittain, A., & Emerson, K. S. 2022, ApJ, 941, 32

Brinch, C., Jørgensen, J. K., & Høgerheijde, M. R. 2009, A&A, 502, 199

Calmonte, U., Altwegg, K., Balsiger, H., et al. 2016, MNRAS, 462, S253

Calmonte, U., Altwegg, K., Balsiger, H., et al. 2017, MNRAS, 469, S787

Caratti o Garatti, A., Ray, T. P., Kavanagh, P. J., et al. 2024, A&A, 691, A134

Carrascosa, H., Muñoz Caro, G. M., Martín-Doménech, R., et al. 2024, MNRAS, 533, 967

Caselli, P., Hartquist, T. W., & Havnes, O. 1997, A&A, 322, 296

Caselli, P., Pineda, J. E., Sipilä, O., et al. 2022, ApJ, 929, 13

Cazaux, S., Carrascosa, H., Muñoz Caro, G. M., et al. 2022, A&A, 657, A100

Ceccarelli, C. 2007, in Molecules in Space and Laboratory, eds. J. L. Lemaire, & F. Combes, 1

Chen, M., & Yu, X. 2019, in Carbonyl Sulfide-Mediated Synthesis of Peptides with Amino Acid Ionic Liquids

Chen, Y., Rocha, W. R. M., van Dishoeck, E. F., et al. 2024, A&A, 690, A205

Chin, Y. N., Henkel, C., Whiteoak, J. B., Langer, N., & Churchwell, E. B. 1996, A&A, 305, 960

Codella, C., Ceccarelli, C., Bianchi, E., et al. 2016, MNRAS, 462, L75

Codella, C., Bianchi, E., Podio, L., et al. 2021, A&A, 654, A52

Collings, M. P., Dever, J. W., Fraser, H. J., McCoustra, M. R. S., & Williams, D. A. 2003, ApJ, 583, 1058

Crockett, N. R., Bergin, E. A., Neill, J. L., et al. 2014, ApJ, 781, 114

Daflon, S., Cunha, K., de la Reza, R., Holtzman, J., & Chiappini, C. 2009, AJ, 138, 1577

Dionatos, O., Nisini, B., Cabrit, S., Kristensen, L., & Pineau Des Forêts, G. 2010, A&A, 521, A7

Dionatos, O., Ray, T., & Güdel, M. 2018, A&A, 616, A84

Drozdovskaya, M. N., van Dishoeck, E. F., Jørgensen, J. K., et al. 2018, MNRAS, 476, 4949

Dutrey, A., Guilloteau, S., & Bachiller, R. 1997, A&A, 325, 758

el Akel, M., Kristensen, L. E., Le Gal, R., et al. 2022, A&A, 659, A100

Enoch, M. L., Evans, Neal J., I., Sargent, A. I., & Glenn, J. 2009, ApJ, 692, 973

Esplugues, G. B., Tercero, B., Cernicharo, J., et al. 2013, A&A, 556, A143

Fiorellino, E., Manara, C. F., Nisini, B., et al. 2021, A&A, 650, A43

Froebrich, D. 2005, ApJS, 156, 169

Fuente, A., Cernicharo, J., Roueff, E., et al. 2016, A&A, 593, A94

Fuente, A., Gerin, M., Pety, J., et al. 2017a, A&A, 606, L3

Fuente, A., Goicoechea, J. R., Pety, J., et al. 2017b, ApJ, 851, L49

Fuente, A., Rivière-Marichalar, P., Beitia-Antero, L., et al. 2023, A&A, 670, A114

Fuente, A., Roueff, E., Le Petit, F., et al. 2024, A&A, 687, A87

Geballe, T. R., Baas, F., Greenberg, J. M., & Schutte, W. 1985, A&A, 146, L6

Gerin, M., Pety, J., Fuente, A., et al. 2015, A&A, 577, L2

Gerin, M., Pety, J., Commerçon, B., et al. 2017, A&A, 606, A35

Gieser, C., Pineda, J. E., Segura-Cox, D. M., et al. 2024, A&A, 692, A55

Goicoechea, J. R., & Cuadrado, S. 2021, A&A, 647, L7

Goldsmith, P. F., & Langer, W. D. 1999, ApJ, 517, 209

Gratier, P., Majumdar, L., Ohishi, M., et al. 2016, ApJS, 225, 25

Grossman, E. N., Masson, C. R., Sargent, A. I., et al. 1987, ApJ, 320, 356

Gueth, F., & Guilloteau, S. 1999, A&A, 343, 571

Hatchell, J., & Dunham, M. M. 2009, A&A, 502, 139

Hatchell, J., & Fuller, G. A. 2008, A&A, 482, 855

Hatchell, J., Wilson, T., Drabek, E., et al. 2013, MNRAS, 429, L10

Hennebelle, P., Commerçon, B., Chabrier, G., & Marchand, P. 2016, ApJ, 830, L8

Herczec, G. J., Karska, A., Bruderer, S., et al. 2012, A&A, 540, A84

HGBS team 2020, Herschel Gould Belt Survey, https://doi.org/10.26131/IRSA72

Hily-Blant, P., Pineau des Forêts, G., Faure, A., & Lique, F. 2022, A&A, 658, A168

Hirano, N., Kamazaki, T., Mikami, H., Ohashi, N., & Umemoto, T. 1999, in Star Formation 1999, ed. T. Nakamoto, 181

Hirano, N., Ho, P. P. T., Liu, S.-Y., et al. 2010, ApJ, 717, 58

Hirano, N., & Liu, F.-c. 2014, ApJ, 789, 50

Hirota, T., Bushimata, T., Choi, Y. K., et al. 2008, PASJ, 60, 37

Holdship, J., Viti, S., Jimenez-Serra, I., et al. 2016, MNRAS, 463, 802

Hsieh, T. H., Segura-Cox, D. M., Pineda, J. E., et al. 2023, A&A, 669, A137

Hsieh, T. H., Pineda, J. E., Segura-Cox, D. M., et al. 2024, A&A, 686, A289

Jennings, R. E., Cameron, D. H. M., Cudlip, W., & Hirst, C. J. 1987, MNRAS, 226, 461

Jiménez-Escobar, A., & Muñoz Caro, G. M. 2011, A&A, 536, A91

Jiménez-Serra, I., Caselli, P., Martín-Pintado, J., & Hartquist, T. W. 2008, A&A, 482, 549

Jiménez-Escobar, A., Muñoz Caro, G. M., Ciaravella, A., et al. 2012, ApJ, 751, L40

Jørgensen, J. K., Schöier, F. L., & van Dishoeck, E. F. 2002, A&A, 389, 908

Jørgensen, J. K., Høgerheijde, M. R., van Dishoeck, E. F., Blake, G. A., & Schöier, F. L. 2004, A&A, 413, 993

Jørgensen, J. K., Bourke, T. L., Myers, P. C., et al. 2008, ApJ, 632, 973

Jørgensen, J. K., Harvey, P. M., Evans, Neal J., I., et al. 2006, ApJ, 645, 1246

Jørgensen, J. K., van Dishoeck, E. F., Visser, R., et al. 2019, A&A, 507, 861

Jørgensen, J. K., Belloche, A., & Garrod, R. T. 2020, ARA&A, 58, 727

Kauffmann, J., Bertoldi, F., Bourke, T. L., Evans, N. J., I., & Lee, C. W. 2008, A&A, 487, 993

Keller, L. P., Hony, S., Bradley, J. P., et al. 2002, Nature, 417, 148

Laas, J. C., & Caselli, P. 2019, A&A, 624, A108

Le Gal, R., Öberg, K. I., Loomis, R. A., Pegues, J., & Bergner, J. B. 2019, ApJ, 876, 72

Le Roy, L., Altwegg, K., Balsiger, H., et al. 2015, A&A, 583, A1

Lee, C.-F., Hirano, N., Palau, A., et al. 2009, ApJ, 699, 1584

Lee, K. I., Dunham, M. M., Myers, P. C., et al. 2015, ApJ, 814, 114

Lefèvre, C., Cabrit, S., Maury, A. J., et al. 2017, A&A, 604, L1

Lefloch, B., Castets, A., Cernicharo, J., Langer, W. D., & Zylka, R. 1998, A&A, 334, 269

Leman, L. J., Orgel, L. E., & Ghadiri, M. R. 2004, Science, 306, 283

Lin, S.-J., Yen, H.-W., & Lai, S.-P. 2024, AJ, 168, 107

Looney, L. W., Mundy, L. G., & Welch, W. J. 2000, ApJ, 529, 477

Machida, M. N., Inutsuka, S.-i., & Matsumoto, T. 2006, ApJ, 647, L151

Marcelino, N., Gerin, M., Cernicharo, J., et al. 2018, A&A, 620, A80

Martín-Doménech, R., Bergner, J. B., Öberg, K. I., et al. 2021, ApJ, 923, 155

Martín-Doménech, R., Öberg, K. I., Muñoz Caro, G. M., et al. 2024, MNRAS, 535, 807

Maury, A. J., Belloche, A., André, P., et al. 2014, A&A, 563, L2

McClure, M. K., Rocha, W. R. M., Pontoppidan, K. M., et al. 2023, Nat. Astron., 7, 431

McGuire, B. A. 2022, ApJS, 259, 30

Mercimek, S., Codella, C., Podio, L., et al. 2022, A&A, 659, A67

Müller, H. S. P., Schlöder, F., Stutzki, J., & Winnewisser, G. 2005, J. Mol. Struct., 742, 215
Navarro-Almaida, D., Le Gal, R., Fuente, A., et al. 2020, A&A, 637, A39

Navarro-Almaida, D., Lebreuilly, U., Hennebelle, P., et al. 2024, A&A, 685, A112

Neufeld, D. A., Godard, B., Gerin, M., et al. 2015, A&A, 577, A49

Olson, K. R., & Straub, K. D. 2016, Physiology, 31, 60

Ossenkopf, V., & Henning, T. 1994, A&A, 291, 943

Palau, A., Zapata, L. A., Rodríguez, L. F., et al. 2014, MNRAS, 444, 833

Palumbo, M. E., Tielens, A. G. G. M., & Tokunaga, A. T. 1995, ApJ, 449, 674

Palumbo, M. E., Geballe, T. R., & Tielens, A. G. G. M. 1997, ApJ, 479, 839

Perrero, J., Beitia-Antero, L., Fuente, A., Ugliengo, P., & Rimola, A. 2024, ApJ, 971, 36

Pineda, J. E., Segura-Cox, D., Caselli, P., et al. 2020, Nat. Astron., 4, 1158

Pineda, J. E., Harju, J., Caselli, P., et al. 2022, AJ, 163, 294

Ray, T. P., McCaughrean, M. J., Caratti o Garatti, A., et al. 2023, Nature, 622, 48

Reipurth, B., Heathcote, S., Roth, M., Noriega-Crespo, A., & Raga, A. C. 1993, ApJ, 408, L49

Reynolds, N. K., Tobin, J. J., Sheehan, P., et al. 2021, ApJ, 907, L10

Rivière-Marichalar, P., Fuente, A., Le Gal, R., et al. 2020, A&A, 642, A32

Rocha, W. R. M., van Dishoeck, E. F., Ressler, M. E., et al. 2024, A&A, 683, A124

Rodríguez, L. F., Zapata, L. A., & Palau, A. 2014, ApJ, 790, 80

Rodríguez-Baras, M., Fuente, A., Riviére-Marichalar, P., et al. 2021, A&A, 648, A120

Ruaud, M., Wakelam, V., & Hersant, F. 2016, MNRAS, 459, 3756

Ruffle, D. P., Hartquist, T. W., Caselli, P., & Williams, D. A. 1999, MNRAS, 306, 691

Sadavoy, S. I., Di Francesco, J., André, P., et al. 2014, ApJ, 787, L18

Sandell, G., Aspin, C., Duncan, W. D., Russell, A. P. G., & Robson, E. I. 1991, ApJ, 376, L17

Santos, J. C., van Gelder, M. L., Nazari, P., Ahmadi, A., & van Dishoeck, E. F. 2024, A&A, 689, A248

Schilke, P., Walmsley, C. M., Pineau des Forets, G., & Flower, D. R. 1997, A&A, 321, 293

Semenov, D., Favre, C., Fedele, D., et al. 2018, A&A, 617, A28

Shingledecker, C. N., Lamberts, T., Laas, J. C., et al. 2020, ApJ, 888, 52

Stephens, I. W., Bourke, T. L., Dunham, M. M., et al. 2019, ApJS, 245, 21

Taillard, A., Martín-Doménech, R., Carrascosa, H., et al. 2025, A&A, 694, A263

Taniguchi, K., Pineda, J. E., Caselli, P., et al. 2024, ApJ, 965, 162

Tobin, J. J., & Sheehan, P. D. 2024, ARA&A, 62, 203

Tobin, J. J., Looney, L. W., Mundy, L. G., Kwon, W., & Hamidouche, M. 2007, ApJ, 659, 1404

Tobin, J. J., Dunham, M. M., Looney, L. W., et al. 2015a, ApJ, 798, 61

Tobin, J. J., Looney, L. W., Wilner, D. J., et al. 2015b, ApJ, 805, 125

Tobin, J. J., Kratter, K. M., Persson, M. V., et al. 2016a, Nature, 538, 483

Tobin, J. J., Looney, L. W., Li, Z.-Y., et al. 2016b, ApJ, 818, 73

Tobin, J. J., Looney, L. W., Li, Z.-Y., et al. 2018, ApJ, 867, 43

Valdivia-Mena, M. T., Pineda, J. E., Segura-Cox, D. M., et al. 2022, A&A, 667, A12

Valdivia-Mena, M. T., Pineda, J. E., Segura-Cox, D. M., et al. 2023, A&A, 677, A92

van Gelder, M. L., Tabone, B., Tychoniec, Ł., et al. 2020, A&A, 639, A87

van Gelder, M. L., Nazari, P., Tabone, B., et al. 2022, A&A, 662, A67

Vastel, C., Quénard, D., Le Gal, R., et al. 2018, MNRAS, 478, 5514

Vidal, T. H. G., & Wakelam, V. 2018, MNRAS, 474, 5575

Vitorino, J., Loison, J. C., Wakelam, V., Congiu, E., & Dulieu, F. 2024, MNRAS, 533, 52

Wakelam, V., Hersant, F., & Herpin, F. 2011, A&A, 529, A112

Wakelam, V., Loison, J. C., Mereau, R., & Ruaud, M. 2017, Mol. Astrophys., 6, 22

Wakelam, V., Gratier, P., Loison, J. C., et al. 2024, A&A, 689, A63

Walmsley, C. M., Flower, D. R., & Pineau des Forêts, G. 2004, A&A, 418, 1035

Yang, Y.-L., Sakai, N., Zhang, Y., et al. 2021, ApJ, 910, 20

Zapata, L. A., Arce, H. G., Brassfield, E., et al. 2014, MNRAS, 441, 3696

Zhang, Z. E., Yang, Y.-L., Zhang, Y., et al. 2023, ApJ, 946, 113
## Appendix A: The H<sub>2</sub><sup>33</sup>S isotopologue

H<sub>2</sub>S is a very abundant molecule in the core of YSOs, showing high column densities too. Observations of the lower excitation temperature lines can result in optically thick lines that do not represent the actual H<sub>2</sub>S budget, and, sometimes, the observation of the less abundant isotopologues is needed in order to achieve an accurate column density estimation.

The H<sub>2</sub><sup>33</sup>S molecule presents two important structural behaviors. On the one hand, the molecule has a para and an orto variants, with a distribution of 1:3 respectively, as in the case of H<sub>2</sub>S. On the other hand, this species exhibits hyperfine structure, which is discernible at lower excitation levels. The transition observed in this work is actually a low-level transition – H<sub>2</sub><sup>33</sup>S: 2<sub>2,0</sub>–2<sub>1,1</sub>. Therefore, we took into account the line splitting when we calculated the column density. Table A.1 displays the hyperfine structure of the molecule for the 2<sub>2,0</sub>–2<sub>1,1</sub> transition. Furthermore, there is an additional factor to take into account, as Acetaldehyde – CH<sub>3</sub>CHO – has a relatively bright transition (~33% of the maximum H<sub>2</sub><sup>33</sup>S intensity) in the same range of frequencies. To deal with this part, we masked the CH<sub>3</sub>CHO peak before the fitting with CLASS by fittin only the not contaminated hyperfine components, in order to avoid misrepresentations of the total emission or the opacity of the line. In Fig. A.1 we show the spectrum of IRAS4B in the spectral window of H<sub>2</sub><sup>33</sup>S and CH<sub>3</sub>CHO. When fitting a line with hyperfine splitting with CLASS, the program also estimates the opacity for the whole transition. The least value possible for the opacity is 0.1, so when the transition is not optically thick, it will return this result. The opacities estimated with CLASS were very low in all cases (0.1) except for IRAS2A and SVS13A, for which we obtained higher opacities. In these two cases, our estimates should be understood only as a lower limit to the real column density.

<table><caption>Table A.1: Hyperfine structure of H<sub>2</sub><sup>33</sup>S 2<sub>2,0</sub>–2<sub>1,1</sub> and nearby CH<sub>3</sub>CHO transition. Data were retrieved or calculated from the CDMS<sup>6</sup> database.</caption><thead><tr><th>Molecule</th><th>ν (GHz)</th><th>Transition</th><th>I (log<sub>10</sub>)</th><th>Norm. I</th></tr></thead><tbody><tr><td rowspan="8">H<sub>2</sub><sup>33</sup>S</td><td>215.4945</td><td>2<sub>2,0,1</sub>–2<sub>1,1,2</sub></td><td>-3.7205</td><td>0.1458</td></tr><tr><td>215.4967</td><td>2<sub>2,0,1</sub>–2<sub>1,1,1</sub></td><td>-3.7205</td><td>0.1458</td></tr><tr><td>215.5008</td><td>2<sub>2,0,4</sub>–2<sub>1,1,3</sub></td><td>-3.6625</td><td>0.1666</td></tr><tr><td>215.5029</td><td>2<sub>2,0,4</sub>–2<sub>1,1,4</sub></td><td>-2.8843</td><td>1.0</td></tr><tr><td>215.5037</td><td>2<sub>2,0,2</sub>–2<sub>1,1,3</sub></td><td>-3.5743</td><td>0.2042</td></tr><tr><td>215.5054</td><td>2<sub>2,0,2</sub>–2<sub>1,1,2</sub></td><td>-3.5163</td><td>0.2333</td></tr><tr><td>215.5076</td><td>2<sub>2,0,2</sub>–2<sub>1,1,1</sub></td><td>-3.7204</td><td>0.1458</td></tr><tr><td>215.5116</td><td>2<sub>2,0,3</sub>–2<sub>1,1,3</sub></td><td>-3.1817</td><td>0.5042</td></tr><tr><td>215.5132</td><td>2<sub>2,0,3</sub>–2<sub>1,1,2</sub></td><td>-3.5743</td><td>0.2042</td></tr><tr><td>215.5136</td><td>2<sub>2,0,3</sub>–2<sub>1,1,4</sub></td><td>-3.6624</td><td>0.1666</td></tr><tr><td>CH<sub>3</sub>CHO</td><td>215.5118</td><td>11(2 9)6–10(2 8)6</td><td>-3.9744</td><td>0.3251</td></tr></tbody></table>

**Notes.** The normalised intensity is the result of dividing each individual intensity (without the log<sub>10</sub>) by the more intense of the lines of the hyperfine splitting, which in this case is 2<sub>2,0,4</sub>–2<sub>1,1,4</sub>.

<figure><img src="image_13.png" /><figcaption>Fig. A.1: Spectrum of IRAS4B around the 215 GHz H<sub>2</sub><sup>33</sup>S line. In black, the integrated emission in a 1.5" circle around the protostar's position. In red, the fit to the H<sub>2</sub><sup>33</sup>S emission. In yellow, the emission of CH<sub>3</sub>CHO, assuming its column density is three times greater than that of H<sub>2</sub><sup>33</sup>S. The purple dashed line represents the sum of the emission of both species. The different lines from Table A.1 appear marked in the bottom side of the figure with small vertical lines. For reference, the most intense H<sub>2</sub><sup>33</sup>S line (the fourth starting from the left) is at a velocity of -3.1 km/s.</figcaption></figure>

## Appendix B: IRAS2A and IRAS4B extended emission integrated maps

The following figures show the extended emission of IRAS2A and IRAS4B as a complement of Fig. 3. The extended emission is found in the north-south outflows direction, and is specially noticeable in IRAS4B.

<figure><img src="image_14.png" /><figcaption>Fig. B.1: Integrated emission of H<sub>2</sub>S and OCS in IRAS2A in a 8"×16" region. (See Fig. 3 caption for more details.)</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Fig. B.2: Integrated emission of H<sub>2</sub>S and OCS in IRAS4B in a 8"×16" region. (See Fig. 3 caption for more details.)</figcaption></figure>

<footnote n="6">6 https://cdms.astro.uni-koeln.de/</footnote>
## Appendix C: H₂S opacity worst-case scenario

Figure C.1 shows the comparison between the measured H₂S/OCS ratios in the sources from our sample, and the corresponding H₂S/OCS ratios in the worst case scenario. As explained at the end of Section 5, this scenario considers that, in the OCS-rich protostars were OC³⁴S was detected but H₂³³S was not, the correction for the H₂S opacity is of a factor of 30. We find that, with these assumptions, the two-family segregation into OCS-poor and OCS-rich sources remains valid.

<figure><img src="image_16.png" /></figure>

Fig. C.1: Ratio of H₂S/OCS in the 24 protostars of our sample. The class of each protostar has been represented using the color code in the legend. Sources are sorted by class, and then ordered by increasing bolometric temperature within each class. Panel a) shows the ratios calculated with the estimated column densities of H₂S and OCS. Panel b) shows the same data but assuming a worst-case scenario where all protostars with OC³⁴S detection but no H₂³³S detection have their H₂S column density scaled by a factor 37, to take into account possible optically thick emission in the H₂S line. The scaled ratios have been represented with a downward triangle, which is a superior limit to the maximum expected H₂S/OCS ratio after accounting for the opacity. In this case, the two differentiated trends still appear.

## Appendix D: Complementary tables

In Tables D.1, D.2, D.3 and D.4, we present crucial information about the 24 sources that form our sample, including tabulated data from previous works and results from our observations. All values have been calculated with their corresponding uncertainties, but they have not been included in these tables for the sake of conciseness and clarity. The complete D.3 and D.4 tables, as well as the spectra and an extra table regarding additional information about the moment-0 maps (see Section 3), can be consulted at https://github.com/JMiranzo/PRODIGE-MissingSulfurProblem.
<table><caption>Table D.1: General source information of the 24 protostars in our sample.</caption><thead><tr><th>Source</th><th>Other Names</th><th>Subregion</th><th>RA (J2000)<br>(hh:mm:ss)</th><th>DEC (J2000)<br>(dd:mm:ss)</th><th>Class</th><th>L<sub>bol</sub><sup>(a)</sup><br>(L<sub>☉</sub>)</th><th>T<sub>bol</sub><sup>(a)</sup><br>(K)</th><th>T<sub>kin</sub><sup>(b)</sup><br>(K)</th><th>Hot<br>Corino</th></tr></thead><tbody><tr><td>B1bN</td><td></td><td>Barnard 1</td><td>03:33:21.21</td><td>31:07:43.6</td><td>0</td><td>0.32±0.10</td><td>14.7±1.0</td><td>100.0</td><td>N</td></tr><tr><td>B1bS</td><td></td><td>Barnard 1</td><td>03:33:21.36</td><td>03:07:26.3</td><td>0</td><td>0.70±0.01</td><td>17.7±1.0</td><td>100.0</td><td>Y</td></tr><tr><td>B5-IRS1</td><td>Per-emb-53</td><td>Barnard 5</td><td>03:47:41.59</td><td>32:51:43.6</td><td>I</td><td>4.70±0.90</td><td>287.0±8.0</td><td>287.0</td><td>N</td></tr><tr><td>HH211MMS</td><td>Per-emb-01</td><td>IC348</td><td>03:43:56.81</td><td>32:00:50</td><td>0</td><td>1.80±0.10</td><td>27.0±1.0</td><td>100.0</td><td>Y</td></tr><tr><td>IC348MMS</td><td>Per-emb-11</td><td>IC348</td><td>03:43:57.07</td><td>32:03:04.8</td><td>0</td><td>1.50±0.10</td><td>30.0±2.0</td><td>100.0</td><td>Y</td></tr><tr><td>IRAS2A</td><td>Per-emb-27</td><td>NGC1333</td><td>03:28:55.57</td><td>31:14:37.0</td><td>0/I</td><td>19.00±0.40</td><td>69.0±1.0</td><td>100.0</td><td>Y</td></tr><tr><td>IRAS2B</td><td>Per-emb-36</td><td>NGC1333</td><td>03:28:57.38</td><td>31:14:15.8</td><td>I</td><td>5.30±1.00</td><td>106.0±12.0</td><td>106.0</td><td>N</td></tr><tr><td>IRAS4B</td><td>Per-emb-13</td><td>NGC1333</td><td>03:29:12.02</td><td>31:13:08.0</td><td>0</td><td>4.00±0.30</td><td>28.0±1.0</td><td>100.0</td><td>Y</td></tr><tr><td>IRAS4C</td><td>Per-emb-14</td><td>NGC1333</td><td>03:29:13.55</td><td>31:13:58.1</td><td>0</td><td>0.70±0.08</td><td>31.0±2.0</td><td>100.0</td><td>N</td></tr><tr><td>L1448-IRS3A</td><td></td><td>L1448</td><td>03:25:36.50</td><td>30:45:21.9</td><td>0</td><td>9.20±1.30</td><td>47.0±2.0</td><td>100.0</td><td>N</td></tr><tr><td>L1448-IRS3B</td><td>Per-emb-33</td><td>L1448</td><td>03:25:36.38</td><td>30:45:14.7</td><td>0</td><td>8.30±0.80</td><td>57.0±3.0</td><td>100.0</td><td>N</td></tr><tr><td>L1448C</td><td>Per-emb-26</td><td>L1448</td><td>03:25:38.88</td><td>30:44:05.3</td><td>0</td><td>8.40±1.50</td><td>47.0±7.0</td><td>100.0</td><td>Y</td></tr><tr><td>L1448NW</td><td></td><td>L1448</td><td>03:25:35.67</td><td>30:45:34.2</td><td>0</td><td>1.40±0.10</td><td>22.0±1.0</td><td>100.0</td><td>N</td></tr><tr><td>Per-emb-02</td><td></td><td>Barnard 1</td><td>03:32:17.92</td><td>30:49:47.8</td><td>0</td><td>0.90±0.07</td><td>27.0±1.0</td><td>100.0</td><td>N</td></tr><tr><td>Per-emb-05</td><td></td><td>...</td><td>03:31:20.94</td><td>30:45:30.2</td><td>0</td><td>1.30±0.10</td><td>32.0±2.0</td><td>100.0</td><td>N</td></tr><tr><td>Per-emb-08</td><td></td><td>IC348</td><td>03:44:43.98</td><td>32:01:35.2</td><td>0</td><td>2.60±0.50</td><td>43.0±6.0</td><td>100.0</td><td>N</td></tr><tr><td>Per-emb-18</td><td>IRAS7</td><td>NGC1333</td><td>03:29:11.27</td><td>31:18:31.1</td><td>0</td><td>2.80±1.70</td><td>59.0±12.0</td><td>100.0</td><td>N</td></tr><tr><td>Per-emb-22</td><td>L1448-IRS2</td><td>L1448</td><td>03:25:22.41</td><td>30:45:13.3</td><td>0</td><td>3.60±0.50</td><td>43.0±2.0</td><td>100.0</td><td>N</td></tr><tr><td>Per-emb-29</td><td>B1c</td><td>Barnard 1</td><td>03:33:17.88</td><td>31:09:31.7</td><td>0</td><td>3.70±0.40</td><td>48.0±1.0</td><td>100.0</td><td>N</td></tr><tr><td>Per-emb-30</td><td></td><td>Barnard 1</td><td>03:33:27.28</td><td>31:07:10.2</td><td>0/I</td><td>1.70±0.01</td><td>78.0±6.0</td><td>100.0</td><td>N</td></tr><tr><td>Per-emb-50</td><td></td><td>NCG1333</td><td>03:29:07.77</td><td>31:21:57.1</td><td>I</td><td>23.20±3.00</td><td>128.0±23.0</td><td>128.0</td><td>N</td></tr><tr><td>Per-emb-62</td><td></td><td>IC348</td><td>03:44:12.98</td><td>32:01:35.4</td><td>I</td><td>1.80±0.40</td><td>378.0±29.0</td><td>378.0</td><td>N</td></tr><tr><td>SVS13A</td><td>Per-emb-44</td><td>NGC1333</td><td>03:29:03.76</td><td>31:16:03.7</td><td>0/I</td><td>32.50±7.10</td><td>188.0±9.0</td><td>188.0</td><td>Y</td></tr><tr><td>SVS13B</td><td></td><td>NGC1333</td><td>03:29:03.08</td><td>31:15:51.7</td><td>0</td><td>1.00±1.00</td><td>20.0±20.0</td><td>100.0</td><td>N</td></tr></tbody></table>

Notes. (a) Data extracted from Tobin et al. (2016b). (b) We have assumed the $T_{\text{kin}}$ of the warm inner core to be the maximum between 100K and $T_{\text{kin}}$ (which implies, $T \ge 100$K). (c) Close to the Barnard 1 region, but isolated.
<table><caption>Table D.2: Detected and undetected transitions in the warm inner cores of the 24 observed protostars.</caption><thead><tr><th>Protostar</th><th>Class</th><th>H<sub>2</sub>S</th><th>H<sub>2</sub><sup>33</sup>S</th><th>OCS</th><th>OCS</th><th>OC<sup>33</sup>S</th><th>OC<sup>34</sup>S</th><th>SO</th><th>SO<sub>2</sub><sup>(a)</sup></th></tr><tr><td></td><td></td><td>2<sub>2,0</sub>-2<sub>1,1</sub></td><td>2(2,0,4)-2(1,1,4)</td><td>J=18-17</td><td>J=19-18</td><td>J=18-17</td><td>J=20-19</td><td>14(0,14)-14(1,13)</td><td></td></tr></thead><tbody><tr><td>B1bN</td><td>0</td><td>-</td><td>✓</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>B1bS</td><td>0</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>-</td><td>-</td></tr><tr><td>B5-IRS1</td><td>I</td><td>✓</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>✓</td><td>-</td></tr><tr><td>HH211MMS</td><td>0</td><td>✓</td><td>-</td><td>✓</td><td>✓</td><td>-</td><td>✓</td><td>✓</td><td>✓</td></tr><tr><td>IC348MMS</td><td>0</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td></tr><tr><td>IRAS2A</td><td>0/I</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td></tr><tr><td>IRAS2B</td><td>I</td><td>✓</td><td>-</td><td>✓</td><td>✓</td><td>-</td><td>-</td><td>✓</td><td>✓</td></tr><tr><td>IRAS4B</td><td>0</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td></tr><tr><td>IRAS4C</td><td>0</td><td>✓</td><td>-</td><td>✓</td><td>✓</td><td>-</td><td>-</td><td>✓</td><td>-</td></tr><tr><td>L1448-IRS3A</td><td>0</td><td>✓</td><td>-</td><td>✓</td><td>✓</td><td>-</td><td>✓</td><td>✓</td><td>✓</td></tr><tr><td>L1448-IRS3B</td><td>0</td><td>✓</td><td>-</td><td>✓</td><td>✓</td><td>-</td><td>✓</td><td>✓</td><td>-</td></tr><tr><td>L1448C</td><td>0</td><td>✓</td><td>-</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td></tr><tr><td>L1448NW</td><td>0</td><td>✓</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>✓</td><td>✓</td></tr><tr><td>Per-emb-2</td><td>0</td><td>✓</td><td>-</td><td>✓</td><td>✓</td><td>-</td><td>-</td><td>✓</td><td>-</td></tr><tr><td>Per-emb-5</td><td>0</td><td>✓</td><td>-</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>-</td></tr><tr><td>Per-emb-8</td><td>0</td><td>✓</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>✓</td><td>-</td></tr><tr><td>Per-emb-18</td><td>0</td><td>✓</td><td>-</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td></tr><tr><td>Per-emb-22</td><td>0</td><td>✓</td><td>-</td><td>✓</td><td>✓</td><td>-</td><td>-</td><td>✓</td><td>✓</td></tr><tr><td>Per-emb-29</td><td>0</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td></tr><tr><td>Per-emb-30</td><td>0/I</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-(b)</td><td>-(b)</td></tr><tr><td>Per-emb-50</td><td>I</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>✓</td><td>✓</td></tr><tr><td>Per-emb-62</td><td>I</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-(b)</td><td>-(b)</td></tr><tr><td>SVS13A</td><td>0/I</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td></tr><tr><td>SVS13B</td><td>0</td><td>✓</td><td>-</td><td>✓</td><td>✓</td><td>-</td><td>✓</td><td>✓</td><td>-</td></tr></tbody></table>

Notes. (a) Retrieved from Artur de la Villarmois et al. (2023). (b) Objects not observed in Artur de la Villarmois et al. (2023).
<table><caption>Table D.3: Column densities and specific ratios in Class 0/I protostars in Perseus.</caption><thead><tr><th>Source</th><th>N(H<sub>2</sub>S) (cm<sup>-2</sup>)×10<sup>13</sup> (cm<sup>-2</sup>)×10<sup>13</sup></th><th>N(H<sub>2</sub><sup>33</sup>S) (cm<sup>-2</sup>)×10<sup>13</sup> (cm<sup>-2</sup>)×10<sup>13</sup></th><th>N(OCS) (cm<sup>-2</sup>)×10<sup>13</sup> (cm<sup>-2</sup>)×10<sup>13</sup></th><th>Calculated H<sub>2</sub>S<sup>(a)</sup> (cm<sup>-2</sup>)×10<sup>13</sup> (cm<sup>-2</sup>)×10<sup>13</sup></th><th>Calculated OCS<sup>(b)</sup> (cm<sup>-2</sup>)×10<sup>13</sup> (cm<sup>-2</sup>)×10<sup>13</sup></th><th>N(SO)<sup>(c)</sup> (cm<sup>-2</sup>)×10<sup>13</sup></th><th>N(SO<sub>2</sub>)<sup>(c)</sup> (cm<sup>-2</sup>)×10<sup>13</sup></th><th><sup>34</sup>S/<sup>33</sup>S H<sub>2</sub>S<sub>13</sub>/H<sub>2</sub>S<sup>(d)</sup> OCS<sub>13</sub>/OCS<sup>(d)</sup></th><th>H<sub>2</sub>S/OCS</th></tr></thead><tbody><tr><td>B1bN</td><td>0.0</td><td>2.890 (N(H<sub>2</sub>S)=154.0)</td><td>0.0</td><td>0.0</td><td>0.0</td><td>&lt;7</td><td>&lt;28</td><td>-</td><td>-</td><td>-</td></tr><tr><td>B1bS</td><td>18.45</td><td>6.055 (N(H<sub>2</sub>S)=2659)</td><td>13.37</td><td>5.099</td><td>12.83 (N(OCS)=288.8)</td><td>&lt;7</td><td>&lt;28</td><td>2.52</td><td>144.1</td><td>21.6</td><td>9.21</td></tr><tr><td>B5-IRS1</td><td>50.40</td><td>0.0</td><td>&lt;7.154</td><td>0.0</td><td>0.0</td><td>[41.3-61.2]</td><td>&lt;38</td><td>-</td><td>-</td><td>-</td><td>&gt;7.05</td></tr><tr><td>HH211MMS</td><td>21.01</td><td>0.0</td><td>9.871</td><td>1.831</td><td>7.095 (N(OCS)=159.6)</td><td>[48.3-71.6]</td><td>[50-200]</td><td>3.88</td><td>-</td><td>16.2</td><td>0.13</td></tr><tr><td>IC348MMS</td><td>33.96</td><td>4.092 (N(H<sub>2</sub>S)=2936)</td><td>17.01</td><td>1.738</td><td>8.735 (N(OCS)=196.5)</td><td>&gt;40.9</td><td>[20-60]</td><td>5.02</td><td>86.5</td><td>11.6</td><td>14.9</td></tr><tr><td>IRAS2A</td><td>446.3</td><td>&gt;20.05 (N(H<sub>2</sub>S)=5687)</td><td>565.7</td><td>13.59</td><td>81.98 (N(OCS)=1845)</td><td>&gt;188.1</td><td>&gt;300</td><td>6.03</td><td>&gt;12.7</td><td>3.26</td><td>&gt;3.08</td></tr><tr><td>IRAS2B</td><td>10.29</td><td>0.0</td><td>10.58</td><td>0.0</td><td>0.0</td><td>&gt;52.7</td><td>[10-40]</td><td>-</td><td>-</td><td>-</td><td>0.97</td></tr><tr><td>IRAS4B</td><td>266.6</td><td>31.93 (N(H<sub>2</sub>S)=20600)</td><td>145.9</td><td>15.20</td><td>63.89 (N(OCS)=1437)</td><td>&gt;11.9</td><td>[20-70]</td><td>4.20</td><td>77.3</td><td>9.85</td><td>14.3</td></tr><tr><td>IRAS4C</td><td>11.54</td><td>0.0</td><td>3.886</td><td>0.0</td><td>0.0</td><td>[33.2-49.2]</td><td>&lt;18</td><td>-</td><td>-</td><td>-</td><td>2.97</td></tr><tr><td>L1448-IRS3A</td><td>21.83</td><td>0.0</td><td>37.36</td><td>2.360</td><td>9.640 (N(OCS)=216.9)</td><td>&gt;35.5</td><td>[20-60]</td><td>4.08</td><td>-</td><td>5.81</td><td>0.10</td></tr><tr><td>L1448-IRS3B</td><td>27.05</td><td>0.0</td><td>11.00</td><td>0.804</td><td>7.152 (N(OCS)=160.9)</td><td>[10.8-16.0]</td><td>&lt;15</td><td>8.90</td><td>-</td><td>14.6</td><td>0.17</td></tr><tr><td>L1448C</td><td>130.9</td><td>0.0</td><td>222.0</td><td>2.965</td><td>16.99 (N(OCS)=382.3)</td><td>&gt;59.0</td><td>[60-400]</td><td>5.73</td><td>-</td><td>1.70</td><td>0.34</td></tr><tr><td>L1448NW</td><td>9.612</td><td>0.0</td><td>&lt;3.426</td><td>0.0</td><td>0.0</td><td>&gt;33.4</td><td>[10-50]</td><td>-</td><td>-</td><td>-</td><td>&gt;2.81</td></tr><tr><td>Per-emb-02</td><td>30.04</td><td>0.0</td><td>2.785</td><td>0.0</td><td>0.0</td><td>[19.4-28.8]</td><td>&lt;49</td><td>-</td><td>-</td><td>-</td><td>10.8</td></tr><tr><td>Per-emb-05</td><td>14.83</td><td>0.0</td><td>12.84</td><td>0.695</td><td>8.509 (N(OCS)=191.5)</td><td>[13.7-20.2]</td><td>&lt;20</td><td>12.2</td><td>-</td><td>14.9</td><td>0.08</td></tr><tr><td>Per-emb-08</td><td>11.64</td><td>0.0</td><td>&lt;2.923</td><td>0.0</td><td>0.0</td><td>[17.6-26.1]</td><td>&lt;40</td><td>-</td><td>-</td><td>-</td><td>&gt;3.98</td></tr><tr><td>Per-emb-18</td><td>27.68</td><td>0.0</td><td>74.03</td><td>2.348</td><td>16.20 (N(OCS)=364.5)</td><td>&gt;65.2</td><td>[30-100]</td><td>6.90</td><td>-</td><td>4.92</td><td>0.08</td></tr><tr><td>Per-emb-22</td><td>37.87</td><td>0.0</td><td>23.66</td><td>0.0</td><td>0.0</td><td>[16.8-51.8]</td><td>[10-40]</td><td>-</td><td>-</td><td>-</td><td>1.60</td></tr><tr><td>Per-emb-29</td><td>98.54</td><td>4.68 (N(H<sub>2</sub>S)=5364)</td><td>36.83</td><td>5.763</td><td>20.04 (N(OCS)=451.0)</td><td>&gt;56.3</td><td>[30-100]</td><td>3.48</td><td>54.4</td><td>12.2</td><td>11.9</td></tr><tr><td>Per-emb-30</td><td>0.0</td><td>0.0</td><td>0.0</td><td>0.0</td><td>0.0</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>Per-emb-50</td><td>0.0</td><td>0.0</td><td>0.0</td><td>0.0</td><td>0.0</td><td>&gt;127.6</td><td>[50-300]</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>Per-emb-62</td><td>0.0</td><td>0.0</td><td>0.0</td><td>0.0</td><td>0.0</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>SVS13A</td><td>974.7</td><td>&gt;67.76 (N(H<sub>2</sub>S)=14820)</td><td>946.6</td><td>39.10</td><td>192.4 (N(OCS)=4329)</td><td>&gt;210.5</td><td>&gt;300</td><td>4.92</td><td>&gt;15.2</td><td>4.57</td><td>&gt;3.42</td></tr><tr><td>SVS13B</td><td>26.76</td><td>0.0</td><td>31.59</td><td>1.964</td><td>6.163 (N(OCS)=138.7)</td><td>[5.0-7.4]</td><td>&lt;18</td><td>3.14</td><td>-</td><td>4.39</td><td>0.19</td></tr></tbody></table>

**Notes.** (a) Calculated using the solar isotopic ratio, $^{32}$S/$^{33}$S=126 (Anders & Grevesse 1989). (b) Calculated using the solar isotopic ratio, $^{32}$S/$^{34}$S=22.5 (Anders & Grevesse 1989). (c) Retrieved from Artur de la Villarmois et al. (2023); no data from Per-emb-30 and Per-emb-62 (not observed). (d) Al, A notation stands for the ratio between the A species column density calculated with the corresponding $^{34}$S isotopologue, and the column density of A estimated without the isotopologue.
<table><caption>Table D.4: Abundances in Class 0/I protostars from Perseus.</caption><thead><tr><th>Source</th><th>H<sub>2</sub>S<br>×10<sup>-9</sup></th><th>H<sub>2</sub>S<br>×10<sup>-9</sup></th><th>OCS<br>×10<sup>-9</sup></th><th>OC<sup>33</sup>S<br>×10<sup>-9</sup></th><th>OC<sup>34</sup>S<br>×10<sup>-9</sup></th><th>SO<sup>(a)</sup><br>×10<sup>-9</sup></th><th>SO<sub>2</sub><sup>(a)</sup><br>×10<sup>-9</sup></th><th>TOTAL<sup>(b)</sup><br>×10<sup>-9</sup></th><th>D<sub>S</sub></th></tr></thead><tbody><tr><td>B1bN</td><td>1.442</td><td>0.025</td><td>0.0</td><td>0.0</td><td>0.0</td><td>&lt;0.007</td><td>&lt;0.029</td><td>1.464</td><td>1×10<sup>4</sup></td></tr><tr><td>B1bS</td><td>16.03</td><td>0.036</td><td>1.741</td><td>0.031</td><td>0.077</td><td>&lt;0.005</td><td>&lt;0.019</td><td>17.78</td><td>8×10<sup>2</sup></td></tr><tr><td>B5-IRS1</td><td>89.98</td><td>0.0</td><td>12.77</td><td>0.0</td><td>0.0</td><td>[0.307-0.950]</td><td>&lt;0.590</td><td>103.8</td><td>1×10<sup>2</sup></td></tr><tr><td>HF211MMS</td><td>0.586</td><td>0.0</td><td>4.453</td><td>0.051</td><td>0.198</td><td>[0.034-0.100]</td><td>[0.035-0.279]</td><td>5.263</td><td>3×10<sup>3</sup></td></tr><tr><td>IC348MMS</td><td>41.88</td><td>0.083</td><td>2.803</td><td>0.025</td><td>0.125</td><td>&gt;0.011</td><td>[0.005-0.030]</td><td>44.72</td><td>3×10<sup>2</sup></td></tr><tr><td>IRAS2A</td><td>&gt;207.1</td><td>&gt;0.411</td><td>67.16</td><td>0.495</td><td>2.985</td><td>&gt;0.608</td><td>&gt;0.970</td><td>&gt;276.3</td><td>&lt;5×10<sup>1</sup></td></tr><tr><td>IRAS2B</td><td>1.142</td><td>0.0</td><td>1.173</td><td>0.0</td><td>0.0</td><td>&gt;0.367</td><td>[0.070-0.828]</td><td>3.642</td><td>4×10<sup>3</sup></td></tr><tr><td>IRAS4B</td><td>84.13</td><td>0.167</td><td>5.871</td><td>0.062</td><td>0.261</td><td>&gt;0.013</td><td>[0.021-0.196]</td><td>90.12</td><td>2×10<sup>2</sup></td></tr><tr><td>IRAS4C</td><td>0.379</td><td>0.0</td><td>&lt;0.128</td><td>0.0</td><td>0.0</td><td>[0.059-0.179]</td><td>&lt;0.066</td><td>0.688</td><td>2×10<sup>4</sup></td></tr><tr><td>L1448-IRS3A</td><td>1.383</td><td>0.0</td><td>13.74</td><td>0.149</td><td>0.611</td><td>&gt;0.198</td><td>[0.111-1.003]</td><td>16.12</td><td>9×10<sup>2</sup></td></tr><tr><td>L1448-IRS3B</td><td>0.537</td><td>0.0</td><td>3.195</td><td>0.016</td><td>0.142</td><td>[0.015-0.048]</td><td>&lt;0.045</td><td>3.793</td><td>4×10<sup>3</sup></td></tr><tr><td>L1448C</td><td>4.113</td><td>0.0</td><td>12.01</td><td>0.093</td><td>0.534</td><td>&gt;0.198</td><td>[0.201-3.133]</td><td>18.12</td><td>8×10<sup>2</sup></td></tr><tr><td>L1448NW</td><td>0.275</td><td>0.0</td><td>&lt;0.097</td><td>0.0</td><td>0.0</td><td>&gt;0.045</td><td>[0.013-0.120]</td><td>0.518</td><td>3×10<sup>4</sup></td></tr><tr><td>Per-emb-02</td><td>0.240</td><td>0.0</td><td>0.022</td><td>0.0</td><td>0.0</td><td>[0.015-0.040]</td><td>&lt;0.068</td><td>0.330</td><td>5×10<sup>4</sup></td></tr><tr><td>Per-emb-05</td><td>0.206</td><td>0.0</td><td>2.658</td><td>0.010</td><td>0.118</td><td>[0.031-0.091]</td><td>&lt;0.090</td><td>2.986</td><td>5×10<sup>3</sup></td></tr><tr><td>Per-emb-08</td><td>0.501</td><td>0.0</td><td>&lt;0.126</td><td>0.0</td><td>0.0</td><td>[0.016-0.042]</td><td>&lt;0.064</td><td>0.697</td><td>2×10<sup>4</sup></td></tr><tr><td>Per-emb-18</td><td>1.271</td><td>0.0</td><td>16.73</td><td>0.108</td><td>0.744</td><td>&gt;0.201</td><td>[0.093-0.681]</td><td>18.76</td><td>8×10<sup>2</sup></td></tr><tr><td>Per-emb-22</td><td>2.683</td><td>0.0</td><td>1.676</td><td>0.0</td><td>0.0</td><td>[0.024-0.164]</td><td>[0.014-0.173]</td><td>4.527</td><td>3×10<sup>3</sup></td></tr><tr><td>Per-emb-29</td><td>65.55</td><td>0.209</td><td>20.18</td><td>0.258</td><td>0.897</td><td>&gt;0.063</td><td>[0.033-0.744]</td><td>85.93</td><td>6×10<sup>1</sup></td></tr><tr><td>Per-emb-30</td><td>0.0</td><td>0.0</td><td>0.0</td><td>0.0</td><td>0.0</td><td>-</td><td>-</td><td>0.0</td><td>-</td></tr><tr><td>Per-emb-50</td><td>0.0</td><td>0.0</td><td>0.0</td><td>0.0</td><td>0.0</td><td>&gt;1.108</td><td>[0.434-7.437]</td><td>&gt;6.408</td><td>&lt;2×10<sup>3</sup></td></tr><tr><td>Per-emb-62</td><td>0.0</td><td>0.0</td><td>0.0</td><td>0.0</td><td>0.0</td><td>-</td><td>-</td><td>0.0</td><td>-</td></tr><tr><td>SVS13A</td><td>&gt;1343</td><td>&gt;6.139</td><td>392.2</td><td>3.542</td><td>17.43</td><td>&gt;1.885</td><td>&gt;2.686</td><td>&gt;1742</td><td>&lt;8.6</td></tr><tr><td>SVS13B</td><td>0.232</td><td>0.0</td><td>1.204</td><td>0.017</td><td>0.053</td><td>[0.019-0.078]</td><td>&lt;0.190</td><td>1.589</td><td>9×10<sup>3</sup></td></tr></tbody></table>

Notes. (a) Retrieved from Artur de la Villarmois et al. (2023); no data from Per-emb-30 and Per-emb-62 (not observed). (b) Total abundance has been calculated using all lower and upper limits, and using the mean value of SO and SO₂ in the cases a range is given. Lower and upper limits to the total sulfur abundance have only been explicit when a limit was given for the most abundant species in the source.
## Appendix E: Complementary figures

In Figs. E.1 and E.2, we show the moment-0 maps and the spectra, respectively, of the 20 sources not shown in the main text. The moment-0 maps cover the H₂S and both OCS lines, while in the spectra we display the emission of the H₂S line, one of the OCS lines and the OC³S line.

<figure><img src="image_17.png" /></figure>

Fig. E.1: Emission of the main species (H₂S, and both OCS lines) in the warm inner core of the 20 protostars from our sample not shown in Section 3. The color map represents the >3σ emission integrated images in a 4"×4" square region. The color scale, shown at the right of each map, is the brightness temperature in K. The white contours represent 10%, 30%, 50%, 70% and 90% of the peak temperature. The pink star shows the position of the protostar, determined by the position of the maximum emission in the continuum (Tobin et al. 2016b). The green cross marks the point with maximum emission of the line.
<figure><img src="image_18.png" /><figcaption>Fig. E.1: (Continuation)</figcaption></figure>
<figure><img src="image_19.png" /></figure>

Fig. E.2: Spectra of the H~2~S (2~2~,0~−~2~1,1), OCS (18~−~17), and OC~34~S (20~−~19) lines in the warm inner core of the 20 protostars from our sample not shown in Section 3. The red dotted lines represent the values of a Gaussian fit to each of the detected lines. In some cases, the combination of two Gaussian fits was necessary for the fit. The OCS (18-17) line was observed with a spectral resolution of 62.5 kHz while the H~2~S and OC~34~S lines were observed with a spectral resolution of 2 MHz. Gaussian fits have been plotted for the >3σ detections, except for the OC~34~S were we also show the <3σ fits when the OC~33~S counterpart was detected with >3σ.
<figure><img src="image_20.png" /><figcaption>Fig. E.2: (Continuation)</figcaption></figure>