Gas phase Elemental abundances in Molecular cloudS
(GEMS)

A. Fuente, P. Rivière-Marichalar, L. Beitia-Antero, P. Caselli, V. Wakelam, G. Esplugues, M. Rodríguez-Baras, D. Navarro-Almaidá, M. Gerin, C. Kramer, et al.

► To cite this version:

A. Fuente, P. Rivière-Marichalar, L. Beitia-Antero, P. Caselli, V. Wakelam, et al. Gas phase Elemental abundances in Molecular cloudS (GEMS). Astronomy & Astrophysics - A&A, 2023, 670, pp.A114. 10.1051/0004-6361/20224843. hal-04292075

HAL Id: hal-04292075

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

Submitted on 2 Dec 2023

**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 CC BY 4.0 - Attribution - International License
# Gas phase Elemental abundances in Molecular cloudS (GEMS)

## VII. Sulfur elemental abundance

A. Fuente<fnref n="1" />, P. Rivière-Marichalar<fnref n="1" />, L. Beitia-Antero<fnref n="1" />,<fnref n="2" />, P. Caselli<fnref n="3" />, V. Wakelam<fnref n="4" />, G. Esplugues<fnref n="1" />,
M. Rodríguez-Baras<fnref n="1" />, D. Navarro-Almaidá<fnref n="5" />, M. Gerin<fnref n="6" />,<fnref n="7" />, C. Kramer<fnref n="8" />,<fnref n="9" />, R. Bachiller<fnref n="1" />, J. R. Goicoechea<fnref n="10" />,
I. Jiménez-Serra<fnref n="11" />, J. C. Loison<fnref n="12" />, A. Ivlev<fnref n="3" />, R. Martín-Doménech<fnref n="13" />, S. Spezzano<fnref n="3" />, O. Roncero<fnref n="10" />,
G. Muñoz-Caro<fnref n="11" />, S. Cazaux<fnref n="14" />,<fnref n="15" />, and N. Marcelino<fnref n="1" />

Received 30 August 2022 / Accepted 2 December 2022

## ABSTRACT

**Context.** Gas phase Elemental abundances in molecular CloudS (GEMS) is an IRAM 30-m Large Program aimed at determining the elemental abundances of carbon (C), oxygen (O), nitrogen (N), and sulfur (S) in a selected set of prototypical star-forming filaments. In particular, the elemental abundance of S remains uncertain by several orders of magnitude, and its determination is one of the most challenging goals of this program.

**Aims.** This paper aims to constrain the sulfur elemental abundance in Taurus, Perseus, and Orion A based on the GEMS molecular database. The selected regions are prototypes of low-mass, intermediate-mass, and high-mass star-forming regions, respectively, providing useful templates for the study of interstellar chemistry.

**Methods.** We have carried out an extensive chemical modeling of the fractional abundances of CO, HCO+, HCN, HNC, CS, SO, H₂S, OCS, and HCS⁺ to determine the sulfur depletion toward the 244 positions in the GEMS database. These positions sample visual extinctions from $A_V$ ~ 3 mag to >50 mag, molecular hydrogen densities ranging from a few × 10³ cm⁻<fnref n="3" /> to 3 × 10⁶ cm⁻<fnref n="3" />, and $T_k$ ~ 10–35 K. We investigate the possible relationship between sulfur depletion and the grain charge distribution in different environments.

**Results.** Most of the positions in Taurus and Perseus are best fitted assuming early-time chemistry, $t = 0.1$ Myr, $\zeta_{H_2}$ ~ (0.5–1) × 10⁻¹⁶ s⁻<fnref n="1" />, and [S/H] ~ 1.5 × 10⁻<fnref n="6" />. On the contrary, most of the positions in Orion are fitted with $t = 1$ Myr and $\zeta_{H_2}$ ~ 10⁻¹⁷ s⁻<fnref n="1" />. Moreover, ~40% of the positions in Orion are best fitted assuming the undepleted sulfur abundance, [S/H] ~ 1.5 × 10⁻<fnref n="5" />. We find a tentative trend of sulfur depletion increasing with density.

**Conclusions.** Our results suggest that sulfur depletion depends on the environment. While the abundances of sulfur-bearing species are consistent with undepleted sulfur in Orion, a depletion factor of ~20 is required to explain those observed in Taurus and Perseus. We propose that differences in the grain charge distribution might explain these variations. Grains become negatively charged at a visual extinction of $A_V$ ~ 3.5 mag in Taurus and Perseus. At this low visual extinction, the S⁺ abundance is high, $X(S⁺) > 10⁻<fnref n="6" />$, and the electrostatic attraction between S⁺ and negatively charged grains could contribute to enhance sulfur depletion. In Orion, the net charge of grains remains approximately zero until higher visual extinctions ($A_V$ ~ 5.5 mag), where the abundance of S⁺ is already low because of the higher densities, thus reducing sulfur accretion. The shocks associated with past and ongoing star formation could also contribute to enhance [S/H].

**Key words.** astrochemistry – ISM: abundances – ISM: kinematics and dynamics – ISM: molecules – stars: formation – stars: low-mass

<footnote n="1">¹ Observatorio Astronómico Nacional (OAN), Alfonso XII 3, 28014 Madrid, Spain
e-mail: a.fuente@oan.es</footnote>

<footnote n="2">² Departamento de Física de la Tierra y Astrofísica, Facultad de CC. Matemáticas, Universidad Complutense de Madrid, Plaza de las Ciencias 3, 28040 Madrid, Spain</footnote>

<footnote n="3">³ Centre for Astrochemical Studies, Max-Planck-Institute for Extraterrestrial Physics, Giessenbachstrasse 1, 85748 Garching, Germany</footnote>

<footnote n="4">⁴ Laboratoire d'astrophysique de Bordeaux, Univ. Bordeaux, CNRS, B18N, allée Geoffroy Saint-Hilaire, 33615 Pessac, France</footnote>

<footnote n="5">⁵ Université Paris-Saclay, CEA, AIM, Département d'Astrophysique (DAp), 91191 Gif-sur-Yvette, France</footnote>

<footnote n="6">⁶ Observatoire de Paris, PSL Research University, CNRS, École Normale Supérieure, Sorbonne Université, UPMC Univ. Paris 6, 75005 Paris, France</footnote>

<footnote n="7">⁷ LERMA, Observatoire de Paris, PSL Research University, CNRS, UMR8112, Place Janssen, 92190, Meudon Cedex, France</footnote>

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

<footnote n="9">⁹ Instituto de Radioastronomía Milimétrica (IRAM), Av. Divina Pastora 7, Nucleo Central, 18012 Granada, Spain</footnote>

<footnote n="10">¹⁰ Instituto de Física Fundamental (CSIC), Calle Serrano 121-123, 28006 Madrid, Spain</footnote>

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

<footnote n="12">¹² Institut des Sciences Moléculaires (ISM), CNRS, Univ. Bordeaux, 351 cours de la Libération, 33400 Talence, France</footnote>

<footnote n="13">¹³ Center for Astrophysics, Harvard & Smithsonian, 60 Garden St., Cambridge, MA 02138, USA</footnote>

<footnote n="14">¹⁴ Faculty of Aerospace Engineering, Delft University of Technology, Delft, The Netherlands</footnote>

<footnote n="15">¹⁵ University of Leiden, PO Box 9513, 2300 RA Leiden, The Netherlands</footnote>
# 1. Introduction

In recent years, space telescopes such as Spitzer and Herschel have revolutionized our view of star-forming regions. Images of giant molecular clouds and dark cloud complexes have revealed spectacular networks of filamentary structures where stars are born (André et al. 2010). Interstellar filaments are almost everywhere in the Milky Way. Now we believe that filaments precede the onset of star formation, funneling interstellar gas and dust into increasingly denser concentrations that will contract and fragment leading to gravitationally bound prestellar cores that will eventually form stars.

Gas chemistry has a key role in the star formation process by determining aspects such as the gas cooling and the gas ionization degree. Molecular filaments can fragment to prestellar cores to a large extent because molecules cool the gas, thus diminishing the thermal support relative to self-gravity. The ionization fraction controls the coupling of magnetic fields with the gas driving the dissipation of turbulence and angular momentum transfer. Therefore, chemistry plays a crucial role in the cloud collapse (isolated vs. clustered star formation) and the dynamics of accretion disks (see Zhao et al. 2016; Padovani et al. 2013). In the absence of other ionization agents (X-rays, UV photons, and J-type shocks), the ionization fraction is proportional to $\sqrt{\zeta_{H_2}}$, where $\zeta_{H_2}$ is the cosmic-ray ionization rate for H$_2$ molecules, which becomes a key parameter in the molecular cloud evolution (McKee 1989; Caselli et al. 2002). In addition to $\zeta_{H_2}$, the gas ionization fraction, $X(e^-)$, depends on the elemental depletion factors (Caselli et al. 1998). In particular, carbon (C) is the main donor of electrons in the cloud surface ($A_V$ < 4 mag). As long as it is not heavily depleted, sulfur (S) is the main donor in the range of ~3.7–7 magnitudes, which encompasses a large fraction of the molecular cloud mass. Depletions of C and O determine the cooling gas rate since CO, [CII], and [OI] are main coolants in molecular clouds.

Gas phase Elemental abundances in Molecular CloudS (GEMS) is an IRAM 30-m Large Program aimed at determining the S, C, N, and O depletions, and $X(e^-)$ in a set of selected prototypical star-forming filaments. The observations corresponding to this program and the resulting molecular database were presented by Rodríguez-Baras et al. (2021). In that paper, we explored the relationship between the abundances of the different molecular species and the local physical parameters, concluding that density is the main parameter determining the chemical abundances in dark clouds. The determination of the elemental abundances is a challenging task, especially in the case of sulfur for which the main potential sulfur reservoirs (atomic S, and H$_2$S ice) cannot be easily observed. The fine-structure line emission of atomic sulfur at ~25 µm has been detected in bipolar outflows using the Spitzer Space Telescope (Anderson et al. 2013). However, the high excitation conditions of this line ($E_u$ = 570 K, $n_{crit}$ = $1.5 \times 10^6 \text{ cm}^{-3}$) prohibits its detection in dark clouds. Alternatively, the S$^+$ abundance can be derived from the emission of the sulfur radio recombination lines (Pankonin & Walmley 1978; Goicoechea & Cuadrado 2021). However, these lines are very weak, thus hindering their use as routine S$^+$ tracers in molecular clouds. The detection of H$_2$S in the ice is hampered by the strong overlap between the 2558 cm$^{-1}$ band and the methanol bands at 2530 and 2610 cm$^{-1}$, and upper limits have been derived thus far (Smith 1991; Jiménez-Escobar & Muñoz Caro 2011). Therefore, the estimate of sulfur elemental abundance in dark clouds must rely on the observation of minor sulfur-bearing species (e.g., CS, and SO) and the predictive power of chemical models. The GEMS project relies on a collaborative effort

among astronomers, modelers, theoretical chemists, and experimentalists, working in a coordinated way to improve the sulfur chemical network and, therefore, the reliability of model predictions. Theoretical ab initio calculations have been carried out to obtain more accurate estimations of the rates of important reactions associated with the SO and CS chemistry (Fuente et al. 2019; Bulut et al. 2021, and in prep.). An important upgrade of the sulfur surface chemistry was carried out by Laas & Caselli (2019). Navarro-Almaida et al. (2020, 2021) used a progressively updated chemical network to investigate the sulfur chemistry in the prototypes of cold dense cores, TMC 1 (CP), TMC 1 (C), and Barnard 1b. Goicoechea et al. (2021) also used an updated network to account for the chemistry of sulfur hydrides in the Orion Bar. Based on the GEMS database, Spezzano et al. (2022) and Esplugues et al. (2022) carried out a detailed study of the CH$_3$OH and H$_2$CS chemistries, respectively. These two species are crucial to constrain the formation of complex organic molecules (COMs) and organo-sulfur compounds in starless cores.

Elemental abundance, as referred in this paper, is the amount of a given atom in volatiles, that is to say in the gas and in the icy grain mantles. The elemental abundances are given as initial condition to most chemical models and remain constant during the calculations, since the atoms are not expected to incorporate to the refractory part of grains in molecular clouds. Usually, the adopted initial abundances of C, N, and O correspond to the values derived in the diffuse cloud $\zeta$ Oph by Jenkins (2009), which result in a C/O abundance ratio of 0.55, giving reasonable predictions for dark clouds (Agúndez & Wakelam 2013). However, the elemental abundance of S is uncertain by several orders of magnitude. While the observed gaseous sulfur accounts for its total solar abundance ([S/H] ~ $1.5 \times 10^{-5}$, Asplund et al. 2009) in diffuse clouds (Neufeld et al. 2015) and highly irradiated photon-dominated regions (Goicoechea & Cuadrado 2021), a sulfur depletion of a factor of >100 is usually needed to account for the observed abundances of sulfur-bearing species in starless cores (Ruffle et al. 1999; Vastel et al. 2018) and proto-planetary disks (Dutrey et al. 2011; Le Gal et al. 2019; Rivière-Marichalar et al. 2020; Le Gal et al. 2021). Sulfur depletions of a factor of ~10 have been estimated in the translucent part of molecular clouds (Fuente et al. 2019), hot cores (Blake et al. 1996; Wang et al. 2013; Fuente et al. 2021; Bouscasse et al. 2022), low irradiated photon-dominated regions (Goicoechea et al. 2006), and bipolar outflows (Anderson et al. 2013; Feng et al. 2020; Taquet et al. 2020). Vidal et al. (2017) were able to reproduce the abundances of the sulfur bearing species detected in TMC 1-CP with undepleted sulfur abundance. Thus far, there is no unified scheme accounting for the observed abundances of sulfur species in different regions of the interstellar medium. The determination of the sulfur depletion is the objective of this paper. The large molecular GEMS database allows us to carry out a uniform and systematic study in regions located in different environments, avoiding the uncertainties coming from the use of different observations and methodology.

# 2. Observations and molecular database

## 2.1. Observations

GEMS observations were carried out using the IRAM 30-m telescope at Pico Veleta (Spain) during several observing periods from July 2017 to December 2020. In addition, complementary observations of the CS 1→0 line were carried
<table><caption>**Table 1.** Cores included in the GEMS sample and observation cuts associated with them.</caption><thead><tr><th>Cloud complex</th><th>Cloud ID</th><th>Core RA (J2000)</th><th>Coordinates Dec (J2000)</th><th>Other names</th><th>Cut</th><th>N<sub>o</sub></th><th>N<sub>v</sub></th></tr></thead><tbody><tr><td rowspan="14">Taurus</td><td rowspan="3">TMC 1</td><td>CP</td><td>04:41:41.90</td><td>+25:41:27.1</td><td></td><td>C1</td><td>6</td><td>6</td></tr><tr><td>NH3</td><td>04:41:21.30</td><td>+25:48:07.0</td><td></td><td>C2</td><td>6</td><td>6</td></tr><tr><td>C</td><td>04:41:38.80</td><td>+25:59:42.0</td><td></td><td>C3</td><td>6</td><td>6</td></tr><tr><td rowspan="5">B 213-N<sup>(1)</sup></td><td>#1</td><td>04:17:41.80</td><td>+28:08:47.0</td><td>5</td><td>C1</td><td>9</td><td>9</td></tr><tr><td>#2</td><td>04:17:50.60</td><td>+27:56:01.0</td><td>—</td><td>C2</td><td>9</td><td>7</td></tr><tr><td>#5</td><td>04:18:03.80</td><td>+28:23:06.0</td><td>7</td><td>C5</td><td>9</td><td>9</td></tr><tr><td>#6</td><td>04:18:08.40</td><td>+28:05:12.0</td><td>8</td><td>C6</td><td>9</td><td>5</td></tr><tr><td>#7</td><td>04:18:11.50</td><td>+27:35:15.0</td><td>9</td><td>C7</td><td>9</td><td>8</td></tr><tr><td rowspan="4">B 213-S<sup>(1)</sup></td><td>#10</td><td>04:19:37.60</td><td>+27:15:31.0</td><td>13a</td><td>C10</td><td>9</td><td>8</td></tr><tr><td>#12</td><td>04:19:51.70</td><td>+27:11:33.0</td><td>—</td><td>C12</td><td>9</td><td>6</td></tr><tr><td>#16</td><td>04:21:21.00</td><td>+27:00:09.0</td><td>—</td><td>C16</td><td>9</td><td>9</td></tr><tr><td>#17</td><td>04:27:54.00</td><td>+26:17:50.0</td><td>26b</td><td>C17</td><td>9</td><td>5</td></tr><tr><td rowspan="16">Perseus</td><td rowspan="2">L1448<sup>(2)</sup></td><td>#32</td><td>03:25:49.00</td><td>+30:42:24.6</td><td></td><td>C1</td><td>8</td><td>7</td></tr><tr><td></td><td>03:29:18.26</td><td>+31:28:02.0</td><td></td><td>C1</td><td>21</td><td>18</td></tr><tr><td rowspan="2">NGC 1333<sup>(2)</sup></td><td></td><td>03:28:41.60</td><td>+31:06:02.0</td><td></td><td>C2</td><td>17</td><td>13</td></tr><tr><td></td><td>03:29:11.00</td><td>+31:18:27.4</td><td>SK20<sup>(3)</sup></td><td>C3</td><td>17</td><td>11</td></tr><tr><td></td><td>03:28:39.40</td><td>+31:18:27.4</td><td></td><td></td><td></td><td></td></tr><tr><td></td><td>03:29:08.80</td><td>+31:15:18.1</td><td>SK16</td><td>C4</td><td>16</td><td>13</td></tr><tr><td></td><td>03:29:04.50</td><td>+31:20:59.1</td><td>SK26</td><td>C5</td><td>11</td><td>11</td></tr><tr><td></td><td>03:29:18.20</td><td>+31:25:10.8</td><td>SK33</td><td>C6</td><td>9</td><td>9</td></tr><tr><td></td><td>03:29:25.50</td><td>+31:28:18.1</td><td></td><td>C7</td><td>1</td><td>1</td></tr><tr><td rowspan="2">Barnard 1<sup>(2)</sup></td><td>1b</td><td>03:33:20.80</td><td>+31:07:34.0</td><td></td><td>C1</td><td>18</td><td>18</td></tr><tr><td></td><td>03:33:01.90</td><td>+31:04:23.2</td><td></td><td>C2</td><td>8</td><td>6</td></tr><tr><td rowspan="3">Barnard 5<sup>(2)</sup></td><td>#79</td><td>03:47:38.99</td><td>+32:52:15.0</td><td></td><td>C1</td><td>13</td><td>9</td></tr><tr><td>#1</td><td>03:44:01.00</td><td>+32:01:54.8</td><td></td><td>C1</td><td>14</td><td>11</td></tr><tr><td>#10</td><td>03:44:05.74</td><td>+32:01:53.5</td><td></td><td></td><td></td><td></td></tr><tr><td rowspan="3">Orion</td><td rowspan="3">Orion A</td><td>OMC 3</td><td>05:35:19.54</td><td>-05:00:41.5</td><td></td><td>C1</td><td>20</td><td>10</td></tr><tr><td>OMC 4</td><td>05:35:08.15</td><td>-05:35:41.5</td><td></td><td>C2</td><td>20</td><td>14</td></tr><tr><td>OMC 2</td><td>05:35:23.68</td><td>-05:12:31.8</td><td></td><td>C3</td><td>13</td><td>9</td></tr><tr><td colspan="4">TOTAL</td><td></td><td>305</td><td>244</td><td></td></tr></tbody></table>

