# On the reservoir of sulphur in dark clouds: chemistry and elemental abundance reconciled

Thomas Vidal, Jean-Christophe Loison, Adam Yassin Jaziri, M. Ruaud, P.

Gratier, Valentine Wakelam

► To cite this version:

Thomas Vidal, Jean-Christophe Loison, Adam Yassin Jaziri, M. Ruaud, P. Gratier, et al. On the reservoir of sulphur in dark clouds: chemistry and elemental abundance reconciled. Monthly Notices of the Royal Astronomical Society: Letters, 2017, 469, pp.435-447. 10.1093/mnras/stx828. hal-01509533

HAL Id: hal-01509533

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

Submitted on 27 Jun 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# On the reservoir of sulphur in dark clouds: chemistry and elemental abundance reconciled

Thomas H. G. Vidal,<fnref n="1" />★ Jean-Christophe Loison,<fnref n="2" />,<fnref n="3" /> Adam Yassin Jaziri,<fnref n="1" /> Maxime Ruaud,<fnref n="1" /> Pierre Gratier<fnref n="1" /> and Valentine Wakelam<fnref n="1" />

Accepted 2017 March 31. Received 2017 March 31; in original form 2016 June 28

## ABSTRACT

Sulphur-bearing species are often used to probe the physical structure of star-forming regions of the interstellar medium, but the chemistry of sulphur in these regions is still poorly understood. In dark clouds, sulphur is supposed to be depleted under a form that is still unknown despite numerous observations and chemical modelling studies that have been performed. In order to improve the modelling of sulphur chemistry, we propose an enhancement of the sulphur chemical network using experimental and theoretical literature. We test the effect of the updated network on the outputs of a three-phases gas–grain chemical model for dark cloud conditions using different elemental sulphur abundances. More particularly, we focus our study on the main sulphur reservoirs as well as on the agreement between model predictions and the abundances observed in the dark cloud TMC-1 (CP). Our results show that depending on the age of the observed cloud, the reservoir of sulphur could either be atomic sulphur in the gas phase or HS and H₂S in icy grain bulks. We also report the first chemical model able to reproduce the abundances of observed S-bearing species in TMC-1 (CP) using as elemental abundance of sulphur its cosmic value.

**Key words:** astrochemistry – methods: numerical – stars: formation – ISM: abundances – ISM: clouds – ISM: molecules.

## 1 INTRODUCTION

Since the first detection of an S-bearing molecule in the interstellar medium (CS; Penzias et al. 1971), sulphur has become a subject of interest for many astrochemists. S-bearing molecules are indeed often used to probe the physical structure of star-forming regions (Lada, Bally & Stark 1991; Plume et al. 1997; Viti et al. 2001; Sakai et al. 2014; Podio et al. 2015) and have been proposed as chemical clocks for hot cores (Charnley 1997; Hatchell et al. 1998; Wakelam et al. 2004). But sulphur chemistry in the dense interstellar medium has also been at the centre of a puzzling issue for many years, namely the sulphur depletion problem (Ruffle et al. 1999). Unlike most of the other elements, in the diffuse medium, the gas-phase abundance of atomic sulphur is observed to be constant with cloud density, around its cosmic value of a few 10⁻⁵ (see for instance Jenkins 2009). However, in dark clouds the total abundance of detected S-bearing molecules only accounts for 0.1 per cent of the cosmic abundance of atomic sulphur (Tieftrunk et al. 1994; Charnley 1997). Therefore, the main reservoirs of sulphur in dark clouds are still un-

known. One could argue that sulphur is in its atomic form, that is, non-observable. However, assuming an elemental abundance of sulphur in chemical models as high as its cosmic value produces predicted abundances of observable S-bearing molecules much higher than the observed one. Consequently, modellers usually assume that the elements (sulphur but other heavy elements as well) are depleted compared to their cosmic reference values and they use the so-called depleted abundance of sulphur of a few 10⁻⁸ (Wakelam & Herbst 2008).

The main hypothesis to explain this missing sulphur is that it is depleted on to interstellar grains. In cometary ices, which are thought to present chemical similarities with the ices processed during hot core formation (Irvine et al. 2000), H₂S is the most abundant S-bearing molecule, at the level of 1.5 per cent compared to water (Bockelée-Morvan et al. 2000). More recently, Holdship et al. (2016) studied the properties of H₂S in the low-mass protostar L1157-B1 and found that a significant fraction of the sulphur is likely to be locked into the form of H₂S on the grains. Chemical models also predict that, in the dense interstellar medium, atomic sulphur would stick on grains and be mostly hydrogenated to form H₂S. To this day, however, only OCS (Palumbo, Geballe & Tielens 1997) and SO₂ (Boogert et al. 1997) have been likely identified

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

<footnote n="2">²Univ. Bordeaux, ISM, UMR 5255, F-33400 Talence, France</footnote>

<footnote n="3">³CNRS, ISM, UMR 5255, F-33400 Talence, France</footnote>

<footnote>* E-mail: thomas.vidal@u-bordeaux.fr</footnote>
in icy grain bulks in dense molecular clouds surrounding high-mass protostars and their estimated total abundance does not account for the missing sulphur. Upper limits for the column density of H₂S in icy grain bulks have been derived by Smith (1991), notably towards the line of sight of three late-type field stars lying behind the Taurus dark cloud, but these are also too low for H₂S to be the reservoir of sulphur in dark clouds.

Laboratory experiments coupled with chemical models have recently brought new insight into the problem by studying the irradiation of H₂S interstellar ice analogues by energetic protons and UV photons by Garozzo et al. (2010) and Jiménez-Escobar & Muñoz Caro (2011). Both studies found that solid H₂S was easily destroyed to form other species such as OCS, SO₂, CS₂ and H₂S₂. Subsequently, a corresponding network was proposed by Druard & Wakelam (2012) to reproduce these experiments but the authors showed that sulphur was converted into other forms, mostly H₂S₂ and H₂S₃, at higher temperatures (>20 K) than those found in dark clouds (around 10 K).

In this paper, we propose an enhancement of the sulphur chemical network for dark cloud simulations motivated by the recent observations of HNCS and HSCN in TMC-1 (CP) (Adande et al. 2010) and of CH₃SH in IRAS 16293–2422 (Majumdar et al. 2016). Our network also includes the network proposed by Druard & Wakelam (2012). We test the effect of the updated network on the outputs of a gas–grain chemical model for dark cloud conditions using different sulphur elemental abundances. More particularly, we focus our study on the main sulphur reservoirs as well as on the agreement between model predictions and the abundances observed in the dark cloud TMC-1 (CP). The gas–grain model we used as well as the modifications of the chemical network are presented in Sections 2 and 3. Sections 4 and 5 describe our new modelling results and include comparisons with the previous version of the network as well as comparisons with observations in the dark cloud TMC-1 (CP). We discuss and conclude on our results in the last section.

## 2 MODEL

The model we use to assess the impact of our new network on the predicted S-bearing species abundances of dark clouds, is the latest version of the NAUTILUS chemical model, described in Ruaud, Wakelam & Hersant (2016). This model allows us to compute the evolution of chemical abundances for a given set of physical and chemical parameters. Recently updated, it simulates a three-phases chemistry including gas-phase, grain surface and grain bulk chemistries, along with the possible exchanges between the phases. These exchanges are: the adsorption of gas-phase species on to grain surfaces, the thermal and non-thermal desorption of species from the grain surface into the gas phase and the surface-bulk and bulk-surface exchange of species. Our reference model for this paper uses the formulation of chemical desorption processes depicted in Garrod, Wakelam & Herbst (2007). They consider that for each surface reactions leading to a single product, a part of the energy released by the reaction can contribute to the desorption of the product in the gas phase using the Rice–Ramsperger–Kessel theory. The fraction of the product desorbed in the gas phase depends on the binding energy of the product to the surface, the enthalpy of the reaction and the fraction of the released energy that is lost to the surface. In our case, we use an *a* parameter of 0.001, which produces approximately a 1 per cent efficiency evaporation for all species. Moreover, the grain chemistry takes into account the standard direct photodissociation by UV photons along with the photodissociation induced by secondary UV photons introduced by Prasad & Tarafdar

<table><caption>**Table 1.** Initial abundances.</caption><thead><tr><th>Element</th><th>n<sub>i</sub>/n<sub>H</sub><sup>a</sup></th><th>References</th></tr></thead><tbody><tr><td>H<sub>2</sub></td><td>0.5</td><td></td></tr><tr><td>He</td><td>0.09</td><td>1</td></tr><tr><td>N</td><td>6.2(−5)</td><td>2</td></tr><tr><td>O</td><td>2.4(−4)</td><td>3</td></tr><tr><td>C<sup>+</sup></td><td>1.7(−4)</td><td>2</td></tr><tr><td>S<sup>+</sup></td><td>8.0(−8)</td><td>4</td></tr><tr><td>Si<sup>+</sup></td><td>8.0(−9)</td><td>4</td></tr><tr><td>Fe<sup>+</sup></td><td>3.0(−9)</td><td>4</td></tr><tr><td>Na<sup>+</sup></td><td>2.0(−9)</td><td>4</td></tr><tr><td>Mg<sup>+</sup></td><td>7.0(−9)</td><td>4</td></tr><tr><td>P<sup>+</sup></td><td>2.0(−10)</td><td>4</td></tr><tr><td>Cl<sup>+</sup></td><td>1.0(−9)</td><td>4</td></tr><tr><td>F</td><td>6.7(−9)</td><td>5</td></tr></tbody></table>

