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

# SiS Formation in the Interstellar Medium through Si+SH Gas-phase Reactions

V. C. Mota¹, A. J. C. Varandas¹,²,³, E. Mendoza⁴, V. Wakelam⁵, and B. R. L. Galvão⁶

¹Departamento de Física, Universidade Federal do Espírito Santo, 29075-910 Vitória, Brazil

²School of Physics and Physical Engineering, Qufu Normal University, 273165, People's Republic of China

³Coimbra Chemistry Centre and Chemistry Department, University of Coimbra, 3004-535 Coimbra, Portugal

⁴Observatório do Valongo, Universidade Federal do Rio de Janeiro, Ladeira Pedro Antônio, 43, Rio de Janeiro—RJ 20.080-090, Brazil

⁵Laboratoire d'astrophysique de Bordeaux, Univ. Bordeaux, CNRS, B18N, allée Geoffroy Saint-Hilaire, F-33615 Pessac, France; valentine.wakelam@u-bordeaux.fr

⁶Centro Federal de Educação Tecnológica de Minas Gerais, CEFET-MG, Av. Amazonas 5253, 30421-169, Belo Horizonte, Minas Gerais, Brazil

brenogalvao@gmail.com

Received 2021 June 10; revised 2021 July 19; accepted 2021 July 27; published 2021 October 11

## Abstract

Silicon monosulfide is an important silicon-bearing molecule detected in circumstellar envelopes and star-forming regions. Its formation and destruction routes are not well understood, partially due to the lack of detailed knowledge on the involved reactions and their rate coefficients. In this work we have calculated and modeled the potential energy surface (PES) of the HSiS system employing highly accurate multireference electronic structure methods. After obtaining an accurate analytic representation of the PES, which includes long-range energy terms in a realistic way via the DMBE method, we have calculated rate coefficients for the Si+SH → SiS+H reaction over the temperature range of 25–1000 K. This reaction is predicted to be fast, with a rate coefficient of ~1 × 10⁻¹⁰ cm³ s⁻¹ at 200 K, which substantially increases for lower temperatures (the temperature dependence can be described by a modified Arrhenius equation with α = 0.770 × 10⁻¹⁰ cm³ s⁻¹, β = -0.756, and γ = 9.873 K). An astrochemical gas-grain model of a shock region similar to L157-B1 shows that the inclusion of the Si+SH reaction increases the SiS gas-phase abundance relative to H₂ from 5 × 10⁻¹⁰ to 1.4 × 10⁻⁸, which perfectly matches the observed abundance of ~2 × 10⁻⁸.

Unified Astronomy Thesaurus concepts: Astrochemistry (75); Reaction rates (2081); Chemical abundances (224)

## 1. Introduction

Silicon monosulfide is the simplest molecule to contain a silicon-sulfur chemical bond and may play an important role in the formation of sulfide grains in the interstellar medium (ISM; Zhukovska & Gail 2008; Massalkhi et al. 2019). After its first detection in the envelope of an evolved carbon star (Morris et al. 1975), it has been observed in several other environments (Dickinson & Kuiper 1981; Danilovich et al. 2018; Lefloch et al. 2018), together with its isovalent counterpart, SiO. From these observations, SiS is consistently found to be less abundant than SiO, and although the origin of silicon monoxide is satisfactorily understood, this is not the case for SiS. Evidence for a strong gradient of [SiO/SiS] ratio across a shock region has been obtained (Podio et al. 2017), indicating different chemical pathways. As a result of a systematic search for SiS in low-mass star-forming regions (Lefloch et al. 2018), it was concluded that SiS cannot be a simple product of gas-phase chemistry in molecular clouds, but rather requires specific conditions, which are best found in shocks. In circumstellar envelopes, on the other hand, SiS is believed to be formed in the densest and hottest parts of the inner envelope, where it is thought to be formed under thermochemical conditions (Schöier et al. 2007; Prieto et al. 2015).

Given the importance of the SiS molecule as a potential gas-phase precursor to sulfide grains, its unknown chemistry, and its poorly explained abundance, several neutral-neutral formation and destruction reactions have been explored recently, both theoretically and experimentally. For example, Zanchet et al. (2018) proposed that one of the main SiS formation reactions in outflows is Si+SO. Further theoretical calculations showed that Si+SH, SiH+S, and SiH+S₂ collisions (Rosi et al. 2018, 2019; Paiva et al. 2020) are also potential candidates. As for destruction