**Notes.** $N_o$ indicates the total number of points observed in the corresponding cut. $N_v$ indicates the number of points where the molecular hydrogen density could be derived. $^{(1)}$B 213 core IDs are from Hacar et al. (2013). IDs indicated in “Other names” column are from Onishi et al. (2002). $^{(2)}$Perseus core IDs (L 1448, NGC 1333, Barnard 1, Barnard 5, IC348) are from Hatchell et al. (2007b). $^{(3)}$NGC 1333 core IDs indicated in “Other names” column are from Sandell & Knee (2001).

out with the Yebes 40-m radiotelescope (Tercero et al. 2021) toward TMC 1 and Barnard 1b. All these observations have been described in detail in the first GEMS paper (Fuente et al. 2019). The complete database is public and available through the IRAM Large Program webpage<fnref n="1" />, or directly from the GEMS webpage<fnref n="2" />.

## 2.2. Sample

Our project focuses on the nearby star-forming regions Taurus, Perseus, and Orion, which are considered as prototypes of low-mass, intermediate-mass, and massive star-forming regions. Moreover, these molecular cloud complexes have been observed with Herschel and SCUBA as part of the Gould Belt Survey (André et al. 2010), and accurate visual extinction ($A_V$) and dust temperature ($T_d$) maps are available (Malinen et al. 2012; Palmeirim et al. 2013; Lombardi et al. 2014; Zari et al. 2016). The angular resolution of the $A_V$ and $T_d$ maps (~36′′) is similar

to that provided by the 30-m telescope at 3 mm allowing a direct comparison of continuum and spectroscopic data.

Chemical calculations show that the determination of elemental abundances requires the observation of more than a dozen of molecular compounds, some of them presenting weak emission (Fuente et al. 2016, 2019). Full-sampling mapping of these regions in such a large number of species would have been unrealistic in terms of telescope observing time. An alternative strategy was adopted in GEMS. Our project observed toward 305 positions distributed in 29 cuts roughly perpendicular to a set of selected filaments. These cuts were defined to intersect the filament along one of the embedded starless cores, avoiding the position of protostars, HII regions and bipolar outflows. The final set of observed positions probes visual extinctions from $A_V$ ~3 mag up to $A_V$ ~ 200 mag. The separation between one position and another in a given cut was selected to sample regular bins of $A_V$. This is, however, not possible close to the visual extinction peaks where the surface density gradient is steep. In Table 1, we show a summary of the cuts forming our sample. A more detailed description of the cuts and their location within the clouds was given by Rodríguez-Baras et al. (2021).

<footnote n="1">¹ https://legacy.iram-institute.org/EN/</footnote>

<footnote n="2">² https://www.oan.es/gems/</footnote>
## 2.3. Molecular database

Carbon sulfide (CS) has been largely used as density and column density tracer in the interstellar medium because it is an abundant species, with a simple rotational spectrum, and its excitation is well-known (Linke & Goldsmith 1980; Zhou et al. 1989; Tatematsu et al. 1993; Zinchenko et al. 1995; Anglada et al. 1996; Bronfman et al. 1996; Launhardt et al. 1998; Shirley et al. 2003; Bayet et al. 2009; Wu et al. 2010; Zhang et al. 2014; Scourfield et al. 2020). In order to calculate the gas density toward the observed positions, we fit the rotational lines of CS and its isotopologues using the non-LTE molecular excitation and radiative transfer code RADEX (van der Tak et al. 2007), and the collisional coefficients calculated by Denis-Alpizar et al. (2018).

Several assumptions have been made in these calculations. First, we adopted fixed isotopic ratios, $^{12}\text{C}/^{13}\text{C} = 60$, $^{32}\text{S}/^{34}\text{S} = 22.5$, and $^{16}\text{O}/^{18}\text{O} = 550$ (Wilson & Rood 1994; Savage et al. 2002; Gratier et al. 2016). The assumption of fixed isotopic ratios is justified in the case of CS because isotopic fractionation is not expected to be important for sulfur. The CS molecule is enriched in $^{34}\text{S}$ at early time because of the $^{34}\text{S}^+$ + CS reaction but not at the characteristic chemical ages of dense clouds (Loison et al. 2019). More controversial could be the $^{12}\text{CS}/^{13}\text{CS}$ ratio; the adopted value is consistent with the results of Gratier et al. (2016) in TMC 1 and Agúndez et al. (2019) in L 483. This value is also consistent with chemical predictions for typical conditions and chemical ages in dark clouds (Colzi et al. 2020; Loison et al. 2020). Isotopic chemical fractionation might be important for HCN (Loison et al. 2020), but the HCN/H$^{13}\text{CN}$ does not differ from the assumed value in more than a factor of ~2. Oxygen isotopic fractionation is expected in some species such as NO, SO, $\text{O}_2$, and $\text{SO}_2$ (Loison et al. 2019). However, it is not expected to be relevant for $\text{C}^{18}\text{O}$ and $\text{HC}^{18}\text{O}^+$, the two oxygen isotopologues in the GEMS database. We assumed a beam filling factor of 1 for all transitions, that is to say that the emission is more extended than the beam size, which is reasonable for most positions.

The observed CS, $\text{C}^{34}\text{S}$, and $^{13}\text{CS}$ line intensities were fitted by varying only two parameters, $n(\text{H}_2)$ and $N(\text{CS})$. The gas kinetic temperature was assumed to be equal to the dust temperature as derived from the Herschel data. The parameter space was explored following the Monte Carlo Markov chain (MCMC) methodology with a Bayesian inference approach. In particular, we use the emcee (Foreman-Mackey et al. 2013) implementation of the Invariant MCMC Ensemble sampler methods by Goodman & Weare (2010). At the very low visual extinctions, not all the observed CS lines were detected, and we could not determine the gas density (see Table 1). Following this methodology, we estimated the gas density toward 244 positions out of the whole sample. The abundances observed toward these positions constitutes the observational basis of this work.

Once the density was determined from the CS line fitting, we assumed this value to estimate the column densities and abundances of the rest of species: $^{13}\text{CO}$, $\text{C}^{18}\text{O}$, $\text{HCO}^+$, $\text{H}^{13}\text{CO}^+$, $\text{HC}^{18}\text{O}^+$, $\text{H}^{13}\text{CN}$, HNC, $\text{HCS}^+$, SO, $^{34}\text{SO}$, $\text{H}_2\text{S}$, and OCS. To do so, we used the RADEX code and the collisional coefficients from the references listed in Table A.1. The lines of the most abundant species, CO, $\text{HCO}^+$, HCN and HNC are expected to be optically thick and the abundances derived from them are indeed lower limits to the real ones. In order to minimize the effect of high opacities, in this paper we use rare isotopologues to calculate the abundance of the main isotopologue as follows, $X(\text{CO}) = X(\text{C}^{18}\text{O}) \times 550$, $X(\text{HCO}^+) = X(\text{H}^{13}\text{CO}^+) \times 60$, and $X(\text{HCN}) = X(\text{H}^{13}\text{CN}) \times 60$. In the case of HNC, we did not

observe $\text{HN}^{13}\text{C}$ in GEMS, and the estimated abundance might be underestimated in positions with $A_V > 8$ mag. A more detailed description of the molecular database and a first statistical analysis was published by Rodríguez-Baras et al. (2021).

## 3. Chemical network

A great effort has been done in the last decade to improve the gas-phase (Fuente et al. 2016, 2017, 2019; Vidal et al. 2017), and surface (Laas & Caselli 2019) sulfur chemical network. Navarro-Almaida et al. (2020) incorporated all these novelties to build the chemical network used in their paper. In addition, we have implemented the reaction rate for the $S + \text{CH} \rightarrow \text{CS} + \text{H}$ recently estimated by Bulut et al. (in prep.). Using this updated chemical network, we performed the abundance calculations using NAUTILUS 1.1 (Ruaud et al. 2016), which is a numerical model suited to study chemistry in astrophysical environments. It is a three-phase model, in which gas, grain surface and grain movable phases, and their interactions, are considered. It solves the kinetic equations for the gas-phase and the solid species at the surface of interstellar dust grains.

In dark clouds, where the temperature of grain particles is below the sublimation temperature of most species, non-thermal desorption processes are needed to maintain significant abundances of molecules in gas phase. In Nautilus, desorption into the gas phase is only allowed for the surface species, considering both thermal and non-thermal mechanisms. The latter include desorption induced by cosmic rays (Hasegawa & Herbst 1993), direct (UV field) and indirect (secondary UV field induced by the cosmic-ray flux) photo-desorption, and reactive chemical desorption (Garrod et al. 2007; Minissale et al. 2016). In the interior of dark clouds where the UV field is highly attenuated by dust grains and photo-desorption is not efficient, desorption induced by cosmic rays and reactive chemical desorption become the main desorbing mechanisms. Since we are considering molecular gas shielded from the UV radiation, we use the prescription proposed by Minissale et al. (2016) for ice coated grains to calculate the reactive chemical desorption.

## 4. Exploration of the parameter space

In this section, we aim to select the best strategy to determine sulfur depletion ($D_S$) based on GEMS database. We ran a grid of models with typical conditions in dark clouds, in particular the physical conditions listed for Taurus in Table 2. The parameter [S/H] was allowed to take three discrete values that represent the cases of no depletion, $D_S = 1$ ([S/H] = $1.5 \times 10^{-5}$), moderate depletion, $D_S = 10$ ([S/H] = $1.5 \times 10^{-6}$), and high depletion, $D_S \sim 187$ ([S/H] = $8.0 \times 10^{-8}$). These values correspond to different estimates of $D_S$ in star-forming regions (Vastel et al. 2018; Navarro-Almaida et al. 2020, 2021; Hily-Blant et al. 2022). The value of $\zeta_{\text{H}_2}$ varied from $10^{-17} \text{ s}^{-1}$ which is typically found in dense and evolved cores (Caselli et al. 2002) to $\zeta_{\text{H}_2} = 5 \times 10^{-16} \text{ s}^{-1}$, expected in low visual extinction regions (Neufeld & Wolfire 2017). Finally, we considered three chemical ages, $t$, that represent early chemistry (0.1 Myr), late chemistry (1 Myr), and steady-state chemistry (10 Myr). For this exploratory work, we only considered models with $A_V = 11.5$ mag, thus neglecting the effect of the UV field. Using these chemical calculations, we performed corner diagrams of the abundances of the S-bearing species included in GEMS database, CS, SO, $\text{H}_2\text{S}$, and OCS, for different values of $D_S$ (see Figs. A.1–A.4). In each diagram, one of the input parameters is kept fixed while the others are allowed to vary in the whole
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Model clouds formed by selecting all the models with $T=10$ K and $A_V = 11.5$ mag from the output of the grid listed in Table 2 for Taurus. These models are used to illustrate the calculation of the parameter $R$ (upper panels) and $R^{30\%}$ (bottom panels) (see Sect. 4 and Table 3). Upper panels: red and blue lines show the polygons obtained by connecting the extreme values of each [S/H] model cloud. The parameters $R$ is then calculated as the ratio of the red and blue contours intersection area over the total area. Bottom panels: grey lines are iso-density contours corresponding to 30% the point density peak in each model cloud. We show in orange the intersection between the grey contours corresponding to the two values of [S/H]. The parameters $R^{30\%}$ is defined as the ratio of the intersection area (area enclosed within the orange contour) over the total area.</figcaption></figure>

<table><caption>**Table 2.** Chemical model parameters.</caption><thead><tr><th></th><th></th></tr></thead><tbody><tr><td>Taurus</td><td></td></tr><tr><td>n<sub>H</sub></td><td>3.16 &times; 10<sup>3</sup> to 3.16 &times; 10<sup>6</sup> cm<sup>-3</sup> in steps of &times; 3.16</td></tr><tr><td>A<sub>V</sub></td><td>1.5, 3.5, 5.5, 7.5, 9.5, 11.5 mag</td></tr><tr><td>T</td><td>10, 15, 20, 25, 30, 35, 45 K</td></tr><tr><td>CRIR<sup>(1)</sup></td><td>(1.0, 5.0, 10, 50) &times; 10<sup>-17</sup> s<sup>-1</sup></td></tr><tr><td>&chi;<sub>UV</sub> (Draine)<br>[S/H]</td><td>5.0<sup>(2)</sup><br>(1.5, 0.15, 0.008) &times; 10<sup>-5</sup></td></tr><tr><td>Perseus</td><td></td></tr><tr><td>n<sub>H</sub></td><td>3.16 &times; 10<sup>3</sup> to 3.16 &times; 10<sup>6</sup> cm<sup>-3</sup> in steps of &times; 3.16</td></tr><tr><td>A<sub>V</sub></td><td>1.5, 3.5, 5.5, 7.5, 9.5, 11.5, 13.5 mag</td></tr><tr><td>T</td><td>10, 15, 20, 25, 30, 35, 45 K</td></tr><tr><td>CRIR<sup>(1)</sup></td><td>(1.0, 5.0, 10, 50) &times; 10<sup>-17</sup> s<sup>-1</sup></td></tr><tr><td>&chi;<sub>UV</sub> (Draine)<br>[S/H]</td><td>25.0<sup>(2)</sup><br>(1.5, 0.15, 0.008) &times; 10<sup>-5</sup></td></tr><tr><td>Orion</td><td></td></tr><tr><td>n<sub>H</sub></td><td>3.16 &times; 10<sup>3</sup> to 3.16 &times; 10<sup>6</sup> cm<sup>-3</sup> in steps of &times; 3.16</td></tr><tr><td>A<sub>V</sub></td><td>1.5, 3.5, 5.5, 7.5, 9.5, 11.5, 15.5, 20 mag</td></tr><tr><td>T</td><td>10, 15, 20, 25, 30, 35, 45 K</td></tr><tr><td>CRIR<sup>(1)</sup></td><td>(1.0, 5.0, 10, 50) &times; 10<sup>-17</sup> s<sup>-1</sup></td></tr><tr><td>&chi;<sub>UV</sub> (Draine)<br>[S/H]</td><td>50<sup>(3)</sup><br>(1.5, 0.15, 0.008) &times; 10<sup>-5</sup></td></tr></tbody></table>