<sup>a</sup>*a*(*b*) stands for *a* × 10<sup>*b*</sup>.

References: (1) Wakelam & Herbst (2008), (2) Jenkins (2009), (3) Hincelin et al. (2011), (4) low metal abundances from Graedel, Langer & Frerking (1982) and (5) depleted value from Neufeld, Wolfire & Schilke (2005).

(1983). These processes are effective on the surface as well as in the bulk of the grains. The model also takes into account the newly implemented competition between reaction, diffusion and evaporation as suggested by Chang, Cuppen & Herbst (2007) and Garrod & Pauly (2011). The diffusion energies of each species are computed as a fraction of their binding energies. We take for the surface a value of this ratio of 0.4 as suggested by experiments and theoretical work made on H (see Ruaud et al. 2016, and reference therein), CO and CO₂ (see Karssemeijer & Cuppen 2014). This value is then extrapolated to every species on the surface. For the bulk, we take a value of 0.8 (see also Ruaud et al. 2016). The reference (gas phase and grains) network is *kida.uva.2014* (see Wakelam et al. 2015a), the same one as in Ruaud et al. (2016).

The set of physical parameters used throughout this paper is the commonly used dark cloud parameter configuration, namely a gas and dust temperature of 10 K, an atomic hydrogen total density of 2 × 10⁴ cm⁻³, a cosmic ionization rate of 1.3 × 10⁻¹⁷ s⁻¹ and a visual extinction of 15 mag. All abundances are expressed with respect to the total H density. Our set of initial abundances is summarized in Table 1. We start with all species in their atomic (or ionized) form, except for hydrogen which is assumed to be entirely in its molecular form.

## 3 MODIFICATION OF THE SULPHUR NETWORK

To update the sulphur chemistry network, we first examine the existing KIDA network (*kida.uva.2014*; Wakelam et al. 2015a) looking systematically at the possible reactions between S and S⁺ with the most abundant species in dense molecular clouds (CO, CH₄, C₂H₂, c-C₃H₈) as well as the potential reactions between sulphur compounds and the most abundant reactive species in molecular clouds (C, C⁺, H, N, O, OH, CN). We found that various neutral reactions were missing from KIDA (or other data bases) such as C + H₂S, C + H₂CS and S + I, c-C₃H₈. When previous experimental or theoretical studies exist, we use them to update the KIDA data base when necessary. However, there are relatively few kinetics data for the reactions of sulphur compounds. To estimate the unknown rate
constants to be used in the network, we use a methodology developed in appendix A (see supplementary material online) including new DFT and ab-initio calculations on the reactions H + CS, H + H₂CS, H + C₂S, H + HNCS, H + HSCN, NH + CS, NH₂ + CS, C₂H₃ + CS, O + C₃S, S + c-C₃H₃⁺, S + c-C₃H₂, N + HCS. All calculations have been performed for gas-phase reactions. The hydrogenation reaction barriers are used in the surface chemistry without changes. The rate constant of the NH + CS, NH₂ + CS, C₂H₃ + CS, O + C₃S reactions are presented in appendix B (see supplementary material online). It should be noted that there are often large uncertainties on rate constant values and branching ratios for the sulphur chemical network. This is particularly true for the branching ratios of dissociative recombination (DR). Most of them are deduced from general rules, developed in appendix A, by analogy with similar oxygenated compounds, but DR of HCS⁺ leads mainly to S + CH (Montaigne et al. 2005), while DR of HCO⁺ leads mainly to H + CO (Hamberg et al. 2014). This may be critical for H₂CS (and also HC₃S) production if DR of CH₃SH⁺ leads mainly to C-S bond breaking in contrast to DR of CH₂OH⁺ which preserves the C-O bond (Hamberg et al. 2007).

In this work, we introduce three isomers for HNCS (HNCS, HSCN, HCNS) and five for the protonated forms. We use the work of Gronowski & Kolos (2014) to describe the ionic chemistry but, unlike them, we do not consider the CSH⁺ isomer but only the HCS⁺ one. Indeed, the CSH⁺ → HCS⁺ isomerization barrier is only 5.4 kJ mol⁻¹ (Puzzarini 2005) and as the reaction producing CSH⁺ (CS + H₃⁺) is exothermic by 78 kJ mol⁻¹, most of the CSH⁺ should isomerize into HCS⁺. The comparison with HOC⁺ is not relevant, as the COH⁺/HCO⁺ isomerization barrier is equal to 150 kJ mol⁻¹ (Nobes & Radom 1981). We carefully searched for other reactions producing HNCS isomers, performing various theoretical calculations on the H-H-C-N-S system. For the H₂ + CS reaction, we found a small barrier (+3.0 kJ mol⁻¹) at the CCSD(T)/cc-pVQZ//MP2/cc-pVTZ level (this barrier being submerged by -2.4 kJ mol⁻¹ at the M06-2X/cc-pVTZ level but this method is thought to be less precise and usually slightly underestimates barrier values). We calculate the rate constant for this reaction using the CCSD(T) values leading to a negligible value at low temperature. We also include various reactions on the grain surface. Indeed, it appears that the s-N + s-HCS reaction is, alongside the following reaction mechanism:

$$ \mathrm{HCS}^{+} + \mathrm{NH}_{2} \rightarrow \mathrm{H} + \mathrm{H}_{2}\mathrm{NCS}^{+} \xrightarrow{e^{-}} \mathrm{HNCS}/\mathrm{HSCN} + \mathrm{H}, \quad (1) $$

the main source of HNCS isomers. The first step of this reaction is very likely N-C bond formation without a barrier through the pairing up of free electrons of the N atom and the HCS radical leading to HC(N)S formation. This species is indeed the transition state linking HNCS and HSCN. The branching ratio between HNCS and HSCN formation can be deduced from statistical theory. HNCS/HSCN production is estimated to be proportional to the density of vibrational states of each isomer near the effective barrier to isomerization in a similar manner to HCN/HNC formation in DR of HCNH⁺ (Herbst, Terzieva & Talbi 2000). Using the MESMER program (Glowacki et al. 2012) for the calculations of the ro-vibrational density of states near the isomerization barrier (see appendix A).

For the loss reactions, we consider no barrier for the reactions of HNCS isomers with carbon atoms taking into account the high reactivity of carbon atoms with unsaturated closed shell molecules (see table 1 in Loison et al. 2014). We neglect the reactions of HNCS isomers with H, N and O atoms in the gas phase either from M06-2X/cc-pVTZ calculations (see appendix A) or by com-

<table><caption>**Table 2.** Table of the 46 S-bearing species added to the network.</caption><thead><tr><th colspan="2">Neutrals</th></tr></thead><tbody><tr><td>CS<sub>2</sub></td><td>HNCS</td></tr><tr><td>H<sub>2</sub>S<sub>3</sub><sup>a</sup></td><td>HSCN</td></tr><tr><td>S<sub>3</sub><sup>a</sup></td><td>H<sub>2</sub>C<sub>3</sub>S</td></tr><tr><td>S<sub>4</sub><sup>a</sup></td><td>HNCHS</td></tr><tr><td>S<sub>5</sub><sup>a</sup></td><td>HSCHN</td></tr><tr><td>S<sub>6</sub><sup>a</sup></td><td>HNCSH</td></tr><tr><td>S<sub>7</sub><sup>a</sup></td><td>NH<sub>2</sub>CHS</td></tr><tr><td>S<sub>8</sub><sup>a</sup></td><td>H<sub>2</sub>C<sub>2</sub>S</td></tr><tr><td>HC<sub>3</sub>S</td><td>NH<sub>2</sub>CS</td></tr><tr><td>HC<sub>2</sub>S</td><td>NH<sub>2</sub>CH<sub>2</sub>S</td></tr><tr><td>CH<sub>3</sub>S</td><td>NH<sub>2</sub>CH<sub>2</sub>SH</td></tr><tr><td>CH<sub>3</sub>SH</td><td>HSO</td></tr><tr><td>CH<sub>2</sub>SH</td><td>HNCHSH</td></tr><tr><td>HCNS</td><td>NH<sub>2</sub>CHSH</td></tr></tbody><thead><tr><th colspan="2">Ions</th></tr></thead><tbody><tr><td>CS<sub>2</sub><sup>+</sup></td><td>HSCN<sup>+</sup></td></tr><tr><td>HCS<sub>2</sub><sup>+</sup></td><td>H<sub>2</sub>CNS<sup>+</sup></td></tr><tr><td>H<sub>2</sub>C<sub>2</sub>S<sup>+</sup></td><td>H<sub>2</sub>NCS<sup>+</sup></td></tr><tr><td>H<sub>2</sub>C<sub>3</sub>S<sup>+</sup></td><td>H<sub>2</sub>SCN<sup>+</sup></td></tr><tr><td>CH<sub>3</sub>S<sup>+</sup></td><td>HCNSH<sup>+</sup></td></tr><tr><td>CH<sub>3</sub>SH<sup>+</sup></td><td>HNCSH<sup>+</sup></td></tr><tr><td>CH<sub>3</sub>SH<sub>2</sub><sup>+</sup></td><td>CH<sub>3</sub>CS<sup>+</sup></td></tr><tr><td>HCNS<sup>+</sup></td><td>NH<sub>2</sub>CHSH<sup>+</sup></td></tr><tr><td>HNCS<sup>+</sup></td><td>NH<sub>3</sub>CH<sub>2</sub>SH<sup>+</sup></td></tr></tbody></table>

