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

The product channel SiS(¹Σ⁺)+H(² S) correlates with the PES of a single electronic state of the HSiS molecule, namely, ${}^{2}A^{\prime}$. 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 ${}^{2}A^{\prime}$ 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)

$$\begin{eqnarray}&&{E}_{X}^{\mathrm{HF}}={E}_{\infty }^{\mathrm{HF}}+B\exp (-\beta X),\end{eqnarray} \tag{ 1 }$$

with B and ${E}_{\infty }^{\mathrm{HF}}$ being parameters to be determined from a fit to a pair of CASSCF/aug-cc-pVXZ calculations. Following Varandas (2019), β 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

$$\begin{eqnarray}&&{E}_{X}^{\mathrm{dc}}={E}_{\infty }^{\mathrm{dc}}+\displaystyle \frac{{A}_{3}}{{\left(X-3/8\right)}^{3}}+\displaystyle \frac{{A}_{5}^{(0)}+{{cA}}_{3}^{n}}{{\left(X-3/8\right)}^{5}},\end{eqnarray} \tag{ 2 }$$

which provides ${E}_{\infty }^{\mathrm{dc}}$ and A₃. 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).

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

$$\begin{eqnarray}&&V({\boldsymbol{R}})=\sum _{i=1}^{3}{V}_{i}^{(2)}({R}_{i})+{V}^{(3)}({R}_{\mathrm{SiS}},{R}_{\mathrm{SiH}},{R}_{\mathrm{SH}}),\end{eqnarray} \tag{ 3 }$$

where ${V}_{i}^{(2)}({R}_{i})$ are two-body energy terms and V⁽³⁾ a three-body one. Both are further partitioned to describe separately short-range (${V}_{\mathrm{EHF}}^{(n)}$) and long-range (${V}_{\mathrm{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⁻⁴ and R⁻⁵) and dispersion energy terms (R⁻⁶, R⁻⁸, and R⁻¹⁰) 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}_{\mathrm{EHF}}^{(n)}$ terms such that the final PES (V( 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}^{2}A^{\prime}$ and ${2}^{2}A^{\prime}$ 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.

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⁴ 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 × 10⁻⁴ kJ mol⁻¹, with the worst trajectory among all 110,000 deviating by only 1.8 kJ mol⁻¹.

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

$$\begin{eqnarray}&&k(T)={g}_{e}(T){\left(\displaystyle \frac{8{k}_{B}T}{\pi \mu }\right)}^{1/2}\pi {b}_{\max }^{2}\displaystyle \frac{{N}^{r}}{N},\end{eqnarray} \tag{ 4 }$$

where k_B is the Boltzmann constant, μ 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(³ P) + SH(X²Π) collisions assumes the form

$$\begin{eqnarray}&&\begin{array}{l}{g}_{e}(T)=2\times {\left[2+2\exp (-542.36/T)\right]}^{-1}\\ \quad \times {\left[1+3\exp (-110.95/T)+5\exp (-321.07/T)\right]}^{-1}.\end{array}\end{eqnarray} \tag{ 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 ${\rm{\Delta }}k=k{\left(\tfrac{N-{N}^{r}}{{{NN}}^{r}}\right)}^{1/2}$.

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 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⁵ 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 × 10⁴ cm⁻³. 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³ 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⁷ cm⁻³. 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 [SiS/H₂] ∼ 2 × 10⁻⁸ 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 → SiS + H.

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 kcal mol⁻¹. 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 kJ mol⁻¹ 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.

Zoom In Zoom Out Reset image size

Figure 2 provides an overview of the abstraction route toward the products as Si + SH → SiSH → SiS + H. It shows a contour plot for a reaction occurring at a fixed SiSH angle of 973, 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.

Zoom In Zoom Out Reset image size

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 is seen as the sharp contours on the lower part of the right-hand side of Figure 3.

Zoom In Zoom Out Reset image size

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.

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 be fitted to the modified Arrhenius expression:

$$\begin{eqnarray}&&k(T)=\alpha {\left(\displaystyle \frac{T}{300}\right)}^{\beta }\exp (-\gamma /T),\end{eqnarray} \tag{ 6 }$$

where α = 0.770 × 10⁻¹⁰ cm³ s⁻¹, β = −0.756, and γ = 9.873 K.

Zoom In Zoom Out Reset image size

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 × 10⁻⁹ cm³ s⁻¹. 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.

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⁵ yr. Part of the silicon depletes on the dust surfaces, forming SiO and mostly SiH₄. Sulfur is also first neutralized and then depletes, forming almost equally HS and H₂S on the grains (see also Vidal et al. 2017).

Zoom In Zoom Out Reset image size

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 ∼2 × 10⁻⁶ (Figure 6). HS being desorbed from the grains produces an abundance to more than 1 × 10⁻⁷ (see Figure 6).

Zoom In Zoom Out Reset image size

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−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, ∼28 times at the early ages (<10² yr) of the shock model, with abundance values (relative to atomic hydrogen) rising from ∼2.5 × 10⁻¹⁰ to 7 × 10⁻⁹. After this time, the silicon depletes again on the grains. The new abundance relative to molecular hydrogen is therefore 1.4 × 10⁻⁸, which matches very well the observed abundance of [SiS/H₂] ∼ 2 × 10⁻⁸ (Podio et al. 2017).

Zoom In Zoom Out Reset image size

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 × 10⁻⁸. 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.