**Notes.** (1)Molecular hydrogen cosmic-ray ionization rate ($\zeta_{H_2}$).

(2)Estimated by Navarro-Almaida et al. (2020) based on the dust temperature and the analytical expression derived by Hocuk et al. (2017).

(3)Assumed value based on UV field estimated toward the photon-dominated region in the Horesehead nebula (Pety et al. 2005).

range. The fixed parameter is the one indicated in Figs. A.1–A.4. Models with different sulfur depletion are differentiated using a color code. When the model clouds corresponding to different values of $D_S$ appear grouped in bands with little overlap among them, we can conclude that the considered molecules are robust tracers of $D_S$. If model clouds largely overlap then the solution is degenerate and we need to fix additional parameters to determine $D_S$.

In order to base our discussion on quantitative grounds, we define the parameters $R$ and $R^{30\%}$, which measure the overlap between model clouds corresponding to different values of $D_S$. Fig. 1 illustrates how these parameters are calculated. In the upper panels of Fig. 1, the extreme points of models with the same $D_S$ are connected to construct polygons. The $R$ parameter is calculated as, $R = (Intersection)/(Total)$, where Intersection refers to the intersection between the areas of the polygons corresponding to two values of $D_S$ (intersection between the areas defined by the blue and red contours in Fig. 1), and Total indicates the area defined by the whole cloud of points, which mathematically would be Area1 + Area2 – Intersection. With this definition, $R=1$ indicates that the two clouds of points fully overlap and it is not possible to differentiate between the two values of $D_S$. On the contrary, if $R=0$, the intersection is null and we have a clear diagnostic. This parameter gives a good measure of the overlapping when the points are uniformly distributed. However, it is not a good descriptor when the points are not evenly distributed, with some located very far from the bulk. In order to account for these cases, we have defined another parameter, $R^{30\%}$, that uses density curves to define the area corresponding to a given value of $D_S$. In particular, Area1$^{30\%}$ and Area2$^{30\%}$ are defined as the areas enclosed by the 30% point
<table><caption>**Table 3.** R and R³⁰% factors.</caption><thead><tr><th>Fixed parameter</th><th>D<sub>S</sub><sup>(1)</sup></th><th>R<sub>min</sub></th><th>R<sub>max</sub></th><th>R<sub>mean</sub></th><th>R<sup>30%</sup><sub>min</sub></th><th>R<sup>30%</sup><sub>max</sub></th><th>R<sup>30%</sup><sub>mean</sub></th></tr></thead><tbody><tr><td rowspan="3">T = 10 K</td><td>187/10</td><td>0.14(H<sub>2</sub>S-CS)</td><td>0.46(CS-SO)</td><td>0.31</td><td>0.05(H<sub>2</sub>S-CS)</td><td>0.51(OCS-SO)</td><td>0.25</td></tr><tr><td>10/1</td><td>0.43(H<sub>2</sub>S-CS)</td><td>0.62(CS-SO)</td><td>0.53</td><td>0.22(H<sub>2</sub>S-CS)</td><td>0.64(OCS-SO)</td><td>0.44</td></tr><tr><td>187/1</td><td>0.06(H<sub>2</sub>S-CS)</td><td>0.29(OCS-SO)</td><td>0.18</td><td>0.00(H<sub>2</sub>S-CS)</td><td>0.41(OCS-SO)</td><td>0.13</td></tr><tr><td rowspan="3">T = 25 K</td><td>187/10</td><td>0.39(H<sub>2</sub>S-SO)</td><td>0.55(OCS-SO)</td><td>0.49</td><td>0.24(H<sub>2</sub>S-CS)</td><td>0.42(OCS-SO)</td><td>0.31</td></tr><tr><td>10/1</td><td>0.66(H<sub>2</sub>S-CS)</td><td>0.83(OCS-SO)</td><td>0.71</td><td>0.58(H<sub>2</sub>S-OCS)</td><td>0.66(H<sub>2</sub>S-CS)</td><td>0.62</td></tr><tr><td>187/1</td><td>0.27(H<sub>2</sub>S-SO)</td><td>0.54(OCS-SO)</td><td>0.37</td><td>0.21(H<sub>2</sub>S-CS)</td><td>0.46(OCS-SO)</td><td>0.35</td></tr><tr><td rowspan="3">T = 35 K</td><td>187/10</td><td>0.32(OCS-H<sub>2</sub>S)</td><td>0.58(OCS-CS)</td><td>0.47</td><td>0.04(OCS-CS)</td><td>0.69(OCS-H<sub>2</sub>S)</td><td>0.31</td></tr><tr><td>10/1</td><td>0.56(OCS-H<sub>2</sub>S)</td><td>0.65(OCS-SO)</td><td>0.60</td><td>0.21(OCS-SO)</td><td>0.46(H<sub>2</sub>S-CS)</td><td>0.35</td></tr><tr><td>187/1</td><td>0.28(OCS-SO)</td><td>0.43(OCS-CS)</td><td>0.34</td><td>0.03(OCS-CS)</td><td>0.40(OCS-H<sub>2</sub>S)</td><td>0.20</td></tr><tr><td rowspan="3">n<sub>H</sub> = 10<sup>4</sup> cm<sup>-3</sup></td><td>187/10</td><td>0.37(OCS-CS)</td><td>0.58(H<sub>2</sub>S-CS)</td><td>0.46</td><td>0.01(H<sub>2</sub>S-CS)</td><td>0.25(OCS-SO)</td><td>0.14</td></tr><tr><td>10/1</td><td>0.48(H<sub>2</sub>S-SO)</td><td>0.62(OCS-SO)</td><td>0.56</td><td>0.25(OCS-CS)</td><td>0.59(OCS-SO)</td><td>0.40</td></tr><tr><td>187/1</td><td>0.18(H<sub>2</sub>S-CS)</td><td>0.30(OCS-SO)</td><td>0.25</td><td>0.00(OCS-CS)</td><td>0.21(OCS-SO)</td><td>0.07</td></tr><tr><td rowspan="3">n<sub>H</sub> = 10<sup>5</sup> cm<sup>-3</sup></td><td>187/10</td><td>0.34(OCS-H<sub>2</sub>S)</td><td>0.57(OCS-SO)</td><td>0.45</td><td>0.32(OCS-CS)</td><td>0.58(OCS-SO)</td><td>0.43</td></tr><tr><td>10/1</td><td>0.33(OCS-H<sub>2</sub>S)</td><td>0.60(CS-SO)</td><td>0.41</td><td>0.29(H<sub>2</sub>S-CS)</td><td>0.63(OCS-SO)</td><td>0.43</td></tr><tr><td>187/1</td><td>0.09(OCS-H<sub>2</sub>S)</td><td>0.32(CS-SO)</td><td>0.19</td><td>0.00(OCS-CS)</td><td>0.34(OCS-SO)</td><td>0.16</td></tr><tr><td rowspan="3">n<sub>H</sub> = 10<sup>6</sup> cm<sup>-3</sup></td><td>187/10</td><td>0.40(H<sub>2</sub>S-SO)</td><td>0.65(OCS-SO)</td><td>0.53</td><td>0.38(H<sub>2</sub>S-SO)</td><td>0.63(CS-SO)</td><td>0.51</td></tr><tr><td>10/1</td><td>0.49(H<sub>2</sub>S-SO)</td><td>0.71(CS-SO)</td><td>0.58</td><td>0.27(OCS-SO)</td><td>0.71(CS-SO)</td><td>0.48</td></tr><tr><td>187/1</td><td>0.19(H<sub>2</sub>S-SO)</td><td>0.50(CS-SO)</td><td>0.35</td><td>0.23(OCS-CS)</td><td>0.53(CS-SO)</td><td>0.35</td></tr><tr><td rowspan="3">t = 0.1 Myr</td><td>187/10</td><td>0.48(H<sub>2</sub>S-CS)</td><td>0.63(CS-SO)</td><td>0.53</td><td>0.28(SO-CS)</td><td>1.0 (OCS-H<sub>2</sub>S)</td><td>0.47</td></tr><tr><td>10/1</td><td>0.55(OCS-CS)</td><td>0.68(CS-SO)</td><td>0.64</td><td>0.43(OCS-CS)</td><td>0.63(OCS-H<sub>2</sub>S)</td><td>0.55</td></tr><tr><td>187/1</td><td>0.32(H<sub>2</sub>S-CS)</td><td>0.42(CS-SO)</td><td>0.37</td><td>0.06(SO-CS)</td><td>1.0 (OCS-H<sub>2</sub>S)</td><td>0.40</td></tr><tr><td rowspan="3">t = 1 Myr</td><td>187/10</td><td>0.41(H<sub>2</sub>S-SO)</td><td>0.59(CS-SO)</td><td>0.51</td><td>0.29(OCS-CS)</td><td>1.0(CS-SO)</td><td>0.50</td></tr><tr><td>10/1</td><td>0.54(H<sub>2</sub>S-SO)</td><td>0.75(CS-SO)</td><td>0.64</td><td>0.42(H<sub>2</sub>S-SO)</td><td>0.61(CS-SO)</td><td>0.52</td></tr><tr><td>187/1</td><td>0.26(H<sub>2</sub>S-SO)</td><td>0.52(OCS-H<sub>2</sub>S)</td><td>0.41</td><td>0.04(OCS-CS)</td><td>0.34(OCS-SO)</td><td>0.24</td></tr><tr><td rowspan="3">t = 10 Myr</td><td>187/10</td><td>0.51(OCS-CS)</td><td>0.61(CS-SO)</td><td>0.58</td><td>0.24(OCS-CS)</td><td>0.61(OCS-SO)</td><td>0.47</td></tr><tr><td>10/1</td><td>0.56(H<sub>2</sub>S-SO)</td><td>0.70(OCS-SO)</td><td>0.65</td><td>0.41(OCS-H<sub>2</sub>S)</td><td>0.64(CS-SO)</td><td>0.53</td></tr><tr><td>187/1</td><td>0.32(H<sub>2</sub>S-SO)</td><td>0.48(SO-CS)</td><td>0.41</td><td>0.16(OCS-CS)</td><td>0.42(OCS-SO)</td><td>0.24</td></tr><tr><td rowspan="3">ζ<sub>H<sub>2</sub></sub> = 1 × 10<sup>-17</sup> s<sup>-1</sup></td><td>187/10</td><td>0.48(H<sub>2</sub>S-CS)</td><td>0.61(CS-SO)</td><td>0.55</td><td>0.47(H<sub>2</sub>S-CS)</td><td>0.61(CS-SO)</td><td>0.54</td></tr><tr><td>10/1</td><td>0.57(H<sub>2</sub>S-CS)</td><td>0.72(CS-SO)</td><td>0.65</td><td>0.58(OCS-SO)</td><td>0.75(CS-SO)</td><td>0.62</td></tr><tr><td>187/1</td><td>0.27(H<sub>2</sub>S-CS)</td><td>0.44(CS-SO)</td><td>0.37</td><td>0.31(H<sub>2</sub>S-CS)</td><td>0.50(CS-SO)</td><td>0.37</td></tr><tr><td rowspan="3">ζ<sub>H<sub>2</sub></sub> = 5 × 10<sup>-17</sup> s<sup>-1</sup></td><td>187/10</td><td>0.45(H<sub>2</sub>S-CS)</td><td>0.67(OCS-SO)</td><td>0.57</td><td>0.36(OCS-CS)</td><td>0.55(OCS-SO)</td><td>0.46</td></tr><tr><td>10/1</td><td>0.50(H<sub>2</sub>S-SO)</td><td>0.70(CS-SO)</td><td>0.59</td><td>0.49(H<sub>2</sub>S-SO)</td><td>0.62(OCS-SO)</td><td>0.54</td></tr><tr><td>187/1</td><td>0.27(H<sub>2</sub>S-SO)</td><td>0.46(CS-SO)</td><td>0.35</td><td>0.26(H<sub>2</sub>S-SO)</td><td>0.29(OCS-CS)</td><td>0.29</td></tr><tr><td rowspan="3">ζ<sub>H<sub>2</sub></sub> = 1 × 10<sup>-16</sup> s<sup>-1</sup></td><td>187/10</td><td>0.41(H<sub>2</sub>S-CS)</td><td>0.64(OCS-SO)</td><td>0.52</td><td>0.31(OCS-CS)</td><td>0.54(OCS-SO)</td><td>0.42</td></tr><tr><td>10/1</td><td>0.51(H<sub>2</sub>S-SO)</td><td>0.65(CS-SO)</td><td>0.57</td><td>0.48(H<sub>2</sub>S-SO)</td><td>0.63(OCS-SO)</td><td>0.56</td></tr><tr><td>187/1</td><td>0.24(H<sub>2</sub>S-SO)</td><td>0.40(OCS-SO)</td><td>0.31</td><td>0.19(OCS-CS)</td><td>0.33(OCS-SO)</td><td>0.26</td></tr></tbody></table>

Notes. $^{(1)}$D<sub>S</sub> = [S/H]<sub>solar</sub>/[S/H]<sub>solar</sub>, where [S/H]<sub>solar</sub> = 1.5 × 10<sup>-5</sup>.

density peak contour, and R³⁰% = (Intersection³⁰%)/(Total³⁰%) (see bottom panels of Fig. 1). As expected, we find R³⁰% < R in the case of a peaky distribution of points.

Table 3 shows the minimum, maximum, and mean values of R and R³⁰% for the panels forming the corner diagrams shown in Figs. A.1–A.4. We remind that we have only used the pairs of molecules, OCS-H<sub>2</sub>S, OCS-SO, OCS-CS, H<sub>2</sub>S-CS, H<sub>2</sub>S-SO, and CS-SO, in our analysis. Thus, the minimum, maximum, and mean R and R³⁰% values are estimated taking into account only these six combination of sulfur species. Although the exact value of R and R³⁰% are slightly different, we find the same trends in these parameters which demonstrates the robustness of our analysis. One first result is that we need to fix temperature to reach values, R³⁰% mean < 0.3. Moreover the lowest values of R³⁰% mean

are found for T = 10 K. In Fig. A.1, we show the corner diagrams when fixing the gas temperature. At low temperature, T ~ 10 K, CS and H<sub>2</sub>S are the best proxies of D<sub>S</sub>. In the panels including these species, the points corresponding to different values of D<sub>S</sub> appear grouped in clearly distinguishable bands. Moreover, the bands are almost orthogonal to the X(H<sub>2</sub>S) axis suggesting that the abundance of H<sub>2</sub>S is an excellent tracer of D<sub>S</sub> with little dependence on the other parameters. This can be easily interpreted in terms of sulfur chemistry. At such a low temperature, solid H<sub>2</sub>S is the main sulfur reservoir, and its abundance in gas-phase is only dependent on the amount of sulfur atoms in the ice and the efficiency of non-thermal desorption mechanisms (Navarro-Almaida et al. 2020). The ubiquitous molecule CS is not such a good tracer of D<sub>S</sub>. Values of X(CS) ≥ 10<sup>-6</sup>, can only
be achieved with low sulfur depletion (see Fig. A.1). However, the situation is less clear for X(CS) < 10⁻⁸ that can be explained assuming the three values of Ds. These low values of X(CS) are the most common in the cold interstellar medium (Agúndez & Wakelam 2013), which hinders the determination of Ds on the basis of only this compound.

Fixing the values of density is also useful, especially for nH < 10⁵ cm⁻³ with Rmean < 0.4. The corner diagrams in Fig. A.2 show that, fixing density, model clouds are distinguishable, although there is always some overlap among them. The overlap increases with increasing density, which implies that the estimated value of Ds is less reliable at high densities.

As mentioned above, the cosmic-ray molecular hydrogen ionization rate, ζH₂, is a key parameter for the chemistry in dark clouds since it is governing the gas ionization degree. Figure A.3 shows the corner diagrams for fixing the values of ζH₂. There is an important overlap between the models corresponding to the different values of Ds in all cases, challenging its estimation if we only fix this parameter. The same behavior is found when we fix the chemical time (see Fig. A.4). This is consistent with the high values of R and R³⁰% shown in Table 3. Therefore, fixing ζH₂ and/or t does not help to determine Ds.

According to these results, it would be desirable to constrain the gas temperature and density in order to have a good estimate of Ds in our sample. As explained in Sect. 5, we assume that gas and dust are thermalized for our fittings, and use the dust temperature derived from Herschel maps to constrain T. We also constrain density using the values obtained by Rodríguez-Baras et al. (2021) through a multitransitional study of CS and its isotopologues. The parameters ζH₂ and t will be fitted together with Ds using our molecular database. It is well known that the HCO⁺/CO and HCS⁺/CS abundance ratios are good tracers of ζH₂ (see e.g., Fuente et al. 2016). The chemical age is better constrained with the CO abundance. This molecule is frozen on dust grains for high densities (nH > 10⁴ cm⁻³) and dust temperatures <16 K (Wakelam et al. 2021). Since CO depletion depends on the local volume density, we can obtain a reliable value of the chemical age as long as the density is known.