parison with similar reactions with HNCO isomers. For He⁺ and H⁺ reactions with HNCS isomers we consider, by comparison with similar reactions (Anicich 2003), a charge transfer mechanism as the main process leading to highly excited species followed by dissociation. For C⁺ reactions with HNCS isomers, the charge exchange is usually a minor channel (Anicich 2003). We consider then that these reactions lead mainly to the most exothermic exit channel involving the fewest steps. For protonation reactions, we use the work of Gronowski & Kolos (2014) leading to five H₂NCS⁺ isomers. The DR rate constant of the five H₂NCS⁺ isomers is assumed to be equal to the ion-electron collision rate estimated to be equal to $3 \times 10^{-7} \times (T/300)^{-0.5}$ considering the size of the cations (Florescu-Mitchell & Mitchell 2006; Fournier et al. 2013). For the branching ratio of the DR of H₂NCS⁺ and HNCSH⁺, a critical parameter for the HNCS/HSCN ratio, we follow the same procedure as used for the DR of HCNH⁺ (Herbst et al. 2000) using the MESMER program for the calculations of the ro-vibrational density of states near the isomerization barrier (see appendix A). Then, the DR reactions of H₂NCS⁺ and HNCSH⁺ leads to similar amounts of HNCS and HSCN.

Overall, we added (or reviewed in the case of reactions) 46 S-bearing species to the network along with 478 reactions in the gas phase, 305 reactions on the grain surface and 147 reactions in the grain bulk (see the table in appendix B). The newly introduced species are listed in Table 2. Among these species, H₂S₃ and Sₙ (n = 3 to 8) are only present on the grains, that is, they are formed on the surfaces and not allowed to evaporate (see also Druard & Wakelam 2012). At 10 K, these species are not efficiently formed so we do not want to add any gas-phase routes for these species as this chemistry is not known. It should be noted that our enhanced network includes the chemical schemes for carbon chains proposed

<footnote><sup>a</sup>These species are not allowed to evaporate in the gas phase (see text).</footnote>
<figure><img src="image_2.png" /><figcaption>Figure 1. Abundances of main sulphur-bearing species relative to H as a function of time for dark cloud physical conditions: (a) in the gas phase, (b) on the grain surface (solid line) and bulk (dashed line). The dotted line represents the elemental abundance of sulphur (here, $X_{\text{ini}} = 8 \times 10^{-8}$).</figcaption></figure>

in Wakelam et al. (2015b), Loison et al. (2016), Hickson, Wakelam & Loison (2016b) and Loison et al. (in preparation).

# 4 IMPACT OF THE NEW NETWORK ON THE CHEMICAL MODEL OF DARK CLOUDS

In this section, we highlight the impact of our enhanced sulphur network on the outputs of the NAUTILUS chemical code, configured for dark cloud physical conditions. We first present the abundance evolution of the most abundant sulphur-bearing species we obtain from the model in this configuration. Then we highlight our results on the newly implemented sulphur-bearing species. Finally, we compare the outputs obtained with our new network with those obtained with the nominal network for the same model configuration.

## 4.1 Abundances of the main sulphur-bearing species

In the following, we study the abundances of the main sulphur-bearing species obtained with our enhanced network, both in the gas phase and on the grains. We choose these species since they contain more than 5 per cent of the elemental sulphur at one point in our simulation.

### 4.1.1 Main gas-phase species: S, CS and SO

Fig. 1(a) shows the time evolution of the main sulphur-bearing species in the gas phase: atomic sulphur (S), carbon monosulphide (CS) and sulphur monoxide (SO). In our simulation, sulphur is initially in the form of $S^+$ and between $7.7 \times 10^3$ and $4.6 \times 10^5$ yr, S becomes the main sulphur reservoir (including species on grains). During this period, its abundance reaches a maximum at $10^5$ yr when it contains up to 73 per cent of the elemental sulphur. The growth of the atomic S abundance up to this maximum seems to be

mainly caused by the electronic recombination mechanisms of $S^+$, as well as of $CS^+$ and $HCS^+$. The last two are related to $S^+$ by the following reaction mechanism:

$$S^+ \xrightarrow{CH} CS^+ \xrightarrow{H_2} HCS^+ \quad (2)$$

This mechanism is therefore very efficient because CH is rapidly formed by electronic recombination of $CH_2^+$ which is itself formed by the ion-neutral reaction between elemental species $C^+$ and $H_2$. Later, S is consumed by reactions with $H_3^+$, $CH_3$, OH and $O_2$, which efficiently produce SO causing the peak in its abundance at $3.6 \times 10^5$ yr. The CS abundance follows closely the atomic sulphur one, first growing from numerous reactions, including notably the electronic recombination of $HCS^+$, $C_2S^+$ and $HC_3S^+$ and destruction of $C_2S$ and by atomic carbon. It reaches a maximum between $10^4$ and $10^5$ yr, becoming the second main reservoir of sulphur, containing up to 29 per cent of the initial sulphur. However, CS does not accumulate, in contrast to CO, as the hydrogenation of CS followed by DR of $HCS^+$ produce much more $S + CH$ than $H + CS$ (Montaigne et al. 2005).

It should be noted that we cannot reproduce the observed $CS/HCS^+$ ratio for TMC-1 of 26 (Gratier et al. 2016) at early time but the agreement become better, if not perfect, for longer age. Although $HCS^+$ and CS are linked through protonation reaction and DR, other reactions are involved in the chemistry of these species. Before $4.6 \times 10^3$ yr, CS is mainly formed by DR of $HCS^+$ and $HCS^+$ is formed by reaction 2. At longer time neutral reactions such as $S + CH \rightarrow CS + H$ plays also a role. Comparing our model results with the observed abundance of CS of $6.5 \times 10^{-9}$ (Gratier et al. 2016), the problem is very likely an overestimation of the CS abundance, which can come from wrong branching ratios of the $HCS^+$ DR (we currently use the measured values from Montaigne et al. 2005), wrong observational estimations of CS (for example, underestimation of CS due to saturation of CS lines), or an additional loss not producing $HCS^+$. A potentially new sink for CS may be the $O + CS$ reaction. We currently use the rate from Lilenfeld & Richardson (1977) determined in the 150–300 K range leading to a negligible rate at 10 K. However, the work of Gonzalez et al. (1996) may suggest a non-negligible rate at low temperature.

### 4.1.2 Main species on grains: HS and $H_2S$

On the grains, it appears that the most abundant species are the hydrosulphide radical (HS) and hydrogen sulphide ($H_2S$). The abundances of those species, both on the surface (solid line) and in the bulk (dashed line), are depicted in Fig. 1(b). At times prior to 100 years, both species are formed on the surface by the successive hydrogenation of atomic sulphur physisorbed on the grains. Once the species are formed, they begin to sink into the bulk. After a hundred years, both species continue to form on the surface and eventually enter into a loop created by the two following reactions:

$$s-HS + s-H \rightarrow s-H_{2S} \quad (3)$$

$$s-H_{2S} + s-H \rightarrow s-H_2 + s-HS. \quad (4)$$

Reaction (4) has a barrier and should not be efficient under dark cloud physical conditions. However, it becomes efficient in our model because of the reaction–diffusion competition whose main effect is to increase the rates of reactions with activation barrier through tunnelling (Ruaud et al. 2016). Therefore, this cycle enabled by the reaction–diffusion competition and supplied in HS by the hydrogenation of S on the grains, keeps increasing at equal rates
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Abundances of HNCS (black) and HSCN (grey) relative to H as a function of time for dark cloud physical conditions in the gas phase (solid line), on the grain surface (dotted line) and in the grain bulk (dashed line).</figcaption></figure>

the abundances of the two species during the rest of the simulation. For this model, we have assumed a barrier of $Ea = 860$ K but there remains uncertainties in this value. The gas-phase reaction $H + H_2S$ has been measured in a wide range of temperature from 190 up to 2237 K (Kurylo, Peterson & Braun 1971; Peng, Hu & Marshall 1999; Yoshimura et al. 1992). The rate constant shows strong non-Arhenius behaviour at low temperature likely due to the importance of tunnelling. Peng et al. (1999) performed theoretical calculations of the barrier of the reaction equal to 1330 at DFT level and 1930 K at QCISD(T) level. Considering these uncertainties, we have used the experimental value of Kurylo et al. (1971), which is the one measured at the lower temperature (in the 190–464 K range). We also tested the higher values by Peng et al. (1999) but it did not change significantly the results shown in this paper.