routes, it has been computationally shown that SiS is stable toward collisions with both atomic and molecular hydrogen (Paiva et al. 2020) but is efficiently destroyed in reaction with atomic oxygen (Paiva et al. 2018; Zanchet et al. 2018). Recent laboratory studies (Doddipatla et al. 2021), together with astrochemical models, have also highlighted the importance of the reaction Si+SH₂ → SiS + H₂ in star-forming regions, in spite of its nonadiabatic character.

Recent astrochemical models predict that the Si+OH reaction is one of the major sources of interstellar SiO in the gas phase (Gusdorf et al. 2008; Rivero Santamaría et al. 2017). It should also be interesting to verify whether the analogous Si+SH reaction could also be a major source of SiS, but accurate values for its rate coefficients are necessary for this verification. For this reason, we develop in this work an accurate potential energy surface (PES) for the HSiS triatomic system, which allows for the calculation of accurate rate coefficients for SiS formation via S+SiH and Si+SH reactions, both of which are exothermic and barrierless. Quasi-classical trajectories are employed to obtain accurate rate constants for the title reaction, which are employed in an astrochemical model to study its role in shock regions.

## 2. Theoretical Methods

### 2.1. Electronic Structure Calculations

The product channel SiS(¹Σ⁺)+H(²S) correlates with the PES of a single electronic state of the HSiS molecule, namely, ²A'. In turn, the ²A'' and quartet states correlate with the excited SiS(a³Π) molecule (Paiva et al. 2020), which lies much higher in energy and does not contribute to formation of SiS. For this reason, only the ground ²A' PES has been calculated and modeled in this work.
The calculations have been performed with the MOLPRO package (Werner et al. 2015) employing the internally contracted version of the MRCI (Werner & Knowles 1988) method including single and double excitations and the Davidson correction (MRCI +Q). The full valence CASSCF wave function (Werner & Knowles 1985) was used as a reference. Specifically, the 1s, 2s, and 2p orbitals of Si and S were considered as the core, and the remaining nine orbitals of the system were treated as active with 11 electrons (11e/9o). A large number of single-point energies covering the whole configuration space of the molecule were calculated with the aug-cc-pV (Q + d)Z and aug-cc-pV (5 + d)Z basis set (Dunning 1989; Kendall et al. 1992), which have been subsequently extrapolated to the complete basis set (CBS; Varandas 2018) limit.

The CASSCF energies were extrapolated with the physically motivated exponential scheme (Feller 1993; Varandas 2019)

$$E_X^{HF} = E_\infty^{HF} + B \exp(-\beta X), \quad (1)$$

with B and $E_\infty^{HF}$ being parameters to be determined from a fit to a pair of CASSCF/aug-cc-pVXZ calculations. Following Varandas (2019), $\beta$ is set to 2.284, and a re-hierarchization procedure has been employed that converts $X=Q$, 5 into $x=4.41$, 5.34 to allow for an accurate dual-level CBS extrapolation of the CASSCF energies.

The dynamical correlation (dc) energy is CBS extrapolated according to the USTE method (Varandas 2007b), in which the dc energy of the $X=Q$, 5 pair of basis sets is fitted to

$$E_X^{dc} = E_\infty^{dc} + \frac{A_3}{(X - 3/8)^3} + \frac{A_5^{(0)} + cA_3^n}{(X - 3/8)^5}, \quad (2)$$

which provides $E_\infty^{dc}$ and $A_3$. The parameters $A_5^{(0)}$, c, and n are system independent and constant for a given post-HF method (Varandas 2007b). Their values for MRCI calculations are $A_5^{(0)} = 0.0037685459 E_h$, $c = -1.17847713 E_h^{1-(n)}$, and $n = 1.25$ (Varandas 2007b).

## 2.2. Potential Energy Surface Modeling

The CBS extrapolated ab initio points are then fitted to an analytic function using the DMBE formalism (Varandas 1985, 1988, 1992, 2000). Accordingly, the PES assumes the form

$$V(\mathbf{R}) = \sum_{i=1}^{3} V_i^{(2)}(R_i) + V_i^{(3)}(R_{SiS}, R_{SiH}, R_{SH}), \quad (3)$$