Our exploratory study also shows that it is always easier to discern between Ds = 187 and Ds = 10, than between Ds = 10 and Ds = 1. This fact should be taken into account when interpreting our fits.

## 5. Methodology

Our goal is to determine Ds in each of the positions observed within GEMS. Our large molecular database with 244 positions allows us to explore possible trends of Ds with the local physical parameters (n, T) and/or environment. We fit the abundances of CO, HCO⁺, HCN, HNC, CS, HCS⁺, H₂S, SO, and OCS as derived by Rodríguez-Baras et al. (2021) with the output obtained with the grid of chemical models shown in Table 2. We adopted χUV = 5 for Taurus and χUV = 25 for Perseus. These values were derived by Fuente et al. (2019) and Navarro-Almaida et al. (2020) using the analytic expression obtained by Hocuk et al. (2017) and the Herschel dust temperature images. In Orion, we adopted χUV = 50, similar to that derived in the photon-dominated region associated with the Horsehead nebula (Pety et al. 2005; Goicoechea et al. 2006; Gerin et al. 2009; Rivière-Marichalar et al. 2019). It should be noticed that some compounds are not detected toward the whole sample of 244 positions. In these cases, we perform the fitting using the detected species as long as the number of species is higher than 3, and contains at least one sulfur-bearing species.

During the fitting process, Ds, ζ(H₂), and t are allowed to freely vary among the range of discrete values adopted in the grid. However, we impose some restrictions to the values of Aν, n(H₂) and T. The visual extinction is assumed to be in the range Aν × 0.5 and Aν × 0.5 + 2, where Aν is the value obtained from Herschel data (Malinen et al. 2012; Palmeirim et al. 2013; Lombardi et al. 2014; Zari et al. 2016), and the factor 0.5 accounts for the fact that the cloud is expected to be illuminated from the back and the front. Therefore, at the center of the molecular cloud, the effective visual extinction would be half of the total along the line of sight. The density is assumed to be in the range from 0.5 × n(H₂) to 5 × n(H₂), being the values of n(H₂), those obtained by Rodríguez-Baras et al. (2021). The gas kinetic temperature is allowed to vary between Td and Td + 5 K.

Finally, we need to select a parameter to describe the goodness of each model in describing the observational results. The standard χ² is not adequate to estimate errors when the molecular abundances included in the fit differ by several orders of magnitude. In order to find the “best-fitting” model at each position, we use the parameter Ddiff that is defined by,

$$D_{\text{diff}} = 1/n_{\text{obs}} \times \sum_i [\log_{10}(X_{\text{mod}}^i) - \log_{10}(X_{\text{obs}}^i)]^2 \quad (1)$$

where nobs is the number of species detected at each position, Xmod is the model predicted abundance for the species i, and Xobs is the abundance of the species i derived from GEMS observations. This parameter is the square of the Disagreement parameter that has been previously used in astrochemistry by different authors (see e.g., Wakelam et al. 2006; Vastel et al. 2018).

## 6. Results

### 6.1. Taurus

The Taurus molecular cloud is considered as a prototype of low-mass star-forming region. At a distance of 145 pc (Yan et al. 2019), it has been extensively studied at infrared and millimeter wavelengths (Cernicharo & Guelin 1987; Onishi et al. 1996; Narayanan et al. 2008; Cambrésy 1999; Padoan et al. 2002; André et al. 2010; Schmalzl et al. 2010). Within GEMS, we have studied two well-known filaments: TMC 1 and B 213/L1495.

In TMC 1, we observed three cuts along the cores TMC 1-CP, TMC 1-NH3, and TMC 1-C, respectively. The positions TMC 1-CP and TMC 1-NH3 (the cyanopolyynes and ammonia emission peaks) are generally adopted as templates of carbon- and oxygen-rich starless cores to compare with chemical codes (e.g., Fehér et al. 2016; Gratier et al. 2016; Agúndez & Wakelam 2013). Less studied from the chemical point of view, TMC 1-C has been identified as an accreting starless core (Schnee et al. 2007, 2010). Fuente et al. (2019) carried out a complete analysis of the data associated with these cores to derive the gas ionization degree and elemental abundances in the TMC 1 translucent cloud. They concluded that the chemistry of the translucent filament is described assuming ζH₂ = (0.5–1.8) × 10⁻¹⁶ s⁻¹ and [S/H] = (0.4–2.2) × 10⁻⁶.

B 213/L 1495 is a prominent filament in Taurus that has been widely studied in the (sub)mm range (Palmeirim et al. 2013; Hacar et al. 2013; Marsh et al. 2014; Tafalla & Hacar 2015; Bracco et al. 2017; Shimajiri et al. 2019). The morphology of the map with striations perpendicular to the filament suggests that the filament is accreting material from its surroundings (Goldsmith et al. 2008; Palmeirim et al. 2013). Shimajiri et al. (2019) proposed that this active star-forming filament was
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Statistics of our molecular fitting in Taurus. The histograms show the fraction of positions that have been fitted with a given value of Age (left), $\zeta_{H_2}$ (center), and [S/H] (right). The positions are grouped by cut in the upper panels (see Table 1), and by cloud in the lower panels. It should be noticed that none of the Taurus positions have been fitted with $[S/H] = 1.5 \times 10^{-5}$.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Same as Fig. 2, but in the upper row, the Taurus positions have been separated in bins of density. In the bottom panels, the positions have been distributed in bins of visual extinction.</figcaption></figure>

initially formed by large-scale compression of HI gas, and is now growing in mass due to the gravitational accretion of ambi-ent molecular gas. A chain of starless cores (Hacar et al. 2013) and protostars (Luhman et al. 2009; Rebull et al. 2010; Davis et al. 2010) is observed along the filament. Within GEMS, we observed 9 cuts along clumps #1, # 2, #5, #6, #7, #10, #12, #16, and #17 (core numbers from the catalog of Hacar et al. 2013). Cuts #1, # 2, #5, #6, and #7 are located in the most active star-forming region (hereafter, B 213-N) where the density of protostars is higher. Cuts #10, #12, #16, and #17 are located in the quiescent part (hereafter, B 213-S). Previous works based on CH$_3$OH data suggest the existence of chemical differentiation between B 213-N and B 213-S (Spezzano et al. 2022; Punanova et al. 2022). In the cores located in B 213-N, the methanol peak is spatially coincident with a visual extinction peak, while the contrary behavior is found toward the south. It has been pro-posed that this behavior could be due to the irradiation on the cores due to nearby protostars which accelerate energetic parti-cles along their outflows. A recent study of the H$_2$CS deuterated compounds (Esplugues et al. 2022), showed that the chemical age of the cores in B 213-N is higher than that of the cores in B 213-S. A combination of environmental and dynamical effects are therefore needed to account for the observed chemical variations.

The GEMS data provide new insights into the variations of chemical age, cosmic-ray ionization rate, and sulfur depletion among the different regions of Taurus. Our approach consists of fitting the abundances of nine species, CO, HCO$^+$, HCN, HNC,

CS, HCS$^+$, H$_2$S, SO, and OCS toward 84 positions of Taurus distributed as shown in Table 1. The results are shown in Figs. 2 and 3. Essentially all the positions in B 213 are best fitted assum-ing early time chemistry ($t = 0.1$ Myr) which is consistent with the idea that this filament is still accreting material from the sur-roundings, keeping the gas chemistry far from the steady state. In TMC 1, we have some dispersion in the values of the chemical age with the positions located at $A_V < 8$ mag being better fitted with $t = 1$ Myr while those at higher visual extinction being bet-ter fitted with $t = 0.1$ Myr. We also find a large dispersion in the values of $\zeta_{H_2}$. In spite of this, we can find some trends. Most positions in TMC 1 and B 213-N are fitted with $\zeta_{H_2} > 10^{-16}$ s$^{-1}$, while B 213-S is better described with $\zeta_{H_2} < 10^{-16}$ s$^{-1}$, consis-tent with the idea of a harsh environment in B 213-N due to the presence of young (proto-)stars. A high value of $\zeta_{H_2}$ has also been argued to explain the richness of TMC 1 in complex carbon chains (Agúndez & Wakelam 2013).

Regarding the sulfur elemental abundance, we find that $D_S = 10$ in TMC 1 and B 213-N, while $D_S = 187$ in B 213-S. Based on observations of NS, Hily-Blant et al. (2022) proposed that the sulfur depletion is higher in the more evolved starless cores such as L1544, compared with others at an earlier evo-lutionary stage. In Fig. 3 we explore the possible dependence of sulfur depletion with density and visual extinction. There is no clear trend with the visual extinction. However, we detect a trend of $D_S$ with density, with the densest cores showing higher sulfur depletion. The density is expected to increase along the starless core evolution to form a collapsing core. This would
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Statistics of our molecular fitting in Perseus. The upper row shows the histograms with the positions separated by cloud. In the middle row, the positions have been distributed in bins of density, and in the bottom row, in bins of visual extinction.</figcaption></figure>

indicate that more evolved cores would present higher values of sulfur depletion in agreement with the results of Hily-Blant et al. (2022). Esplugues et al. (2022) found that together with density, higher values of sulfur depletion favors the formation of singly and doubly deuterated H₂CS.

NAUTILUS 1.1 gas-grain code includes the freeze out of molecules onto grain surfaces, which is known to increase with density when the grain temperature is below the evaporation temperature of a given species. This means that an additional process which removes sulfur from the volatile chemistry is needed to explain the increase of the depletion with density. On the other hand, as explained in Sect. 4, the reliability of our calculation of D~S~ decreases in dense regions. Therefore, we need to be cautious with this result.

## 6.2. Perseus

Located at a distance of 310 pc (Ortiz-León et al. 2018), the Perseus molecular cloud is considered the prototype of intermediate-mass star-forming region. This large molecular cloud complex is associated with two clusters containing pre-main-sequence stars: IC 348, with an estimated age of 2 Myr (Luhman et al. 2003), and NGC 1333, which is thought to be younger than <1 Myr (Lada et al. 1996; Wilking et al. 2004). The cloud has been extensively studied using a large variety of techniques and wavelengths (Bachiller & Cernicharo 1984; Warin et al. 1996; Hatchell et al. 2005, 2007a,b; Ridge et al. 2006; Kirk et al. 2006; Enoch et al. 2006; Curtis & Richer 2011; Pineda et al. 2008, 2010, 2015; Zari et al. 2016; Friesen et al. 2017; Hacar et al. 2017b). In particular, the content in dense cores was described in a series of papers based on continuum maps at 850 and 450 µm obtained with SCUBA at the JCMT (Hatchell et al. 2005, 2007a,b). Hatchell et al. (2007b) classified the 91 dense cores detected using their Spectral Energy Distribution (SED), resulting in 47 starless cores, 34 Class 0, and 22 Class I protostars. Later, Hatchell et al. (2007a) surveyed the outflow activity in the region to gain a deeper insight into the evolutionary stage of the protostars. In

contrast to Taurus, a significant fraction of these protostars are forming proto-clusters.

We have observed 11 cuts along starless cores distributed in Barnard 1, IC 348, L 1448, NGC 1333, and B5. The group of cores in IC 348 and NGC 1333 are close to the star clusters and therefore immersed in a harsh environment. The cuts associated with B5 and L 1448 are located in more quiescent regions. The Barnard 1 dark cloud contains several dense cores at different evolutionary stages of the star formation process. While B1-a and B1-c are known to host Class 0 sources associated with high velocity outflows (Hatchell et al. 2007a), the B1-b core appears to be more quiet, although its western edge is interacting with an outflow that is possibly arising from sources B1-a or B1-d, which are both located 1' SW of B1-b. Recently, the B1-b core has been associated with an extremely young Class 0 object, B1b-S, and a First Hydrostatic Core (FSHC) candidate, B1b-N (Gerin et al. 2015, 2017; Fuente et al. 2017; Marcelino et al. 2018).

We have performed the fitting of the observed molecular abundances using the grid shown in Table 2. Essentially, all positions are better fitted assuming early time chemistry (~0.1 Myr) (see Fig. 4). Similarly to the case of Taurus, the regions at low visual extinction are better fitted with *t* = 1 Myr. Moreover, there is a large dispersion in the values of ζ~H~₂ without any clear trend with visual extinction and density (see Fig. 4). Regarding sulfur depletion, most of the positions are fitted with D~S~ = 10. Contrary to Taurus, we find some positions toward which the best fit corresponds to undepleted sulfur. These positions are found in IC 348, NGC 1333, and B5. It should be noticed that IC 348 and NGC 1333 are the nearest regions to star clusters. Moreover, NGC 1333 hosts a large population of Class 0 and I protostars. The high values of the sulfur elemental abundance in these regions is very likely related to the star formation activity in the vicinity. Finally, we also show the behavior of the sulfur depletion as a function of the density and visual extinction in Fig. 4. We do not detect any clear trend with these parameters. This is not unexpected since the highest density and highest visual extinction regions are located in NGC 1333 which is hosting numerous protostars and bipolar outflows that are modifying the molecular chemistry.
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Statistics of our molecular fitting in Orion. The upper row shows the histograms with the positions separated by cloud. In the middle row, the positions have been distributed in bins of density, and in the bottom row, in bins of visual extinction.</figcaption></figure>

### 6.3. Orion

Orion is the nearest massive star-forming cluster to Earth (e.g., Hillenbrand 1997; Lada et al. 2000; Muench et al. 2002; Da Rio et al. 2012; Robberto et al. 2013; Zari et al. 2019) and it is located at a distance of ~428 pc (Zucker et al. 2019). Traditionally, the “Orion nebula” refers to the visible part of the region, the HII region, powered by the ionizing radiation of the Trapezium OB association. However, this nebula is part of a much larger complex, referred to as the Orion Molecular Cloud (OMC), that is formed by two giant molecular clouds: Orion A hosting the M42 nebula, and the more quiescent Orion B (see, e.g., Pety et al. 2017).

Orion A has been extensively mapped in molecular lines using single-dish telescopes and large millimeter arrays (Hacar et al. 2017a, 2018, 2020; Goicoechea et al. 2019; Nakamura et al. 2019; Tanabe et al. 2019; Ishii et al. 2019; Kirk et al. 2017; Monsch et al. 2018; Suri et al. 2019; Kong et al. 2019). Different clouds can be identified within Orion A based on millimeter, submillimeter, and infrared observations. Orion molecular cloud 1 (OMC 1) was identified as the dense gas directly associated with Orion KL (Wilson et al. 1970; Zuckerman 1973; Liszt et al. 1974). The molecular cloud associated with the HII region M 43 is referred to as OMC 2 (Gatley et al. 1974). OMC 3 is composed by a series of clumps detected in CO emission that are located about 25' to the north of OMC 1 (Kutner et al. 1976). The ¹³CO (*J*=1→0) observations by Bally et al. (1987) revealed that all these clouds form the so-called integral-shaped filament (ISF) of molecular gas, itself part of a larger filamentary structure extending from north to south over 4°. After that, the SCUBA maps at 450 µm and 850 µm presented concentrations of submillimeter continuum emission in the southern part of the ISF, which are now referred to as OMC 4 (Johnstone & Bally 1999) and OMC 5 (Johnstone & Bally 2006).

We have observed three cuts: Ori-C1 in OMC 3, Ori-C2 in OMC 4, and Ori-C3 in OMC 2. These cuts are selected to avoid the (proto-)stars, probing quiescent environments far from the Orion nebula. We have fitted the detected molecular abundances using the grid shown in Table 2, and the results are shown in

Fig. 5. This region presents clear differences relative to Taurus and Perseus. Although the cut Ori-C1 in OMC 3 is fitted with early time chemistry (*t* = 0.1 Myr), the cut Ori-C2 in OMC 4 is fitted with a chemical age of *t*=10 Myr, and the cut Ori-C3 in OMC 2, with *t* = 1 Myr. Also contrary to Taurus and Perseus, we obtain a low value of cosmic-ray molecular hydrogen ionization rate, ζ<sub>H₂</sub> < 10<sup>-16</sup> s<sup>-1</sup>, in most positions of the three cuts. This is an unexpected result taking into account the high star formation rate in the Orion A. We remind that the cuts are selected in quiescent regions far from the massive protostar cluster. On the other hand, the higher temperatures and UV flux in Orion could reduce the sensitivity of the studied molecular abundances to the value of ζ<sub>H₂</sub>. Finally, a significant fraction of the observations are best fitted assuming D<sub>S</sub> = 1. This is a notorious difference with Taurus, and with most of the clouds in Perseus, and supports the role of the environment on the sulfur elemental abundance. As in previous cases, we have explored possible trends of the sulfur elemental abundance with the gas density and the visual extinction. We do not detect any trend with density (see Fig. 5). However, there is a weak trend with visual extinction where the higher values of D<sub>S</sub> are found at high visual extinctions. It is reasonable that the visual extinction plays a more important role in this giant molecular cloud, which is bathed by a more intense UV radiation field.