As HS and $H_2S$ keep on sinking into the bulk, their abundances in the latter eventually becomes higher than at the surface. At approximately $6 \times 10^5$ yr and until $10^7$ yr, HS and $H_2S$ (in the bulk) are the main reservoirs of sulphur, sharing more than 80 per cent of the initial sulphur nearly equally divided between the two species, with a slight excess for HS. It should be noted that during this period of time, the 20 per cent of the initial sulphur remaining is mostly divided between atomic S in the gas phase and in the bulk, as well as NS in the bulk.

## 4.2 Abundances of the newly implemented S-bearing species

Our enhanced sulphur network includes the addition of 45 new S-bearing species. Among those species, we now present the results of our calculations for HNCS, the species around which most of our new network has been built, and CH$_3$SH, which was recently detected in the envelope of the low-mass protostar IRAS 16293–2422 (Majumdar et al. 2016).

### 4.2.1 HNCS and HSCN

Isothiocyanic acid HNCS is the most stable among possible CHNS isomers. It has been very well studied since its first detection in Sgr B2 by Frerking, Linke & Thaddeus (1979), and has since been detected in the dark cloud TMC-1 (CP) (Adande et al. 2010). In contrast, thiocyanic acid HSCN is highly unstable and is about a factor of 3 less abundant in Sgr B2 than HNCS (Halfen et al. 2009), and has been also detected in TMC-1 (CP) with a similar abundance to HNCS (Adande et al. 2010). Fig. 2 shows the time evolution of their abundances in the three phases of the model (in the gas phase, on the grain surface and in the grain bulk) for dark cloud conditions. It should be noted that our network for HNCS and

HSCN is an improved version of the one proposed in Gronowski & Kolos (2014), which itself is a revised version of the network proposed in Adande et al. (2010). It appears from our simulations that both species are essentially produced at the surface of the grains where their main formation reaction is

$$s-N + s-HCS \rightarrow s-HNCS \text{ or } s-HSCN. \quad (5)$$

HNCS and HSCN abundances on the grain surface are also regulated by hydrogenation loops similar to the one described in Section 4.1.2: for HNCS with HNCSH and $NH_2CS$, and for HSCN with HNCSH and HSCHN. It should be noted that the hydrogenation of HSCN is calculated to be much less efficient than the hydrogenation of HNCS (see appendix A). At times $>2.8 \times 10^5$ yr, as less and less atomic N is available, abundances of both species decrease on the surface. For HNCS, it happens through the chemical desorption of the products of its hydrogenation, HNCSH and $NH_2CS$. For HSCN, it happens through its successive hydrogenations. Moreover, as the chemistry of both species is not efficient in the bulk, its abundance increase in this phase is primarily due to sinking from the surface.

In the gas phase, HNCS and HSCN abundances increase efficiently from chemical desorption from the grains following reaction (5), as well as from the following reaction mechanism:

$$HCS^{+} + NH_{2} \rightarrow H + H_{2}NCS^{+} \xrightarrow{e^{-}} HNCS/HSCN + H. \quad (6)$$

Consequently, after $2.8 \times 10^5$ yr, as reaction (5) is less efficient because there is less N available, reaction (6) then becomes the main formation pathway of HNCS and HSCN in the gas phase. They are also predicted to be efficiently destroyed by atomic C via:

$$C + HNCS \rightarrow HNC + CS \quad (7)$$

$$\rightarrow HCN + CS \quad (8)$$

$$C + HSCN \rightarrow HCN + CS. \quad (9)$$

The depletion of atomic C on grains after $10^5$ yr then explains the increase in the gas-phase abundance gradients for both species at this time.

### 4.2.2 CH$_3$SH

Numerous previous detections of methyl mercaptan CH$_3$SH in hot cores such as Sgr B2 (Linke, Frerking & Thaddeus 1979), G327.3–0.6 (Gibb et al. 2000) and Orion KL (Kolesniková et al. 2014) suggest that this species initially forms in the ices and then evaporates in hot cores. More recently, Majumdar et al. (2016) presented the detection of CH$_3$SH in IRAS 16293–2422 and proposed an associated chemical network to explain its observed abundance. This work is included in our enhanced sulphur network and has been completed, notably regarding grain chemistry. Fig. 3 shows the time evolution of the CH$_3$SH abundance. It appears from our results that CH$_3$SH is mainly formed on the grain surface and is released into the gas phase by chemical desorption in dark cloud conditions. The main formation processes on the grains are:

$$s-H + s-CH_3S \rightarrow s-CH_3SH \quad (10)$$

$$s-H + s-CH_2SH \rightarrow s-CH_3SH. \quad (11)$$

Moreover, similarly to HNCS, it seems that CH$_3$SH chemistry in the grain bulk is relatively inert, which throughout the simulation explains its accumulation in this phase from sinking from the grain surface.
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Abundance of CH₃SH relative to H as a function of time for dark cloud physical conditions in the gas phase (solid line), on the grain surface (dotted line) and in the grain bulk (dashed line).</figcaption></figure>

As said above, the gas-phase chemistry of CH₃SH is essentially ruled by its chemical desorption from the grain surface during the entire simulation. As for HNCS, it should be noted that it is consumed by atomic C at times <2.8 × 10⁵ yr following:

$$C + \mathrm{CH_3SH} \rightarrow \mathrm{CH_3} + \mathrm{HCS}. \qquad (12)$$

However, after this time, in contrast to HNCS, gas-phase CH₃SH is still produced from grain surface chemistry and chemical desorption, which can be explained by the fact that the CS abundance in the gas phase remains high, and gradually depletes to form CH₃SH through successive hydrogenations. Moreover, as the hydrogen atom abundance remains relatively constant on the grain surface during the simulation, due to cosmic ray interactions with H₂, CH₃SH is efficiently produced in contrast to HNCS and HSCN as nitrogen atomic strongly decreases.

## 4.3 Comparison with the previous version of the network

We now highlight the differences between the nominal model (hereafter Model A) and Model B, which takes into account our enhanced sulphur network. It should be noted that Model A already includes the reactions between HS and H₂S studied in Section 4.1.2, as well as the reaction-diffusion competition mechanism, two key parameters for sulphur chemistry. As our calculations show that numerous S-bearing species are impacted by the changes we made to the network, we present here only a selection of these species (C₃S, H₂CS, HCS, H₂S, HS and SO). We selected these molecules because they present a difference between Models A and B of more than one order of magnitude in the time period between 10⁵ and 10⁶ yr and because their abundances rise above 10⁻¹² during this period. Indeed, the estimated ages of well-studied dark clouds such as TMC-1 or L134N are thought to correspond to this time period (see discussion in Agúndez & Wakelam 2013). Hence, an abundance of 10⁻¹² is considered as the detection limit for species in the gas phase.

### 4.3.1 C₃S

Fig. 4 shows the comparison of the abundances of C₃S calculated for Models A and B, both in the gas phase and in the grain bulk. It can be seen that our new network induces differences for this species that can reach up to two orders of magnitude in the gas phase and four orders of magnitude in the grain bulk at times around 10⁵ yr. In both models, C₃S is mainly formed in the gas phase by the reaction:

$$\mathrm{HC_3S^+} + \mathrm{e^-} \rightarrow \mathrm{H} + \mathrm{C_3S}. \qquad (13)$$

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Comparison of the abundances of C₃S relative to H as a function of time for dark cloud physical conditions in the gas phase (solid line) and in the grain bulk (long dashed line).</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Comparison of the abundances of H₂CS relative to H as a function of time for dark cloud physical conditions in the gas phase (solid line) and on the grain surface (dotted line).</figcaption></figure>

HC₃S⁺ being efficiently produced through reactions of atomic sulphur with c,l-C₃H₂⁺ and c,l-C₃H₃⁺. We thus consider that the decrease of its abundance in the grain bulk in Model B (as compared to Model A) is a direct result of its decrease in the gas phase. This decrease is due to the following added reaction:

$$C + \mathrm{C_3S} \rightarrow \mathrm{C_3} + \mathrm{CS}. \qquad (14)$$

This reaction becomes efficient in Model B at times >10³ yr when the abundance of atomic C in the gas phase reaches 10⁻⁴. Finally, after 10⁴ yr, while in Model A the abundance of C₃S continues to increase because of reaction (13), in Model B the small peak at 2.8 × 10⁵ yr is due to the reaction:

$$\mathrm{O} + \mathrm{HC_3S} \rightarrow \mathrm{C_3S} + \mathrm{OH}. \qquad (15)$$

Reaction (15) becomes more efficient than reaction (14) at this time because depletion is more efficient for atomic carbon (~10⁵ yr) than for atomic oxygen (~8 × 10⁵ yr). It should be noted that branching ratios for the DR of HC₃S⁺ are unknown and our network may overestimate C₃S production leading to an overestimation of its abundance.

### 4.3.2 H₂CS

The abundance of H₂CS in the gas phase presents a similar change in its evolution between Models A and B (see Fig. 5). Its main reaction of formation in both models is the electronic recombination of H₃CS⁺:

$$\mathrm{H_3CS^+} + \mathrm{e^-} \rightarrow \mathrm{H} + \mathrm{H_2CS}. \qquad (16)$$
<figure><img src="image_7.png" /><figcaption>Figure 6. Same as Fig. 5 but for HCS.</figcaption></figure>