where $V_i^{(2)}(R_i)$ are two-body energy terms and $V_i^{(3)}$ a three-body one. Both are further partitioned to describe separately short-range ($V_{EHF}^{(n)}$) and long-range ($V_{dc}^{(n)}$) interactions. The latter play a major role in barrier-free reactions, especially at low temperatures, which makes them a very important aspect for the present study. Specifically, we incorporate electrostatic ($R^{-4}$ and $R^{-5}$) and dispersion energy terms ($R^{-6}$, $R^{-8}$, and $R^{-10}$) suitably modeled from accurately calculated dipolar and quadrupolar moments, as well as polarizabilities.

After the long-range interactions are modeled, the ab initio energies are fitted to the $V_{EHF}^{(n)}$ terms such that the final PES ($V(\mathbf{R})$) accurately represents the CBS energies. For this, a polynomial times a range decaying function is employed. It should be noted that the ab initio calculations show the existence of a conical intersection (CI) between $1^2A'$ and $2^2A'$ electronic states for linear SiS-H configurations, which causes a cusp in the ground adiabatic

PES that is being modeled. Because such cusps cannot be modeled with smooth functions, we improve their description by using a method that allows for a realistic description of CIs on adiabatic PESs (Galvão et al. 2015, 2016); for further details, the reader is referred to previous work (Galvão et al. 2015, 2016; Rocha & Varandas 2016; Gonçalves et al. 2018), where the method has been successfully applied in a variety of triatomic systems.

## 2.3. Rate Coefficients

To compute the rate coefficients for the title reaction, we employed the quasiclassical trajectory (QTC; Hase et al. 1996; Peslherbe et al. 1999) method using a version of the Venus96 (Hase et al. 1996) code costumized by the Coimbra Theoretical and Computational Group. For each temperature, a batch of $10^4$ trajectories were integrated, starting with reactants separated by 25 Å and with the initial conditions of every trajectory sampled such that the whole batch mimics the correct energy distributions at the fixed temperature. In this spirit, the atom-diatom translational energy follows the Maxwell-Boltzmann distribution, the rovibrational quantum states of SH follow Boltzmann's distribution, and the impact parameter (b) is uniformly distributed from 0 to $b_{max}$, where $b_{max}$ was independently tuned for each temperature, ranging from 20 Å for 10 K to 9 Å for 1000 K. A time step of 0.2 fs was used to ensure a satisfactory conservation of the total energy in all trajectories, with an average error of $5 \times 10^{-4}$ kJ mol$^{-1}$, with the worst trajectory among all 110,000 deviating by only 1.8 kJ mol$^{-1}$.

The Monte Carlo integrated rate coefficient for the reaction assumes the form

$$k(T) = g_e(T) \left( \frac{8k_B T}{\pi \mu} \right)^{1/2} \pi b_{\max}^2 \frac{N^r}{N}, \quad (4)$$

where $k_B$ is the Boltzmann constant, $\mu$ is the reduced mass of the reactants, and $N^r/N$ is the fraction of reactive trajectories. In turn, $g_e(T)$ is the electronic degeneracy factor, which for the $Si(^3P) + SH(X^2\Pi)$ collisions assumes the form

$$g_e(T) = 2 \times [2 + 2 \exp(-542.36/T)]^{-1} \times [1 + 3 \exp(-110.95/T) + 5 \exp(-321.07/T)]^{-1}. \quad (5)$$

To estimate an error margin for the rate coefficients caused by the trajectory sampling, we plot $k(T)$ with error bars corresponding to the statistical 68% confidence interval, which is given by $\Delta k = k \left( \frac{N - N^r}{NN^r} \right)^{1/2}$.

## 2.4. Astrochemical Modeling

To study the impact of the new proposed reaction for interstellar chemistry, we have used the Nautilus gas-grain code (Ruaud et al. 2016) and simulated a shock similar to L1157-B1 (Podio et al. 2017). Time-dependent models with and without the newly proposed reaction were performed to assess its impact on the SiS abundances.

Nautilus is a full gas-grain astrochemical model that can compute the ice and gas chemical composition under interstellar conditions. The gas-phase chemistry is based on the kida.uva.2014 network (Wakelam et al. 2015). In addition, for all appropriate gas-phase processes, the model computes the
<table><caption>**Table 1**