## 7. The complete sample: Determining D<sub>S</sub>

In Sect. 6, we have discussed the results of our model fitting considering the three regions, Taurus, Perseus, and Orion, separately. We can also merge the complete set of data in order to obtain an overall view. Our first result is that environment is a key parameter to determine sulfur depletion (see Fig. 6a). More than 50% of the positions in Taurus and Perseus are best fitted with D<sub>S</sub> ~ 10. Basically, all the other positions need a higher depletion, D<sub>S</sub> ~ 187. A different behavior is observed in Orion, in which a significant fraction of the positions (~40%) can be fitted with D<sub>S</sub> = 1. This suggests that sulfur depletion is dependent on the star formation activity in the environment. It
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Statistics performed with all the positions (Taurus, Perseus, and Orion) of the molecular database. (a) Histograms with the positions ordered by cloud complex (upper panel), in bins of density (middle panel), and in bins of visual extinction (bottom panel). (b) Mean value of the sulfur elemental abundance derived in each cut. The value of [S/H] toward the high extinction peak is indicated with the dashed area.</figcaption></figure>

is also interesting to explore possible trends of $D_S$ with density and visual extinction. However, we need to put a word of caution in this analysis. Since only ~13% of the considered positions belong to Orion (see Table 1), the results are biased to the physical conditions prevailing in Taurus and Perseus. There is a trend of increasing sulfur depletion with density, which suggests a progressive incorporation of sulfur into grains. However, we cannot say that there is a clear correlation between depletion and density. This trend is dominated by the positions with densities $n_H \le 10^4$ cm$^{-3}$ that present lower values of $D_S$. Since these low-density positions are located in Taurus and Perseus, the statistics is biased. In addition, as discussed in Sect. 4, the uncertainties in the estimated sulfur elemental abundance also increase with density, hindering to establish firm conclusions. No trend is observed in the histograms of $D_S$ as a function of visual extinction. As commented above, only in Orion we can find some correlation, and this correlation is missed when merging all positions.

In our first grid, we have only selected three values of $D_S$ in order to keep calculation time within reasonable terms. This coarse grid allowed us to consider a wide range of physical parameters but prevented us from determining the value of $D_S$ with an accuracy better than a factor of ~10. In a second step, we fixed the values of $t$ and $\zeta_{H_2}$ according to the results obtained in our first grid, and then ran a second one with finer steps of $D_S$. In particular, we used $\zeta_{H_2} = 10^{-16}$ s$^{-1}$ for Taurus, $\zeta_{H_2} = 5 \times 10^{-17}$ s$^{-1}$ in Perseus, and $\zeta_{H_2} = 10^{-17}$ s$^{-1}$ for Orion. The chemical age was fixed to 0.1 Myr in Taurus and Perseus. In the case of Orion where the $t$ is not well determined, we repeated the fitting with an age of 0.1 Myr and 1 Myr, and then selected the chemical age providing the lowest value of $D_{\text{diff}}$ for each cut. According to the results, we selected $t=0.1$ Myr for Ori-C1 (OMC 3), and 1 Myr for Ori-C2 (OMC 4) and Ori-C3

(OMC 2). The value of $D_S$ for each cut was varied from 1 to 187 in multiplicative steps of a factor of ~2 ($D_S = 1, 2, 4, 8, 16, 32, 65, 187$).

In Table A.2, we show a summary of the results thus obtained for the 29 cuts of our sample. The comparison between the model-predicted abundances and those derived from observations toward the visual extinction peak in each cut are shown in Figs. A.5 and A.6. As shown in Table A.2, the values of $D_{\text{diff}}$ vary between 0.1 and 0.5. This means that the abundances of most of the observed molecules are predicted within a factor of 2–5. However, there are outliers with errors of a factor of ~10 toward some positions (see Figs. A.5 and A.6). These outliers are not necessarily sulfur-bearing species. Indeed we have large errors when fitting the CO, HCN, and HNC abundances in a few locations. These errors cannot be attributed to opacity effects in the case of CO and HCN because, as explained in Sect. 5, we are using the less abundant isotopologues $C^{18}\text{O}$ and $H^{13}\text{CN}$ to estimate their abundances. As commented above, isotopic chemical fractionation might be important for HCN (Loison et al. 2020), but the HCN/$H^{13}\text{CN}$ does not differ from the assumed value in more than a factor of ~2. The discrepancy between the observations and model predictions are due to uncertainties in the physical structure and chemistry. We remind that we are using a 0D model that neglects the temperature and density gradients along the line of sight.

Within the group of sulfur-bearing species, we find some systems in the errors, with the SO abundance being usually over-predicted by the chemical model, while the $H\text{CS}^+$ abundance is under-predicted. This reflects the problem already discussed by Navarro-Almaida et al. (2020) and Bulut et al. (2021) on the difficulty of fitting the abundances of all sulfur-bearing species using the same input parameters. This also explains the discrepancy between the results obtained in this work and previous
papers. With a slightly different chemical network, Vidal et al. (2017) concluded that the abundances of sulfur bearing species in TMC 1-CP can be fitted assuming undepleted or low depleted (a factor of ~3) sulfur abundance. With a chemical network very similar to that used in this paper, Navarro-Almaida et al. (2020) fitted the H₂S and emission in B1b and TMC 1 with undepleted sulfur abundance. This fit, however, over-estimated the abundance of CS by factor ≥10 along the cores. In a later paper, Bulut et al. (2021) needed to assume a depletion of a factor of 20 to fit the CS and HCS⁺ abundances in TMC 1-CP, TMC1-NH3, and TMC1-C. In this work, we follow a uniform and systematic method to derive the sulfur depletion in all the GEMS positions. Although our study is limited by the current knowledge of the sulfur chemistry, we are providing a uniform database that allows a reliable comparison between different regions. Moreover, we can take advantage of the large number of positions measured in GEMS to obtain average values, and hence minimize possible observational errors.

In Table A.2, we list the values of [S/H] obtained toward the visual extinction peak, and the mean value of [S/H] along each of the GEMS cuts. These two values intend to be representative of the sulfur elemental abundance in the dense core and the envelope, respectively. This same information is graphically represented in Fig. 6b. The envelopes in all the cuts located Taurus and Perseus, except TMC1-NH3 and B1b, are well reproduced with a depletion of ≥20. In Orion, the estimates of [S/H] are compatible with undepleted sulfur within the uncertainties. The mean value and that toward the extinction peak differ only by a factor of 2–3 in most cuts. This is expected behavior taking into account that our chemical model takes into account the freeze-out of molecules on the grain surfaces in cold and dense regions. However, there are some exceptions to this rule: B 213-C6, B 213-C7, TMC 1-NH3, IC 348-C10, NGC1333-C5, B1b, ORI-C1, ORI-C2, and ORI-C3. In these cases, the sulfur elemental abundance in the envelope is a factor of >4 higher than that toward the extinction peak. Interestingly, the value of [S/H] toward the extinction peak is similar to those measured in the rest of regions. Although the correspondence is not perfect, all these cuts are located in regions where molecular chemistry is known to present some peculiarities. The cuts B 213-C6 and B 213-C7 are located in the northern part of the L1495/B 213 region, which is associated with a cluster of (proto-)stars (Hacar et al. 2013). Recent observations of methanol in this region shows that, contrary to the behavior in the more quiescent starless cores in the south, the maximum of the N(CH₃OH)/N(C¹⁸O) ratio is found toward the extinction peak in these cuts (Spezzano et al. 2022). These authors suggest that this behavior is related to the existence of low-mass protostars in the neighborhood. Low-mass stars can accelerate energetic particles, thus increasing the local cosmic-ray flux and the methanol abundance in the densest inner region of the starless core. A high cosmic-ray ionization rate flux has been also proposed to explain the active chemistry in TMC 1-CP and TMC 1-NH3 by several authors (Fuente et al. 2019). Relative to NGC 1333-C5, the cut associated with this core is crossing a cluster of protostars associated with energetic bipolar outflows (Hatchell et al. 2007a; Hatchell & Dunham 2009). A different explanation needs to be found for the cuts studied in Orion since our fittings suggest that the cosmic-ray ionization rate is low in these regions. The values derived in Orion are consistent with the scenario of sulfur being undepleted in the envelope of this giant molecular cloud. Goicoechea & Cuadrado (2021) determined that sulfur is undepleted in the photon-dominated surface of OMC1, with a gas-phase carbon to sulfur abundance ratio of ~10. Observations of

the 158 µm line of C⁺ shows that this giant molecular cloud is surrounded by a partially ionized envelope that extends to the north (OMC 3) and south (OMC 4), probing that a high UV flux is still found in these suburb regions (Pabst et al. 2020), in line with our results.

In Table 4, we list the 20 most abundant sulfur-bearing species, in order of decreasing abundance, toward three regions that are representative of Taurus, Perseus, and Orion. In regions with $n_H > 10^5$ cm⁻³, less than 1% of sulfur is in the gas where the most abundant species is atomic sulfur with a fractional abundance ~4–10 times larger than the rest of the gaseous compounds. Contrary to carbon and oxygen, sulfur remains ionized until visual extinctions of ~4 mag. Atomic sulfur is then formed by recombination of S⁺ and/or dissociative recombinations of HCS⁺ and CS⁺. CS is the most abundant sulfur-bearing molecule at early time but does not become an important sulfur reservoir as it is CO for carbon. This is because the hydrogenation of CS followed by dissociative recombination of HCS⁺ produces much more S + CH than H + CS (see discussion in Vidal et al. 2017). At later times, in dense clouds, SO is expected to become the most abundant molecule. The amount of SO hence depends a lot on the amount of O₂ and OH in molecular clouds. If oxygen is heavily depleted, atomic sulfur will be the most dominant sulfur bearing species even at late times (Fuente et al. 2016). In the ice, m-H₂S, m-HS and m-S (m- indicates ice bulk molecule) are the most abundant ones. A wealth of organo-sulfur species are also found in solid phase as proposed by Laas & Caselli (2019) with lower abundances. At later times, ~1 Myr, the allotrope m-S₈ can form in the ice. As commented above, this stable allotrope can be considered as semi-refractory material. However, current chemical models predict the formation of S₈ mainly in dense regions and with very low abundances, clearly insufficient to explain the observed sulfur depletion. In the following section, we propose an alternative scenario to efficiently produce allotropes in the interface between the diffuse and translucent phase in dark clouds where Dₛ ~ 10 is already observed.

## 8. Discussion: The influence of grain charge distribution in sulfur depletion

In these dense and cold regions, the sum of the observed gas-phase abundances of S-species (the most abundant are SO, SO₂, H₂S, CS, HCS⁺, H₂CS, C₂S and C₃S) only accounts for <1% of the cosmic sulfur abundance (Vastel et al. 2018; Fuente et al. 2019). One could think that most of the sulfur is locked on the icy grain mantles, but surprisingly a similar trend is encountered within the solid phase, where s-OCS (Geballe et al. 1985; Palumbo et al. 1995) and s-SO₂ (Boogert et al. 1997) are the only sulfur-bearing species detected thus far, and only upper limits to the s-H₂S abundance have been derived (Smith 1991; Jiménez-Escobar & Muñoz Caro 2011). According to these data, the abundances of the observed icy species account for <5% of the total expected sulfur abundance. This means that 90–95% of the sulfur is missing in our counting. It has been suggested that this so-called depleted sulfur may be locked in hitherto undetected reservoirs in gas and icy grain mantles, or as refractory material. In particular, laboratory experiments and theoretical work show that sulfur allotropes, such as S₈, could be an important reservoir (Wakelam et al. 2004; Jiménez-Escobar et al. 2012; Shingledecker et al. 2020; Cazaux et al. 2022). The sublimation temperature of this allotrope is >500 K and can be considered as semi-refractory material.

We have determined the sulfur elemental abundance, that is, the amount of sulfur atoms in volatiles in starless cores located
<table><caption>**Table 4.** Most abundant volatile sulfur-bearing species<sup>(1)</sup>.</caption><thead><tr><th>B213-C12-8</th><th>B213-C12-1</th></tr></thead><tbody><tr><td><table><tr><td>n<sub>H</sub> = 3.16 &times; 10<sup>4</sup> cm<sup>-3</sup></td></tr><tr><td>D<sub>S</sub> ~ 65</td></tr><tr><td>97.7% of S atoms in refractories</td></tr><tr><td>1.21% of S atoms in gas</td></tr><tr><td>0.32% of S atoms in ice</td></tr></table></td><td><table><tr><td>n<sub>H</sub> = 3.16 &times; 10<sup>5</sup> cm<sup>-3</sup></td></tr><tr><td>D<sub>S</sub> ~ 65</td></tr><tr><td>97.7% of S atoms refractories</td></tr><tr><td>0.06% of S atoms gas</td></tr><tr><td>1.47% of S atoms ice</td></tr></table></td></tr><tr><td>Top-20: g-S, m-HS, m-H<sub>2</sub>S, m-NS, g-CS, m-SO, g-SO, m-OCS, g-S<sup>+</sup>, g-H<sub>2</sub>S, s-HS, s-H<sub>2</sub>S, s-NS, g-OCS, s-SO, g-H<sub>2</sub>CS, s-OCS, g-C<sub>3</sub>S, g-SO<sub>2</sub>, g-C<sub>2</sub>S</td><td>Top-20: m-HS, m-H<sub>2</sub>S, m-NS, m-SO, m-S, g-S, s-HS, s-H<sub>2</sub>S, s-NS, g-H<sub>2</sub>CS, s-SO, g-CS, m-CH<sub>3</sub>SH, m-OCS, m-CH<sub>3</sub>S, m-H<sub>2</sub>CS, g-H<sub>2</sub>S, s-S, m-NH<sub>2</sub>CH<sub>2</sub>SH, m-CH<sub>2</sub>SH</td></tr><tr><td>Barnard 1b (+140,0)</td><td>Barnard 1b (0,0)</td></tr><tr><td><table><tr><td>n<sub>H</sub> = 3.16 &times; 10<sup>4</sup> cm<sup>-3</sup></td></tr><tr><td>D<sub>S</sub> ~ 8</td></tr><tr><td>87.3% of S atoms in refractories</td></tr><tr><td>9.60% of S atoms in gas</td></tr><tr><td>3.06% of S atoms in ice</td></tr></table></td><td><table><tr><td>n<sub>H</sub> = 3.16 &times; 10<sup>5</sup> cm<sup>-3</sup></td></tr><tr><td>D<sub>S</sub> ~ 14</td></tr><tr><td>92.7% of S atoms in refractories</td></tr><tr><td>0.22% of S atoms in gas</td></tr><tr><td>7.11% of S atoms in ice</td></tr></table></td></tr><tr><td>Top-20: g-S, g-CS, m-HS, m-H<sub>2</sub>S, m-NS, m-SO, m-OCS, g-S<sup>+</sup>, g-SO, s-HS, s-H<sub>2</sub>S, s-NS, s-SO, g-OCS, s-OCS, m-NH<sub>2</sub>CH<sub>2</sub>SH, m-SO<sub>2</sub>, g-C<sub>3</sub>S, g-H<sub>2</sub>CS, g-C<sub>2</sub>S</td><td>Top-20: m-HS, m-H<sub>2</sub>S, m-NS, m-SO, m-S, g-S, s-HS, s-H<sub>2</sub>S, s-NS, m-CH<sub>3</sub>SH, m-OCS, m-CH<sub>3</sub>S, s-SO, g-H<sub>2</sub>CS, m-NH<sub>2</sub>CH<sub>2</sub>SH, m-H<sub>2</sub>CS, g-SO, g-CS, g-H<sub>2</sub>S, m-SO<sub>2</sub></td></tr><tr><td>Ori-C3-11</td><td>Ori-C3-1</td></tr><tr><td><table><tr><td>n<sub>H</sub> = 1.0 &times; 10<sup>5</sup> cm<sup>-3</sup></td></tr><tr><td>D<sub>S</sub> ~ 2</td></tr><tr><td>50.0% of S atoms in refractories</td></tr><tr><td>0.29% of S atoms in gas</td></tr><tr><td>49.71% of S atoms in ice</td></tr></table></td><td><table><tr><td>n<sub>H</sub> = 3.16 &times; 10<sup>5</sup> cm<sup>-3</sup></td></tr><tr><td>D<sub>S</sub> ~ 16</td></tr><tr><td>93.72% of S atoms in refractories</td></tr><tr><td>0.04% of S atoms in gas</td></tr><tr><td>6.24% of S atoms in ice</td></tr></table></td></tr><tr><td>Top-20: m-H<sub>2</sub>S, m-OCS, s-H<sub>2</sub>S, m-H<sub>2</sub>S<sub>3</sub>, m-SO, m-SO<sub>2</sub>, s-OCS, m-HS, s-SO<sub>2</sub>, s-H<sub>2</sub>S<sub>3</sub>, g-S, m-HSSH, m-CS<sub>2</sub>, m-HSS, s-SO, g-OCS, m-S, s-HS, g-H<sub>2</sub>S, s-HSSH</td><td>Top-20: m-H<sub>2</sub>S, m-SO, s-H<sub>2</sub>S, m-HS, m-HSSH, m-NS, m-CS<sub>2</sub>, m-HSS, g-S, m-OCS, m-H<sub>2</sub>S<sub>3</sub>, m-CH<sub>3</sub>S, g-H<sub>2</sub>S, m-SO<sub>2</sub>, s-SO, s-HS, m-H<sub>2</sub>CCS, m-H<sub>2</sub>C<sub>3</sub>S, m-S<sub>8</sub>, s-HSSH</td></tr></tbody></table>