$H_3CS^+$ being efficiently produced through $S^+ + CH_4$ reaction. $H_2CS$ is less efficiently formed in Model B after $10^3$ yr because, like $C_3S$, it is consumed by atomic C via the following two reactions:

$$C + H_2CS \rightarrow H + HC_2S \quad (17)$$

$$C + H_2CS \rightarrow CH_2 + CS. \quad (18)$$

These reactions induce differences in the $H_2CS$ abundance of more than two orders of magnitude between $10^4$ and $10^5$ yr. As atomic carbon is depleted on to the grains, these reactions become less efficient after $10^5$ yr and the abundance of $H_2CS$ increases again through reaction (16).

In contrast to $C_3S$, $H_2CS$ is also efficiently formed on the grain surface in both models by the hydrogenation of HCS. However, in Model B, two hydrogenation processes of $H_2CS$ have been added, which cause its abundance at the grain surface to drop compared to Model A by more than two orders of magnitude. Those processes are:

$$s-H_2CS + s-H \rightarrow s-CH_2SH \quad (19)$$

$$\rightarrow s-CH_3S, \quad (20)$$

which also lead to the formation of $CH_3SH$. Moreover, as already noted, the branching ratios for the DR of $H_3CS^+$ are unknown and may lead mainly to breaking of the C-S bond which will strongly limit the $H_2CS$ abundance.

### 4.3.3 HCS

Fig. 6 shows the time evolution of the HCS abundance in the gas phase in both Models A and B. We can see in this figure that the overall effect of our enhanced network on this species is to increase its abundance in the gas phase above the detection level of $10^{-12}$. This potential observability of HCS in dark clouds will be discussed later. In Model A, HCS is mainly formed by the two following reactions:

$$S + CH_2 \rightarrow H + HCS \quad (21)$$

$$H_2CS^+ + e^- \rightarrow H + HCS. \quad (22)$$

However, in Model B two new reactions show a higher HCS production efficiency. Those reactions are:

$$H_3CS^+ + e^- \rightarrow H + H + HCS \quad (23)$$

$$C + H_2S \rightarrow H + HCS. \quad (24)$$

<figure><img src="image_8.png" /><figcaption>Figure 7. Comparison of the abundances of HS, $H_2S$ and SO, relative to H as a function of time for dark cloud physical conditions in the gas phase (solid line) and on the grain surface (dotted line).</figcaption></figure>

It should be noted that HCS is also formed at the grain surface by the hydrogenation of CS, and that its abundance in this phase is decreased in Model B because of numerous added reactions with atomic O and N, the latter forming the newly implemented species HNCS on the grains.

### 4.3.4 HS, $H_2S$ and SO

We choose to describe the changes brought by our new model to HS, $H_2S$ and SO abundances together, notably on the grain surface, because they are correlated.

On the one hand, in the gas phase, among those three species, the only one that presents a significant difference in Model B (compared to Model A) is $H_2S$ (see Fig. 7). This difference (that reaches more than one order of magnitude) appears between $5 \times 10^3$ and $3 \times 10^5$ yr and is mainly due to the fact that reaction (24) slows down the increase of $H_2S$ during this period of time. Besides this change, gas-phase $H_2S$ is efficiently formed in both models on the grain surface by the hydrogenation of HS (reaction 4) enabled by the reaction-diffusion competition, followed by desorption.

On the other hand, our modifications of the network bring about a new reaction mechanism on the grain surface that causes several
<table><caption>**Table 3.** Observed abundances of S-bearing species in TMC-1 (CP).</caption><thead><tr><th>Species</th><th>n<sub>i</sub>/n<sub>H</sub><sup>a</sup></th><th>References</th></tr></thead><tbody><tr><td>OCS</td><td>1.1(−9)</td><td>1</td></tr><tr><td>NS</td><td>4.0(−10)</td><td>2</td></tr><tr><td>CS</td><td>6.5(−9)</td><td>3</td></tr><tr><td>HCS<sup>+</sup></td><td>2.5(−10)</td><td>3</td></tr><tr><td>H<sub>2</sub>CS</td><td>2.6(−9)</td><td>3</td></tr><tr><td>C<sub>2</sub>S</td><td>5.0(−9)</td><td>3</td></tr><tr><td>C<sub>3</sub>S</td><td>7.0(−10)</td><td>3</td></tr><tr><td>SO</td><td>1.0(−8)</td><td>4</td></tr><tr><td>SO<sub>2</sub></td><td>1.5(−10)</td><td>5</td></tr><tr><td>HNCS</td><td>4.15(−12)</td><td>6</td></tr><tr><td>HSCN</td><td>6.3(−12)</td><td>6</td></tr></tbody></table>

<sup>a</sup>*a(b)* stands for *a* × 10<sup>*b*</sup>.

References (1) Matthews et al. (1987), (2) McGonagle, Irvine & Ohishi (1994), (3) Gratier et al. (2016), (4) Lique, Cernicharo & Cox (2006), (5) Cernicharo et al. (2011) and (6) Adande et al. (2010).

changes to the abundances of HS, H<sub>2</sub>S and SO:

$$s-S \xrightarrow{s-H} s-HS \qquad (25)$$

$$s-O \xrightarrow{s-SO} s-SO \xrightarrow{s-H} s-HSO \xrightarrow{s-H} s-HS. \qquad (26)$$

From Model A to Model B, the main effects of this reaction mechanism are:

(i) to globally increase the grain surface abundance of HS because of the added path (26),

(ii) to globally increase the grain surface abundance of H<sub>2</sub>S because of the loop described in Section 4.1.2, and hence its abundance in the gas phase (except during the time period described above),

(iii) to globally decrease the grain surface abundance of SO which is, in our enhanced network, successively hydrogenated to form HSO then HS.

## 5 COMPARISON WITH OBSERVATIONS

In this section, we compare our modelling results with observations of sulphur bearing species in dense molecular clouds. We use the method of the so-called distance of disagreement described in Wakelam, Herbst & Selsis (2006) and which is computed as follows:

$$D(t) = \frac{1}{N_{\text{obs}}} \sum_{i} | \log[n(X)_{\text{obs},i}] - \log[n(X)_i(t)] |, \qquad (27)$$

where *n(X)obs,i* and *n(X)i(t)* are the observed abundances (relative to H) of the *i*th species and the modelled ones at time *t*, respectively. *Nobs* is the total number of observed species considered. The minimum of function *D* (hereafter noted *Dmin*) is then the minimum average difference (in magnitude) between modelled and observed abundances. We call the time for which this minimum is obtained the ‘best fit’ time (hereafter noted *tBF*), which consequently is the estimated age of the observed object according to the model.

We detail the observed abundances we use for S-bearing species in Table 3, and for other species we use table 4 of Agúndez & Wakelam (2013) and reference therein.

<table><caption>**Table 4.** Results of the comparison of Models A and B with observations in the dark cloud TMC-1 (CP).</caption><thead><tr><th>Model</th><th>D<sub>min</sub></th><th>t<sub>BF</sub> (yr)</th><th>Fraction of reproduced molecules</th></tr></thead><tbody><tr><td>A</td><td>0.908</td><td>2.8 × 10<sup>5</sup></td><td>64 per cent</td></tr><tr><td>B</td><td>0.766</td><td>2.8 × 10<sup>5</sup></td><td>67 per cent</td></tr></tbody></table>

### 5.1 Comparison with Models A and B

In order to validate our enhanced sulphur network, we compare the results of both Models A and B with the observations of the well-studied TMC-1 (CP) dark cloud for which numerous observational constraints are available. Indeed, more than 60 gas-phase species have been detected in this source, and upper limits on the abundance of 7 more are available (see a review of these species in Adande et al. 2010; Agúndez & Wakelam 2013). For the calculation of the corresponding distance of disagreement for both models, we take into account 58 of the 64 detected species because six of them are not implemented in one or either of the two networks (nominal and enhanced). Using equation (27), we find a similar level of agreement between Models A and B, with however a slightly better agreement with observations for Model B with *Dmin* = 0.766 (against *Dmin* = 0.908 for Model A), giving a ‘best fit’ time of 2.8 × 10<sup>5</sup> yr, same as Model A. At this time, and assuming that the observed abundances were reproduced by the model when the difference between the two was smaller than one order of magnitude, the fractions of reproduced molecules are 64 and 67 per cent, respectively, for Model A and Model B. A summary of these results is shown in Table 4.

For both models, five of the nine S-bearing species taken into account in the comparison are reproduced: OCS, CS, C<sub>3</sub>S, SO and SO<sub>2</sub>, and Model B gives a slightly better agreement than Model A for these species. The similarity of the results comes from the fact that Model A includes the reaction-diffusion competition as well as the hydrogenation loop described in Section 4.1.2, both key parameters for sulphur chemistry. Moreover, considering only the 58 species taken into account in Model A, Model B reproduces only two more species than Model A which are carbon chains. However, the goal of this paper is to evaluate the efficiency of our enhanced network (Model B) to reproduce S-bearing species abundances, which should be done by considering all the species it takes into account. Therefore, we hereafter choose to exclude Model A from our study and focus only on Model B.

### 5.2 Variation of the elemental sulphur abundance