Elemental Abundances Used for the Model (Jenkins 2009; Ruaud et al. 2018)</caption><thead><tr><th>Species</th><th>Abundances</th></tr></thead><tbody><tr><td>H</td><td>1.0</td></tr><tr><td>He</td><td>0.09</td></tr><tr><td>N</td><td>6.2 × 10<sup>-5</sup></td></tr><tr><td>O</td><td>3.3 × 10<sup>-4</sup></td></tr><tr><td>C<sup>+</sup></td><td>1.8 × 10<sup>-4</sup></td></tr><tr><td>S<sup>+</sup></td><td>1.5 × 10<sup>-5</sup></td></tr><tr><td>Si<sup>+</sup></td><td>1.8 × 10<sup>-6</sup></td></tr><tr><td>Fe<sup>+</sup></td><td>2.0 × 10<sup>-7</sup></td></tr><tr><td>Na<sup>+</sup></td><td>2.3 × 10<sup>-7</sup></td></tr><tr><td>Mg<sup>+</sup></td><td>2.3 × 10<sup>-6</sup></td></tr><tr><td>P<sup>+</sup></td><td>7.8 × 10<sup>-8</sup></td></tr><tr><td>Cl<sup>+</sup></td><td>3.4 × 10<sup>-8</sup></td></tr><tr><td>F</td><td>1.8 × 10<sup>-8</sup></td></tr></tbody></table>

adsorption and desorption of the species onto and from the surfaces, as well as the reactions on the surfaces. The model is a three-phase model, which means that, in addition to the gas phase, the species on top of the grains are differentiated as “surface” or “bulk.” The surface represents the most external two layers of molecules, while the rest below is considered as bulk. Only species from the surface can desorb (through thermal desorption, chemical desorption, cosmic-ray-induced desorption, and photodesorption). Diffusion of the species is faster than in the bulk, but both are allowed. Photodissociations (due to direct or indirect photons) can occur both on the surface and in the bulk. Only “typical” silicate grains (of 0.1 µm in radius) are considered. All chemical parameters are the same as in Ruaud et al. (2016). We refer to this paper for a full description of the code.

Since the SiS molecule has been observed in shock regions, we have used a three-step method. As a first step, we simulated a cold molecular cloud chemical composition by computing the chemical evolution for $10^5$ yr, starting from an initial atomic abundance except for hydrogen, which is entirely molecular. The list of elemental abundances is given in Table 1, and they are the atomic abundances observed in the ζ Oph diffuse medium from Jenkins (2009; see also Ruaud et al. 2018, for discussion). The gas and dust temperatures were set to 10 K and the total H density to $2 \times 10^4$ cm$^{-3}$. The integration time is reasonable for the molecular cloud “chemical age” derived from the comparing observations and modeling (see, e.g., Agúndez & Wakelam 2013). The value we have used for the silicon elemental abundance is rather high with respect to the abundance very often considered in the “low-metal” case (Graedel et al. 1982), but it is the one observed in the diffuse medium and that can be reasonably assumed to be available for volatile interstellar chemistry. The sulfur abundance is also high, but the recent chemical models have been able to reproduce cold core observations with such high elemental abundance of sulfur (Vidal et al. 2017). The second step of our model simulates the shock itself. Starting from the chemical composition obtained at the end of step 1, we then set the temperature to 100 K, keeping the density the same. We run the model for $10^3$ yr, with the final ice and gas compositions then used as an initial composition for the third step, during which the gas and dust cool down to 30 K and the density then increases to $10^7$ cm$^{-3}$. These physical conditions are in agreement with those observed in sources like L1157-B1,

which is a protostellar shock driven by the low-mass Class 0 source L1157 mm.

Several S- and Si-bearing molecules have been cataloged and analyzed toward the L1157-B1 shock region (e.g., Bachiller & Pérez Gutiérrez 1997; Lefloch et al. 2012; Podio et al. 2017; Holdship et al. 2019; Feng et al. 2020; Spezzano et al. 2020); regarding SiS, Podio et al. (2017) reported its detection based on single-dish and interferometric observations toward the shock. They observed various lines of SiS and used the rotational diagram method to obtain an estimate of the gas temperature and SiS average column densities over the cavity of L1157-B1. Their results suggest a gas temperature around 24 K and an abundance value of $[\text{SiS}/\text{H}_2] \sim 2 \times 10^{-8}$ for the head of the cavity. In comparison with SiO, they also found that SiO/SiS ratios might vary in the shock, estimating values of, e.g., 25 and over 180 in the head of the cavity B1c and at the jet impact region B1a, respectively.

