Open Access
Issue
A&A
Volume 670, February 2023
Article Number A114
Number of page(s) 23
Section Interstellar and circumstellar matter
DOI https://doi.org/10.1051/0004-6361/202244843
Published online 15 February 2023

© The Authors 2023

Licence Creative CommonsOpen Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.

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 ζH2$\sqrt {{\zeta _{{{\rm{H}}_2}}}} $, where ζH2${\zeta _{{{\rm{H}}_2}}}$ is the cosmic-ray ionization rate for H2 molecules, which becomes a key parameter in the molecular cloud evolution (McKee 1989; Caselli et al. 2002). In addition to ζH2${\zeta _{{{\rm{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 (AV < 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 H2S 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 (Eu = 570 K, ncrit = 1.5 × 106 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 & Walmsley 1978; Goicoechea & Cuadrado 2021). However, these lines are very weak, thus hindering their use as routine S+ tracers in molecular clouds. The detection of H2S 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 CH3OH and H2CS 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 ζ 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 × 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 photondominated 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 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 webpage1, or directly from the GEMS webpage2.

Table 1

Cores included in the GEMS sample and observation cuts associated with them.

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 (AV) and dust temperature (Td) maps are available (Malinen et al. 2012; Palmeirim et al. 2013; Lombardi et al. 2014; Zari et al. 2016). The angular resolution of the AV and Td 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 AV ~3 mag up to AV ~ 200 mag. The separation between one position and another in a given cut was selected to sample regular bins of AV. 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).

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, 12C/13C = 60, 32S/34S = 22.5, and 16O/18O = 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 34S at early time because of the 34S+ + CS reaction but not at the characteristic chemical ages of dense clouds (Loison et al. 2019). More controversial could be the 12CS/13CS 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/H13CN 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, O2, and SO2 (Loison et al. 2019). However, it is not expected to be relevant for C18O and HC18O+, 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, C34S, and 13CS line intensities were fitted by varying only two parameters, n(H2) and N(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: 13CO, C18O, HCO+, H13CO+, HC18O+, H13CN, HNC, HCS+, SO, 34SO, H2S, 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, 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(CO) = X(C18O) × 550, X(HCO+) = X(H13CO+) × 60, and X(HCN) = X(H13CN) × 60. In the case of HNC, we did not observe HN13C in GEMS, and the estimated abundance might be underestimated in positions with AV > 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 + CH → CS + 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 mantle 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 (DS) 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, DS = 1 ([S/H] = 1.5 × 10−5), moderate depletion, Ds = 10 ([S/H] = 1.5 × 10−6), and high depletion, Ds ~ 187 ([S/H] = 8.0 × 10−8). These values correspond to different estimates of DS in star-forming regions (Vastel et al. 2018; Navarro-Almaida et al. 2020, 2021; Hily-Blant et al. 2022). The value of ζH2${\zeta _{{{\rm{H}}_2}}}$ varied from ζH2=1017s1${\zeta _{{{\rm{H}}_2}}} = {10^{ - 17}}{{\rm{s}}^{ - 1}}$ which is typically found in dense and evolved cores (Caselli et al. 2002) to ζH2=5×1016s1${\zeta _{{{\rm{H}}_2}}} = 5 \times {10^{ - 16}}{{\rm{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 (1Myr), and steady-state chemistry (10 Myr). For this exploratory work, we only considered models with AV = 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, H2S, and OCS, for different values of DS (see Figs. A.1A.4). In each diagram, one of the input parameters is kept fixed while the others are allowed to vary in the whole range. The fixed parameter is the one indicated in Figs. A.1A.4. Models with different sulfur depletion are differentiated using a color code. When the model clouds corresponding to different values of DS appear grouped in bands with little overlap among them, we can conclude that the considered molecules are robust tracers of DS – If model clouds largely overlap then the solution is degenerate and we need to fix additional parameters to determine DS.

In order to base our discussion on quantitative grounds, we define the parameters R and R30%, which measure the overlap between model clouds corresponding to different values of DS. Fig. 1 illustrates how these parameters are calculated. In the upper panels of Fig. 1, the extreme points of models with the same DS are connected to construct polygons. The R parameter is calculated as, R = (Inter section)/ (Total), where Intersection refers to the intersection between the areas of the polygons corresponding to two values of DS (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 Areal + 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 DS. 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, R30%, that uses density curves to define the area corresponding to a given value of DS. In particular, Area130% and Area230% are defined as the areas enclosed by the 30% point density peak contour, and R30% = (Intersection30%)/(Total30%) (see bottom panels of Fig. 1). As expected, we find R30% < R in the case of a peaky distribution of points.

Table 3 shows the minimum, maximum, and mean values of R and R30% for the panels forming the corner diagrams shown in Figs. A.1A.4. We remind that we have only used the pairs of molecules, OCS-H2S, OCS-SO, OCS-CS, H2S-CS, H2S-SO, and CS-SO, in our analysis. Thus, the minimum, maximum, and mean R and R30% values are estimated taking into account only these six combination of sulfur species. Although the exact value of R and R30% 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, Rmean30%<0.3$R_{{\rm{mean}}}^{30\% } < 0.3$. Moreover the lowest values of Rmean30%$R_{{\rm{mean}}}^{30\% }$ 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 H2S are the best proxies of DS. In the panels including these species, the points corresponding to different values of DS appear grouped in clearly distinguishable bands. Moreover, the bands are almost orthogonal to the X(H2S) axis suggesting that the abundance of H2S is an excellent tracer of DS with little dependence on the other parameters. This can be easily interpreted in terms of sulfur chemistry. At such a low temperature, solid H2S is the main sulfur reservoir, and its abundance in gasphase 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 DS. Values of X(CS) ≥ 10−6, can only be achieved with low sulfur depletion (see Fig. A.1). However, the situation is less clear for X(CS) < 10−8 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 < 105 cm−3 with Rmean30%<0.4$R_{{\rm{mean}}}^{30\% } < 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, ζH2${\zeta _{{{\rm{H}}_2}}}$, 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 ζH2${\zeta _{{{\rm{H}}_2}}}$. 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 R30% shown in Table 3. Therefore, fixing ζH2${\zeta _{{{\rm{H}}_2}}}$ 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 ζH2${\zeta _{{{\rm{H}}_2}}}$ 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 ζH2${\zeta _{{{\rm{H}}_2}}}$ (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 > 104 cm−3) 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.

thumbnail Fig. 1

Model clouds formed by selecting all the models with T = 10 Κ and AV = 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 R30% (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 R30% is defined as the ratio of the intersection area (area enclosed within the orange contour) over the total area.

Table 2

Chemical model parameters.

Table 3

R and R30% factors.

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+, H2S, 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, ζ(H2), 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 AV, n(H2) and T. The visual extinction is assumed to be in the range AV × 0.5 and AV × 0.5 + 2, where AV 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(H2) to 5 × n(H2), being the values of n(H2), 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 χ2 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, Ddiff=1/nobs×Σi[ log10(Xmodi)log10(Xobsi) ]2$ {D_{{\rm{diff}}}} = 1/{n_{{\rm{obs}}}} \times {\Sigma _i}{\left[ {{{\log }_{10}}\left( {X_{\bmod }^i} \right) - {{\log }_{10}}\left( {X_{{\rm{obs}}}^i} \right)} \right]^2} $(1)

where nobs is the number of species detected at each position, Xmodi$X_{\bmod }^i$ is the model predicted abundance for the species i, and Xobsi$X_{{\rm{obs}}}^i$ 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 ζH2=(0.51.8)×1016s1${\zeta _{{{\rm{H}}_2}}} = \left( {0.5-1.8} \right) \times {10^{ - 16}}{{\rm{s}}^{ - 1}}$ and [S/H] = (0.4–2.2) × 10−6.

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 initially formed by large-scale compression of HI gas, and is now growing in mass due to the gravitational accretion of ambient 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, Β 213-N) where the density of protostars is higher. Cuts #10, #12, #16, and #17 are located in the quiescent part (hereafter, Β 213-S). Previous works based on CH3OH data suggest the existence of chemical differentiation between Β 213-N and Β 213-S (Spezzano et al. 2022; Punanova et al. 2022). In the cores located in Β 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 proposed that this behavior could be due to the irradiation on the cores due to nearby protostars which accelerate energetic particles along their outflows. A recent study of the H2CS deuterated compounds (Esplugues et al. 2022), showed that the chemical age of the cores in Β 213-N is higher than that of the cores in Β 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+, H2S, 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 Β 213 are best fitted assuming early time chemistry (t = 0.1 Myr) which is consistent with the idea that this filament is still accreting material from the surroundings, 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 AV < 8 mag being better fitted with t = 1 Myr while those at higher visual extinction being better fitted with t = 0.1 Myr. We also find a large dispersion in the values of ζH2${\zeta _{{{\rm{H}}_2}}}$. In spite of this, we can find some trends. Most positions in TMC 1 and Β 213-N are fitted with ζH2>1016s1${\zeta _{{{\rm{H}}_2}}} > {10^{ - 16}}{{\rm{s}}^{ - 1}}$, while Β 213-S is better described with ζH2<1016s1${\zeta _{{{\rm{H}}_2}}} < {10^{ - 16}}{{\rm{s}}^{ - 1}}$, consistent with the idea of a harsh environment in Β 213-N due to the presence of young (proto-)stars. A high value of ζH2${\zeta _{{{\rm{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 DS = 10 in TMC 1 and and Β 213-N, while DS = 187 in Β 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 LI544, compared with others at an earlier evolutionary 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 DS 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 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 H2CS.

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 DS decreases in dense regions. Therefore, we need to be cautious with this result.

thumbnail 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), ζH2${\zeta _{{{\rm{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 × 10−5.

thumbnail 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.

thumbnail 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.

6.2 Perseus

Located at a distance of 310 pc (Ortiz-Leon et al. 2018), the Perseus molecular cloud is considered the archetype 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 Β5 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 Bl-a and Β1-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 Β1-a or Β1-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, Β1b-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 ζH2${\zeta _{{{\rm{H}}_2}}}$ without any clear trend with visual extinction and density (see Fig. 4). Regarding sulfur depletion, most of the positions are fitted with DS = 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.

thumbnail 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.

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 13CO (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, ζH2<1016s1${\zeta _{{{\rm{H}}_2}}} < {10^{ - 16}}{{\rm{s}}^{ - 1}}$, 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 ζH2${\zeta _{{{\rm{H}}_2}}}$. Finally, a significant fraction of the observations are best fitted assuming DS = 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 DS 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 DS

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 DS ~ 10. Basically, all the other positions need a higher depletion, DS ~ 187. A different behavior is observed in Orion, in which a significant fraction of the positions (~40%) can be fitted with DS = 1. This suggests that sulfur depletion is dependent on the star formation activity in the environment. It is also interesting to explore possible trends of DS 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 nH ≤ 104 cm−3 that present lower values of DS. 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 DS 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 DS 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 DS with an accuracy better than a factor of ~10. In a second step, we fixed the values of t and ζH2${\zeta _{{{\rm{H}}_2}}}$ according to the results obtained in our first grid, and then ran a second one with finer steps of DS. In particular, we used ζH2=1016s1${\zeta _{{{\rm{H}}_2}}} = {10^{ - 16}}{{\rm{s}}^{ - 1}}$ for Taurus, ζH2=5×1017s1${\zeta _{{{\rm{H}}_2}}} = 5 \times {10^{ - 17}}{{\rm{s}}^{ - 1}}$ in Perseus, and ζH2=1017s1${\zeta _{{{\rm{H}}_2}}} = {10^{ - 17}}{{\rm{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 Ddiff 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 DS for each cut was varied from 1 to 187 in multiplicative steps of a factor of ~2 (DS = 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 Ddiff 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 C18O and H13CN to estimate their abundances. As commented above, isotopic chemical fractionation might be important for HCN (Loison et al. 2020), but the HCN/H13CN 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 systematics in the errors, with the SO abundance being usually over-predicted by the chemical model, while the HCS+ 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 H2S 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(CH3OH)/N(C18O) 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 gasphase 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 nH > 105 cm−3, 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 O2 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-H2S, 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-S8 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 S8 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 DS ~ 10 is already observed.

thumbnail 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.

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

In these dense and cold regions, the sum of the observed gasphase abundances of S-species (the most abundant are SO, SO2, H2S, CS, HCS+, H2CS, C2S and C3S) 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-SO2 (Boogert et al. 1997) are the only sulfur-bearing species detected thus far, and only upper limits to the s-H2S 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 S8, 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 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 S8 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+ ions frozen onto grain surfaces could produce large allotropes in a few times 104 but, in order to promote the formation of allotropes instead of sulfur hydrides, one needs to assume that the sticking coefficient of S+ 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+, S+, Na+,…) 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+, sulfur would remain ionized in the translucent medium until Av ~ 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 S8 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+ 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 some simple calculations to test whether there is some correlation 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 accretion 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 envelope. 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 ζH2=1016s1${\zeta _{{{\rm{H}}_2}}} = {10^{ - 16}}{{\rm{s}}^{ - 1}}$, G0 = 9 in unit of Habing field (χ = 5 in unit of Draine field), Τ = 15 Κ, and two values of density, nH = 3.16 × 103 cm−3 and nH = 3.16 × 104 cm−3. 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, AV < 1.5 mag (see Fig. 7). At higher visual extinctions, the S+ abundance decreases but remains >10−6 until AV ~ 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 ζH2=1017s1${\zeta _{{{\rm{H}}_2}}} = {10^{ - 17}}{{\rm{s}}^{ - 1}}$, more than 40% of the positions in OMC 3 are better fitted with ζH2=5×1017s1${\zeta _{{{\rm{H}}_2}}} = 5 \times {10^{ - 17}}{{\rm{s}}^{ - 1}}$ (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 density is high in all the cuts of this massive star-forming region (Rodríguez-Baras et al. 2021). Therefore, we fixed the values of G0 = 90 (χ = 50) and Τ = 25 Κ, and considered four sets of physical conditions in our calculations: 1) ζH2=1017s1,nH=3.16×104cm3;2)ζH2=1017s1,nH=3.16×105cm3;3)ζH2=1016s1,nH=3.16×104cm3${\zeta _{{{\rm{H}}_2}}} = {10^{ - 17}}{{\rm{s}}^{ - 1}},{n_{\rm{H}}} = 3.16 \times {10^4}{\rm{c}}{{\rm{m}}^{ - 3}};2){\zeta _{{{\rm{H}}_2}}} = {10^{ - 17}}{{\rm{s}}^{ - 1}},{n_{\rm{H}}} = 3.16 \times {10^5}{\rm{c}}{{\rm{m}}^{ - 3}};3){\zeta _{{{\rm{H}}_2}}} = {10^{ - 16}}{{\rm{s}}^{ - 1}},{n_{\rm{H}}} = 3.16 \times {10^4}{\rm{c}}{{\rm{m}}^{ - 3}};$; and 4) ζH2=1016s1,nH=3.16×105cm3${\zeta _{{{\rm{H}}_2}}} = {10^{ - 16}}{{\rm{s}}^{ - 1}},{n_{\rm{H}}} = 3.16 \times {10^5}{\rm{c}}{{\rm{m}}^{ - 3}}$. Moreover, we repeated the calculations for t = 0.1 Myr and t = 1 Myr. Assuming high density, nH = 3.16 × 105 cm−3, the S+ abundance is low, X(S+)<10−7, even at a low extinction of AV ~ 3.5 mag, regardless of the other parameters (see Figs. 8 and 9). In the case of moderate density, nH = 3.16 × 104 cm−3, we obtain high S+ abundance only in the case of ζH2=1016s1${\zeta _{{{\rm{H}}_2}}} = {10^{ - 16}}{{\rm{s}}^{ - 1}}$ and early-time chemistry.

The other essential ingredient in our interpretation is the grain charge. The charging of dust grains in the interstellar medium is governed by the equilibrium among different processes. On the one hand, a dust grain may acquire a positive charge due to photoelectric emission of electrons induced either by UV radiation (Draine 1978) or by fluorescence of H2 (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-Mejí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. 79. The grain charge distributions have been calculated for AV = 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 (AV ~ 1.5 mag), and changes to negative at higher extinction. In the case of Taurus, the net charge of grains is negative for AV ~ 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 AV = 5.5 mag where the abundance of S+ is already very low, especially for n(H2) ~ 105 cm−3. 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.

Table 4

Most abundant volatile sulfur-bearing species

thumbnail 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 (χUV = 5, ζH2=1016s1${\zeta _{{{\rm{H}}_2}}} = {10^{ - 16}}{{\rm{s}}^{ - 1}}$, Τ = 15 Κ, t = 0.1 Myr). Vertical lines indicate AV= 1.5, 3.5, and 5.5 mag. The horizontal blue line shows X(S+) = 10−7. Two values of density are considered, nH = 3.16 × 103 cm−3 (black) and nH = 3.16 × 104 cm−3 (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 Av = 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.

thumbnail 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 (χUV = 50, T = 25 Κ) and t = 0.1 Myr. We consider two values of density, nH = 3.16×104 cm−3 (black) and nH = 3.16×105 cm−3 (red). Continuous and dashedlines show the results for ζH2=1017s1${\zeta _{{{\rm{H}}_2}}} = {10^{ - 17}}{{\rm{s}}^{ - 1}}$ and ζH2=1016s1${\zeta _{{{\rm{H}}_2}}} = {10^{ - 16}}{{\rm{s}}^{ - 1}}$, respectively. Vertical lines indicate AV = 1.5, 3.5, and 5.5 mag. The horizontal blue line shows X(S+)=10−7. Right: Grain charge distribution for silicates with a size of 0.1 μm at AV = 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 ζH2=1017s1${\zeta _{{{\rm{H}}_2}}} = {10^{ - 17}}{{\rm{s}}^{ - 1}}$ and ζH2=1016s1${\zeta _{{{\rm{H}}_2}}} = {10^{ - 16}}{{\rm{s}}^{ - 1}}$, respectively.

thumbnail Fig. 9

Same as Fig. 8, but for t = 1 Myr.

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 × 10−5 (DS = 1), 1.5 × 10−6 (DS = 10), and 8 × 10−8 (DS = 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 ζH2${\zeta _{{{\rm{H}}_2}}}$ with a distribution peaking at ζH2=1×1016s1${\zeta _{{{\rm{H}}_2}}} = 1 \times {10^{ - 16}}{{\rm{s}}^{ - 1}}$, 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 Β 213-N are fitted with ζH2>1016s1${\zeta _{{{\rm{H}}_2}}} > {10^{ - 16}}{{\rm{s}}^{ - 1}}$, while the positions in Β 213-S are fitted with ζH2<1016s1${\zeta _{{{\rm{H}}_2}}} < {10^{ - 16}}{{\rm{s}}^{ - 1}}$, consistent with the idea of a harsh environment in Β 213-N due to the presence of young (proto-)stars. Regarding sulfur elemental abundance, we find DS = 10 in TMC 1 and and Β 213-N, and DS = 187 in Β 213-S. In addition to environment, we find some correlation of DS 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 ζH2${\zeta _{{{\rm{H}}_2}}}$ with a distribution peaking at ζH2~(510)×1017s1${\zeta _{{{\rm{H}}_2}}}\~\left( {5-10} \right) \times {10^{ - 17}}{{\rm{s}}^{ - 1}}$, without any clear trend with visual extinction and density. Regarding the sulfur elemental abundance, most of the positions are fitted with DS = 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 ζH2${\zeta _{{{\rm{H}}_2}}}$. While the majority of the positions in OMC 4 are fitted with ζH2=1017s1${\zeta _{{{\rm{H}}_2}}} = {10^{ - 17}}{{\rm{s}}^{ - 1}}$, more than 40% of the positions in OMC 3 are better fitted with ζH2=5×1017s1${\zeta _{{{\rm{H}}_2}}} = 5 \times {10^{ - 17}}{{\rm{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 ζH2${\zeta _{{{\rm{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 AV ~ 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'Études 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 Fundacion Ramon Areces for funding support through their international postdoc grant program. The authors also thank the anonymous referee for their interesting suggestions.

Appendix A Additional figures and tables

Table A.1

References for the collisional rate coefficients used for the volume density estimates.

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 Ddiff is the uncertainty as defined in Sect. 5.

thumbnail Fig. A.1

Corner diagrams performed by selecting the models with Av = 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 Κ (left), 25 Κ (center), and 35 Κ (right).

thumbnail Fig. A.2

Same as Fig. A.1, but in each panel we show all the models sharing the same value of density: nH = 104 cm−3 (left), 105 cm−3 (center), and 106 cm−3 (right).

thumbnail 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: ζH2=1017s1${\zeta _{{{\rm{H}}_2}}} = {10^{ - 17}}{{\rm{s}}^{ - 1}}$ (left), 5 × 10−17 s−1 (center), and 10−16 s−1 (right).

thumbnail 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).

thumbnail 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.

thumbnail 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.

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 (Jion), and by radiative charging arising from the interstellar radiation field (Jpe) or from the cosmic-ray induced H2 fluorescence (Jpe,CR, Prasad & Tarafdar 1983). On the other hand, negative charging of dust grains mainly arise from collisions with Maxwellian electrons in the plasma (Je), 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: nH = 3.16 × 103 cm−3 and nH = 3.16 × 104 cm−3 for Taurus, and nH = 3.16 × 104 cm−3 and nH = 3.16 × 105 cm−3 for Orion; at these densities the main ionic species are H+, He+, H+, H+, 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 (ζH2)$\left( {{\zeta _{{{\rm{H}}_2}}}} \right)$ and the slope of the extinction curve at visible wavelengths (RV). These terms influence the final far-ultraviolet photon flux based on the following relationship derived by Cecchi-Pestellini & Aiello (1992): FUV960(11ω)(ζ1017s1)(NH2/AV1021cm2mag1)(RV3.2)1.5$ {F_{{\rm{UV}}}} \simeq 960\left( {{1 \over {1 - \omega }}} \right)\left( {{\zeta \over {{{10}^{ - 17}}{{\rm{s}}^{ - 1}}}}} \right)\left( {{{{N_{{{\rm{H}}_2}/{A_V}}}} \over {{{10}^{21}}{\rm{c}}{{\rm{m}}^{ - 2}}{\rm{ma}}{{\rm{g}}^{ - 1}}}}} \right){\left( {{{{R_V}} \over {3.2}}} \right)^{1.5}} $(B.1)

For the dust albedo ω and the gas-to-extinction ratio NH2/AV${N_{{{\rm{H}}_2}}}/{A_V}$ we have assumed the same values than Ibáñez-Mejía et al. (2019) (0.5 and 1.87 × 1021 cm−2 mag−1 respectively). However, based on the best-fit of the chemical models we have assumed ζH2=1016s1${\zeta _{{{\rm{H}}_2}}} = {10^{ - 16}}{{\rm{s}}^{ - 1}}$ for Taurus, while for Orion we have explored two possible scenarios: the low-ζ case suggested by our models (ζH2=1017s1)$\left( {{\zeta _{{{\rm{H}}_2}}} = {{10}^{ - 17}}{{\rm{s}}^{ - 1}}} \right)$ and a high-ζ model (ζH2=1016s1)$\left( {{\zeta _{{{\rm{H}}_2}}} = {{10}^{ - 16}}{{\rm{s}}^{ - 1}}} \right)$ 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 RV 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 RV = 3.46. For Orion, we find one star for each pointing in a search radius between 10’ and 15’: HD 294264 (RV = 5.48), HD 36982 (RV = 5.60) and HD 37061 (RV = 4.55); as a compromise, we have taken a characteristic value of RV = 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)[ Jpe(Z)+Jpe,CR+Jion(Z) ]=f(Z+1)Je(Z+1)$ f\left( Z \right)\left[ {{J_{{\rm{pe}}}}\left( Z \right) + {J_{{\rm{pe,CR}}}} + {J_{{\rm{ion}}}}\left( Z \right)} \right] = f\left( {Z + 1} \right){J_{\rm{e}}}\left( {Z + 1} \right) $(B.2)

The probability distribution is computed between two extreme values Zmin and Zmax set by the user, for integer values of Z, and the final curve is normalized such that: Z=ZminZmaxf(Z)=1$ \sum\limits_{Z = {Z_{\min }}}^{{Z_{\max }}} {f\left( Z \right) = 1} $(B.3)

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

References

  1. Agúndez, M., & Wakelam, V. 2013, Chem. Rev., 113, 8710 [Google Scholar]
  2. Agúndez, M., Marcelino, N., Cernicharo, J., Roueff, E., & Tafalla, M. 2019, A&A, 625, A147 [Google Scholar]
  3. Aharonian, F., Yang, R., & de Oña Wilhelmi, E. 2019, Nat. Astron., 3, 561 [Google Scholar]
  4. Anderson, D. E., Bergin, E. A., Maret, S., & Wakelam, V. 2013, ApJ, 779, 141 [Google Scholar]
  5. André, P., Men’shchikov, A., Bontemps, S., et al. 2010, A&A, 518, L102 [Google Scholar]
  6. Anglada, G., Estalella, R., Pastor, J., Rodriguez, L. F., & Haschick, A. D. 1996, ApJ, 463, 205 [Google Scholar]
  7. Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481 [NASA ADS] [CrossRef] [Google Scholar]
  8. Bachiller, R., & Cernicharo, J. 1984, A&A, 140, 414 [NASA ADS] [Google Scholar]
  9. Bally, J., Langer, W. D., Stark, A. A., & Wilson, R. W. 1987, ApJ, 312, L45 [Google Scholar]
  10. Bayet, E., Aladro, R., Martín, S., Viti, S., & Martín-Pintado, J. 2009, ApJ, 707, 126 [Google Scholar]
  11. Berné, O., Marcelino, N., & Cernicharo, J. 2014, ApJ, 795, 13 [CrossRef] [Google Scholar]
  12. Blake, G. A., Mundy, L. G., Carlstrom, J. E., et al. 1996, ApJ, 472, L49 [NASA ADS] [CrossRef] [Google Scholar]
  13. Boogert, A. C. A., Schutte, W. A., Helmich, F. P., Tielens, A. G. G. M., & Wooden, D. H. 1997, A&A, 317, 929 [NASA ADS] [Google Scholar]
  14. Bouscasse, L., Csengeri, T., Belloche, A., et al. 2022, A&A, 662, A32 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  15. Bracco, A., Palmeirim, P., André, P., et al. 2017, A&A, 604, A52 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  16. Bronfman, L., Nyman, L. A., & May, J. 1996, A&AS, 115, 81 [Google Scholar]
  17. Bulut, N., Roncero, O., Aguado, A., et al. 2021, A&A, 646, A5 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  18. Cambrésy, L. 1999, A&A, 345, 965 [NASA ADS] [Google Scholar]
  19. Caselli, P., Walmsley, C. M., Terzieva, R., & Herbst, E. 1998, ApJ, 499, 234 [Google Scholar]
  20. Caselli, P., Walmsley, C. M., Zucconi, A., et al. 2002, ApJ, 565, 344 [Google Scholar]
  21. Cazaux, S., Carrascosa, H., Muñoz Caro, G. M., et al. 2022, A&A, 657, A100 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  22. Cecchi-Pestellini, C., & Aiello, S. 1992, MNRAS, 258, 125 [NASA ADS] [CrossRef] [Google Scholar]
  23. Cernicharo, J., & Guelin, M. 1987, A&A, 176, 299 [NASA ADS] [Google Scholar]
  24. Colgan, S. W. J., Schultz, A. S. B., Kaufman, M. J., Erickson, E. F., & Hollenbach, D. J. 2007, ApJ, 671, 536 [NASA ADS] [CrossRef] [Google Scholar]
  25. Colzi, L., Sipilä, O., Roueff, E., Caselli, P., & Fontani, F. 2020, A&A, 640, A51 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  26. Curtis, E. I., & Richer, J. S. 2011, MNRAS, 410, 75 [Google Scholar]
  27. Da Rio, N., Robberto, M., Hillenbrand, L. A., Henning, T., & Stassun, K. G. 2012, ApJ, 748, 14 [Google Scholar]
  28. Dagdigian, P. J. 2020, MNRAS, 494, 5239 [NASA ADS] [CrossRef] [Google Scholar]
  29. Davis, C. J., Chrysostomou, A., Hatchell, J., et al. 2010, MNRAS, 405, 759 [NASA ADS] [Google Scholar]
  30. Denis-Alpizar, O., Stoecklin, T., Guilloteau, S., & Dutrey, A. 2018, MNRAS, 478, 1811 [Google Scholar]
  31. Draine, B. T. 1978, ApJS, 36, 595 [Google Scholar]
  32. Draine, B. T., & Salpeter, E. E. 1979, ApJ, 231, 438 [Google Scholar]
  33. Draine, B. T., & Sutin, B. 1987, ApJ, 320, 803 [NASA ADS] [CrossRef] [Google Scholar]
  34. Dumouchel, F., Kłos, J., & Lique, F. 2011, Phys. Chem. Chem. Phys. (Incorp. Faraday Trans.), 13, 8204 [NASA ADS] [CrossRef] [Google Scholar]
  35. Dutrey, A., Wakelam, V., Boehler, Y., et al. 2011, A&A, 535, A104 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  36. Enoch, M. L., Young, K. E., Glenn, J., et al. 2006, ApJ, 638, 293 [Google Scholar]
  37. Esplugues, G., Fuente, A., Navarro-Almaida, D., et al. 2022, A&A, 662, A52 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  38. Feddersen, J. R., Arce, H. G., Kong, S., et al. 2020, ApJ, 896, 11 [CrossRef] [Google Scholar]
  39. Fehér, O., Tóth, L. V., Ward-Thompson, D., et al. 2016, A&A, 590, A75 [Google Scholar]
  40. Feng, S., Codella, C., Ceccarelli, C., et al. 2020, ApJ, 896, 37 [NASA ADS] [CrossRef] [Google Scholar]
  41. Fitzpatrick, E. L., & Massa, D. 2007, ApJ, 663, 320 [Google Scholar]
  42. Flower, D. R. 1999, MNRAS, 305, 651 [Google Scholar]
  43. Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306 [Google Scholar]
  44. Friesen, R. K., Pineda, J. E., co-PIs, et al. 2017, ApJ, 843, 63 [NASA ADS] [CrossRef] [Google Scholar]
  45. Fuente, A., Cernicharo, J., Roueff, E., et al. 2016, A&A, 593, A94 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  46. Fuente, A., Goicoechea, J. R., Pety, J., et al. 2017, ApJ, 851, L49 [NASA ADS] [CrossRef] [Google Scholar]
  47. Fuente, A., Navarro, D. G., Caselli, P., et al. 2019, A&A, 624, A105 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  48. Fuente, A., Treviño-Morales, S. P., Alonso-Albi, T., et al. 2021, MNRAS, 507, 1886 [NASA ADS] [CrossRef] [Google Scholar]
  49. Garrod, R. T., Wakelam, V., & Herbst, E. 2007, A&A, 467, 1103 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  50. Gatley, I., Becklin, E. E., Mattews, K., et al. 1974, ApJ, 191, L121 [Google Scholar]
  51. Geballe, T. R., Baas, F., Greenberg, J. M., & Schutte, W. 1985, A&A, 146, L6 [NASA ADS] [Google Scholar]
  52. Gerin, M., Pety, J., & Goicoechea, J. R. 2009, ASP Conf. Ser., 417, 165 [NASA ADS] [Google Scholar]
  53. Gerin, M., Ruaud, M., Goicoechea, J. R., et al. 2015, A&A, 573, A30 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  54. Gerin, M., Pety, J., Commerçon, B., et al. 2017, A&A, 606, A35 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  55. Goicoechea, J. R., & Cuadrado, S. 2021, A&A, 647, L7 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  56. Goicoechea, J. R., Pety, J., Gerin, M., et al. 2006, A&A, 456, 565 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  57. Goicoechea, J. R., Santa-Maria, M. G., Bron, E., et al. 2019, A&A, 622, A91 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  58. Goicoechea, J. R., Aguado, A., Cuadrado, S., et al. 2021, A&A, 647, A10 [EDP Sciences] [Google Scholar]
  59. Goldsmith, P. F., Heyer, M., Narayanan, G., et al. 2008, ApJ, 680, 428 [Google Scholar]
  60. Goodman, J., & Weare, J. 2010, Commun. Appl. Math. Comput. Sci., 5, 65 [Google Scholar]
  61. Gratier, P., Majumdar, L., Ohishi, M., et al. 2016, ApJS, 225, 25 [Google Scholar]
  62. Green, S., & Chapman, S. 1978, ApJS, 37, 169 [Google Scholar]
  63. Gustafsson, M., Kristensen, L. E., Clénet, Y., et al. 2003, A&A, 411, 437 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  64. Hacar, A., Tafalla, M., Kauffmann, J., & Kovács, A. 2013, A&A, 554, A55 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  65. Hacar, A., Alves, J., Tafalla, M., & Goicoechea, J. R. 2017a, A&A, 602, L2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  66. Hacar, A., Tafalla, M., & Alves, J. 2017b, A&A, 606, A123 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  67. Hacar, A., Tafalla, M., Forbrich, J., et al. 2018, A&A, 610, A77 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  68. Hacar, A., Bosman, A. D., & van Dishoeck, E. F. 2020, A&A, 635, A4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  69. Hasegawa, T. I., & Herbst, E. 1993, MNRAS, 261, 83 [NASA ADS] [CrossRef] [Google Scholar]
  70. Hatchell, J., & Dunham, M. M. 2009, A&A, 502, 139 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  71. Hatchell, J., Richer, J. S., Fuller, G. A., et al. 2005, A&A, 440, 151 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  72. Hatchell, J., Fuller, G. A., & Richer, J. S. 2007a, A&A, 472, 187 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  73. Hatchell, J., Fuller, G. A., Richer, J. S., Harries, T. J., & Ladd, E. F. 2007b, A&A, 468, 1009 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  74. Hernández Vera, M., Lique, F., Dumouchel, F., Hily-Blant, P., & Faure, A. 2017, MNRAS, 468, 1084 [Google Scholar]
  75. Hillenbrand, L. A. 1997, AJ, 113, 1733 [Google Scholar]
  76. Hily-Blant, P., Pineau des Forêts, G., Faure, A., & Lique, F. 2022, A&A, 658, A168 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  77. Hocuk, S., Szűcs, L., Caselli, P., et al. 2017, A&A, 604, A58 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  78. Holdship, J., Jimenez-Serra, I., Viti, S., et al. 2019, ApJ, 878, 64 [NASA ADS] [CrossRef] [Google Scholar]
  79. Ibáñez-Mejía, J. C., Walch, S., Ivlev, A. V., et al. 2019, MNRAS, 485, 1220 [CrossRef] [Google Scholar]
  80. Ishii, S., Nakamura, F., Shimajiri, Y., et al. 2019, PASJ, 71, S9 [Google Scholar]
  81. Ivlev, A. V., Padovani, M., Galli, D., & Caselli, P. 2015, ApJ, 812, 135 [Google Scholar]
  82. Jenkins, E. B. 2009, ApJ, 700, 1299 [Google Scholar]
  83. Jiménez-Escobar, A., & Muñoz Caro, G. M. 2011, A&A, 536, A91 [Google Scholar]
  84. Jiménez-Escobar, A., Muñoz Caro, G. M., Ciaravella, A., et al. 2012, ApJ, 751, L40 [Google Scholar]
  85. Johnstone, D., & Bally, J. 1999, ApJ, 510, L49 [Google Scholar]
  86. Johnstone, D., & Bally, J. 2006, ApJ, 653, 383 [NASA ADS] [CrossRef] [Google Scholar]
  87. Kirk, H., Johnstone, D., & Di Francesco, J. 2006, ApJ, 646, 1009 [Google Scholar]
  88. Kirk, H., Friesen, R. K., Pineda, J. E., et al. 2017, ApJ, 846, 144 [Google Scholar]
  89. Kong, S., Arce, H. G., Sargent, A. I., et al. 2019, ApJ, 882, 45 [NASA ADS] [CrossRef] [Google Scholar]
  90. Kutner, M. L., Evans, N. J. I., & Tucker, K. D. 1976, ApJ, 209, 452 [NASA ADS] [CrossRef] [Google Scholar]
  91. Laas, J. C., & Caselli, P. 2019, A&A, 624, A108 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  92. Lada, C. J., Alves, J., & Lada, E. A. 1996, AJ, 111, 1964 [Google Scholar]
  93. Lada, C. J., Muench, A. A., Haisch, Karl E. J., et al. 2000, AJ, 120, 3162 [NASA ADS] [CrossRef] [Google Scholar]
  94. Launhardt, R., Evans, I., Neal, J., Wang, Y., et al. 1998, ApJS, 119, 59 [NASA ADS] [CrossRef] [Google Scholar]
  95. Le Gal, R., Öberg, K. I., Loomis, R. A., Pegues, J., & Bergner, J. B. 2019, ApJ, 876, 72 [Google Scholar]
  96. Le Gal, R., Öberg, K. I., Teague, R., et al. 2021, ApJS, 257, 12 [NASA ADS] [CrossRef] [Google Scholar]
  97. Linke, R. A., & Goldsmith, P. F. 1980, ApJ, 235, 437 [Google Scholar]
  98. Lique, F., & Spielfiedel, A. 2007, A&A, 462, 1179 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  99. Liszt, H. S., Wilson, R. W., Penzias, A. A., et al. 1974, ApJ, 190, 557 [NASA ADS] [CrossRef] [Google Scholar]
  100. Loison, J.-C., Wakelam, V., Gratier, P., et al. 2019, MNRAS, 485, 5777 [Google Scholar]
  101. Loison, J.-C., Wakelam, V., Gratier, P., & Hickson, K. M. 2020, MNRAS, 498, 4663 [Google Scholar]
  102. Lombardi, M., Bouy, H., Alves, J., & Lada, C. J. 2014, A&A, 566, A45 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  103. Luhman, K. L., Stauffer, J. R., Muench, A. A., et al. 2003, ApJ, 593, 1093 [Google Scholar]
  104. Luhman, K. L., Mamajek, E. E., Allen, P. R., & Cruz, K. L. 2009, ApJ, 703, 399 [NASA ADS] [CrossRef] [Google Scholar]
  105. Malinen, J., Juvela, M., Rawlings, M. G., et al. 2012, A&A, 544, A50 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  106. Marcelino, N., Gerin, M., Cernicharo, J., et al. 2018, A&A, 620, A80 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  107. Marsh, K. A., Griffin, M. J., Palmeirim, P., et al. 2014, MNRAS, 439, 3683 [NASA ADS] [CrossRef] [Google Scholar]
  108. McKee, C. F. 1989, ApJ, 345, 782 [NASA ADS] [CrossRef] [Google Scholar]
  109. Minissale, M., Dulieu, F., Cazaux, S., & Hocuk, S. 2016, A&A, 585, A24 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  110. Monsch, K., Pineda, J. E., Liu, H. B., et al. 2018, ApJ, 861, 77 [NASA ADS] [CrossRef] [Google Scholar]
  111. Muench, A. A., Lada, E. A., Lada, C. J., & Alves, J. 2002, ApJ, 573, 366 [NASA ADS] [CrossRef] [Google Scholar]
  112. Nakamura, F., Ishii, S., Dobashi, K., et al. 2019, PASJ, 71, S3 [NASA ADS] [CrossRef] [Google Scholar]
  113. Narayanan, G., Heyer, M. H., Brunt, C., et al. 2008, ApJS, 177, 341 [NASA ADS] [CrossRef] [Google Scholar]
  114. Navarro-Almaida, D., Le Gal, R., Fuente, A., et al. 2020, A&A, 637, A39 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  115. Navarro-Almaida, D., Fuente, A., Majumdar, L., et al. 2021, A&A, 653, A15 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  116. Neufeld, D. A., & Wolfire, M. G. 2017, ApJ, 845, 163 [Google Scholar]
  117. Neufeld, D. A., Godard, B., Gerin, M., et al. 2015, A&A, 577, A49 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  118. Onishi, T., Mizuno, A., Kawamura, A., Ogawa, H., & Fukui, Y. 1996, ApJ, 465, 815 [NASA ADS] [CrossRef] [Google Scholar]
  119. Onishi, T., Mizuno, A., Kawamura, A., Tachihara, K., & Fukui, Y. 2002, ApJ, 575, 950 [NASA ADS] [CrossRef] [Google Scholar]
  120. Ortiz-León, G. N., Loinard, L., Dzib, S. A., et al. 2018, ApJ, 865, 73 [Google Scholar]
  121. Pabst, C. H. M., Goicoechea, J. R., Teyssier, D., et al. 2020, A&A, 639, A2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  122. Padoan, P., Cambrésy, L., & Langer, W. 2002, ApJ, 580, L57 [NASA ADS] [CrossRef] [Google Scholar]
  123. Padovani, M., Hennebelle, P., & Galli, D. 2013, A&A, 560, A114 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  124. Padovani, M., Marcowith, A., Sánchez-Monge, Á., Meng, F., & Schilke, P. 2019, A&A, 630, A72 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  125. Palmeirim, P., André, P., Kirk, J., et al. 2013, A&A, 550, A38 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  126. Palumbo, M. E., Tielens, A. G. G. M., & Tokunaga, A. T. 1995, ApJ, 449, 674 [NASA ADS] [CrossRef] [Google Scholar]
  127. Pankonin, V., & Walmsley, C. M. 1978, A&A, 64, 333 [NASA ADS] [Google Scholar]
  128. Pety, J., Teyssier, D., Fossé, D., et al. 2005, A&A, 435, 885 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  129. Pety, J., Guzmán, V. V., Orkisz, J. H., et al. 2017, A&A, 599, A98 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  130. Pineda, J. E., Caselli, P., & Goodman, A. A. 2008, ApJ, 679, 481 [Google Scholar]
  131. Pineda, J. E., Goodman, A. A., Arce, H. G., et al. 2010, ApJ, 712, L116 [Google Scholar]
  132. Pineda, J. E., Offner, S. S. R., Parker, R. J., et al. 2015, Nature, 518, 213 [NASA ADS] [CrossRef] [Google Scholar]
  133. Prasad, S. S., & Tarafdar, S. P. 1983, ApJ, 267, 603 [Google Scholar]
  134. Punanova, A., Vasyunin, A., Caselli, P., et al. 2022, ApJ, 927, 213 [CrossRef] [Google Scholar]
  135. Rebull, L. M., Padgett, D. L., McCabe, C. E., et al. 2010, ApJS, 186, 259 [NASA ADS] [CrossRef] [Google Scholar]
  136. Ridge, N. A., Di Francesco, J., Kirk, H., et al. 2006, AJ, 131, 2921 [Google Scholar]
  137. Rivière-Marichalar, P., Fuente, A., Goicoechea, J. R., et al. 2019, A&A, 628, A16 [Google Scholar]
  138. Rivière-Marichalar, P., Fuente, A., Le Gal, R., et al. 2020, A&A, 642, A32 [EDP Sciences] [Google Scholar]
  139. Robberto, M., Soderblom, D. R., Bergeron, E., et al. 2013, ApJS, 207, 10 [NASA ADS] [CrossRef] [Google Scholar]
  140. Rodríguez-Baras, M., Fuente, A., Riviére-Marichalar, P., et al. 2021, A&A, 648, A120 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  141. Ruaud, M., Wakelam, V., & Hersant, F. 2016, MNRAS, 459, 3756 [Google Scholar]
  142. Ruffle, D. P., Hartquist, T.W., Caselli, P., & Williams, D. A. 1999, MNRAS, 306, 691 [NASA ADS] [CrossRef] [Google Scholar]
  143. Sandell, G., & Knee, L. B. G. 2001, ApJ, 546, L49 [NASA ADS] [CrossRef] [Google Scholar]
  144. Savage, C., Apponi, A. J., Ziurys, L. M., & Wyckoff, S. 2002, ApJ, 578, 211 [NASA ADS] [CrossRef] [Google Scholar]
  145. Schmalzl, M., Kainulainen, J., Quanz, S. P., et al. 2010, ApJ, 725, 1327 [NASA ADS] [CrossRef] [Google Scholar]
  146. Schnee, S., Caselli, P., Goodman, A., et al. 2007, ApJ, 671, 1839 [NASA ADS] [CrossRef] [Google Scholar]
  147. Schnee, S., Enoch, M., Noriega-Crespo, A., et al. 2010, ApJ, 708, 127 [NASA ADS] [CrossRef] [Google Scholar]
  148. Scourfield, M., Viti, S., García-Burillo, S., et al. 2020, MNRAS, 496, 5308 [NASA ADS] [CrossRef] [Google Scholar]
  149. Shimajiri, Y., André, P., Palmeirim, P., et al. 2019, A&A, 623, A16 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  150. Shingledecker, C. N., Lamberts, T., Laas, J. C., et al. 2020, ApJ, 888, 52 [Google Scholar]
  151. Shirley, Y. L., Evans, Neal J. I., Young, K. E., Knez, C., & Jaffe, D. T. 2003, ApJS, 149, 375 [NASA ADS] [CrossRef] [Google Scholar]
  152. Smith, R. G. 1991, MNRAS, 249, 172 [NASA ADS] [CrossRef] [Google Scholar]
  153. Spezzano, S., Fuente, A., Caselli, P., et al. 2022, A&A, 657, A10 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  154. Spitzer, L. Jr 1941, ApJ, 93, 369 [NASA ADS] [CrossRef] [Google Scholar]
  155. Suri, S., Sánchez-Monge, Á., Schilke, P., et al. 2019, A&A, 623, A142 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  156. Tafalla, M., & Hacar, A. 2015, A&A, 574, A104 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  157. Tanabe, Y., Nakamura, F., Tsukagoshi, T., et al. 2019, PASJ, 71, S8 [NASA ADS] [Google Scholar]
  158. Taquet, V., Codella, C., De Simone, M., et al. 2020, A&A, 637, A63 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  159. Tatematsu, K., Umemoto, T., Kameya, O., et al. 1993, ApJ, 404, 643 [NASA ADS] [CrossRef] [Google Scholar]
  160. Tercero, F., López-Pérez, J. A., Gallego, J. D., et al. 2021, A&A, 645, A37 [EDP Sciences] [Google Scholar]
  161. Umebayashi, T., & Nakano, T. 1980, PASJ, 32, 405 [NASA ADS] [Google Scholar]
  162. 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 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  163. Vastel, C., Quénard, D., Le Gal, R., et al. 2018, MNRAS, 478, 5514 [Google Scholar]
  164. Vera, M. H., Kalugina, Y., Denis-Alpizar, O., Stoecklin, T., & Lique, F. 2014, J. Chem. Phys., 140, 224302 [CrossRef] [Google Scholar]
  165. Vidal, T. H. G., Loison, J.-C., Jaziri, A. Y., et al. 2017, MNRAS, 469, 435 [Google Scholar]
  166. Wakelam, V., Castets, A., Ceccarelli, C., et al. 2004, A&A, 413, 609 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  167. Wakelam, V., Herbst, E., Selsis, F., & Massacrier, G. 2006, A&A, 459, 813 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  168. Wakelam, V., Dartois, E., Chabot, M., et al. 2021, A&A, 652, A63 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  169. Wang, K. S., Bourke, T. L., Hogerheijde, M. R., et al. 2013, A&A, 558, A69 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  170. Warin, S., Castets, A., Langer, W. D., Wilson, R. W., & Pagani, L. 1996, A&A, 306, 935 [NASA ADS] [Google Scholar]
  171. Watson, W. D., & Salpeter, E. E. 1972, ApJ, 174, 321 [NASA ADS] [CrossRef] [Google Scholar]
  172. Wilking, B. A., Meyer, M. R., Greene, T. P., Mikhail, A., & Carlson, G. 2004, AJ, 127, 1131 [NASA ADS] [CrossRef] [Google Scholar]
  173. Wilson, T. L., & Rood, R. 1994, ARA&A, 32, 191 [Google Scholar]
  174. Wilson, R. W., Jefferts, K. B., & Penzias, A. A. 1970, ApJ, 161, L43 [NASA ADS] [CrossRef] [Google Scholar]
  175. Wu, J., Evans, Neal J.I., Shirley, Y. L., & Knez, C. 2010, ApJS, 188, 313 [NASA ADS] [CrossRef] [Google Scholar]
  176. Yang, B., Stancil, P. C., Balakrishnan, N., & Forrey, R. C. 2010, ApJ, 718, 1062 [Google Scholar]
  177. Yan, Q.-Z., Zhang, B., Xu, Y., et al. 2019, A&A, 624, A6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  178. Yazidi, O., Ben Abdallah, D., & Lique, F. 2014, MNRAS, 441, 664 [NASA ADS] [CrossRef] [Google Scholar]
  179. Zari, E., Lombardi, M., Alves, J., Lada, C. J., & Bouy, H. 2016, A&A, 587, A106 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  180. Zari, E., Brown, A. G. A., & de Zeeuw, P. T. 2019, A&A, 628, A123 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  181. Zhang, Z.-Y., Gao, Y., Henkel, C., et al. 2014, ApJ, 784, L31 [NASA ADS] [CrossRef] [Google Scholar]
  182. Zhao, B., Caselli, P., Li, Z.-Y., et al. 2016, MNRAS, 460, 2050 [Google Scholar]
  183. Zhou, S., Wu, Y., Evans, Neal J. I., Fuller, G. A., & Myers, P. C. 1989, ApJ, 346, 168 [NASA ADS] [CrossRef] [Google Scholar]
  184. Zinchenko, I., Mattila, K., & Toriseva, M. 1995, A&AS, 111, 95 [Google Scholar]
  185. Zucker, C., Speagle, J. S., Schlafly, E. F., et al. 2019, ApJ, 879, 125 [Google Scholar]
  186. Zuckerman, B. 1973, ApJ, 183, 863 [NASA ADS] [CrossRef] [Google Scholar]

All Tables

Table 1

Cores included in the GEMS sample and observation cuts associated with them.

Table 2

Chemical model parameters.

Table 3

R and R30% factors.

Table 4

Most abundant volatile sulfur-bearing species

Table A.1

References for the collisional rate coefficients used for the volume density estimates.

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 Ddiff is the uncertainty as defined in Sect. 5.

All Figures

thumbnail Fig. 1

Model clouds formed by selecting all the models with T = 10 Κ and AV = 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 R30% (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 R30% is defined as the ratio of the intersection area (area enclosed within the orange contour) over the total area.

In the text
thumbnail 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), ζH2${\zeta _{{{\rm{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 × 10−5.

In the text
thumbnail 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.

In the text
thumbnail 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.

In the text
thumbnail 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.

In the text
thumbnail 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.

In the text
thumbnail 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 (χUV = 5, ζH2=1016s1${\zeta _{{{\rm{H}}_2}}} = {10^{ - 16}}{{\rm{s}}^{ - 1}}$, Τ = 15 Κ, t = 0.1 Myr). Vertical lines indicate AV= 1.5, 3.5, and 5.5 mag. The horizontal blue line shows X(S+) = 10−7. Two values of density are considered, nH = 3.16 × 103 cm−3 (black) and nH = 3.16 × 104 cm−3 (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 Av = 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.

In the text
thumbnail 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 (χUV = 50, T = 25 Κ) and t = 0.1 Myr. We consider two values of density, nH = 3.16×104 cm−3 (black) and nH = 3.16×105 cm−3 (red). Continuous and dashedlines show the results for ζH2=1017s1${\zeta _{{{\rm{H}}_2}}} = {10^{ - 17}}{{\rm{s}}^{ - 1}}$ and ζH2=1016s1${\zeta _{{{\rm{H}}_2}}} = {10^{ - 16}}{{\rm{s}}^{ - 1}}$, respectively. Vertical lines indicate AV = 1.5, 3.5, and 5.5 mag. The horizontal blue line shows X(S+)=10−7. Right: Grain charge distribution for silicates with a size of 0.1 μm at AV = 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 ζH2=1017s1${\zeta _{{{\rm{H}}_2}}} = {10^{ - 17}}{{\rm{s}}^{ - 1}}$ and ζH2=1016s1${\zeta _{{{\rm{H}}_2}}} = {10^{ - 16}}{{\rm{s}}^{ - 1}}$, respectively.

In the text
thumbnail Fig. 9

Same as Fig. 8, but for t = 1 Myr.

In the text
thumbnail Fig. A.1

Corner diagrams performed by selecting the models with Av = 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 Κ (left), 25 Κ (center), and 35 Κ (right).

In the text
thumbnail Fig. A.2

Same as Fig. A.1, but in each panel we show all the models sharing the same value of density: nH = 104 cm−3 (left), 105 cm−3 (center), and 106 cm−3 (right).

In the text
thumbnail 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: ζH2=1017s1${\zeta _{{{\rm{H}}_2}}} = {10^{ - 17}}{{\rm{s}}^{ - 1}}$ (left), 5 × 10−17 s−1 (center), and 10−16 s−1 (right).

In the text
thumbnail 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).

In the text
thumbnail 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.

In the text
thumbnail 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.

In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.