In order to assess the issue of sulphur depletion in chemical simulations of dark clouds with our enhanced network, we run the model for four different values of the elemental sulphur abundance and compare the results with observations to determine which one leads to better agreement with the observed species abundances. These values vary from the so-called depleted one (8 × 10<sup>-8</sup>), which is commonly used by chemical models to reproduce observations and which we used previously in this paper, to the cosmic one (1.5 × 10<sup>-5</sup>), which is considered to be the real elemental abundance of sulphur for dark cloud formation. The different values of the sulphur elemental abundances used and the labels of their associated models are summarized in Table 5. Apart from sulphur, the other species elemental abundances we use are those reported in Table 1.

In the following, we first compare for each model our calculated abundances of gas-phase species with those determined from the
<table><caption>**Table 5.** Set of elemental sulphur abundances.</caption><thead><tr><th>Model</th><th>n(S)<sub>ini</sub><sup>a</sup></th></tr></thead><tbody><tr><td>1</td><td>8(−8)</td></tr><tr><td>2</td><td>5(−7)</td></tr><tr><td>3</td><td>5(−6)</td></tr><tr><td>4</td><td>1.5(−5)</td></tr></tbody><tfoot><tr><td colspan="2"><sup>a</sup> a(b) stands for a × 10<sup>b</sup>.</td></tr></tfoot></table>

<figure><img src="image_9.png" /><figcaption>**Figure 8.** Comparison of the distances of disagreement for Models 1, 2, 3 and 4 considering 61 detected species in TMC-1 (CP).</figcaption></figure>

well-studied dark cloud TMC-1 (CP). The only sulphur-bearing species observed on ice being OCS and SO₂ in the molecular cloud surrounding the deeply embedded protostar W33A, we then compare the modelled abundances of these species on the grains with the observed ones.

### 5.2.1 Comparison with all observed gas-phase species in TMC-1 (CP)

To assess the overall agreement of our enhanced network with observations, we compute the distance of disagreement for Models 1, 2, 3 and 4 with the 61 observed species that are taken into account in this network. Fig. 8 displays these distances of disagreement between 10⁵ and 10⁷ yr, the time period when the minimum values are reached. On the one hand, it appears that Model 3 is the best model in that it presents the distance of disagreement with the lowest minimum of $D_3(t_{BF}) = 0.736$ with $t_{BF} = 10^6$ yr which is still a reasonable estimated age for TMC-1 (CP). Moreover at that time, 44 of the 61 species considered (72 per cent) are reproduced by the model within a factor of 10. This result means that the elemental sulphur abundance which our new model theoretically needs to best reproduce the observations is around $X(S)_{ini} = 5 \times 10^{-6}$. Model 4 also reproduce 44 of the 61 species at its 'best fit' time but its minima of distance of disagreement is slightly higher than Model 3 (0.772).

On the other hand, Fig. 8 shows that the maximum difference between the distances of disagreement of the four models is lower than 0.3. Hence, for a given time, the average of the differences of abundances between two models (for a same species) is lower than a factor of 2. The distance of disagreement is thus not very sensitive to the elemental sulphur abundance. This result can be explained by the fact that among the 61 observed species we consider for the calculations, only 11 of them are S-bearing species.

<figure><img src="image_10.png" /><figcaption>**Figure 9.** Comparison of the distances of disagreement for Models 1, 2, 3 and 4 considering only the 11 detected S-bearing species in TMC-1 (CP).</figcaption></figure>

### 5.2.2 Comparison with observed S-bearing gas-phase species in TMC-1 (CP)

Considering the low sensitivity of the distance of disagreement (61 observed species) to the elemental sulphur abundance and because our improvement of the network concerns mainly the sulphur chemistry, we also assess the efficiency of our enhanced network by checking its agreement with only the observations of the S-bearing species detected in TMC-1 (CP). Fig. 9 is the same as Fig. 8 but considering only the 11 detected S-bearing gas-phase species in TMC-1 (CP). Those species are OCS, NS, HCS⁺, CS, H₂CS, CCS, C₃S, SO, SO₂ (see Ágúndez & Wakelam 2013, for a review), HNCS and HSCN (Adande et al. 2010). As expected, the distances of disagreement are much more sensitive to the elemental sulphur abundance than in the previous case, with differences going up to an average order of magnitude of 1.6 at some time steps. Moreover, the figure shows that the distance of disagreement for Models 1 and 2 never goes down to 0.7, whereas the other two present their minima under 0.4. Models 1 and 2 are therefore the 'worst' models, which means the depleted elemental sulphur abundance of $8 \times 10^{-8}$ and $5 \times 10^{-7}$ do not allow the model to reproduce the observations of S-bearing species with our enhanced network.

However, at the times of their respective minimums, Models 3 and 4 reproduce all the S-bearing molecules considered. Model 4 is the one that presents the distance of disagreement with the lowest minimum of $D_4(t_{BF}) = 0.371$ with $t_{BF} = 1.3 \times 10^6$ yr. Furthermore, Model 3 allows us to obtain a minimum of the distance of disagreement that is nearly equal to the one found for Model 4, namely $D_3(t_{BF}) = 0.376$ with $t_{BF} = 10^6$ yr. In both cases, the time of 'best fit' is in the acceptable range for the age of TMC-1 (CP). These results suggest that our model is able to reproduce the observations of S-bearing species in TMC-1 (CP) using the cosmic elemental abundance of sulphur $X(S)_{ini} = 1.5 \times 10^{-5}$ or an abundance three times lower.

### 5.2.3 The case of HNCS and HSCN

As our enhanced network was mainly built around HNCS and its isomers, it is important to check that our results fit well with the observations of these species. Adande et al. (2010) identified the $7_{0,7} \rightarrow 6_{0,6}$ and $8_{0,8} \rightarrow 7_{0,7}$ transitions of HNCS and HSCN towards TMC-1 (CP) and were able to derive a corresponding ratio $[HNCS]/[HSCN] = 1.4 \pm 0.7$ and abundances for HNCS and HSCN of, respectively, $4.15 \times 10^{-12}$ and $3 \times 10^{-12}$. Fig. 10 shows the simulated abundances of HNCS and HSCN in the case of the cosmic elemental abundance for sulphur (Model 4) as well as their observed abundances in TMC-1 (CP) (horizontal dashed lines). At
<figure><img src="image_11.png" /><figcaption>**Figure 10.** Comparison of the simulated gas-phase abundances of HNCS and HSCN with the observed ones in TMC-1 (CP) (horizontal dashed lines), relative to H.</figcaption></figure>

$t_{\text{BF}} = 1.3 \times 10^6$ yr, we obtain a ratio of 1.68 that seems to fit well with the observed value. Moreover, the figure shows that at this time, the differences between the observed and calculated abundances of HNCS and HSCN are less than a factor 1.2 and 1.1, respectively, which confirms the efficiency of the network to model the chemistry of HNCS. It should be noted that Model 3 also allows us to reproduce well the observed abundance of HNCS and HSCN at $t_{\text{BF}} = 10^6$ yr.

### 5.2.4 Comparison with observed S-bearing species on grains towards W33A

Only two solid S-bearing species have been detected in icy grain bulks: OCS (Palumbo et al. 1997) and $SO_2$ (Boogert et al. 1997). Both species have been detected towards the deeply embedded protostar W33A and we can therefore use the results of our models in the dark cloud configuration to compare with these observations. In order to do so, we have to compare for a given species its observed abundance ($7 \times 10^{-8}$ and $6.2 \times 10^{-7}$ for OCS and $SO_2$, respectively) with the sum of its simulated abundances on the grain surface and in the grain bulk. For both species, Fig. 11 displays the comparison between their simulated total abundances for our four models in icy grain bulks and their observed ones between $10^5$ and $10^7$ yr. The age of the giant molecular cloud W33 is indeed supposed to lie in that time frame (Messineo et al. 2015).

Fig. 11(a) shows that both Models 3 and 4 reproduce the observed abundance of OCS in icy grain bulks, at $7.74 \times 10^5$ and $2.15 \times 10^6$ yr, respectively. Moreover, considering the fact that most of the young stellar objects detected in W33 have an estimated age of a few million years (see e.g. Messineo et al. 2015), the most favoured model is Model 4. Regarding $SO_2$, Fig. 11(b) shows that the model for which its modelled abundance in icy grain bulks comes the closest to its observed one is Model 4. Both comparisons suggest once again that our enhanced network needs an elemental sulphur abundance close to the cosmic one in order to best reproduce the observations. Note that we assume here that the detections of OCS and $SO_2$ ices are real although the two species are not in the list of ‘firmly detected’ species of Boogert, Gerakines & Whittet (2015). The JWST mission may give us more constraints on this matter.

## 5.3 Impacts of the use of the cosmic abundance of sulphur on the reservoirs

In Section 4, we studied the chemistry of the main S-bearing species using the elemental abundances listed in Table 1, more particularly

<figure><img src="image_12.png" /><figcaption>**Figure 11.** Comparison of the calculated abundances in icy grain bulks of OCS (a) and $SO_2$ (b) for Models 1, 2, 3 and 4 (solid lines) with the observed one towards W33A (dashed lines).</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 12.** Abundances of main sulphur-bearing species relative to H as a function of time for dark cloud physical conditions: (a) in the gas phase, (b) on the grain surface (solid line) and bulk (dashed line). The dotted line represents the elemental abundance of sulphur (here, $X_{\text{ini}} = 1.5 \times 10^{-5}$).</figcaption></figure>