We did two sets of models. In the first we have used the original chemical networks, whereas in the second we have included the new reaction $Si + HS \rightarrow SiS + H$.

## 3. Results

### 3.1. Features of the Modeled PES

The final fit employed in the construction of the PES used 212 linear parameters and 56 nonlinear ones, which mimicked all 2267 energy points with an overall rms deviation below $1 \text{ kcal mol}^{-1}$. As first pointed out by Rosi et al. (2018), this triatomic system has two covalently bound potential energy minima. The deepest one corresponds to the thiosilaformyl (HSiS) radical, which can isomerize via a hydrogen atom migration to the SiSH isomer. Both structures present a nonlinear geometry. The isomerization transition state (TS1) lies lower in energy than any dissociation channel. In collisions between Si+SH, the highest-energy isomer (SiSH) is formed in a direct manner, while the other one can be accessed after the isomerization step. The S+SiH asymptote lies $55 \text{ kJ mol}^{-1}$ higher in energy than Si+SH; hence, collisions leading to this channel unlikely contribute to the abundance of silylidyne (SiH) in the ISM, except at very high temperature regions.

Figure 1 summarizes the main attributes of the PES. The zero-point energy (ZPE) corrected values for the stationary structures obtained by the present PES are given in bold. For comparison, also shown are the values obtained by Rosi et al. (2019) at the CCSD(T)/aug-cc-pV(T+d)Z level and those of de Castro et al. (2020). As seen, our results agree well with the CCSD(T) ones. Since our energies are at the CBS limit, and the HSiS minimum shows a relatively high multireference character ($T1 = 0.035$; Goettl et al. 2021), we trust our results as the most accurate. The agreement with the results of de Castro et al. (2020) is acceptable only at the SiH+S asymptotic channel.

Figure 2 provides an overview of the abstraction route toward the products as $Si + SH \rightarrow SiSH \rightarrow SiS + H$. It shows a contour plot for a reaction occurring at a fixed SiSH angle of 97.3, where it can be seen that the interaction between atomic silicon and the mercapto radical leads to the SiSH minimum without a potential energy barrier. From this minimum to SiS + H dissociation, there is a small exit transition state (TS2). This barrier lies much lower in energy than the Si+SH reactants and hence does not prevent reaction to occur at low temperatures, as will be shown later.
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Schematic representation of the PES for the reaction, including ZPE correction. The bold numbers refer to the values obtained by the fit performed in this work, while those of Rosi et al. (2019) and de Castro et al. (2020) are in italics and within parentheses, respectively. The CI is also presented and connected with gray lines, to emphasize that this is not a minimum energy path. Atoms are colored as follows: silicon (purple), sulfur (yellow), hydrogen (white).</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Contour plot for an Si–S–H bond stretching at fixed valence angle of 97°.3 (showing the SiSH minimum). The zero of energy is set at the reactants limit (Si+SH). Solid contours correspond to energies above this limit, while dashed–dotted ones are lower. The contours are incremented by ±20 kJ mol⁻¹.</figcaption></figure>

For an improved visualization of the two isomers of the triatomic molecule, Figure 3 shows a contour plot for a hydrogen atom moving around the SiS diatomic with its bond length fixed at equilibrium position. The SiSH isomer that can be obtained directly from reactants is seen on the right-hand side of this plot, and TS2 can also be seen above the SiS+H products, as mentioned earlier. The CI that was modeled by our fitted PES can also be seen in this figure (see also Figure 1). It occurs when going from the SiSH potential energy basin toward SiS+H but for linear arrangements, which is not the minimum energy path for reaction. The cusp caused by this CI

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Contour plots for H moving around SiS. The zero of energy is set at the products limit (H+SiS). Solid contours correspond to energies above this limit, while dashed ones are lower. The contours are incremented by ±8.5 kJ mol⁻¹.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Calculated rate coefficients for the Si+SH reaction. The error bars correspond to a statistical 68% confidence interval, and for most temperatures they are smaller than the circles. The results of Rosi et al. (2019) using a capture model are shown in the inset for comparison.</figcaption></figure>

is seen as the sharp contours on the lower part of the right-hand side of Figure 3.

