Modeling Sulfur Depletion in Interstellar Clouds

The elemental depletion of interstellar sulfur from the gas phase has been a recurring challenge for astrochemical models. Observations show that sulfur remains relatively non-depleted with respect to its cosmic value throughout the diffuse and translucent stages of an interstellar molecular cloud, but its gas-phase constituents cannot account for this cosmic value towards higher-density environments. We have attempted to address this issue by modeling the evolution of an interstellar cloud from its pristine state as a diffuse atomic cloud to a molecular environment of much higher density, using a gas/grain astrochem. code and an enhanced sulfur reaction network. A common gas/grain reaction network has been systematically updated and greatly extended based on previous lit. and models, with a focus on the grain chemistry and processes. A simple model was used to benchmark the resulting network updates, and the results of the model were compared to typical astronomical observations sourced from the literature. Our new gas/grain model is able to reproduce the elemental depletion of sulfur, whereby sulfur can be depleted from the gas-phase by two orders of magnitude, and this process may occur under dark cloud conditions if the cloud has a chemical age of at least 1 Myrs. The resulting mix of sulfur-bearing species on the grain ranges across all the most common chemical elements (H/C/N/O), not dissimilar to the molecules observed in cometary environments. Notably, this mixture is not dominated simply by H2S, unlike all other current astrochem. models. Despite our relatively simple physical model, most of the known gas-phase S-bearing molecular abundances are accurately reproduced under dense conditions, however they are not expected to be the primary molecular sinks of sulfur. Our model predicts that most of the missing sulfur is in the form of organo-sulfur species trapped on grains.


Introduction
Sulfur poses an interesting challenge to models of interstellar chemistry. Within primitive interstellar environments, it is known that sulfur remains in ionized atomic form and close to the cosmic abundance (Jenkins 2009). However, in molecular clouds and star-forming regions, this cosmic abundance is drastically reduced from the gas-phase molecular inventory. Many simple organo-sulfur species can be detected, but not seemingly enough to fully account for its cosmic abundances.
Modern astrochemical models still today predict that the bulk of sulfur resides as condensed H 2 S, despite upper limits from observations (Smith 1991;van der Tak et al. 2003;Jiménez-Escobar & Muñoz Caro 2011). The effective workaround has been to use a severely depleted value of the elemental abundance ( 1% with respect to the cosmic standard abundance) for the initial sulfur content when modeling dense interstellar environments. This has resulted in a rather poor understanding of sulfur with respect to interstellar matter, and the inability to correctly model sulfur remains a severe shortcoming of astrochemical models. That's not to say that astrochemical models of sulfur are not improving. For example, a recent model (Woods et al. 2015) helps to provide constraints on how much sulfur may be locked in the refractory residue that is known from laboratory Full Tables A.1 and B.4 are only available at the CDS via anonymous ftp to cdsarc.u-strasbg.fr (130.79.128.5) or via http: //cdsarc.u-strasbg.fr/viz-bin/qcat?J/A+A/624/A108 studies of photo-processed sulfur; the reaction network pertaining to certain categories of molecules has also been expanded by the study reported by Vidal et al. (2017); and, a better understanding of dynamical models is provided by Vidal & Wakelam (2018) to help put into perspective the apparent chemical variability of sulfur chemistry across cloud evolution.
Laboratory studies are also continuously helping to shed light on sulfur chemistry, which can be notoriously complex. Sulfur bonds are generally not as strong as those of the firstand second-row elements (i.e., H, C, N, O), and processing of interstellar analogs of ice mixtures containing simple S-bearing molecules can yield a highly heterogeneous mixture of products (e.g., Jiménez-Escobar et al. 2014, and references therein). Qualitatively, these processed ices sometimes even resemble the type of chemistry that has been detected in both cometary ices (Calmonte et al. 2016) and meteoritic material (see, e.g., Ehrenfreund et al. 2002). However, many important laboratory studies have not yet been incorporated into a modern gas-grain astrochemical model.
We have set out to update the sulfur chemistry within a modern astrochemical reaction network in a systematic and more complete way than has been done in the past, and then use this gas-grain reaction network to model the evolution of an interstellar cloud from its pristine diffuse stage to a dense and dark quiescent state to check for clues pertaining to the gas-phase depletion of sulfur. A collection of observational studies has also been used to provide general constraints for benchmarking this updated astrochemical model to ensure reasonable results for all the sulfur-bearing interstellar molecules that are typically observed in interstellar environments. We have found that this updated sulfur network does in fact reproduce observations of interstellar sulfur chemistry in dense environments, and this relatively simple astrochemical model also reproduces the severe gas-phase depletion of sulfur thanks to the production a variety of (mostly-stable) organo-sulfur molecules on grains (see Sect. 3).
In Sect. 2, we first present the astrochemical model in terms of its chemical and physical characteristics. In Sect. 3, we present our results, beginning with a look at the question of elemental depletion of sulfur during the dense stage (m gas ≈ 10 4 −10 6 cm −3 ) of an interstellar cloud, and then followed up with a brief analysis of commonly observed molecules across a range of evolutionary stages and densities of interstellar clouds with our model predictions. In Sect. 4 we provide a summary comparing our model with previous studies and to discuss key results. Finally, Sect. 5 contains some general conclusions.

Computer code
The various model trials presented here have been performed using a new python-based code, which has not yet been described elsewhere 1 . The code is based on the OSU gas-grain astrochemical model presented by Garrod et al. (2008), which is a timedependent, single-point, rate-based model for general application to interstellar cloud environments. Benchmarks were performed during its development to ensure that its operation is consistent with the code it is based on, and it has been further enhanced with a number of internal consistency checks and software tools for the development of a chemical network and for exploring its results.

Chemical network
The core reaction network was based on the OSU gas-grain model, which consists of ca. 8300 reactions and 860 chemical species, and remains in high circulation even a decade after being released. Aside from a major revamping of the sulfur chemistrywhich we address in a separate section below -a number of modifications were made to the chemical inputs to reflect recent literature reports. These changes include: (a) modifications to the thermodynamical properties of many grain surface species; (b) updated photochemical rates as per the report from Heays et al. (2017); and, (c) a number of other updates to various reactions. For completion, all the binding energies and heats of formation for the grain species are listed in Table A.1; except for all the updates from Heays et al. (2017), we also list out new and/or modified reactions in Table B.1. In the particular, we have added routes for the heavy atoms (Fe + , Mg + , Na + , S + , Si + , Cl + , and P + ) to adsorb directly onto negatively-charged grains, as discussed in Ruffle et al. (1999) and references therein.
Whereas a number of photodesorption reactions have been added based on Hollenbach et al. (2009), we have completely disabled/ignored chemical (i.e., reactive) desorption. The new photodesorption entries are crucial for accurate results during Stages 1 and 2 (i.e., the diffuse and translucent cloud phases), however, recent studies suggest that chemical desorption is largely not efficient in ices (Minissale et al. 2016;Chuang et al. 2018) and additionally introduces an added complexity which 1 https://bitbucket.org/notlaast/ggackal/ Notes. (a) We define the fractional abundance X(i) as the ratio of the abundance of species i wrt the total gas-phase hydrogen number density, which is initially equal to X(H) + 2X(H 2 ). (b) The total electron fraction is simply initialized as the sum of the gas-phase cations. References.
we felt detracts too much focus from development on the chemical network and the study of sulfur depletion. We have used, when possible, known cosmic standard elemental abundances for the initial abundances of our model (i.e., for Stage 1), and these values are contained in Table 1. Whereas the details of the stages are presented below, in Sect. 2.3, we note here that the initial fractional abundances of subsequent stages are simply the final abundances of the preceding stages.