the depleted abundance of sulphur of $X(S)_{\text{ini}} = 8 \times 10^{-8}$. Now that we determined that in order to reproduce S-bearing species observations in dark clouds, our model needs an elemental abundance of sulphur close to the cosmic one, we need to assess the impact of the use of such an abundance on the results.

Fig. 12 is the same as Fig. 1 but for an elemental abundance of sulphur of $X_{\text{ini}} = 1.5 \times 10^{-5}$. What appears is that the differences are mainly quantitative and that the chemistry stays the same as
described in Section 4, with the exception of the two following points.

(i) Adsorption of atomic S on the grain is much more efficient with $X_{ini} = 1.5 \times 10^{-5}$, which cause its abundance to drop much sooner than for $X(S)_{ini} = 8 \times 10^{-8}$, near $3 \times 10^3$ yr. However, atomic S still stays the main reservoir of sulphur between $2.8 \times 10^3$ and $4.6 \times 10^5$, containing at its maximum 88 per cent of the initial sulphur.

(ii) With $X_{ini} = 1.5 \times 10^{-5}$, we can no longer consider SO as a main sulphur-bearing species as defined in Section 4. Indeed, at $3.6 \times 10^5$ yr, time when its abundance peaks, it now only contains less than 2 per cent of the initial sulphur.

Hence, our results regarding the reservoirs of sulphur in dark clouds stay the same (with the exception of SO).

# 6 DISCUSSIONS AND CONCLUSIONS

## 6.1 About the elemental abundance of sulphur

The comparison of our enhanced model with the gas-phase observations in TMC-1 (CP) (see Section 5) favours the use of an elemental abundance of sulphur between $X(S)_{ini} = 5 \times 10^{-6}$ and $X(S)_{ini} = 1.5 \times 10^{-5}$. As most of the estimations of the cosmic abundance of sulphur lie between these two values (Federman et al. 1993; Sofia, Cardelli & Savage 1994; Tieftrunk et al. 1994; Ruffle et al. 1999; Shalabiea 2001), our model does not require additional depletion of sulphur. Moreover, our results on the comparison with the observations of OCS and SO₂ in icy grain bulks towards W33A confirm that a cosmic abundance of sulphur as the elemental sulphur abundance in our model is required to reproduce the observed abundances of these two species.

## 6.2 About the reservoirs of sulphur

Two main hypotheses exist to explain the observed depletion of sulphur in dark clouds: either the sulphur is in an as yet undiscovered form in the gas phase, which could be atomic sulphur, or it could be locked in icy grain bulks. The latter hypothesis is supported by the observation of H₂S as the most important S-bearing species in cometary ices (Bockelée-Morvan et al. 2000), which are thought to present chemical similarities with the ices processed during hot cores formation (Irvine et al. 2000). Moreover, as hydrogenation is the most effective process in grains chemistry, models often predict that when sulphur atoms stick on the grain, they will consequently form H₂S (Garrod et al. 2007). However, H₂S has never been detected in dark cloud ices. Additionally, upper limits to its column density towards three high-mass protostars and three late-type stars lying behind TMC-1 (CP) have been derived from observations by Smith (1991). These limits, despite large uncertainties, are believed to be too small to account for the missing sulphur in dark clouds and alternative hypotheses have recently been proposed (Jiménez-Escobar & Muñoz Caro 2011; Druard & Wakelam 2012; Martín-Doménech et al. 2016).

As seen in Section 4.1, the results of simulations with our enhanced sulphur network seem to support both of the main hypotheses about the reservoirs of sulphur, depending on the age of the observed dark cloud. Indeed, according to our results, if the age of the cloud is less than $5 \times 10^5$ yr then atomic sulphur in the gas phase is the main reservoir of sulphur. However, if the age of the cloud is greater, then our results show that because of hydrogenation of sulphur on the grain surface (the cycle between HS and H₂S enabled by

the reaction-diffusion competition and the accumulation of these species in the bulk), the reservoirs of sulphur are nearly equally HS and H₂S in the grain bulk. The upper limit to the abundance of H₂S in ices provided by the observations is between $3 \times 10^{-7}$ and $3 \times 10^{-6}$ (assuming an abundance of H₂O in ices of $10^{-4}$ and based on Smith 1991). It should be noted that these upper limits are derived using the infrared transition strengths obtained by Ferraro, Sill & Fink (1980) which has been obtained in pure H₂S ice. The intensity strengths can be notably different in interstellar ice analogues. As the total abundances of H₂S on the grains in Models 3 and 4 at the time of 'best fit' $10^6$ yr are, respectively, $1.5 \times 10^{-6}$ and $4.4 \times 10^{-6}$, it appears that our enhanced network results are in agreement with the upper limits of Smith (1991) and therefore in favour of a joint HS and H₂S sulphur reservoir in icy grain bulks in dark clouds. This result is in agreement with the results depicted in appendix E of Furuya et al. (2015). Indeed, they also find that HS and H₂S are likely to be the main reservoirs of sulphur, notably due to the same hydrogenation loop described in equation (4). Constraints on the icy grain bulk abundance of HS would be useful to validate this hypothesis. However, the $v = 0 \rightarrow v = 1$ transition of HS lies around $2599 \text{ cm}^{-1}$ in the gas phase (Bernath, Amano & Wong 1983) and has never been observed in ice. The calculated IR strength (at DFT or MP2 level using Gaussian) is three orders of magnitude less efficient than water and four order of magnitude less efficient than OCS. Subsequently, even if solid HS is the reservoir of sulphur in dark clouds, the lack of sensibility of IR detection may prevent its observation.

## 6.3 About the observability of HCS

In order to assess the observability of HCS in the gas phase of TMC-1 (CP) (see Section 4.3.3), we proceed using a local thermodynamic equilibrium calculation of the HCS emission spectrum in the 3- and 2-mm atmospheric windows with $T_{ex} = 10$ K, a column density of HCS of $1 \times 10^{12} \text{ cm}^{-2}$ (corresponding to a gas-phase abundance of HCS of approximately $10^{-10}$ as obtained with Model 4 in the time frame between $10^5$ and $10^7$ yr), assuming a wide source (no beam dilution) and a linewidth of 1 km s⁻¹. We estimate that the brightest line would be the one at 80 553.516 MHz with a main beam temperature of 2 mK. This result infers that for a $5\sigma$ detection of HCS in the gas phase, it would take more than a thousand hours of observations with the IRAM 30 m. We therefore conclude that HCS is unlikely to be detected in the gas phase in dark clouds even if its abundance grows notably above $10^{-12}$ in our simulations.

# ACKNOWLEDGEMENTS

This work has been founded by the European Research Council (Starting Grant 3DICE, grant agreement 336474). The authors are grateful to the CNRS programme 'Physique et Chimie du Milieu Interstellaire' (PCMI) for partial funding of their work.

# REFERENCES

* Adande G. R., Halfen D. T., Ziurys L. M., Quan D., Herbst E., 2010, ApJ, 725, 561

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

* Anicich V. G., 2003, JPL Technical Report, JPL-Publ-03-19, NASA, Washington, DC

* Baulch D. et al., 2005, J. Phys. Chem. Ref. Data, 34, 757

* Benson S., 1976, Thermochemical Kinetics: Methods for the Estimation of Thermochemical Data and Rate Parameters, 2nd edn. Wiley, New York
Bernath P., Amano T., Wong M., 1983, J. Mol. Spectrosc., 98, 20

Bockelée-Morvan D. et al., 2000, A&A, 353, 1101

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

Boogert A. C. A., Gerakines P. A., Whittet D. C. B., 2015, ARA&A, 53, 541

Cernicharo J. et al., 2011, A&A, 531, A103

Chabot M. et al., 2010, A&A, 524, A39

Chang Q., Cuppen H. M., Herbst E., 2007, A&A, 469, 973

Charnley S. B., 1997, ApJ, 481, 396

Druard C., Wakelam V., 2012, MNRAS, 426, 354

Federman S. R., Sheffer Y., Lambert D. L., Gilliland R. L., 1993, ApJ, 413, L51

Ferraro J. R., Sill G., Fink U., 1980, Appl. Spectrosc., 34, 525

Florescu-Mitchell A. I., Mitchell J. B. A., 2006, Phys. Rep., 430, 277

Fournier J. A., Shuman N. S., Melko J. J., Ard S. G., Viggiano A. A., 2013, J. Chem. Phys., 138, 154201

Frerking M. A., Linke R. A., Thaddeus P., 1979, ApJ, 234, L143

Furuya K., Aikawa Y., Hincelin U., Hassel G. E., Bergin E. A., Vasyunin A. I., Herbst E., 2015, A&A, 584, A124

Galland N., Caralp F., Rayez M.-T., Hannachi Y., Loison J.-C., Dorthe G., Bergeat A., 2001, J. Phys. Chem. A, 105, 9893

Galland N., Caralp F., Hannachi Y., Bergeat A., Loison J.-C., 2003, J. Phys. Chem. A, 107, 5419