Figure 3 also shows the possible H atom migration between the two isomers (through TS1). Note that the isomerizing mechanism Si + SH → SiSH → HSiS → SiS + H does not show an exit transition state. This means that the inverse reaction SiS + H → HSiS can occur barrierlessly. Although the SiH+S or Si+SH channels cannot be accessed from SiS+H collisions, the barrierless path to the HSiS minimum should have important consequences for the SiS collisional vibrational energy transfer. This subject will be investigated in the future.

## 3.2. Quasi-classical Calculations

After integrating 10⁴ trajectories for each temperature and using them to calculate the rate coefficients, the final results are displayed in Figure 4. As shown, the calculated rate coefficients below 200 K increase drastically with decreasing temperature, as expected for a barrierless reaction. In fact, they are similar to the ones reported by Rivero-Santamaría et al. (2014) for the analogous Si + OH → SiO + H reaction. Moreover, they can
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Abundances as a function of time obtained from the step corresponding to the cold core phase. The terms in parentheses, “g” and “s,” stand for abundances computed in the gaseous and solid phase, respectively.</figcaption></figure>

be fitted to the modified Arrhenius expression:

$$k(T) = \alpha \left(\frac{T}{300}\right)^{\beta} \exp(-\gamma/T), \quad (6)$$

where $\alpha = 0.770 \times 10^{-10} \text{ cm}^3 \text{ s}^{-1}$, $\beta = -0.756$, and $\gamma = 9.873$ K.

We should point out at this stage that previous attempts to estimate the rate coefficients for the Si+SH reaction have been proposed (Rosi et al. 2019) within a simple capture model, achieving a high-temperature limit of $1 \times 10^{-9} \text{ cm}^3 \text{ s}^{-1}$. This is higher than our results and is also shown for comparison in Figure 4. The capture model does not show the same temperature dependence, with rate coefficients decreasing for lower temperatures. Another QCT study has been performed (de Castro et al. 2020), but the authors predict that the reaction shows a barrier-like behavior, thus contradicting other theoretical work (Rosi et al. 2018, 2019; Paiva et al. 2020) and the possibility of using a capture-type model since this implies a barrier-free reaction. Not surprisingly, their thermal rate coefficients are much smaller than ours and nearly four orders of magnitude smaller than that of the analogous Si+OH reaction. In the absence of evidence that justifies the existence of any barrier, we trust the present work as realistic and appropriate for astrochemical databases.

### 3.3. Astrochemical Modeling Results

During the first step (cold molecular chemistry), given in Figure 5, the silicon, initially in its atomic ionized form, is mostly neutralized, and the neutral Si in the gas phase is still very high at $10^5$ yr. Part of the silicon depletes on the dust surfaces, forming SiO and mostly SiH$_4$. Sulfur is also first neutralized and then depletes, forming almost equally HS and H$_2$S on the grains (see also Vidal et al. 2017).

The results of the second fast step can be seen in Figure 6, where the only effect of temperature is to desorb the Si and S molecules in the gas phase. The Si gas-phase atomic abundance already high at the beginning of this step stays as high as $\sim 2 \times 10^{-6}$ (Figure 6). HS being desorbed from the grains produces an abundance to more than $1 \times 10^{-7}$ (see Figure 6).

Starting the last step (step 3) of the model with such high Si and SH abundances in the gas phase, and with the new neutral

<figure><img src="image_7.png" /><figcaption>**Figure 6.** Abundances as a function of time obtained from the second and fast step. As indicated in Figure 5, the terms “g” and “s” stand for abundances computed in the gaseous and solid phase, respectively.</figcaption></figure>

neutral reaction included in the model, strongly enhanced the SiS abundance predicted by the model, as can be seen in Figure 7. Such an increment is, for instance, $\sim 28$ times at the early ages ($<10^2$ yr) of the shock model, with abundance values (relative to atomic hydrogen) rising from $\sim 2.5 \times 10^{-10}$ to $7 \times 10^{-9}$. After this time, the silicon depletes again on the grains. The new abundance relative to molecular hydrogen is therefore $1.4 \times 10^{-8}$, which matches very well the observed abundance of [SiS/H$_2$] $\sim 2 \times 10^{-8}$ (Podio et al. 2017).