**Notes.** <sup>(1)</sup>Notation: gas (g-), ice surface (s-) and ice bulk (m-)species.

in different environments. An updated chemical network (Laas & Caselli 2019; Navarro-Almaida et al. 2020; Bulut et al. 2021) has been used to determine sulfur depletion in the 244 positions forming the GEMS molecular database. We find that sulfur depletion depends on the star formation activity in the neighborhood. While sulfur depletions of a factor of 10–20 are needed to explain the observations in Taurus and Perseus, the abundances of the sulfur-bearing species in the outer parts of Orion are well explained assuming the cosmic value. Within Taurus and Perseus, we detect an additional trend with sulfur depletion increasing with density. We would like to remind that these conclusions are based on the state-of-the-art knowledge of the sulfur chemistry. So what we have determined is the fraction of sulfur that does not participate to the volatile chemistry as we know it nowadays. Still, there could be volatile species missing from our models and locking the major fraction of sulfur. In the following, we discuss a possible explanation for the lack of sulfur atoms in volatiles, within the state-of-the-art knowledge.

One possibility is that sulfur atoms adsorb on the grain surfaces and form large allotropes in the cloud surface. The evaporation temperature of large allotropes like S<sub>8</sub> is of hundreds of K, therefore these compounds would remain on grain surfaces even in the extreme conditions of a hot core. This scenario is supported by recent theoretical Monte Carlo simulations carried out by Cazaux et al. (2022), which showed that S<sup>+</sup> ions frozen onto grain surfaces could produce large allotropes in a few times 10<sup>4</sup> but, in order to promote the formation of allotropes instead

of sulfur hydrides, one needs to assume that the sticking coefficient of S<sup>+</sup> is higher than that of H. Based on the treatment developed by Umebayashi & Nakano (1980) and Draine & Sutin (1987) for collisions with charged grains, Ruffle et al. (1999) proposed that the sticking coefficient of positive ions (C<sup>+</sup>, S<sup>+</sup>, Na<sup>+</sup>, ...) increases as A = (1 + 167/T) relative to that of neutrals (A = 1) in regions where the grains are negatively charged. Contrary to other elements like C<sup>+</sup>, sulfur would remain ionized in the translucent medium until A<sub>V</sub> ~ 4–7 mag where the grains are expected to be negatively charged, thus increasing its adsorption to the grain surface. With this same approach to describe the adsorption of ions onto grain surfaces, Laas & Caselli (2019) and Shingledecker et al. (2020) were able to explain gas-phase abundances of most compounds in dark clouds, and predicted that organo-sulfur compounds and S<sub>8</sub> could be an important sulfur reservoir in solid phase. This interpretation is, however, controversial since the sticking efficiency of positively charged atoms on a negatively charged grain is not fully understood, yet. It would depend on the details of the recombination process and the possibility of chemisorption of the neutral atom. Watson & Salpeter (1972) concluded for C<sup>+</sup> that the rate of sticking could be as much as ~4 times larger than for neutral species in the most favorable case, but also much smaller, <1, if the energy liberated during the electron recombination is used for the neutral atom to leave the grain surface. The case of sulfur should be more favorable for accretion than carbon since the energy released during the recombination would be smaller. In the following, we make
<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Comparison between the spatial distribution of grain charges and sulfur cation in a molecular cloud with similar physical conditions to Taurus. Left: abundance of S⁺ as a function of the visual extinction assuming the physical parameters derived in this work for Taurus (χ<sub>UV</sub> = 5, ζ<sub>H₂</sub> = 10<sup>-16</sup> s<sup>-1</sup>, T = 15 K, t = 0.1 Myr). Vertical lines indicate A<sub>V</sub> = 1.5, 3.5, and 5.5 mag. The horizontal blue line shows X(S<sup>+</sup>) = 10<sup>-7</sup>. Two values of density are considered, n<sub>H</sub> = 3.16 × 10<sup>3</sup> cm<sup>-3</sup> (black) and n<sub>H</sub> = 3.16 × 10<sup>4</sup> cm<sup>-3</sup> (red), which are representative of the translucent and dense phases (Fuente et al. 2019). Right: grain charge distribution for silicates with a size of 0.1 µm at A<sub>V</sub> = 1.5, 3.5, and 5.5 mag assuming the physical conditions described in the left panel. In these calculations we have assumed undepleted sulfur.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Comparison between the spatial distribution of grain charges and sulfur cation in a molecular cloud with similar physical conditions to Orion. Left: Abundance of S⁺ as a function of the visual extinction assuming the physical parameters derived for Orion (χ<sub>UV</sub> = 50, T = 25 K) and t = 0.1 Myr. We consider two values of density, n<sub>H</sub> = 3.16×10<sup>4</sup> cm<sup>-3</sup> (black) and n<sub>H</sub> = 3.16×10<sup>5</sup> cm<sup>-3</sup> (red). Continuous and dashed lines show the results for ζ<sub>H₂</sub> = 10<sup>-17</sup> s<sup>-1</sup> and ζ<sub>H₂</sub> = 10<sup>-16</sup> s<sup>-1</sup>, respectively. Vertical lines indicate A<sub>V</sub> = 1.5, 3.5, and 5.5 mag. The horizontal blue line shows X(S<sup>+</sup>)=10<sup>-7</sup>. Right: Grain charge distribution for silicates with a size of 0.1 µm at A<sub>V</sub> = 1.5, 3.5, and 5.5 mag assuming the physical conditions described in the left panel. In these calculations we have assumed undepleted sulfur. Continuous black and dashed blue lines show the results for ζ<sub>H₂</sub> = 10<sup>-17</sup> s<sup>-1</sup> and ζ<sub>H₂</sub> = 10<sup>-16</sup> s<sup>-1</sup>, respectively.</figcaption></figure>

some simple calculations to test whether there is some correla-
tion between the grain charge distribution and sulfur depletion
in our sample. Although not conclusive, this correlation would
support the role of the grain charge in sulfur depletion.

Assuming that the main sulfur depletion is due to the accre-
tion of S⁺ on the negatively charged grain surfaces, sulfur
depletion would depend on the amount of S⁺ available, the
grain charge distribution, and the density in the translucent enve-
lope. We have calculated the abundance of S⁺ and the grain
charge distribution as a function of visual extinction in Taurus
and Orion. The S⁺ abundance has been calculated using our
chemical network and the physical parameters derived in this
work and undepleted sulfur abundance. For Taurus, we have
adopted ζ<sub>H₂</sub> = 10<sup>-16</sup> s<sup>-1</sup>, G₀=9 in unit of Habing field (χ=5 in
unit of Draine field), T = 15 K, and two values of density,
n<sub>H</sub> = 3.16 × 10<sup>3</sup> cm<sup>-3</sup> and n<sub>H</sub> = 3.16 × 10<sup>4</sup> cm<sup>-3</sup>. These values are
representative of the translucent and dense gas phases in TMC 1
(Fuente et al. 2019). All the sulfur is in the form of S⁺ at low
visual extinctions, A<sub>V</sub> ≤ 1.5 mag (see Fig. 7). At higher visual
extinctions, the S⁺ abundance decreases but remains >10<sup>-6</sup> until
A<sub>V</sub> ~ 3.5–5.5 mag.

The case of Orion is more complex, since our fitting revealed
different physical conditions in each of the three cuts observed.
While the majority of the positions in OMC 4 are fitted with

ζ<sub>H₂</sub> = 10<sup>-17</sup> s<sup>-1</sup>, more than 40% of the positions in OMC 3
are better fitted with ζ<sub>H₂</sub> = 5 × 10<sup>-17</sup> s<sup>-1</sup> (see Fig. 5). The
same problem appears with the chemical age, while early-time
chemistry provides the best fit in OMC 3, late-time chemistry
better fits the positions in OMC 2 and OMC 4. It seems clear,
however, that the incident UV flux as well as the mean den-
sity is high in all the cuts of this massive star-forming region
(Rodríguez-Baras et al. 2021). Therefore, we fixed the values
of G₀=90 (χ=50) and T = 25 K, and considered four sets
of physical conditions in our calculations: 1) ζ<sub>H₂</sub> = 10<sup>-17</sup> s<sup>-1</sup>,
n<sub>H</sub>=3.16 × 10<sup>4</sup> cm<sup>-3</sup>; 2) ζ<sub>H₂</sub> = 10<sup>-17</sup> s<sup>-1</sup>, n<sub>H</sub>=3.16 × 10<sup>5</sup> cm<sup>-3</sup>;
3) ζ<sub>H₂</sub> = 10<sup>-16</sup> s<sup>-1</sup>, n<sub>H</sub>=3.16 × 10<sup>4</sup> cm<sup>-3</sup>; and 4) ζ<sub>H₂</sub> = 10<sup>-16</sup> s<sup>-1</sup>,
n<sub>H</sub>=3.16 × 10<sup>5</sup> cm<sup>-3</sup>. Moreover, we repeated the calculations
for t = 0.1 Myr and t = 1 Myr. Assuming high density,
n<sub>H</sub>=3.16 × 10<sup>5</sup> cm<sup>-3</sup>, the S<sup>+</sup> abundance is low, X(S<sup>+</sup>) < 10<sup>-7</sup>,
even at a low extinction of A<sub>V</sub> ~ 3.5 mag, regardless of the other
parameters (see Figs. 8 and 9). In the case of moderate density,
n<sub>H</sub>=3.16 × 10<sup>4</sup> cm<sup>-3</sup>, we obtain high S<sup>+</sup> abundance only in the
case of ζ<sub>H₂</sub> = 10<sup>-16</sup> s<sup>-1</sup> and early-time chemistry.

The other essential ingredient in our interpretation is the
grain charge. The charging of dust grains in the interstel-
lar medium is governed by the equilibrium among different
processes. On the one hand, a dust grain may acquire a posi-
tive charge due to photoelectric emission of electrons induced
<figure><img src="image_10.png" /><figcaption>**Fig. 9.** Same as Fig. 8, but for t = 1 Myr.</figcaption></figure>

either by UV radiation (Draine 1978) or by fluorescence of H₂ (Ivlev et al. 2015) depending on the intensity of the local radiation field and on the gas density; additional positive charging of the grain is enhanced by the accretion of free ions in the plasma (Spitzer 1941), although this phenomenon plays a minor role in dense environments where the mobility of the species is drastically reduced. On the other hand, negative grain charging can be produced by the accretion of free electrons in the plasma (Spitzer 1941) independently on their origin, although the contribution of secondary electrons arising from the cosmic-ray population (Draine & Salpeter 1979) is of minor importance (Ivlev et al. 2015). For this work, we have calculated the grain charge distributions in a similar way to Ibáñez-Méjía et al. (2019) for a population of silicate grains of 0.1 µm; more details on the computation of the distribution can be found in Appendix B. Assuming the same gas densities and radiation field intensities as above, and the abundances for ions and electron density predicted by our chemical calculations, we obtain the grain charge distributions shown in Figs. 7–9. The grain charge distributions have been calculated for $A_V = 1.5$ mag, 3.5 mag, and 5.5 mag. In all cases, the grain net charge centroid is close to ~0 or positive in the outer part of the molecular cloud ($A_V \sim 1.5$ mag), and changes to negative at higher extinction. In the case of Taurus, the net charge of grains is negative for $A_V \sim 3.5$ mag where $X(S^+) > 10^{-6}$, thus promoting sulfur depletion and the formation of sulfur chains as described by Cazaux et al. (2022). On the contrary, in Orion the centroid of the grain charge distribution remains close to ~0 until $A_V = 5.5$ mag where the abundance of S⁺ is already very low, especially for $n(H_2) \sim 10^5$ cm⁻³. This low abundance of negatively charged grains in the cloud envelope could contribute to maintain most of the sulfur atoms in gas phase. These calculations are very simple as long as only one grain size and composition are taken into account. One could think that the inclusion of the grain charges in chemical models would also affect other important species as those formed from C⁺. However, the situation could be different for carbon. Because of its lower ionization potential, sulfur remains ionized until visual extinctions of ~4 mag where grains are negatively charged under dark cloud conditions. However, carbon would be neutral in this region, thus avoiding an enhanced sticking probability. Even though, a more rigorous treatment is necessary to fully disentangle the role of the grain charge in the gas chemistry.

Finally, we cannot discard the effect that the shocks produced by the expansion of the HII region and those associated with the bipolar outflows of the population of young stellar objects in the region could produce in the sulfur elemental abundance in very active star-forming regions such as Orion A (Gustafsson et al. 2003; Colgan et al. 2007; Berné et al. 2014;

Feddersen et al. 2020). The sulfur depletion is known to be moderate, of the order of ~10, in the shocks associated with bipolar outflows driven by low-mass stars (Anderson et al. 2013; Holdship et al. 2019; Feng et al. 2020), still larger than the almost null depletion we have found in outer envelope of Orion A. Moreover, our cuts have been selected to avoid the jets and molecular outflows that could affect the molecular chemistry hindering the chemical composition of the pristine gas. However, we cannot discard the possibility that previous shocks have released sulfur from the refractory grain cores, that is expelled and turbulently mixed with the molecular gas, thus enriching the sulfur content in the environment.

## 9. Summary and conclusions

This work uses the GEMS molecular database to derive the sulfur elemental abundance in a wide sample of starless cores located in the nearby star-forming regions Taurus, Perseus, and Orion. These regions have different degrees of star formation activity, and therefore different physical conditions, providing a possibility to explore the effect of environment. In order to derive the sulfur elemental abundance we have modeled the abundances of 9 species using a state-of-the-art chemical code with an updated sulfur network. In a first step, in order to explore a wide range of physical conditions, we only considered three values of $[S/H] = 1.5 \times 10^{-5}$ ($D_S = 1$), $1.5 \times 10^{-6}$ ($D_S = 10$), and $8 \times 10^{-8}$ ($D_S = 187$). Our results can be summarized as follows:

- Most of the positions in Taurus are best fitted assuming early-time chemistry ($t = 0.1$ Myr). We find a large dispersion in the values of $\zeta_{H_2}$ with a distribution peaking at $\zeta_{H_2} \sim 1 \times 10^{-16}$ s⁻¹, and without any clear trend with visual extinction and/or molecular hydrogen density. However, we do find a trend with the environment. Most positions in TMC 1 and B 213-N are fitted with $\zeta_{H_2} > 10^{-16}$ s⁻¹, while the positions in B 213-S are fitted with $\zeta_{H_2} < 10^{-16}$ s⁻¹, consistent with the idea of a harsh environment in B 213-N due to the presence of young (proto-)stars. Regarding sulfur elemental abundance, we find $D_S = 10$ in TMC 1 and B 213-N, and $D_S = 187$ in B 213-S. In addition to environment, we find some correlation of $D_S$ with density, with depletion increasing with increasing density.

- Similarly to Taurus, essentially most of the positions in Perseus are better fitted with a chemical age of $t = 0.1$ Myr. There is also large dispersion in the values of $\zeta_{H_2}$ with a distribution peaking at $\zeta_{H_2} \sim (5-10) \times 10^{-17}$ s⁻¹, without any clear trend with visual extinction and density. Regarding the sulfur elemental abundance, most of the positions are fitted with $D_S = 10$. Contrary to Taurus, we find some positions
toward which the best fit corresponds to undepleted sulfur abundance. These positions are found toward IC 348, NGC 1333, and B5. IC 348 and NGC 1333 are the nearest regions to star clusters. NGC 1333 hosts a large population of Class 0 and I protostars. The low values of sulfur depletion in these regions are very likely related to the high local star formation activity.

- We have observed three cuts in the massive star-forming region Orion A: Ori-C1 in OMC 3, Ori-C2 in OMC 4, and Ori-C3 in OMC 2. Although the number of positions observed in Orion A is low, it seems clear that this region presents clear differences relative to Taurus and Perseus. Although the cut Ori-C1 in OMC 3 is fitted with early time chemistry ($t = 0.1$ Myr), the cut Ori-C2 in OMC 4 is fitted with a chemical age of $t=10$ Myr, and the cut Ori-C3 in OMC 2, with $t = 1$ Myr. We find also dispersion in the values of $\zeta_{H_2}$. While the majority of the positions in OMC 4 are fitted with $\zeta_{H_2} = 10^{-17} s^{-1}$, more than 40% of the positions in OMC 3 are better fitted with $\zeta_{H_2} = 5 \times 10^{-17} s^{-1}$. Regarding sulfur, a significant fraction of the positions (~40%), especially those at low visual extinction, are best fitted assuming undepleted sulfur abundance.