Garozzo M., Fulvio D., Kanuchova Z., Palumbo M. E., Strazzulla G., 2010, A&A, 509, A67

Garrod R., Herbst E., 2006, A&A, 457, 927

Garrod R. T., Pauly T., 2011, ApJ, 735, 15

Garrod R. T., Wakelam V., Herbst E., 2007, A&A, 467, 1103

Georgievskii Y., Klippenstein S. J., 2005, J. Chem. Phys., 122, 194103

Georgievskii Y., Klippenstein S. J., 2007, J. Phys. Chem. A, 111, 3802

Gibb E., Nummelin A., Irvine W. M., Whittet D. C. B., Bergman P., 2000, ApJ, 545, 309

Glowacki D. R., Liang C.-H., Morley C., Pilling M. J., Robertson S. H., 2012, J. Phys. Chem. A, 116, 9545

Gonzalez M., Hijazo J., Novoa J. J., Sayos R., 1996, J. Chem. Phys., 105, 10999

Graedel T. E., Langer W. D., Frerking M. A., 1982, ApJS, 48, 321

Gratier P., Majumdar L., Ohishi M., Roueff E., Loison J. C., Hickson K. M., Wakelam V., 2016, ApJS, 225, 25

Gronowski M., Kolos R., 2014, ApJ, 792, 89

Halfen D. T., Ziyurys L. M., Brünken S., Gottlieb C. A., McCarthy M. C., Thaddeus P., 2009, ApJ, 702, L124

Hama T., Watanabe N., 2013, Chem. Rev., 113, 8783

Hamberg M. et al., 2007, Mol. Phys., 105, 899

Hamberg M. et al., 2014, J. Phys. Chem. A, 118, 6034

Hasegawa T. I., Herbst E., Leung C. M., 1992, ApJS, 82, 167

Hatchell J., Thompson M. A., Millar T. J., MacDonald G. H., 1998, A&A, 338, 713

Herbst E., Terzieva R., Talbi D., 2000, MNRAS, 311, 869

Hickson K., Loison J., Cavalié T., Hébrard E., Dobrijevic M., 2014, A&A, 572, A58

Hickson K. M., Loison J.-C., Nunez-Reyes D., Méreau R., 2016a, J. Phys. Chem. Lett., preprint (arXiv:1608.08877)

Hickson K. M., Wakelam V., Loison J.-C., 2016b, Mol. Astrophys., 3, 1

Hincelin U., Wakelam V., Hersant F., Guilloteau S., Loison J. C., Honvault P., Troe J., 2011, A&A, 530, A61

Hippler H., Viskolcz B., 2002, Phys. Chem. Chem. Phys., 4, 4663

Holdship J. et al., 2016, MNRAS, 463, 802

Irvine W. M., Schloerb F. P., Crovisier J., Fegley B., Jr, Mumma M. J., 2000, Protostars and Planets IV. Univ. Arizona Press, Tucson, AZ, p. 1159

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

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

Karssemeijer L. J., Cuppen H. M., 2014, A&A, 569, A107

Kolesniková L., Tercero B., Cernicharo J., Alonso J. L., Daly A. M., Gordon B. P., Shipman S. T., 2014, ApJ, 784, L7

Korth M., Grimme S., 2009, J. Chem. Theory Comput., 5, 993

Kurylo M. J., Peterson N. C., Braun W., 1971, J. Chem. Phys., 54, 943

Lada E. A., Bally J., Stark A. A., 1991, ApJ, 368, 432

Leonori F., Occhiogrosso A., Balucani N., Bucci A., Petrucci R., Casavecchia P., 2011, J. Phys. Chem. Lett., 3, 75

Leonori F., Balucani N., Nevrly V., Bergeat A., Falcinelli S., Vanuzzo G., Casavecchia P., Cavallotti C., 2015, J. Phys. Chem. C, 119, 14632

Lilenfeld H. V., Richardson R. J., 1977, J. Chem. Phys., 67, 3991

Linke R. A., Frerking M. A., Thaddeus P., 1979, ApJ, 234, L139

Lique F., Cernicharo J., Cox P., 2006, ApJ, 653, 1342

Loison J.-C., Wakelam V., Hickson K. M., Bergeat A., Mereau R., 2014, MNRAS, 437, 930

Loison J.-C. et al., 2016, MNRAS, 456, 4101

McGonagle D., Irvine W. M., Ohishi M., 1994, ApJ, 422, 621

Majumdar L., Gratier P., Vidal T., Wakelam V., Loison J.-C., Hickson K. M., Caux E., 2016, MNRAS, 458, 1859

Martín-Doménech R., Jiménez-Serra I., Muñoz Caro G. M., Müller H. S. P., Occhiogrosso A., Testi L., Woods P. M., Viti S., 2016, A&A, 585, A112

Matthews H. E., MacLeod J. M., Broten N. W., Madden S. C., Friberg P., 1987, ApJ, 315, 646

Mendes M. B. et al., 2012, ApJ, 746, L8

Messineo M. et al., 2015, ApJ, 805, 110

Montaigne H. et al., 2005, ApJ, 631, 653

Neufeld D. A., Wolfire M. G., Schilke P., 2005, ApJ, 628, 260

Nobes R. H., Radom L., 1981, Chem. Phys., 60, 1

Oehlers C., Wagner H. G., Ziemer H., Temps F., Dobe S., 2000, J. Phys. Chem. A, 104, 10500

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

Peng J., Hu X., Marshall P., 1999, J. Phys. Chem. A, 103, 5307

Penzias A. A., Solomon P. M., Wilson R. W., Jefferts K. B., 1971, ApJ, 168, L53

Peters P. S., Duflot D., Wiesenfeld L., Toubin C., 2013, J. Chem. Phys., 139, 164310

Plessis S., Carrasco N., Dobrijevic M., Pernot P., 2012, Icarus, 219, 254

Plume R., Jaffe D. T., Evans N. J., II, Martín-Pintado J., Gómez-González J., 1997, ApJ, 476, 730

Podio L. et al., 2015, A&A, 581, A85

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

Puzzarini C., 2005, J. Chem. Phys., 123, 024313

Reiter D., Janev R., 2010, Contrib. Plasma Phys., 50, 986

Rice B. M., Chabalowski C. F., 1994, J. Phys. Chem., 98, 9488

Rice B. M., Cartland H. E., Chabalowski C. F., 1993, Chem. Phys. Lett., 211, 283

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

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

Sabbah H., Biennier L., Sims I. R., Georgievskii Y., Klippenstein S. J., Smith I. W., 2007, Science, 317, 102

Sakai N. et al., 2014, Nature, 507, 78

Sander S. et al., 2011, JPL Technical Report, JPL-Publ-10-6, NASA, Washington, DC

Scott G. B., Fairley D. A., Freeman C. G., McEwan M. J., Anicich V. G., 1999, J. Phys. Chem. A, 103, 1073

Scott G. B., Milligan D. B., Fairley D. A., Freeman C. G., McEwan M. J., 2000, J. Chem. Phys., 112, 4959

Shalabiea O. M., 2001, A&A, 370, 1044

Shannon R. J., Blitz M. A., Goddard A., Heard D. E., 2013, Nature Chem., 5, 745

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

Smith I. W., 2006, Angew. Chem. Int. Ed., 45, 2842

Smith I. W. M., 2011, ARA&A, 49, 29

Sofia U. J., Cardelli J. A., Savage B. D., 1994, ApJ, 430, 650

Tieftrunk A., Pineau des Forets G., Schilke P., Walmsley C. M., 1994, A&A, 289, 579

Tielens A., Hagen W., 1982, A&A, 114, 245

Viti S., Caselli P., Hartquist T. W., Williams D. A., 2001, A&A, 370, 1017

Wagner A. F., Bowman J. M., 1987, J. Phys. Chem., 91, 5314

Wakelam V., Herbst E., 2008, ApJ, 680, 371

Wakelam V., Caselli P., Ceccarelli C., Herbst E., Castets A., 2004, A&A, 422, 159
Wakelam V., Herbst E., Selsis F., 2006, A&A, 451, 551

Wakelam V. et al., 2010, Space Sci. Rev., 156, 13

Wakelam V. et al., 2012, ApJS, 199, 21

Wakelam V. et al., 2015a, ApJS, 217, 20

Wakelam V., Loison J.-C., Hickson K. M., Ruaud M., 2015b, MNRAS, 453, L48

Wang H., Eyre J., Dorfman L. M., 1973, J. Chem. Phys., 59, 5199

Wierzejewska M., Moc J., 2003, J. Phys. Chem. A, 107, 11209

Woods P. M., Occhiogrosso A., Viti S., Kaňuchová Z., Palumbo M. E., Price S. D., 2015, MNRAS, 450, 1256

Woon D. E., Herbst E., 2009, ApJS, 185, 273

Yoshimura M., Koshi M., Matsui H., Kamiya K., Umeyama H., 1992, Chem. Phys. Lett., 189, 199

Zhao Y., Truhlar D. G., 2008, J. Chem. Theory Comput., 4, 1849

## SUPPORTING INFORMATION

Supplementary data are available at MNRAS online.

**Supplementary material (file for review).pdf**

**Supplementary material (online).pdf**

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the article.

This paper has been typeset from a TeX/Latex file prepared by the author.