To assess the dependency of the results on changes in the rate coefficients, we also employ the values obtained by Rosi et al. (2019) to rerun the same astrochemical model. These rate coefficients are about one order of magnitude higher than ours and approximate our values only below 25 K. As a result, the time evolution of the abundances stays essentially the same as those previously presented in Figures 5, 6, and 7. However, the final abundance of SiS relative to molecular hydrogen is approximately 2.5 times higher than when using our rate coefficients, attaining a value of about $3.5 \times 10^{-8}$. In summary, using the rate coefficients of Rosi et al. (2019) yields an abundance higher than the observed one, while deviating from it more than when using our calculated values. Nevertheless, the predicted abundances of SiS compare reasonably well with the observed values within a given observational uncertainty.

### 4. Concluding Remarks

In this work we have performed accurate multireference calculations based on the MRCI(Q) method, which were subsequently extrapolated to the complete basis set limit. The calculated energies cover the whole configurational space of the HSiS triatomic and were accurately modeled via the DMBE method. This approach further warrants an accurate representation of the long-range interactions, which are of crucial importance for predicting rate coefficients at low temperatures. The fitted PES was then employed in QCT calculations of the rate coefficients for the Si+SH $\rightarrow$ SiS+H reaction over a wide range of temperatures.

The relevance of the title reaction for SiS formation was tested with an astrochemical gas-grain model of the L1157-B1 shock region. The results show that the inclusion of this neutral-neutral
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Abundances as a function of time obtained from the step corresponding to the shock model: comparative results (a) without and (b) with the reaction Si + HS → SiS + H. Similar to Figure 5, the terms “g” and “s” stand for abundances computed in the gaseous and solid phase, respectively.</figcaption></figure>

process enhances by 28-fold the SiS gas-phase abundance relative to H₂, thence from $5 \times 10^{-10}$ to $1.4 \times 10^{-8}$. It is therefore sufficient to explain the observed abundance of $[SiS/H_2] \sim 2 \times 10^{-8}$. This reaction therefore helps elucidate the mechanisms of formation and abundance of silicon monosulfide in the ISM.

The new rate coefficients can be used in other astrochemical models and can also be added to existing databases. Additionally, the novel PES can be used to predict the rate coefficients for other molecular collisions such as S+SiH and SiS+H, which are also relevant for astrochemistry. Additional details regarding the PES, comparison with previous work, treatment of the ZPE leakage (Varandas 2007a, and references therein) in QCT trajectories, and state-specific rate coefficients will be given in a future publication.

This work has been financed in part by the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior—Brasil (CAPES) —Finance Code 001, Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq), grant 305469/2018-5 (B.R.L. G.) and grant 150465/2019-0 (E.M.), and Fundação de Amparo à Pesquisa do estado de Minas Gerais (FAPEMIG). The authors

also acknowledge the CNRS program “Physique et Chimie du Milieu Interstellaire” (PCMI) co-funded by the Centre National d’Etudes Spatiales (CNES). B.R.L.G. is also thankful to Rede Mineira de Química (RQ-MG), while V.C.M. acknowledges the support of Edital 2015 do Programa institucional de Fundo de Apoio à Pesquisa da Universidade Federal do Espírito Santo and also the Edital universal FAPES 2018. This work has also the support of Foundation for Science and Technology, Portugal, and Coimbra Chemistry Centre, Portugal, through the project UID/QUI/00313/2019. A.J.C.V. thanks also China’s Shandong Province “Double-Hundred Talent Plan” (2018).

## ORCID iDs

V. C. Mota https://orcid.org/0000-0001-8368-0803

A. J. C. Varandas https://orcid.org/0000-0003-1501-3317

E. Mendoza https://orcid.org/0000-0001-9381-7826

V. Wakelam https://orcid.org/0000-0001-9676-2605

B. R. L. Galvão https://orcid.org/0000-0002-4184-2437

## References

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

Bachiller, R., & Pérez Gutiérrez, M. 1997, ApJL, 487, L93

Danilovich, T., Ramstedt, S., Gobrecht, D., et al. 2018, A&A, 617, A132

de Castro, D. G., Ballester, M. Y., & Song, Y.-Z. 2020, J. Phys. B: At. Mol. Opt. Phys., 53, 175203

Dickinson, D., & Kuiper, E. R. 1981, ApJ, 247, 112

Doddipatla, S., He, C., Goettl, S., et al. 2021, SciA, 7, eabg7003