In order to have a deeper insight into the sulfur elemental abundance in the 29 studied starless cores, we have run a model with a finer grid with the values of [S/H] varying in steps of a factor of 2. In addition, we fix the values of $t$ and $\zeta_{H_2}$ to diminish possible degeneracies. Furthermore, in each cut we considered separately the value of [S/H] derived in the visual extinction peak from that in the envelope. We found sulfur depletion of a factor >20 is always needed toward the visual extinction peaks. However, lower values of sulfur depletion are measured in the envelopes of regions with enhanced star formation activity, and especially in Orion A.

In addition, we have explored the possibility that the observed variations of sulfur depletion is consequence of the influence of the local physical conditions on the abundance of $S^+$ and the grain charges. In the case of Taurus, we find that grains become negatively charged at a visual extinction of $A_V \sim 3.5$ mag. In this region, the abundance of $S^+$ might be $>10^{-6}$, and the electrostatic attraction between $S^+$ and negatively charged grains increases the accretion rate of $S^+$ on dust, and leads to the formation of S chains (Cazaux et al. 2022), which would enhance sulfur depletion. This could explain that sulfur depletion is already significant, ~20, in the translucent region of Taurus filaments. In the case of Orion, the net charge of grains is close to 0 in the region where the abundance of $S^+$ is high, which could slow down, or even suppress, the depletion of sulfur in the cloud envelope. Therefore, our calculations suggest that grain charge could play an important role to explain the observed differences in the sulfur depletion. However, we have assumed a single composition (silicates) and size (0.1 µm) for grains, which is a simple scenario. A full distribution of grain sizes and compositions, as well as possible differences of the grain sizes along the cloud, should be taken into account for a more rigorous approach to quantify this effect.

Summarizing, our data show that the environment is the driving agent of the sulfur depletion in molecular clouds. The mechanisms responsible for this differentiation are not known in full detail. The influence of the grain charges on the chemistry, not considered in most chemical models, is very likely one of the causes. We cannot discard that the shocks associated with massive star formation could erode the grain cores, thus contributing to enhance the sulfur elemental abundance in the hosting molecular cloud. Additional observations of massive

star-forming regions using large interferometers would be desirable to have a global and reliable picture of the sulfur chemistry in these intriguing regions.

**Acknowledgements.** We thank the Spanish MICIN for funding support from PID2019-106235GB-I00. I.J.-S. acknowledges support from grant No. PID2019-105552RB-C41 by the Spanish Ministry of Science and Innovation/State Agency of Research MCIN/AEI/10.13039/501100011033. J.R.G. acknowledges Spanish MICIN for support from PID2019-106110GB-I00. V.W. acknowledges the CNRS program "Physique et Chimie du Milieu Interstellaire" (PCMI) co-funded by the Centre National d'Etudes Spatiales (CNES). L.B.-A. acknowledges the receipt of a Margarita Salas postdoctoral fellowship from Universidad Complutense de Madrid (CT31/21), funded by "Ministerio de Universidades" with Next Generation EU funds. D.N.-A. acknowledges the Fundación Ramón Areces for funding support through their international postdoc grant program. The authors also thank the anonymous referee for their interesting suggestions.

## References

Agúndez, M., & Wakelam, V. 2013, Chem. Rev., 113, 8710

Agúndez, M., Marcelino, N., Cernicharo, J., Roueff, E., & Tafalla, M. 2019, A&A, 625, A147

Aharonian, F., Yang, R., & de Oña Wilhelmi, E. 2019, Nat. Astron., 3, 561

Anderson, D. E., Bergin, E. A., Maret, S., & Wakelam, V. 2013, ApJ, 779, 141

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

Anglada, G., Estalella, R., Pastor, J., Rodriguez, L. F., & Haschick, A. D. 1996, ApJ, 463, 205

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

Bachiller, R., & Cernicharo, J. 1984, A&A, 140, 414

Bally, J., Langer, W. D., Stark, A. A., & Wilson, R. W. 1987, ApJ, 312, L45

Bayet, E., Aladro, R., Martín, S., Viti, S., & Martín-Pintado, J. 2009, ApJ, 707, 126

Berné, O., Marcelino, N., & Cernicharo, J. 2014, ApJ, 795, 13

Blake, G. A., Mundy, L. G., Carlstrom, J. E., et al. 1996, ApJ, 472, L49

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

Bouscasse, L., Csengeri, T., Belloche, A., et al. 2022, A&A, 662, A32

Bracco, A., Palmeirim, P., André, P., et al. 2017, A&A, 604, A52

Bronfman, L., Nyman, L. A., & May, J. 1996, A&AS, 115, 81

Bulut, N., Roncero, O., Aguado, A., et al. 2021, A&A, 646, A5

Cambrésy, L. 1999, A&A, 345, 965

Caselli, P., Walmsley, C. M., Terzieva, R., & Herbst, E. 1998, ApJ, 499, 234

Caselli, P., Walmsley, C. M., Zucconi, A., et al. 2002, ApJ, 565, 344

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

Cecchi-Pestellini, C., & Aiello, S. 1992, MNRAS, 258, 125

Cernicharo, J., & Guelin, M. 1987, A&A, 176, 299

Colgan, S. W. J., Schultz, A. S. B., Kaufman, M. J., Erickson, E. F., & Hollenbach, D. J. 2007, ApJ, 671, 536

Colzi, L., Sipilä, O., Roueff, E., Caselli, P., & Fontani, F. 2020, A&A, 640, A51

Curtis, E. I., & Richer, J. S. 2011, MNRAS, 410, 75

Da Rio, N., Robberto, M., Hillenbrand, L. A., Henning, T., & Stassun, K. G. 2012, ApJ, 748, 14

Dagdigian, P. J. 2020, MNRAS, 494, 5239

Davis, C. J., Chrysostomou, A., Hatchell, J., et al. 2010, MNRAS, 405, 759

Denis-Alpizar, O., Stoecklin, T., Guilloteau, S., & Dutrey, A. 2018, MNRAS, 478, 1811

Draine, B. T. 1978, ApJS, 36, 595

Draine, B. T., & Salpeter, E. E. 1979, ApJ, 231, 438

Draine, B. T., & Sutin, B. 1987, ApJ, 320, 803

Dumouchel, F., Klos, J., & Lique, F. 2011, Phys. Chem. Chem. Phys. (Incorp. Faraday Trans.), 13, 8204

Dutrey, A., Wakelam, V., Boehler, Y., et al. 2011, A&A, 535, A104

Enoch, M. L., Young, K. E., Glenn, J., et al. 2006, ApJ, 638, 293

Esplugues, G., Fuente, A., Navarro-Almaida, D., et al. 2022, A&A, 662, A52

Feddersen, J. R., Arce, H. G., Kong, S., et al. 2020, ApJ, 896, 11

Fehér, O., Tóth, L. V., Ward-Thompson, D., et al. 2016, A&A, 590, A75

Feng, S., Codella, C., Ceccarelli, C., et al. 2020, ApJ, 896, 37

Fitzpatrick, E. L., & Massa, D. 2007, ApJ, 663, 320

Flower, D. R. 1999, MNRAS, 305, 651

Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306

Friesen, R. K., Pineda, J. E., co-PIs, et al. 2017, ApJ, 843, 63

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

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

Fuente, A., Navarro, D. G., Caselli, P., et al. 2019, A&A, 624, A105

Fuente, A., Treviño-Morales, S. P., Alonso-Albi, T., et al. 2021, MNRAS, 507, 1886
Garrod, R. T., Wakelam, V., & Herbst, E. 2007, A&A, 467, 1103

Gatley, I., Becklin, E. E., Mattews, K., et al. 1974, ApJ, 191, L121

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

Gerin, M., Pety, J., & Goicoechea, J. R. 2009, ASP Conf. Ser., 417, 165

Gerin, M., Ruaud, M., Goicoechea, J. R., et al. 2015, A&A, 573, A30

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

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

Goicoechea, J. R., Pety, J., Gerin, M., et al. 2006, A&A, 456, 565

Goicoechea, J. R., Santa-Maria, M. G., Bron, E., et al. 2019, A&A, 622, A91

Goicoechea, J. R., Aguado, A., Cuadrado, S., et al. 2021, A&A, 647, A10

Goldsmith, P. F., Heyer, M., Narayanan, G., et al. 2008, ApJ, 680, 428

Goodman, J., & Weare, J. 2010, Commun. Appl. Math. Comput. Sci., 5, 65

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

Green, S., & Chapman, S. 1978, ApJS, 37, 169

Gustafsson, M., Kristensen, L. E., Clénet, Y., et al. 2003, A&A, 411, 437

Hacar, A., Tafalla, M., Kauffmann, J., & Kovács, A. 2013, A&A, 554, A55

Hacar, A., Alves, J., Tafalla, M., & Goicoechea, J. R. 2017a, A&A, 602, L2

Hacar, A., Tafalla, M., & Alves, J. 2017b, A&A, 606, A123

Hacar, A., Tafalla, M., Forbrich, J., et al. 2018, A&A, 610, A77

Hacar, A., Bosman, A. D., & van Dishoeck, E. F. 2020, A&A, 635, A4

Hasegawa, T. I., & Herbst, E. 1993, MNRAS, 261, 83

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

Hatchell, J., Richer, J. S., Fuller, G. A., et al. 2005, A&A, 440, 151

Hatchell, J., Fuller, G. A., & Richer, J. S. 2007a, A&A, 472, 187

Hatchell, J., Fuller, G. A., Richer, J. S., Harries, T. J., & Ladd, E. F. 2007b, A&A, 468, 1009

Hernández Vera, M., Lique, F., Dumouchel, F., Hily-Blant, P., & Faure, A. 2017, MNRAS, 468, 1084

Hillenbrand, L. A. 1997, AJ, 113, 1733

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

Hocuk, S., Szűcs, L., Caselli, P., et al. 2017, A&A, 604, A58

Holdship, J., Jimenez-Serra, I., Viti, S., et al. 2019, ApJ, 878, 64

Ibáñez-Mejía, J. C., Walch, S., Ivlev, A. V., et al. 2019, MNRAS, 485, 1220

Ishii, S., Nakamura, F., Shimajiri, Y., et al. 2019, PASJ, 71, S9

Ivlev, A. V., Padovani, M., Galli, D., & Caselli, P. 2015, ApJ, 812, 135

Jenkins, E. B. 2009, ApJ, 700, 1299

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

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

Johnstone, D., & Bally, J. 1999, ApJ, 510, L49

Johnstone, D., & Bally, J. 2006, ApJ, 653, 383

Kirk, H., Johnstone, D., & Di Francesco, J. C. 2006, ApJ, 646, 1009

Kirk, H., Friesen, R. K., Pineda, J. E., et al. 2017, ApJ, 846, 144

Kong, S., Arce, H. G., Sargent, A. I., et al. 2019, ApJ, 882, 45

Kutner, M. L., Evans, N. J., I., & Tucker, K. D. 1976, ApJ, 209, 452

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

Lada, C. J., Alves, J., & Lada, E. A. 1996, AJ, 111, 1964

Lada, C. J., Muench, A. A., Haisch, Karl E., et al. 2000, AJ, 120, 3162

Launhardt, R., Evans, Neal J., I., Wang, Y., et al. 1998, ApJS, 119, 59

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

Le Gal, R., Öberg, K. I., Teague, R., et al. 2021, ApJS, 257, 12

Linke, R. A., & Goldsmith, P. F. 1980, ApJ, 235, 437

Lique, F., & Spielfiedel, A. 2007, A&A, 462, 1179

Liszt, H. S., Wilson, R. W., Penzias, A. A., et al. 1974, ApJ, 190, 557

Loison, J.-C., Wakelam, V., Gratier, P., et al. 2019, MNRAS, 485, 5777

Loison, J.-C., Wakelam, V., Gratier, P., & Hickson, K. M. 2020, MNRAS, 498, 4663

Lombardi, M., Bouy, H., Alves, J., & Lada, C. J. 2014, A&A, 566, A45

Luhman, K. L., Stauffer, J. R., Muench, A. A., et al. 2003, ApJ, 593, 1093

Luhman, K. L., Mamajek, E. E., Allen, P. R., & Cruz, K. L. 2009, ApJ, 703, 399

Malinen, J., Juvela, M., Rawlings, M. G., et al. 2012, A&A, 544, A50

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

Marsh, K. A., Griffin, M. J., Palmeirim, P., et al. 2014, MNRAS, 439, 3683

Mckee, C. F. 1989, ApJ, 345, 782

Minissale, M., Dulieu, F., Cazaux, S., & Hocuk, S. 2016, A&A, 585, A24

Monsch, K., Pineda, J. E., Liu, H. B., et al. 2018, ApJ, 861, 77

Muench, A. A., Lada, E. A., Lada, C. J., & Alves, J. 2002, ApJ, 573, 366

Nakamura, F., Ishii, S., Dobashi, K., et al. 2019, PASJ, 71, S3

Narayanan, G., Heyer, M. H., Brunt, C., et al. 2008, ApJS, 177, 341

Navarro-Almaida, D., Le Gal, R., Fuente, A., et al. 2020, A&A, 637, A39

Navarro-Almaida, D., Fuente, A., Majumdar, L., et al. 2021, A&A, 653, A15

Neufeld, D. A., & Wolfire, M. G. 2017, ApJ, 845, 163

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

Onishi, T., Mizuno, A., Kawamura, A., Ogawa, H., & Fukui, Y. 1996, ApJ, 465, 815

Onishi, T., Mizuno, A., Kawamura, A., Tachihara, K., & Fukui, Y. 2002, ApJ, 575, 950

Ortiz-León, G. N., Loinard, L., Dzib, S. A., et al. 2018, ApJ, 865, 73

Pabst, C. H. M., Goicoechea, J. R., Teyssier, D., et al. 2020, A&A, 639, A2

Padoan, P., Cambrésy, L., & Langer, W. 2002, ApJ, 580, L57

Padovani, M., Hennebelle, P., & Galli, D. 2013, A&A, 560, A114

Padovani, M., Marcowith, A., Sánchez-Monge, Á., Meng, F., & Schilke, P. 2019, A&A, 630, A72

Palmeirim, P., André, P., Kirk, J., et al. 2013, A&A, 550, A38

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

Pankonin, V., & Walmsley, C. M. 1978, A&A, 64, 333

Pety, J., Teyssier, D., Fossé, D., et al. 2005, A&A, 435, 885

Pety, J., Guzmán, V. V., Orkisz, J. H., et al. 2017, A&A, 599, A98

Pineda, J. E., Caselli, P., & Goodman, A. A. 2008, ApJ, 679, 481

Pineda, J. E., Goodman, A. A., Arce, H. G., et al. 2010, ApJ, 712, L116

Pineda, J. E., Offner, S. S. R., Parker, R. J., et al. 2015, Nature, 518, 213

Prasad, S. S., & Tarafdar, S. P. 1983, ApJ, 267, 603

Punanova, A., Vasyunin, A., Caselli, P., et al. 2022, ApJ, 927, 213

Rebull, L. M., Padgett, D. L., McCabe, C. E., et al. 2010, ApJS, 186, 259

Ridge, N. A., Di Francesco, J., Kirk, H., et al. 2006, AJ, 131, 2921

Rivière-Marichalar, P., Fuente, A., Goicoechea, J. R., et al. 2019, A&A, 628, A16

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

Robberto, M., Soderblom, D. R., Bergeron, E., et al. 2013, ApJS, 207, 10

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

Sandell, G., & Knee, L. B. G. 2001, ApJ, 546, L49

Savage, C., Apponi, A. J., Ziurys, L. M., & Wyckoff, S. 2002, ApJ, 578, 211

Schmalzl, M., Kainulainen, J., Quanz, S. P., et al. 2010, ApJ, 725, 1327

Schnee, S., Caselli, P., Goodman, A., et al. 2007, ApJ, 671, 1839

Schnee, S., Enoch, M., Noriega-Crespo, A., et al. 2010, ApJ, 708, 127

Scourfield, M., Viti, S., García-Burillo, S., et al. 2020, MNRAS, 496, 5308

Shimajiri, Y., André, P., Palmeirim, P., et al. 2019, A&A, 623, A16

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

Shirley, Y. L., Evans, Neal J., I., Young, K. E., Knez, C., & Jaffe, D. T. 2003, ApJS, 149, 375

Smith, R. G. 1991, MNRAS, 249, 172