Sulfur chemistry
In order to allow for a number of chemical routes that a high (cosmic) sulfur abundance may take, we have sorted through a wealth of literature data to systematically build a more complete gas-grain sulfur network. These new reactions were based on both experimental studies relating specifically to interstellar ice analogs and also to atmospheric chemistry. Furthermore, sulfur is known for its rich chemistry, and behaves more similarly to carbon than to oxygen, despite being found below the latter in the periodic table of the elements. The sulfur chemistry in our network takes into account oxidation states and a plurality of isomers and bond types.
As with the other heavy atoms and what we have pointed out above, our model includes a route for direct accretion of S + onto negatively-charged grains. To counterbalance this effect, we found that it was important to add also the photodesorption of H 2 S, so that it behaves similarly to a number of other simple species -in particular, CH 4 , CO, and H 2 O.
In Tables B.2-B.4, we list all of the new reactions which have supplemented the base network. In general, many of the included rates have been estimated from scratch, based on similar chemistry (yielding reasonably large uncertainties) or calculated from basic electrostatics data and capture rate theory (Troe 1985). We also provide highly detailed explanations or notes for these collections of reactions, which can be found also within the Appendices. A108, page 2 of 17  (d) 10 7 10 6 10 5 -10 6 Notes. (a) Classification is based on Snow & McCall (2006); see text for details. (b) This the H 2 column density external to our single-point calculations, whereby the total column density used for the shielding factor is the sum of the external and "internal" components; the component of external N(CO) is always set to a factor of 10 −6 below that of the external N(H 2 ), as per general agreement with observations (Snow & McCall 2006, and references therein). (c) Equivalent to that of Rho Oph (Indriolo & McCall 2012

Physical model
The physical model we have employed is defined by three distinct, physically-static stages. Whereas the code is capable of time-dependent densities and temperatures, we have opted to use only static physics for each stage, to avoid significant uncertainties related to the dynamical evolution of molecular clouds. Here we focus on the effects of the new chemistry for comparison with previous models, which are also mainly static; a dynamicalchemical model is deferred to future work. Stage 1 follows the evolution of a diffuse cloud during its atomic-to-molecular phase transition, with the final time of this transition being defined as the moment when less than 50% of the original atomic hydrogen content remains. Stage 2 simulates a translucent cloud, whose defining characteristics are its increased visual extinction compared to diffuse clouds, as well as the conversion of ionized carbon to molecular species. Finally, Stage 3 is meant to simulate the aging of a dense core, and we have included a range of densities in order to explore the effect of density on chemical differentiation and timescales, as this is known to be a varying parameter across astrophysical environments and should play a critical role on the rate of sulfur depletion. We employ the convention of resetting the defined time to t = 0 at the beginning of each stage, and the initial abundances of Stages 2 and 3 are set to those which have been predicted in the preceding stage at the times defined in Table 2.
These three stages adhere well to the categorization of interstellar clouds presented by Snow & McCall (2006, Sect. 2), and their defining characteristics are summarized in Table 2. We also show below ( Fig. 1) that this sequence of stages is able to approximate well the main astrochemical results supported by observations during the early stages of molecular cloud evolution, agreeing particularly well with Fig. 1 from Snow & McCall (2006).
Whereas the model is a single-point calculation, the column densities and radii listed above are used for determining the total effective column densities of H 2 and CO, which are then used to determine on-the-fly their self-shielding factors for photodissociation. For the work presented here, we have used parameterizations based on Draine & Bertoldi (1996, Eq. (37)) and Lee et al. (1996), for H 2 and CO, respectively. Besides the stage-dependent physical parameters listed above, a large number of other parameters are involved behind the scenes, particularly for reactions on the grain. These additional parameters are standard values nowadays, but for completeness' sake we list them in Table 3.

General early-stage features
In the plots shown in Fig. 1, we present the time evolution of key astrochemical species during the early stages of cloud evolution, showing consistent agreement with the defining chemical characteristics summarized by Snow & McCall (2006, Sect. 2, Fig. 1) for the diffuse-to-translucent phase transition of an interstellar cloud. The single significant exception with respect to Snow & McCall (2006) is the ionization fraction, which does not deplete so quickly in our model, continuing to be supported by the heavy atoms Mg + , Si + , S + , and Fe + after carbon has been significantly neutralized. This steadier ionization fraction is somewhat consistent with models of translucent clouds (A v ≤ 2 mag) (e.g., Herbst & Leung 1986;. While none of these results are particularly novel, they clearly benchmark our gas-grain astrochemical model against a general consensus of astronomical observations, providing reasonable context for the comparisons of S-bearing species presented below.

Sulfur depletion and its primary budget
Regarding the early stages of interstellar cloud evolution, we show in Fig. 2 that sulfur remains fully in the atomic ionized form during the diffuse cloud, and it is eventually partially neutralized (ca. 20%) late into the translucent stage. It is only during the late stage of a more dense environment that it finally naturally depletes (>99%) from the gas-phase, into a diverse spread of simple molecules. To illustrate this, we have collected stacked histogram-like plots in Fig. 3 showing the molecules that dominate the sulfur budget during Stage 3, where the horizontal axis tracks the time, and the vertical component represents the abundance of each species as a percentage of the total sulfur composition. We have also overlaid these plots with a semitransparent box to denote the gravitational free-fall time of a cloud at their respective densities. At the left of these plots, it can be seen that sulfur is primarily atomic at early times, but then H 2 S and other grain species become increasingly abundant A108, page 3 of 17 A&A 624, A108 (2019) as freeze-out takes over. At short timescales, S adsorbs onto the grain and is rapidly hydrogenated to H 2 S, however, only to a limited extent of ∼7%. This upper limit agrees very well with the non-detections of solid H 2 S (e.g., Smith 1991) toward a number of interstellar environments, and contrast starkly with all other gas-grain astrochemical models. This trend also agrees well with laboratory studies that have shown that solid H 2 S can be rapidly converted to a variety of other species (Garozzo et al. 2010;Chen et al. 2015); we discuss some of its chemistry in more detail below (Appendix C.1.2). An intriguing observation that should also be noted is the detection of solid OCS toward IRC-L1041-2, reported by Aikawa et al. (2012), yielding a fractional abundance of approximately 1.6 × 10 −6 relative to hydrogen. This is an incredibly high fraction of the total abundance of cosmic elemental sulfur, but is indeed a value within the range of our high-density models between gas densities of 10 4 cm −3 and 10 5 cm −3 . In this case, there is clear evidence that a significant fraction of the elemental sulfur is trapped in the form of a condensed species, and our model clearly suggests the rest may be found across a variety of other S-bearing organics.
We show in our model that simple freeze-out can account for the appropriate timescale of the gas-phase depletion of sulfur. The major result is that our expanded reaction network produces a significant amount of organo-sulfur molecules, which prevents the build-up of H 2 S on the grain. The direct adsorption of positively-charged ions onto negatively-charged grains has been added. This enhances the freeze-out kinetics by a factor of 50-100% (not shown here). We speculate that the organosulfur chemistry together with the freeze-out can explain why it has been so challenging to reconcile interstellar sulfur depletion from an observational perspective. We predict that sulfur gas-phase depletion does not significantly occur at low densities; however, it does occur rapidly ( 10 5 yr) at densities higher than 10 4 cm −3 , producing solid organo-sulfur molecules that are hard to detect.

Comparison to observations
For the three evolutionary stages of interstellar clouds on which we have based our model, there have been a number of molecular surveys that provide general observational constraints of S-bearing species. These observations have all been summarized below, in Table 4, and are also illustrated in various figures below, alongside related model results. In this section we look at these constraints in more detail through comparisons with our model.

Diffuse clouds
A comparison of the model results against observations of diffuse clouds reported by Neufeld et al. (2012Neufeld et al. ( , 2015 is shown in the plots of Fig. 4. Despite a cosmic initial abundance of elemental sulfur (i.e. X 0 (S + ) = 1.66 × 10 −5 ), all gas-phase S-bearing species that we focus on are underestimated to varying degrees with respect to observations. From the figure, it can be seen that the observed relative abundances of the species agree in a qualitative manner with the model results. Quantitatively, these species are all under-produced in the model, by a factor of 10. The significant underproduction of HS + , as well as the underproduction of HCS + and SO + during A108, page 4 of 17       Fig. 3. Primary S-carrying species from the molecule as a function of time and different volume densities (from 10 4 cm −3 , top, to 10 6 cm −3 , bottom) during Stage 3. The vertical width of each component represents how much sulfur is accounted for. We have also used partially opaque boxes to highlight the timescale that is defined by the free fall time of a homogeneous cloud to collapse under ideal gravitational conditions (t ff ≈ 3/2πGρ 0 , i.e., ignoring magnetic fields, angular momentum, and other resistive effects), for reference to other models. The gas number densities, m, and approximate free fall time limits, t ff , are noted for each plot.
A108, page 5 of 17 Transl. (6) 6.5 × 10 −9 Dense (7) Notes. (a) The environments are classified accordingly: "Diffuse", "Transl.", and "Dense" refer to diffuse (A v ≤ 1), translucent (1 A v 2), and dense (A v > 5) interstellar clouds, respectively. the translucent stage (see below) suggests a shortcoming of our simplified physical model for ion chemistry during the early stages of cloud evolution, which we speculate could be due to a number of reasons, such as our simplistic homogeneous medium, our simplistic model of monodisperse dust grains, or our neglect of energy injection via shocks and turbulence dissipation. We have not attempted to improve the model during this early stage, and, instead, refer the reader to previous models of ion chemistry within these environments (e.g., Neufeld & Dalgarno 1989;Scappini et al. 2000;Le Petit et al. 2004;Wakelam & Herbst 2008;Neufeld et al. 2015;Ivlev et al. 2018).

Translucent clouds
A thorough survey of many S-bearing species in translucent cloud environments has been presented in Turner (1995aTurner ( , 1996a and Turner et al. (1998). These reports have provided self-consistent constraints that we can use for comparison with our model, which we show in Fig. 5. In general, the model performs much better at predicting gas-phase abundances of S-bearing species for the translucent stage compared to the dif- fuse stage, however, most species remain underpredicted. Of these species, CS is predicted at the highest abundance within the model, at a level of nearly 10× above the upper range of abundances reported by Turner (1996a) before t = 1 Myr, but then settling within this range afterwards. The two species observed at the highest abundances by Turner (1995aTurner ( , 1996b. and Turner et al. (1998) are H 2 S and SO, which are both underpredicted by only a factor of a few around t = 1 Myr compared to observations. C 2 S is predicted at nearly the value reported by Turner et al. (1998), however it was detected in only a small number of sources with respect to total number, and this value therefore serves as more like an upper limit. The remaining species -H 2 CS, SO + , OCS, HCS + and SO 2 -are all severely underpredicted by the model compared to their ranges of reported abundances. Again, we speculate that our underpredictions might simply be due to our overly-simplistic model of low-density regions that is lacking certain processes important for diatomics and molecular ions.

Dense environments
For our comparisons to dense clouds, it was desired to focus on a single environment, though few systematic studies of large sets of S-bearing species could be identified. It was necessary, therefore, to combine a number of datasets. TMC-1 is a cold, dark cloud that has been very well-characterized regarding both its chemical and physical characteristics, making it one of the best matches for our Stage 3 model at a density of 10 4 cm −3 . These observational constraints are shown alongside our model results through the time plot in Fig. 6.
For the model trial at a density of 10 4 cm −3 , all molecules are predicted in very good agreement with observations at times greater than 10 6 years. Some species are clearly very sensitive to the "chemical age" of the model, particularly SO 2 , which has an enhanced abundance at early times of the dense stage that better agrees with translucent cloud observations (Turner 1995a) rather than toward TMC-1 (Cernicharo et al. 2011). On the other hand, the relative abundances of SO and SO 2 (SO/SO 2 ≈ 5−10) are consistent with the range that is typical for many other interstellar observations (e.g., van der Tak et al. 2003), though their observed abundances do appear to vary significantly from source to source. A108, page 6 of 17   We can also compare the relative abundance ratios of various S-bearing species against more commonly observed molecular species. In Table 5 we compare the gas-phase abundances of several S-bearing species as a ratio with respect to CO, where the values have been chosen at times that give approximately the peak abundance for the species shown in Fig. 6. These abundance ratios vary significantly across three interstellar environments, but each modeled density of our dense stage gives reasonable agreement with these values. Furthermore, OCS has been a well-studied component of interstellar ices via its absorption band near 2040 cm −1 . It is the only S-bearing species detected in interstellar ices, and our model predicts that its formation is driven primarily by S-atom addition to CO (see Appendix C.4.2 for more details). In Table 6 we com- Notes. Values are in terms of percentage with respect to CO gas . Observed values are determined from the following: (a) Observed values for the S-bearing species are based on those referenced in Table 4 and Fig. 6, and scaled with respect to a value of f (CO) = 8 × 10 −5 as per Irvine et al. (1987). (b) All values are taken from Table 1 of Charnley et al. (2001). <0.8 n = 10 4 cm −3 0.7 t = 2 Myr n = 10 5 cm −3 4 t = 0.3 Myr n = 10 6 cm −3 5 t = 0.03 Myr Notes. Observed values are in terms of percentage with respect to CO ice , and are taken from Table 2 of Palumbo et al. (1997). pare a summary of OCS ice observations with respect to solid CO toward a variety of protostars, where our selected points in time match those of the previous table (Table 5).

Discussion
Our gas-grain model incorporates a number of general updates to the parent OSU gas-grain model (Garrod et al. 2008) pertaining to photochemistry so that the model may also be used at earlier stages of cloud evolution, namely low-A v conditions (0.5 A v (mag) 2). These photochemistry updates have allowed us to include an explicit third stage for the aging of a diffuse cloud, and the inclusion of the earlier stages and cosmic abundances provides the opportunity to check for elements of sulfur depletion or unexpected effects caused by our changes. In general terms, however, our model predictions do not perform well for minor gas-phase species at earlier stages, but we are confident that the underprediction of many gas-phase sulfurcontaining species does not significantly affect the major results regarding sulfur depletion during the dense stage. This has been confirmed by alternatively running the later stages with initial abundances of select species fixed to better match observational values of the preceding stages, but their abundances rapidly (within 10 2 yr) converge to the same results at intermediateand late-times (not shown here). We have also adopted a higher binding energy for O surf and a lower (compared to the previous OSU model) binding energy for NH 3 surf under the guidance of recent laboratory studies, and these differences will affect basic ice chemistry of water and ammonia at later stages.
Our model considers a number of novel approaches with regard to sulfur chemistry, to better bring it into alignment with other chemistries, such as carbon, nitrogen, and oxygen. In general terms, our new model provides the most expansive astrochemical network for sulfur to date. Unlike many previous models, we have made broad updates across the entire sulfur reaction network, and have sourced as much as possible all previous laboratory and theoretical studies of sulfur chemistry in an unbiased manner, including the correction of many erroneous rate coefficients and activation barriers that were used in the past prior to these more recent studies.
Regarding sulfur allotropes (sulfur chains and rings), we have made significant updates to the network and molecular parameters, to be consistent with the literature about solid sulfur, so that they are treated just like any other chemical species and may interact with the reaction network through a number of processes. This is unlike other studies which have used a more heuristic approach with only a subset of reactions and processes to test for their importance under specific conditions (e.g., Wakelam et al. 2005;Woods et al. 2015). On the other hand, these heavy species do not appear active in our model, at such a low temperature and quiescent physics.
We have also newly added a number of species related to CS 2 and the C x S carbon-sulfur chains, many of which have been both studied in the laboratory and detected in interstellar environments, but not yet considered in a general astrochemical model. In the case of CS 2 , it appears to be much more important in cometary ices than interstellar environments, though its interstellar significance may be underestimated due to its detection being significantly impeded without a permanent dipole moment. Nevertheless, we have included it and many of its hydrogenated derivatives. As is the case with the oxygen counterparts, the singly-hydrogenated forms of C 2 S and C 3 S are expected to be important molecules relative to their parents. With the exception of ethyl mercaptan (CH 3 CH 2 SH), our model therefore includes all known interstellar sulfur-bearing species to date, and is therefore the most complete astrochemical reaction network to date.
Finally, unlike recent models of interstellar sulfur, we show that sulfur can be depleted from its gas-phase cosmic abundance onto the grain at late-times during the dense molecular phase, and this grain composition is not dominated by H 2 S. To prevent this unrealistically large build-up of solid H 2 S, it was important to include its enhanced photodissociation on the grain (compared to its gas-phase rate), its photodesorption under conditions of low-A V , and the extensive grain reactions that provide more routes to a generally more diverse range of sulfur chemistry.

Conclusion
We have modeled sulfur chemistry using an expanded gas-grain astrochemical model across the evolution of an interstellar cloud. Our model suggests that the vast majority of the sulfur eventually becomes locked up in the condensed phase and across a variety of simple organic molecules, which may explain the ongoing challenge of detecting enough gas-phase sulfur-bearing species to account for its cosmic abundance. Our model demonstrates that an enhanced gas-grain astrochemical model of sulfur chemistry can provide insight into the observed sulfur depletion from the cosmic gas-phase abundance, without breaking agreement with many gas-phase S-bearing species which models have generally already succeeded.
Our model suggests that sulfur depletion from the gas phase does occur under dark cloud conditions if the cloud has an age greater than 10 6 years. This process is caused by rapid freezeout at enhanced densities, which also results in highly nonequilibrium chemistry that appears to be highly sensitive to the gas density and timescale of the dense cloud stage, resulting in a variable molecular inventory that may present a challenge in modeling specific interstellar environments.
Our model provides no evidence of sulfur chains or rings forming a significant refractory component, and alternatively suggests that elemental sulfur is converted to a range of simple, stable organo-sulfur molecules that are primarily formed there. Many of these species resemble those which have been detected in cometary environments, however some differences also exist, such as our negligible reservoir of solid H 2 S. During the preparation of this manuscript, a number of new reports related to interstellar sulfur molecules have showed up in the literature (i.e., HCS, S 2 H, HSO, and NS + ), and their detections and upper limits have also been found to agree very well with our model. This suggests that our work is a significant step forward for astrochemical models of sulfur, and we hope we have provided enough information to serve as a stepping stone for future studies.
We have provided a summary of key formation reactions behind the full range of S-bearing interstellar molecules, many of which are new to a gas-grain astrochemical model. The variety of reactions presented here reflects the diverse chemistry of sulfur, highlighting the ability of sulfur to react with a variety of types of molecules/elements -including itself -and across a range of oxidation states. Unsurprisingly, many reaction rates presented here are uncertain, leaving room for much improvement by way of laboratory and/or theoretical studies.  In Table A.1, we show all the adopted binding energies for all species that may interact/reside on the grain. These values have been the subject of many recent studies, and we have tried to implement a number of suggestions from the literature. In the many cases where the values have only been guessed, we have included a note or formula behind the estimation. In most cases, it has been our experience that similar formulae for many complex molecules are consistent with measured and/or theoretical values.

Appendix B: Reaction tables
In the tables below, we list all the modifications introduced to the reaction network, with respect to the original OSU gasgrain astrochemical network (Garrod et al. 2008). We have separated these changes out into many different tables according to whether they pertain to sulfur or not, as well as the type of modification (removal, modification, or new additions). For the tables listing chemical reactions, we have tried to include as much information as possible to help with their adoption, including the reaction itself (i.e. the reactions and products), the kinetics that define the rate coefficient, the form of the rate law, and finally any notes or references that were used for defining the entry. The chemical reactions include many entries that are entirely based on educated guesses, and were synthesized based on a number of general chemical/physical principles, which we describe here.

B.1. Rate laws
Following the conventions established already decades ago, the reaction network is made up of classes (or types) of reactions that can be categorized by either the types of reactants/products involved or the forms of the rate laws that define their individual reaction rates. In the reaction tables below (Tables B.1-B.4), we have sorted the reactions according to the reaction types, which we designate as their "rtype". The list below can be used to look up the rate laws that define their rate coefficients, k x , according to their "rtype" (in subscript):

B.2. Gas-phase kinetics: Electrostatics
Estimates for gas-phase rates involving electrostatics are based on two general assumptions, depending on the specific case. For the adsorption of atomic cations onto negatively-charged grains, we have simply copied the rates previously defined for the collisional rates involving the neutralization reactions of the atomic ions and negative grains. This was done for all the heavy atoms (Fe, Mg, Na, S, Si, Cl, and P) that have relatively large binding energies, in agreement with the study by Umebayashi & Nakano (1980). On the other hand, capture rates between molecular ions and neutral molecules have been calculated explicitly for each case that lacks an external reference. Capture rates have been estimated using the formulae generalized by Troe (1985), where we include both the classic Langevin rate constant, k L , as well as the dipole-enhanced rate constant, k D . The Langevin rate constant is defined by where q is the charge of the ion, α is the isotropic polarizability of the neutral molecule, and µ is the center of mass of the two species. The dipole-enhanced rate constant is defined by where the permanent dipole moment µ D of the neutral species may also contribute to an enhanced, temperature-dependent collision rate, with k referring to the Boltzmann constant, and T being the kinetic temperature. An additional factor C that is less than one must be used to compensate for the imperfectly-aligned dipole of a rotating gas-phase molecule: At the most rigorous level of this theory, one might consider how factors such as temperature, rotation constants, internal state populations, and partition functions might all affect the factor C.
In practice, it appears that this value ranges somewhere between ca. 0.386 and 1, and only nearer to unity for very light molecules (i.e., not S-bearing molecules) at very low temperatures (≤10 K). We therefore just use a value of C = 0.4. The determination of k L and k D may easily be estimated from only a few molecular properties, most of which are readily available from external references or easily calculable with computational chemistry A108, page 10 of 17 packages. We have used Psi4 2 for the estimation of all proton affinities, the online service Chemicalize 3 for estimates of molecular polarizabilities, and, when possible, the CDMS and JPL databases for molecular spectroscopy 4 for dipole moments, otherwise our own simple Psi4 calculations. The uncertainties for these calculations are expected to be similar to the alreadypresent entries: ±50%. The specific selection of which ion-neutral reactions to include is roughly based on which of those were already present for the O-substituted species, with the additional qualification that proton exchange reactions (i.e. A(H) + + B → A + B(H + )) are allowed only when the relative proton affinities of the reacting species yields an exothermic reaction. Besides this manual selection of reactions that were included, the inclusion was also extended through a number of routines built in to our astrochemical code to ensure that all new neutral species reacted with a minimal set of highly-abundant ions, and that all new ionic species could be neutralized through recombination. 2 http://www.psicode.org/ 3 https://chemicalize.com/ 4 https://cdms.astro.uni-koeln.de/ and https://spec. jpl.nasa.gov/

B.3. Gas-phase kinetics: neutral-neutral rates
There are many reactions that satisfy the basic requirements for inclusion in our model (endothermic, a low or missing activation barrier, at least two products so one may serve as a kinetic/energetic sink), but have not been suitably targeted by theory or experiment to provide accurate gas-phase kinetics. Due to the fact that gas-phase rates for neutral-neutral reactions are notoriously difficult to estimate accurately, we have not attempted to include such reactions. In many of these cases, we have already included notes in the text above as suggestions for future studies.

B.4. Photochemistry
Not all molecular species that were addressed in this study were subjected to previous quantitative laboratory or theoretical photochemical studies. However, it is critical that all species are able to interact with cosmic-rays and/or energetic photons, and this was ensured through a built-in routine that checks for certain inadequacies within the chemical network. Therefore, we had to make a number of guesses on reaction channels and cross sections of cosmic-rays and/or UV photons. As a general rule of thumb, all the molecules considered in this study were at A108, page 11 of 17 Notes. (a) The "rtype" of each reaction can be used to look up their respective rate coefficient using the listing shown in Appendix B.1.

References.
(1) Isoniemi et al. (1999); (2)  least addressed either directly or indirectly by particle or photon bombardment studies, and therefore relative cross sections of known species could at least be estimated.
We also note that a number of photodissociation and photoionization reactions were updated according to the suggestions by Heays et al. (2017). In those cases, we have also included the reaction coefficient with which the dust-attenuation factor may be calculated using a second-order exponential integral in place of the standard exponential function, resulting in a rate law that is instead equal to k 13,19 = αE 2 (−βA v ) (s −1 ). (B.4) This alternative equation improves the fitted rates across a wider range of magnitudes and was found to be advantageous for our multi-stage model.

B.5. Grain reaction rates
The rates for surface reactions are computed in the same manner as that of Garrod et al. (2008), including allowing for tunneling to be considered in case the reaction involves the light species H or H 2 and an activation barrier. In general terms, the thermallycontrolled reaction rate is defined under the scheme presented in Hasegawa et al. (1992). That is, each species has a characteristic diffusion rate based on its mass, the grain temperature, and diffusion-and binding-energies, and the reaction rates proceed according to these classical diffusion rates, activation barriers or endothermic reactions excepted. Endothermic reactions are strictly prohibited in our model, and there is logic coded in to check against all the surface reactions, both old and new. Indeed there are some cases where our tabulated binding energies and A108, page 12 of 17 Cl + + grain − → grain + Cl surf 0 2.90e−17 0.5 0 Mass-scaled from Fe Fe + + grain − → grain + Fe surf 0 2.30e−17 0.5 0 Duplicated from neutral atom grain − + Mg + → grain + Mg surf 0 3.70e−17 0.5 0 Duplicated from neutral atom grain − + Na + → grain + Na surf 0 3.60e−17 0.5 0 Duplicated from neutral atom grain − + P + → grain + P surf 0 3.10e−17 0.5 0 Mass-scaled from Fe grain − + S + → grain + S surf 0 3.00e−17 0.5 0 Duplicated from neutral atom grain − + Si + → grain + Si surf 0 3.30e−17 0.5 0 Duplicated from neutral atom CS 2 + cr → CS + S 1 1.70e+03 0.0 0 Copied from CO2 CSSH + cr → CS 2 + H 1 1.70e+03 0.0 0 Duplicated from CS2 HC 2 S + cr → CH + CS 1 1.50e+03 0.0 0 Duplicated from HC2O Notes. This version of the table only contains the first ten entries. The full version is only available at the CDS. (a) The "rtype" of each reaction can be used to look up their respective rate coefficient using the listing shown in Appendix B.1. References.
(1) Isoniemi et al. (1999); (2)  heats of formation combine to yield an endothermic reaction, and make note of this in the tables below for such cases.

Appendix C: Details of chemistry
C.1. Sulfur allotropes and polysulfanes C.1.1. Pure S One of the most significant aspects of sulfur in the condensed phase is that it exhibits the most allotropes (i.e., different forms of a pure element) known to exist for any single element in the periodic table (e.g., Steudel & Eckert 2003). Whereas discussion of heavy sulfur allotropes is typically limited to ambient pressures and at relatively high temperatures ( 300 K), there are clues supporting their importance in astrophysically-relevant ices (Muñoz Caro 2002;Jiménez-Escobar & Muñoz Caro 2011). There have been instances where such chemistry was considered within the context of astrochemical modeling (e.g. Wakelam et al. 2005;Woods et al. 2015;Vidal et al. 2017), but not yet in an explicit and systematic manner. For all species S x , we have treated them as normal gas-grain species like all others present in the model. The diradical chains (x = 2−4) may form through barrierless atomaddition reactions on the grain, and the rings (x = 5−8) may form only through combination reactions of the smaller chains (i.e. diradical-diradical ring closures). Their binding energies are assumed to scale with the number of sulfur atoms, and therefore these aforementioned grain processes are not likely to be efficient at dust temperatures of 10-15 K. Unfortunately there is also a lack of information regarding photochemistry cross sections. Therefore, we have made the assumption that the photostabilities increases with size of the chains and rings (i.e., the cross sections decrease).
The only species within this category that reaches an appreciable gas-phase abundance is S 2 , reaching a fractional abundance X(S 2 ) ≈ 10 −9 after 0.2 Myr during the translucent cloud stage. Its abundance is nearly constant after this, including during the dense stage, except for a lull at intermediate times under dense conditions, when its formation cannot keep up with its rapid destruction by C and O atoms. On the grain, S 3 reaches the highest abundance among the pure allotropes, to a fractional abundance of ∼10 −11 after 0.2 Myr during the translucent cloud stage. Its rotational spectrum is relatively well-studied (Thorwirth et al. 2005), but its detection would prove challenging at such a low abundance. We had hoped to see more significant abundances of these pure allotropes, but it appears that even the small chains are too heavy and/or too sticky (high binding energies/diffusion barriers) for thermal roaming to react with each other before alternative processes destroy them.

C.1.2. (Poly)sulfanes (S x H y )
Aside from the pure allotropes of sulfur, a number of hydrogenated forms may also readily co-exist and are collectively known as sulfanes or polysulfanes. The simplest species of this class are HS and HS + , but it was only relatively recently that they were detected in the ISM (Neufeld et al. 2012;Menten et al. 2011, respectively), specifically within diffuse interstellar clouds, at fractional abundances of ≤4 × 10 −9 . Our model does not predict an agreeable abundance of HS during the diffuse stage in comparison with Neufeld et al. (2015), however our model does reach X(HS) ≈ 10 −9 late into the translucent cloud stage and it remains there for much of the dense cloud stage. On the other hand, HS + is not predicted at a significant abundance at all during our model's early stages. Savage et al. (2004) have suggested two possible formation routes for HS + : S + + H 2 → HS + + H, H 2 S + + H → HS + + H 2 , however, the former requires an activation barrier of nearly 10 000 K (Millar et al. 1987), and the latter relies on atomic hydrogen and H 2 S + , neither of which are highly abundant during the dense stage and this may explain the low prediction of HS + . This shortcoming of our cold, dense model lends support to the hypothesis that shock chemistry may play a crucial role to form these species (Neufeld et al. 2015).
H 2 S becomes an important molecule in the gas-phase during the translucent cloud stage and for much of the dense cloud A108, page 13 of 17 stage, with a near-constant fractional abundance of ∼10 −9 . In the gas-phase, the main formation route of H 2 S during the translucent cloud stage (i.e. low-A v ) is photodesorption of its grain component, which was missing completely from previous astrochemical models, and then during the dense stage from H 2 S + via non-dissociative recombination with electrons and chargetransfer reactions involving a slew of other heavy atoms. On the grain, H 2 S is formed rapidly through successive hydrogenation of elemental sulfur and is therefore an important primary constituent of sulfur on the grain. During the early stages of sulfur depletion, it begins to account for up to ca. 10% of the total sulfur budget, it never builds up to more significant levels, as cosmic rays keep its abundance in check while a number of recently-adsorbed species -notably atomic carbon -scavenge it and convert it to other species. The limited peak abundance of solid H 2 S agrees very well with interstellar non-detections, but stands in stark contrast with the cometary value reported by Calmonte et al. (2016), who suggest that it accounts for the majority of the sulfur composition in the bulk ice of comet 67P/Churyumov-Gerasimenko. Additionally, we would like to point out that very recent results (Oba et al. 2018) suggest that H 2 S formation in ice may also have a non-negligible reactive desorption, which may drive down its abundance on the grain even lower than our model predicts.
H 2 S + itself can be formed quite efficiently from the radiative association reaction, S + + H 2 → H 2 S + + hν, which was already included in the OSU network, as well as The reaction above is introduced by Vidal et al. (2017) as a second route with a relatively small branching ratio in comparison to the route that yields the products H 2 /HS + . On the grain, H 2 S may form rapidly through association of the HS radical with atomic hydrogen. The ionic (H 2 S + ) and protonated (H 3 S + ) forms are also predicted to reach non-negligible gasphase abundances (∼10 −11 ) at late-times during the dense stage. These ions may prove fruitful targets for astronomical detection, providing better constraints for astrochemical modeling, and in the case of H 2 S + , a good target for additional laboratory spectroscopy.
Two other important (poly)sulfanes are S 2 H and H 2 S 2 (also sometimes written as HSSH), the former of which has now recently been reported in the interstellar medium (Fuente et al. 2017). These species have long been suspected as interstellar candidates, because of the importance of H 2 S and how readily it may be converted to alternative chemistry upon irradiation (Jiménez-Escobar & Muñoz Caro 2011). Whereas our model predicts a near-constant abundance of S 2 after a specific time during the translucent cloud stage, our model predicts only a significant gas-phase abundance of S 2 H and H 2 S 2 during the dense stage, and their abundances are predicted to nearly match each other at a peak fractional abundance of ∼10 −10 , which is comparable the detection reported by Fuente et al. (2017) toward the Horsehead PDR. S 2 H always reaches a higher peak abundance in the gas-phase compared to the grain, as it lacks an efficient formation route on the grain, however H 2 S 2 does reach a higher abundance in the ice by a factor of 10. In the gas-phase, S 2 H and H 2 S 2 are expected to form from a similar, but rather convoluted ionic pathway: H + 3 + S 2 → S 2 H + + H 2 , S 2 H + + H 2 → H 3 S + 2 + hν, where the final dissociative recombination reaction is estimated to have a branching ratio of 1:1 between the two sets of products. This final reaction originates from the OSU gas-grain model, but we could not find any experimental or theoretical studies to support this particular branching ratio; any update to this product branching ratio would have a significant impact on the relative abundances of S 2 H and H 2 S 2 in the gas phase. Under our model conditions, the ions involved in the aforementioned pathway are not predicted at high abundance, but they aren't well-studied in the laboratory and knowledge of their rotational spectra might yield some surprises under certain interstellar conditions. On the grain, S 2 H has no significant formation pathway except for the photo-and CR-induced dissociation reactions of H 2 S 2 . H 2 S 2 primarily forms on the grain from the reaction, which has been shown to be barrierless despite involving a closed-shell species and a number of transient intermediates (Zhou et al. 2008). However, the reaction itself is not particularly efficient because of the reactants' high surface binding energies (ca. 1100 K and 2700 K, respectively), and it benefits only from their high abundances. H 2 S 2 may also form from the reaction of the HS radical with itself, but this reaction must compete with the more rapid hydrogenation of HS to form H 2 S.

C.2.1. CS
CS is one of the most abundant S-bearing molecules under all physical conditions considered. It has significant abundances in translucent and dense stages, and is always the most abundant C-bearing sulfur species in our model. During the translucent and dense stages, our model predicts a fractional abundance of 10 −8 −10 −7 at intermediate times, which then drops at later times, and its abundance at these later times agrees well with observations. On the grain, its abundance during the translucent stage is similar to the gas phase, but it is enhanced an additional factor of ten during the dense stages. Its gas-phase formation is primarily (∼80% during translucent stage, ≤30% during the dense stage) through the series of ion-neutral reactions: CH + S + → CS + + H, CS + + H 2 → H + HCS + , In the dense stage, its gas-phase formation is also enhanced by atomic carbon stripping sulfur from SO and OCS. Its grain-surface formation is largely driven by the reaction between atomic C and the HS radical, however its abundance also significantly relies on the accretion of gas-phase CS, as well as radical-radical reactions that result in fragmentation. The CS + cation is predicted to reach a fractional abundance of ∼10 −11 during intermediate times of the translucent cloud stage, however its lower dipole moment and much lower abundance compared to its neutral counterpart would likely hamper an astronomical detection.

C.2.2. HCS, H 2 CS
At high densities, HCS and H 2 CS are also abundant C-bearing sulfur species in the ice, and the latter is also present at a significant level in the gas-phase, reaching peak fractional abundances ≥10 −9 . HCS is formed in the ice from the well-studied (see Deeyamulla & Husain 2006, and refs therein) insertion reaction between atomic C and H 2 S. HCS may also reach reasonably high abundances in the gas-phase at intermediate densities (∼10 4 ), via two neutral-neutral reactions: The recent detection of HCS toward L483 at a fractional abundance of 2 × 10 −10 (Águndez et al. 2018) agrees perfectly with our dense stage model at a gas-phase density of 10 4 cm −3 . On the grain, H 2 CS relies significantly on the availability of HCS, whereby H 2 CS may form directly from the hydrogenation of HCS. In the gas-phase, it may derive from two possible routes: the ion-neutral/DR reaction combination that is important while the ionization fraction remains high: At late timescales, it may also form directly from neutralneutral displacement reaction,

C.2.3. Methyl mercaptan (CH 3 SH)
Methyl mercaptan is the S-bearing analog of methanol (CH 3 OH), and as such, much of its chemistry is similar in our model. In the gas-phase, CH 3 SH is severely lacking in abundance in our dark cloud model, with a peak abundance a factor of nearly 10 5 below that of CH 3 OH. Its most efficient gas-phase formation route in the model is through the ion-neutral radiative association reaction, CH + 3 + H 2 S → CH 3 SH + 2 + hν, followed by successive dissociative recombination, CH 3 SH + 2 + e − → CH 3 SH + H. This analogous route for CH 3 OH is much less efficient than the ion-neutral reaction, However the analogous sulfur reaction of the latter route, as well as the species HS − , are both missing from our model due to uncertainties behind them.
In the ice, on the other hand, our model shows that CH 3 SH can be efficiently formed already during the translucent stage, at a fractional abundance >10 −10 ; and in the most dense trials, up to 10 −7 . Its grain surface formation is synonymous with that of CH 3 OH: it likely forms through successive hydrogenation of the CH 2 SH and CH 3 S radicals. These precursor radicals are also synonymous with the O-bearing counterparts, whereby they may form from hydrogenation of H 2 CS through high barriers (800 K and 1200 K to form CH 2 SH and CH 3 S, respectively, Vidal et al. 2017), or from the barrierless radical-radical reactions, CH 2 + HS → CH 2 SH, CH 3 + S → CH 3 S.
Our predicted peak abundances of CH 3 SH in the ice are only slightly higher than its gas-phase detections toward dense protostellar environments (ca. 10 −10 −10 −8 , Gibb et al. 2000;Müller et al. 2016;Majumdar et al. 2016), in agreement with previous suggestions that it is predominantly formed on the surface. Without a process to release it into the gas phase, it may not be an important component of gas toward dark clouds.

C.2.4. Carbon-sulfur chains (H x C y S)
A great variety of mixed carbon chains are known to exist in the ISM, including those containing sulfur. The OSU network was already quite developed for this class of molecules, but we have also added a significant number of reactions. As can be seen in the plots within Fig. 6, our model predictions for C 2 S are in good agreement with observations during the dense cloud stage at lower densities. On the other hand, C 2 S appears to be quite sensitive to density, dropping off quickly and below the plotted range at higher densities.
In the gas phase, C 2 S is formed primarily via the somewhat convoluted ionic pathway, S + + C 2 H → C 2 S + + H, C 2 S + + H 2 → HC 2 S + + H, None of these intermediate cations reach appreciable abundances. On the grain, C 2 S mostly forms via the neutral-neutral reaction, CH + CS → H + C 2 S, benefiting from the high abundances of CH and CS.
We also note that our model predicts a significant amount of gas-phase HC 2 S during the dense stage, even more than C 2 S. This result is similar to the report by Águndez et al. (2015) and their observation of more HC 2 O than C 2 O. During the dense stage, HC 2 S is formed solely through the reaction C + H 2 CS ↔ H + HC 2 S; Yamada et al. (2002) suggests a number of other potential neutral-neutral gas-phase routes that may enhance its abundance even further, however the kinetics of these additional reactions remain so highly uncertain that we have left them out of the gas-phase portion of the network.
Larger chains have somewhat different formation routes, relying heavily on rapid neutral-neutral reaction even in the gas phase, in line with the investigations reported by Yamada et al. (2002) and references therein. The gas-phase formation of C 3 S is partially driven by a similar ionic pathway as that of C 2 S, however approximately half its formation is also due to contributions from the following neutral-neutral reactions during the dense stage: On the grain, only the reaction C + C 2 S → C 3 S dominates. C 4 S forms in a much simpler way, whereby the reaction S + C 4 H → C 4 S + H A108, page 15 of 17 accounts for its formation in both the gas phase and on the grain, though this process may not be efficient enough to provide an easy interstellar detection. There is laboratory data for the larger chains (i.e., five or more carbon atoms, Gordon et al. 2001Gordon et al. , 2002, but we have left them out of our network; both our model and observations show that the gas-phase abundance rapidly drops off for larger chains.

C.2.5. Carbon disulfide (CS 2 )
There has been much interest in CS 2 raised by laboratory studies that have observed its formation within irradiated ice mixtures intended to simulate interstellar processes (e.g., Ferrante et al. 2008;Garozzo et al. 2010), and we have therefore introduced this species into the network. In the gas-phase, its presence would be difficult to detect, as its structure lacks a permanent dipole moment. In the condensed phase, it reaches a non-negligible fractional abundance at late times during the translucent cloud stage (>10 −10 ), and becomes a significant sink of elemental sulfur at late times of the dense stage. It can form from the following reactions, where the first entry is important only briefly during early times of the dense stage, and the latter two routes account for the majority of its formation, at late times of both the translucent and dense cloud stages.

C.2.6. HCSSH
We have also introduced the molecule dithioformic acid (HCSSH) to our gas-grain network, due to its focus in a recent laboratory study (Prudenzano et al. 2018). Although not much is known about this species in terms of kinetics and thermodynamics, its similarity to formic acid (HCOOH) provides a number of insights that may be used for generating a minimal reaction network.
In the gas phase, it may be formed either directly through the neutral-neutral reaction, HS + H 2 CS → HCSSH + H, in a similar way to that of HCOOH, or also the ion-neutral reaction, HCSSH + 2 + e − → HCSSH + H. The latter pathway is much less important for HCSSH than it is for HCOOH, as its protonated form does not reach as high abundance as its O-containing counterpart. On the grain, it may form from a series of neutral-neutral addition reactions via radicals: HS + CS → CS 2 H, H + CS 2 H → CS 2 H.
Its grain surface formation appears to be more efficient than the gas-phase route, with its gas-phase fractional abundance reaching ca. 10 −13 at the moment when other sulfur species peak, and its grain abundance being mostly constant at 10 −12 throughout the dense stages. Despite its similarity to HCOOH, HCSSH does not seem to be as important for sulfur chemistry as HCOOH is for oxygen chemistry.

C.3.1. NS and derivatives
The NS radical is the most important interstellar nitrogen/sulfur species, having been detected in a variety of environments, including a PDR (Leurini et al. 2006), a molecular outflow (Sánchez Contreras et al. 2000) and several dense regions (Gottlieb et al. 1975;McGonagle et al. 1994). During the translucent cloud stage, it forms in the gas phase primarily via the reaction N + SH → NS + H, which also drives its formation on the grain surface. At higher densities and after nitrides have been established, the ionic pathway involving the reactions, NH + 2 + S → HNS + + H, HNS + + e − → NS + H, also becomes somewhat (up to 30%) important, although the former reaction is a significant bottleneck because NH + 2 can much more easily react with H 2 , up to 10 4 times more rapidly. On the grain, it becomes an important sink for sulfur in our model, accounting for up to 7-10% of the total sulfur budget.
The ubiquity of NS suggests that other interstellar species may easily derive from it. The closest related species is the NS + cation, which has recently been studied in the laboratory and detected in the ISM toward B1-b as reported by Cernicharo et al. (2018). NS + is a closed-shell ion that may form from ionmolecule reactions between NS and a number of cations serving as electron scavengers due to sulfur contributing to a relatively low electron affinity compared to the higher-row elements. Our model predicts a marginal fractional abundance (<10 −13 ) during the translucent cloud stage, primarily through charge-exchange reactions with atomic cations, C + + NS → NS + + C, H + + NS → NS + + H.
Our model predicts peak fractional abundances of NS + in the range of ca. 10 −12 -10 −11 , which is in very good agreement with its detection toward B1-b. We also predict that NS + might also form directly from photoionization of NS, since NS has an estimated ionization potential of ∼8.87 eV (Dyke et al. 1977), making it an easy target in PDRs. This latter process has, however, yet to be studied in either the laboratory or theory.

C.3.2. HNCS/HSCN/HCNS isomeric family
The isomeric family containing isothiocyanic acid (HNCS), thiocyanic acid (HSCN), and thiofulminic acid (HCNS) has received some focus from the astrochemical community. Of the three isomers listed here, HNCS is the lowest-energy species, and the first to be detected in the ISM (Frerking et al. 1979). The slightly less stable isomer HSCN was finally detected much more recently at an abundance slightly below that of HNCS (Halfen et al. 2009). The two least stable isomers, HCNS and HSNC, have not yet been reported in the ISM, though it was only very recently that their rotational spectra were characterized sufficiently to merit an interstellar search (McGuire et al. 2016).
The chemical network pertaining to these species was adapted solely from that of Vidal et al. (2017), as these molecules A108, page 16 of 17 have not been well-studied in the laboratory. It should also be noted that the least stable isomer, HSNC, remains missing from the astrochemical model. Given these circumstances, there is little we can add to the understanding of these isomers, except to state that they are not well-modeled with our model. Because we do not consider reactive desorption in our model -in contrast to the work reported by Vidal et al. (2017) -this highlights some interesting behavior in the formation of these isomers. In the gas-phase, HCNS is produced most efficiently in comparison to the other isomers, HNCS is predicted at an abundance approximately tenfold below HCNS, and HSCN is nearly another tenfold lower. On the grain, on the other hand, HNCS and HSCN have identical formation routes and (estimated) surface properties, and therefore their abundances match at all times. During the dense stage, HNCS and HSCN are enhanced compared to HCNS at early times (i.e. much earlier than the free-fall time), and this role is reversed at later times. It is conceivable that reactive desorption could significantly alter the gas-phase abundances of these isomers, causing behavior that reproduces that of Vidal et al. (2017). On the other hand, the chemistry of these species are not yet well understood, and they could significantly benefit from further physico-chemical laboratory studies.

C.4.1. SO and SO 2
With the exception of the diffuse cloud stage, SO is predicted to be the most abundant O-bearing sulfur species in the gas phase, however our model does not match many observations, being underpredicted during the early stages and overpredicted at higher densities. During the diffuse stage, its primary formation mechanism in our model is the gas-phase reaction, O + HS → SO + H, and its destruction is driven by photodissociation and photoionization. During the translucent cloud stage, these processes still dominate, but its destruction is now also governed by the reaction with C + , until later times, when carbon has been sufficiently neutralized. During the dense stage, the reaction OH + S → H + SO also becomes an important route. Whereas gas-phase SO is an important S-bearing molecule, the majority of SO can actually be found on the grain, where the aforementioned neutral-neutral reaction may also take place, as well as the simple direct addition reaction of S and O atoms. This latter process is no longer as efficient as in previous models, due to the large binding energy of the oxygen atom. SO 2 is also an abundant O-bearing sulfur species in the gas-phase, albeit only at higher densities. In general, it peaks in gas-phase abundance approximately ten times lower than SO, in general agreement with most astronomical observations. It may form directly from SO, through the radiative association reaction O + SO → SO 2 + hν.
On the grain, SO 2 plays a lesser role at such low temperatures, as this latter reaction is still possible but suffers from the high binding and diffusion energies of both O and SO. O 2 is slightly more reactive, and provides a formation route via O 2 + SO → SO 2 + O, but not at significant levels.
The radical HSO is a derivative of SO and also an unexpectedly important oxygen-containing species in our model despite the recent report of its gas-phase non-detection toward several extraterrestrial environments (Cazzoli et al. 2016). In our model it does not reach appreciable abundance in the gas-phase, however it appears to collect significantly on the grain, as it is modeled to form readily from the hydrogenation of SO. Its peak abundance during the dense stage and at late times accounts for 5-20% of the sulfur budget, which is admittedly quite high for an open-shell species.

C.4.2. Carbonyl sulfide (OCS)
In the dense stage, our model consistently predicts a peak fractional abundance of ∼2 × 10 −9 in the gas-phase, which sits squarely within the range of values observed in interstellar environments. Its formation derives from the following reactions (in order of decreasing contribution): HCS + O → H + OCS, HCO + S → H + OCS, CS + OH → H + OCS.
Its destruction is dominated by carbon atoms that cleave the C−O bond, yielding CO and CS radicals. On the grain, it is also one of the primary sinks of elemental sulfur, forming directly from the association of CO and sulfur atoms. Its destruction is controlled simply by photodissociation from cosmic rays, though there is non-negligible contribution from sulfur atoms scavenging to form S 2 and CO at later times, after atomic sulfur is able to deplete onto the grain.