Dunning, T. H. 1989, JChPh, 90, 1007

Feller, D. 1993, JChPh, 98, 7059

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

Galvão, B. R. L., Mota, V. C., & Varandas, A. J. C. 2015, JPCA, 119, 1415

Galvão, B. R. L., Mota, V. C., & Varandas, A. J. C. 2016, CPL, 660, 55

Goettl, S. J., Doddipatla, S., Yang, Z., et al. 2021, PCCP, 23, 3647

Gonçalves, C. E. M., Galvão, B. R. L., Mota, V. C., Braga, J. P., & Varandas, A. J. C. 2018, JPCA, 122, 4198

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

Gusdorf, A., des Forêts, G. P., Cabrit, S., & Flower, D. R. 2008, A&A, 490, 695

Hase, W. L., Duchovic, R. J., Hu, X., et al. 1996, QCPE Bull, 16, 43

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

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

Kendall, R. A., Dunning, T. H., Jr., & Harrison, R. J. 1992, JChPh, 96, 6796

Lefloch, B., Cabrit, S., Busquet, G., et al. 2012, ApJL, 757, L25

Lefloch, B., Bachiller, R., Ceccarelli, C., et al. 2018, MNRAS, 477, 4792

Massalkhi, S., Agúndez, M., & Cernicharo, J. 2019, A&A, 628, A62

Morris, M., Gilmore, W., Palmer, P., Turner, B., & Zuckerman, B. 1975, ApJ, 199, L47

Paiva, M. A. M., Lefloch, B., & Galvão, B. 2018, MNRAS, 481, 1858

Paiva, M. A. M., Lefloch, B., & Galvão, B. R. L. 2020, MNRAS, 493, 299

Peslherbe, G. H., Wang, H., & Hase, W. L. 1999, AdChP, 105, 171

Podio, L., Codella, C., Lefloch, B., et al. 2017, MNRAS, 470, L16

Prieto, L. V., Cernicharo, J., Quintana-Lacaci, G., et al. 2015, ApJ, 805, L13

Rivero-Santamaría, A., Dayou, F., Rubayo-Soneira, J., & Monnerville, M. 2014, CPL, 610, 335

Rivero Santamaría, A., Dayou, F., Rubayo-Soneira, J., & Monnerville, M. 2017, JPCA, 121, 1675

Rocha, C. M. R., & Varandas, A. J. C. 2016, JChPh, 144, 064309

Rosi, M., Mancini, L., Skouteris, D., et al. 2018, CPL, 695, 87

Rosi, M., Skouteris, D., Balucani, N., et al. 2019, in Int. Conf. on Computational Science and Its Applications (Berlin: Springer), 306

Ruaud, M., Wakelam, V., Gratier, P., & Bonnell, I. A. 2018, A&A, 611, A96

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

Schöier, F. L., Bast, J., Olofsson, H., & Lindqvist, M. 2007, A&A, 473, 871

Spezzano, S., Codella, C., Podio, L., et al. 2020, A&A, 640, A74

Varandas, A. J. C. 1985, J. Mol. Struct. Theochem., 120, 401

Varandas, A. J. C. 1988, Adv. Chem. Phys., 74, 255

Varandas, A. J. C. 1992, CPL, 194, 333

Varandas, A. J. C. 2000, in Lecture Notes in Chemistry, ed. A. Laganà & A. Riganelli, Vol. 75 (Berlin: Springer) p.33
Varandas, A. J. C. 2007a, CPL, 439, 386

Varandas, A. J. C. 2007b, JChPh, 126, 244105

Varandas, A. J. C. 2018, ARPC, 69, 177

Varandas, A. J. C. 2019, PCCP, 21, 8022

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

Wakelam, V., Loison, J.-C., Herbst, E., et al. 2015, ApJS, 217, 20

Werner, H.-J., & Knowles, P. J. 1985, JChPh, 82, 5053

Werner, H.-J., & Knowles, P. J. 1988, JChPh, 89, 5803

Werner, H.-J., Knowles, P. J., Knizia, G., et al. 2015, WIRES Comput. Mol. Sci., 2, 242

Zanchet, A., Roncero, O., Agúndez, M., & Cernicharo, J. 2018, ApJ, 862, 38

Zhukovska, S., & Gail, H.-P. 2008, A&A, 486, 229