Spezzano, S., Fuente, A., Caselli, P., et al. 2022, A&A, 657, A10

Spitzer, L. Jr 1941, ApJ, 93, 369

Suri, S., Sánchez-Monge, Á., Schilke, P., et al. 2019, A&A, 623, A142

Tafalla, M., & Hacar, A. 2015, A&A, 574, A104

Tanabe, Y., Nakamura, F., Tsukagoshi, T., et al. 2019, PASJ, 71, S8

Taquet, V., Codella, C., De Simone, M., et al. 2020, A&A, 637, A63

Tatematsu, K., Umemoto, T., Kameya, O., et al. 1993, ApJ, 404, 643

Tercero, F., López-Pérez, J. A., Gallego, J. D., et al. 2021, A&A, 645, A37

Umebayashi, T., & Nakano, T. 1980, PASJ, 32, 405

van der Tak, F. F. S., Black, J. H., Schöier, F. L., Jansen, D. J., & van Dishoeck, E. F. 2007, A&A, 468, 627

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

Vera, M. H., Kalugina, Y., Denis-Alpizar, O., Stoecklin, T., & Lique, F. 2014, J. Chem. Phys., 140, 224302

Vidal, T. H. G., Loison, J.-C., Jaziri, A. Y., et al. 2017, MNRAS, 469, 435

Wakelam, V., Castets, A., Ceccarelli, C., et al. 2004, A&A, 413, 609

Wakelam, V., Herbst, E., Selsis, F., & Massacrier, G. 2006, A&A, 459, 813

Wakelam, V., Dartois, E., Chabot, M., et al. 2021, A&A, 652, A63

Wang, K. S., Bourke, T. L., Hogerheijde, M. R., et al. 2013, A&A, 558, A69

Warin, S., Castets, A., Langer, W. D., Wilson, R. W., & Pagani, L. 1996, A&A, 306, 935

Watson, W. D., & Salpeter, E. E. 1972, ApJ, 174, 321

Wilking, B. A., Meyer, M. R., Greene, T. P., Mikhail, A., & Carlson, G. 2004, AJ, 127, 1131

Wilson, T. L., & Rood, R. 1994, ARA&A, 32, 191

Wilson, R. W., Jefferts, K. B., & Penzias, A. A. 1970, ApJ, 161, L43

Wu, J., Evans, Neal J., I., Shirley, Y. L., & Knez, C. 2010, ApJS, 188, 313

Yang, B., Stancil, P. C., Balakrishnan, N., & Forrey, R. C. 2010, ApJ, 718, 1062

Yan, Q.-Z., Zhang, B., Xu, Y., et al. 2019, A&A, 624, A6

Yazidi, O., Ben Abdallah, D., & Lique, F. 2014, MNRAS, 441, 664

Zari, E., Lombardi, M., Alves, J., Lada, C. J., & Bouy, H. 2016, A&A, 587, A106

Zari, E., Brown, A. G. A., & de Zeeuw, P. T. 2019, A&A, 628, A123

Zhang, Z.-Y., Gao, Y., Henkel, C., et al. 2014, ApJ, 784, L31

Zhao, B., Caselli, P., Li, Z.-Y., et al. 2016, MNRAS, 460, 2050

Zhou, S., Wu, Y., Evans, Neal J., I., Fuller, G. A., & Myers, P. C. 1989, ApJ, 346, 168

Zinchenko, I., Mattila, K., & Toriseva, M. 1995, A&AS, 111, 95

Zucker, C., Speagle, J. S., Schlaflay, E. F., et al. 2019, ApJ, 879, 125

Zuckerman, B. 1973, ApJ, 183, 863
## Appendix A: Additional figures and tables

<table><caption>**Table A.1.** References for the collisional rate coefficients used for the volume density estimates.</caption><thead><tr><td>Molecule</td><td>Reference</td></tr></thead><tbody><tr><td>CS</td><td>Denis-Alpizar et al. (2018)</td></tr><tr><td><sup>13</sup>CO</td><td>Yang et al. (2010)</td></tr><tr><td>C<sup>18</sup>O</td><td>Yang et al. (2010)</td></tr><tr><td>HCO<sup>+</sup></td><td>Yazidi et al. (2014)</td></tr><tr><td>H<sup>13</sup>CO<sup>+</sup></td><td>Yazidi et al. (2014)</td></tr><tr><td>HC<sup>18</sup>O<sup>+</sup></td><td>Flower (1999)</td></tr><tr><td>HCS<sup>+</sup></td><td>Flower (1999)</td></tr><tr><td>H<sup>13</sup>CN</td><td>Vera et al. (2014); Hernández Vera et al. (2017)</td></tr><tr><td>SO</td><td>Lique & Spielfiedel (2007)</td></tr><tr><td><sup>34</sup>SO</td><td>Lique & Spielfiedel (2007)</td></tr><tr><td>HNC</td><td>Dumouchel et al. (2011)<br>Vera et al. (2014)</td></tr><tr><td>OCS</td><td>Green & Chapman (1978)</td></tr><tr><td>o-H<sub>2</sub>S</td><td>Dagdigian (2020)</td></tr></tbody></table>
<table><caption>**Table A.2.** Comparison of the mean value of [S/H] in each GEMS cut with that derived toward the visual extinction peak (coordinates in Table 1). The parameter *D*<sub>diff</sub> is the uncertainty as defined in Sect. 5.</caption><thead><tr><th>Region</th><th>Cloud</th><th>cut</th><th>Peak</th><th>D<sub>diff</sub><sup>peak</sup></th><th>Mean</th><th>D<sub>diff</sub><sup>mean</sup></th></tr></thead><tbody><tr><td rowspan="14">Taurus</td><td rowspan="3">TMC 1</td><td>CP</td><td>2.3 × 10<sup>-7</sup></td><td>0.39</td><td>7.8 × 10<sup>-7</sup></td><td>0.46</td></tr><tr><td>NH3</td><td>2.3 × 10<sup>-7</sup></td><td>0.16</td><td>3.1 × 10<sup>-6</sup></td><td>0.33</td></tr><tr><td>C</td><td>4.7 × 10<sup>-7</sup></td><td>0.48</td><td>5.6 × 10<sup>-7</sup></td><td>0.41</td></tr><tr><td rowspan="5">B 213-N</td><td>C1</td><td>1.1 × 10<sup>-6</sup></td><td>0.14</td><td>8.7 × 10<sup>-7</sup></td><td>0.29</td></tr><tr><td>C2</td><td>2.3 × 10<sup>-7</sup></td><td>0.34</td><td>3.6 × 10<sup>-7</sup></td><td>0.43</td></tr><tr><td>C5</td><td>4.7 × 10<sup>-7</sup></td><td>0.31</td><td>4.1 × 10<sup>-7</sup></td><td>0.55</td></tr><tr><td>C6</td><td>2.3 × 10<sup>-7</sup></td><td>0.27</td><td>1.2 × 10<sup>-6</sup></td><td>0.52</td></tr><tr><td>C7</td><td>2.3 × 10<sup>-7</sup></td><td>0.25</td><td>8.4 × 10<sup>-7</sup></td><td>0.22</td></tr><tr><td rowspan="4">B 213-S</td><td>C10</td><td>2.3 × 10<sup>-7</sup></td><td>0.32</td><td>2.7 × 10<sup>-7</sup></td><td>0.28</td></tr><tr><td>C12</td><td>2.3 × 10<sup>-7</sup></td><td>0.41</td><td>2.1 × 10<sup>-7</sup></td><td>0.38</td></tr><tr><td>C16</td><td>4.7 × 10<sup>-7</sup></td><td>0.38</td><td>2.2 × 10<sup>-7</sup></td><td>0.48</td></tr><tr><td>C17</td><td>4.7 × 10<sup>-7</sup></td><td>0.14</td><td>4.3 × 10<sup>-7</sup></td><td>0.29</td></tr><tr><td rowspan="16">Perseus</td><td rowspan="3">B5 (#79)</td><td>C1</td><td>4.7 × 10<sup>-7</sup></td><td>0.18</td><td>5.5 × 10<sup>-7</sup></td><td>0.28</td></tr><tr><td>C1</td><td>1.1 × 10<sup>-6</sup></td><td>0.16</td><td>6.0 × 10<sup>-7</sup></td><td>0.17</td></tr><tr><td>C1</td><td>2.3 × 10<sup>-7</sup></td><td>0.15</td><td>4.1 × 10<sup>-7</sup></td><td>0.48</td></tr><tr><td rowspan="5">L1448</td><td>C10</td><td>8.0 × 10<sup>-8</sup></td><td>0.20</td><td>1.1 × 10<sup>-7</sup></td><td>0.35</td></tr><tr><td>C1</td><td>9.4 × 10<sup>-7</sup></td><td>0.48</td><td>3.6 × 10<sup>-7</sup></td><td>0.48</td></tr><tr><td>C2</td><td>2.3 × 10<sup>-7</sup></td><td>0.28</td><td>2.9 × 10<sup>-7</sup></td><td>0.44</td></tr><tr><td>C3</td><td>4.7 × 10<sup>-7</sup></td><td>0.16</td><td>8.2 × 10<sup>-7</sup></td><td>0.44</td></tr><tr><td>C4</td><td>1.1 × 10<sup>-6</sup></td><td>0.09</td><td>7.7 × 10<sup>-7</sup></td><td>0.35</td></tr><tr><td rowspan="4">IC348</td><td>C5</td><td>2.3 × 10<sup>-7</sup></td><td>0.23</td><td>1.8 × 10<sup>-6</sup></td><td>0.94</td></tr><tr><td>C6</td><td>2.3 × 10<sup>-7</sup></td><td>0.23</td><td>2.7 × 10<sup>-7</sup></td><td>0.66</td></tr><tr><td>C7</td><td>4.7 × 10<sup>-7</sup></td><td>0.28</td><td>4.7 × 10<sup>-7</sup></td><td>0.28</td></tr><tr><td>C1</td><td>1.1 × 10<sup>-6</sup></td><td>0.43</td><td>4.3 × 10<sup>-6</sup></td><td>0.34</td></tr><tr><td rowspan="3">NGC 1333</td><td>C2</td><td>4.7 × 10<sup>-7</sup></td><td>0.28</td><td>5.4 × 10<sup>-7</sup></td><td>0.30</td></tr><tr><td>C3</td><td>9.4 × 10<sup>-7</sup></td><td>0.22</td><td>6.5 × 10<sup>-7</sup></td><td>0.42</td></tr><tr><td>B1b</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td rowspan="3">Orion</td><td rowspan="3">A</td><td>C1</td><td>2.3 × 10<sup>-7</sup></td><td>0.51</td><td>4.8 × 10<sup>-6</sup></td><td>0.51</td></tr><tr><td>C2</td><td>7.5 × 10<sup>-6</sup></td><td>0.06</td><td>5.1 × 10<sup>-6</sup></td><td>0.28</td></tr><tr><td>C3</td><td>9.4 × 10<sup>-7</sup></td><td>0.10</td><td>3.0 × 10<sup>-6</sup></td><td>0.23</td></tr></tbody></table>
<figure><img src="image_11.png" /><figcaption>**Fig. A.1.** Corner diagrams performed by selecting the models with $A_V = 11.5$ mag from the output of the grid described in Table 2 for Taurus. In each panel we show all the models corresponding to a single value of the temperature: T = 10 K (left), 25 K (center), and 35 K (right).</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. A.2.** Same as Fig. A.1, but in each panel we show all the models sharing the same value of density: $n_H = 10^4$ cm$^{-3}$ (left), $10^5$ cm$^{-3}$ (center), and $10^6$ cm$^{-3}$ (right).</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Fig. A.3.** Same as Fig. A.1, but in each panel we show all the models sharing the same value of molecular hydrogen cosmic-ray ionization rate: $\zeta_{H_2} = 10^{-17} \text{ s}^{-1}$ (left), $5 \times 10^{-17} \text{ s}^{-1}$ (center), and $10^{-16} \text{ s}^{-1}$ (right).</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Fig. A.4.** Same as Fig. A.1, but in each panel we show all the models sharing the same value of chemical age: $t = 0.1$ Myr (left), 1 Myr (center), and 10 Myr (right).</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>**Fig. A.5.** Comparison between model predictions and observed abundances toward the position with the highest visual extinction in each of the observed cuts in Taurus.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Fig. A.6.** Same as Fig A.5, but toward the position with the highest visual extinction in each of the observed cuts in Perseus and Orion.</figcaption></figure>
## Appendix B: Computation of the grain charge distribution

For the computation of the grain charges we rely on the statistical equilibrium between positive and negative charging. A detailed explanation on how to compute the grain charge is out of the scope of this paper, but the recent work by Ibáñez-Mejía et al. (2019) can be consulted for more details; we have used this work as a baseline for our implementation, and in this appendix we present the main differences between our procedure and that presented in Ibáñez-Mejía et al. (2019).

Positive charging of dust grains is mainly governed by collisions with the ions in the plasma ($J_{ion}$), and by radiative charging arising from the interstellar radiation field ($J_{pe}$) or from the cosmic-ray induced H$_2$ fluorescence ($J_{pe,CR}$, Prasad & Tarafdar 1983). On the other hand, negative charging of dust grains mainly arise from collisions with Maxwellian electrons in the plasma ($J_e$), since the electron flux induced from cosmic rays is negligible at these column densities (Ivlev et al. 2015).

The main singularity of our computation of grain charges is that we have tuned the input quantities to be representative of the two main molecular complexes studied in this paper: Taurus and Orion. With that purpose, the particle number density of the charged species in gas phase has been taken from the models built with Nautilus. We have considered two densities for each cloud: $n_H = 3.16 \times 10^3$ cm$^{-3}$ and $n_H = 3.16 \times 10^4$ cm$^{-3}$ for Taurus, and $n_H = 3.16 \times 10^4$ cm$^{-3}$ and $n_H = 3.16 \times 10^5$ cm$^{-3}$ for Orion; at these densities the main ionic species are H$^+$, He$^+$, H$_2^+$, H$_3^+$, N$^+$, O$^+$, HCO$^+$, C$^+$, and S$^+$. In addition, there are two parameters related to the cosmic-ray induced ultraviolet flux that we have set based on the particular traits of Orion and Taurus: the cosmic-ray ionization rate ($\zeta_{H_2}$) and the slope of the extinction curve at visible wavelengths ($R_V$). These terms influence the final far-ultraviolet photon flux based on the following relationship derived by Cecchi-Pestellini & Aiello (1992):

$$F_{UV} \approx 960 \left(\frac{1}{1-\omega}\right) \left(\frac{\zeta}{10^{-17} \text{s}^{-1}}\right) \left(\frac{N_{H_2}/A_V}{10^{21} \text{cm}^{-2} \text{mag}^{-1}}\right) \left(\frac{R_V}{3.2}\right)^{1.5} \quad (B.1)$$

For the dust albedo $\omega$ and the gas-to-extinction ratio $N_{H_2}/A_V$ we have assumed the same values than Ibáñez-Mejía et al. (2019) (0.5 and $1.87 \times 10^{21}$ cm$^{-2}$ mag$^{-1}$ respectively). However, based on the best-fit of the chemical models we have assumed $\zeta_{H_2} = 10^{-16}$ s$^{-1}$ for Taurus, while for Orion we have explored two possible scenarios: the low-$\zeta$ case suggested by our models ($\zeta_{H_2} = 10^{-17}$ s$^{-1}$) and a high-$\zeta$ model ($\zeta_{H_2} = 10^{-16}$ s$^{-1}$) based on observational evidence of some massive star-forming regions (Aharonian et al. 2019; Padovani et al. 2019). Finally, we have relied on the values for $R_V$ computed by Fitzpatrick & Massa (2007) based on the full extinction curve of hot stars toward these complexes. For the observed cuts explored in this paper, we have found that in Taurus the nearest star (search radius of 15') with measurements by Fitzpatrick & Massa (2007) is HD 29647 with $R_V = 3.46$. For Orion, we find one star for each pointing in a search radius between 10' and 15': HD 294264 ($R_V = 5.48$), HD 36982 ($R_V = 5.60$) and HD 37061 ($R_V = 4.55$); as a compromise, we have taken a characteristic value of $R_V = 5.5$ for the whole region, and we have checked that there are not any appreciable differences between the chosen value and that reported by the extinction curves.

The dust charge distribution is then computed by solving the equilibrium relation:

$$f(Z)[J_{pe}(Z) + J_{pe,CR} + J_{ion}(Z)] = f(Z+1)J_e(Z+1) \quad (B.2)$$

The probability distribution is computed between two extreme values $Z_{min}$ and $Z_{max}$ set by the user, for integer values of Z, and the final curve is normalized such that:

$$\sum_{Z=Z_{min}}^{Z_{max}} f(Z) = 1 \quad (B.3)$$

The python source code is available for download in the public github repository https://github.com/lbeitia/dust_distribution.