Formation and Destruction of SiS in Space

The molecule SiS was first observed in the molecular envelope around the carbon-rich AGB star IRC +10216 and in the molecular cloud toward the Galactic center Sgr B2 (Morris et al. 1975). This species is known to be present in circumstellar envelopes around evolved stars, where it is thought to be formed under thermochemical equilibrium conditions at the stellar surface and later ejected into the expanding wind (Bujarrabal et al. 1994; Cernicharo et al. 2000; Schöier et al. 2007; Agúndez et al. 2012; Fonfría et al. 2015; Velilla Prieto et al. 2015). The molecule is observed in both carbon-rich and oxygen-rich envelopes, but SiS maser emission is only observed in the carbon-rich object IRC +10216 (Fonfría Expósito et al. 2006). Apart from the envelopes of evolved stars, SiS has only been found in a handful of star-forming regions with evidence of outflows, like Sgr B2 and Orion KL (Dickinson & Kuiper 1981; Ziurys 1988, 1991; Tercero et al. 2011). Recently, Podio et al. (2017) reported the detection of SiS in L1157-B1, a shocked region associated with an outflow driven by a low-mass protostar.

Despite the fact that SiS is a stable molecule involving two elements with high cosmic abundances, it is much less commonly observed than its oxygen chemical analogue SiO. While SiO is prevalent in almost every outflow affected by shocks, and as such it is considered a good tracer of shock activity (Martin-Pintado et al. 1992), the number of SiS detections in shocked outflows is scarce. There are various obvious questions regarding the different prevalence of SiO and SiS in the interstellar medium, in general, and in shocked outflows, in particular. Is SiS an elusive molecule that is only formed under very particular conditions? Could it be that lines of SiS have not been targeted as often as those of SiO, meaning that there has been an observational bias and that both SiO and SiS are ubiquitous in shocked gas, SiS being merely less abundant than SiO? Answers to these questions may come from more sensitive searches for SiS in regions with intense SiO emission. Also chemical models may hold the clue for the apparent different behavior of SiS and SiO. However, the kinetics of chemical reactions involving SiS is very poorly known, which hampers a good understanding of the chemistry of this molecule in interstellar space.

In this article we study various potentially important processes that can drive the formation and destruction of SiS in interstellar and circumstellar media. In a recent work, Rosi et al. (2018) proposed the formation of this molecule through the reactions SiH + S and SiH + S₂ and carried out electronic structure calculations of the stationary points along the reactive potential energy surfaces (PESs). In this work we study a different chemical route, the formation of SiS via the reaction Si + SO and its destruction through the SiS + O reaction. For this purpose, we calculate the ab initio PESs of the SiOS system in the two lower electronic states and perform quasiclassical trajectory (QCT) calculations of the rate coefficients for several vibrational states of the reactants in the 50–2500 K temperature range. These reaction rate coefficients are then used in a chemical model to evaluate their impact on the formation and destruction of SiS in outflows. The manuscript is organized as follows. First, in Section 2 we describe the calculations carried out. In Section 3 we discuss the impact of the studied processes on the chemistry of SiS in outflows. Finally, in Section 4 we present our conclusions.

2.1. SiOS PESs

The SiOS system correlates with the open shell fragments SO(³Σ⁻) + Si(³P), SiS(¹Σ⁺) + O(³P), and SiO(¹Σ⁺) + S(³P), as shown in Figure 1. All of the atomic fragments are open shell ³P states, and there are several electronic states correlating to them. In the present study we neglect spin–orbit couplings, and only three states in each asymptote will be considered here (correlating to the P states of Si, S, and O atoms, respectively). These three states are classified according to their symmetry by reflection through the plane of the molecule: two of ³A'' symmetry and one of ³A' symmetry. The excited 2³A'' state presents a very high barrier for the reaction, and will not be considered in this work. The minimum energy paths connecting the three rearrangement channels for the 1³A'' and 1³A' states are shown in Figure 1, and are clearly connected by rather deep wells. These paths are obtained from the full-dimension PESs described below.

Zoom In Zoom Out Reset image size

The PESs of the 1³A'' and 1³A' states have been calculated with an internally contracted multireference configuration interaction method with simple and double excitations (icMRCI; Werner et al. 1988), including the Davidson correction (Davidson 1975). Aug-cc-pVTZ basis sets have been considered for all of the atoms, including spdf basis functions, and the calculations have been done with the MOLPRO suite of programs (Werner et al. 2012). In order to get an accurate and homogeneous description of the PESs, the molecular orbitals and reference configurations were determined with a state-average complete active space (SA-CASSCF) method (Werner & Knowles 1985) for the calculation of the first two ³A' and the first two ³A'' electronic states. An active space of 12 electrons in 10 orbitals (12–18a' and 3–5a'') was used in order to include all of the relevant valence orbitals with affordable calculation time. The 1s orbital of oxygen was kept frozen as well as the 1s, 2s, and 2p orbitals of Si and S.

Ab initio icMRCI+Q calculations were performed for 3042 geometries for the 1³A' and 1³A'' states. These points have been fitted using a method based on the reproducing kernel Hilbert space for each adiabatic state independently. This procedure is known to be well adapted for potentials of triatomic systems with complex topologies consisting of several wells and saddle points (Zanchet et al. 2006, 2009). In addition, a fast evaluation algorithm can be implemented when points are sampled on a regular grid (Hollebeek et al. 1997) like in this work, making this kind of analytical potentials very appropriate to perform dynamics calculation.

The two PESs are shown in Figure 2 and are very similar; they present two deep wells associated to SSiO and OSSi isomers, and a shallow van der Waals well in the entrance channel. The saddle point between the van der Waals and the SSiO wells gives rise to a submerged barrier in the entrance channel (see Figure 2). At these submerged barriers, the variation of the potential with the SiS stretching distance becomes narrower than for free SiO(¹Σ⁺). This will have interesting effects in the dynamics described below. A more detailed description of the fits and the PESs will be given elsewhere (A. Zanchet et al. 2018, in preparation).

Zoom In Zoom Out Reset image size

In this work we focus on two reactions. The formation of SiS(¹Σ⁺) from Si(³P) + SO(³Σ⁻), which competes with the more exothermic channel leading to SiO(¹Σ⁺) + S(³P),

$$\begin{eqnarray}\mathrm{Si}{(}^{3}P)+\mathrm{SO}{(}^{3}{{\rm{\Sigma }}}^{-}) & \to & \mathrm{SiS}{(}^{1}{{\rm{\Sigma }}}^{+})+{\rm{O}}{(}^{3}P)+1\,\mathrm{eV}\\ & \to & \mathrm{SiO}{(}^{1}{{\rm{\Sigma }}}^{+})+{\rm{S}}{(}^{3}P)+2.9\,\mathrm{eV},\end{eqnarray} \tag{ 1 }$$

and the destruction of SiS(¹Σ⁺) by the reaction

$$\begin{eqnarray}&&\mathrm{SiS}{(}^{1}{{\rm{\Sigma }}}^{+})+{\rm{O}}{(}^{3}P)\to \mathrm{SiO}{(}^{1}{{\rm{\Sigma }}}^{+})+{\rm{S}}{(}^{3}P)+1.9\,\mathrm{eV}.\end{eqnarray} \tag{ 2 }$$

These reactions are exothermic and are discussed below. They correspond to the first triplet states on the three arrangement channels. It is worth mentioning that the Si(³P)+SO(³Σ⁻) correlates not only with these triplet states, but also with one singlet and five quintuplet states. The quintuplets correlate with very high states in the SiS and SiO asymptotes and their contribution to the reactions studied here can be neglected. The singlet correlates with the other two arrangement channels as

$$\begin{eqnarray}\mathrm{Si}{(}^{3}P)+\mathrm{SO}{(}^{3}{{\rm{\Sigma }}}^{-}) & \to & \mathrm{SiS}{(}^{1}{{\rm{\Sigma }}}^{+})+{\rm{O}}{(}^{1}D)-1\,\mathrm{eV}\\ & \to & \mathrm{SiO}{(}^{1}{{\rm{\Sigma }}}^{+})+{\rm{S}}{(}^{1}D)+1.75\,\mathrm{eV}.\end{eqnarray} \tag{ 3 }$$

SiS(¹Σ⁺) + O(${}^{1}D$) is endothermic by ≈1 eV, and it is closed for the energies considered. The SiO(¹Σ⁺) + S(${}^{1}D$) channel is open and can be formed, but since the statistical weight of the singlet is one-third with respect to that of the triplets, its contribution is expected to be small, and will not be considered in the present work.

2.2. SiS Destruction: SiS + O $\to$ SiO+ S

The use of quantum methods in systems with such heavy atoms and deep wells is still a challenge nowadays. However, classical mechanics are expected to give rather accurate results, and for this reason in this work we use the QCT method (Karplus et al. 1965) as implemented in the miQCT code (Dorta-Urra et al. 2015; Zanchet et al. 2016). For each vibrational state of SiS, $5\times {10}^{5}$ trajectories are run changing the initial conditions, consistent with a Boltzmann distribution of translation and the rotation energy at a given temperature, T. The initial distance between reactants is set to 50 au, with a maximum impact parameter of 35 au The trajectories are stopped when any distance becomes larger than 55 au A conservation of energy better than 0.01 meV is imposed. From these calculations, the vibrational state-selected rate coefficients, K_v, are calculated and are shown in Figure 3. The rate for each electronic state, ³A' and ³A'', has been multiplied by one-third to account for the orbital angular momentum degeneracy, while the spin degeneracy is absent because all of the triplet states lead to the same reactivity.

Zoom In Zoom Out Reset image size

The rate coefficients for the 1³A'' state are 3-4 times larger than those obtained for the 1³A'. This is due to the larger cone of acceptance (i.e., the amplitude of the configuration space between reactants and products) of the 1³A'' state simply because the submerged barrier in the entrance channel is lower in energy. In all cases, the rate coefficient increases with increasing translational temperature, also as a consequence of the increase of cone of acceptance. It is interesting to note that the rates coefficient decreases with increasing the initial vibrational excitation of SiS for the two electronic states. This behavior is due to the narrowing of the PES as a function of SiS internuclear distance in the entrance channel. As can be seen in Figure 2, it occurs for SiO distances of ≈5 au, prior to reach the van der Waals well. Finally, in the case of the 1³A'' state the rate coefficients increase for $T\lt 300$ K, which is interpreted as a manifestation of a complex forming mechanism occurring at these low temperatures due to the van der Waals well in the reactant channel. A similar rise is expected for the ${1}^{1}A^{\prime}$ state but at lower temperatures. A more careful description of the reaction dynamics is now in progress (A. Zanchet et al. 2018, in preparation).

2.3. SiS Formation: Si + SO $\to$ SiS + O

The formation of SiS(¹Σ⁺) has also been analysed by performing a QCT study of the reactive collision between Si and SO. The method and analysis are similar to those discussed in Section 2.2. The difference is that in this case there are two different products: SiS and SiO. The rate coefficients calculated for each of these two products are shown in Figure 4. This asymptote correlates with one singlet, three triplets, and five quintuplets, and the three P states of Si, making a total of 27 states. Neglecting the contributions of the singlet and quintuplets, and considering the triple degeneracy of triplets, the total electronic partition function is 3/27 for ³A' and ³A'', included in Figure 4. The formation rate of SiO is always larger than that of SiS. The SiO/SiS ratio between these two rate coefficients changes with temperature, being about 7 at 500 K and ≈5 at 2000 K. This reduction of the ratio with increasing temperature can be explained by statistical arguments: the number of accessible states for SiS increases faster than for SiO because SiS vibrational frequencies and rotational constants are lower than those of SiO.

Zoom In Zoom Out Reset image size

The total reaction rate coefficients obtained for the formation and destruction rates (summed over electronic and vibrational states) have been fitted to the usual expression

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

and the α, β, and γ parameters are listed in Table 1, together with those corresponding to other reactions used in the astrochemical model described below.

Table 1.  Selected Reactions and Associated Rate Coefficients

Reaction	α	β	γ	References
	(cm3 s−1)		(K)
SiS + O $\to$ SiO + S	$9.53\times {10}^{-11}$	0.29	−32	(1)
Si + SO $\to$ SiO + S	$1.53\times {10}^{-10}$	−0.11	32	(1)
Si + SO $\to$ SiS + O	$1.77\times {10}^{-11}$	0.16	−20	(1)
Si + SO2 $\to$ SiO + SO	$1.00\times {10}^{-10}$	0.00	0	(2)
Si + SO2 $\to$ SiS + O2	$1.00\times {10}^{-10}$	0.00	0	(2)
SiH + S $\to$ SiS + H	$1.00\times {10}^{-10}$	0.00	0	(3)
SiH + S2 $\to$ SiS + SH	$1.00\times {10}^{-10}$	0.00	0	(3)

Note. The rate coefficient is given by α (T/300)${}^{\beta }$ $\exp (-\gamma /T)$.

References. (1) This study. The temperature range of validity is 50–2500 K for SiS + O and 100–2500 K for Si + SO. (2) Assumed value for exothermic and spin-allowed reactions. (3) Assumed value for barrierless reaction studied by Rosi et al. (2018).

Download table as:  ASCIITypeset image

The kinetics of chemical reactions involving SiS is very poorly known and the processes studied here are potentially important in the formation and destruction of SiS in space. Because we have quantitatively evaluated the kinetics of the studied processes, we can draw some conclusions about their role in the chemistry of SiS in space. As already mentioned in Section 1, SiS has been detected in envelopes around evolved stars and in outflows driven by protostars.

In envelopes around evolved stars, SiS is thought to be formed under thermochemical equilibrium at the high-density and high-temperature conditions of the stellar photosphere. SiS is then injected into the expanding envelope, where it may vary its abundance due to different types of processes such as shocks induced by the periodic pulsation of the star, interaction with dust grains, and photoprocesses driven by the UV radiation emitted by nearby hot stars. In this sense, the reaction of the destruction with O atoms studied here may play an important role in regulating the SiS abundance. In circumstellar envelopes, oxygen atoms can be released by the photodissociation of CO, SiO, and H₂O. Although their relative abundances vary depending on the oxygen- or carbon-rich character, these molecules are the main oxygen reservoirs. Photodissociation of these molecules does certainly takes place in the outer regions of the envelope, which are no longer shielded against UV photons by the circumstellar dust, but it may also take place in the inner regions of clumpy envelopes (Agúndez et al. 2010). Since the reaction between SiS and O is rapid, it may provide an efficient way to destroy SiS and recycle the trapped silicon into SiO. In this line, it is worth noting that interferometric maps of SiO and SiS in the C-star envelope IRC +10216 indicate that SiS emission disappears at shorter distances from the star than SiO (Velilla Prieto et al. 2015). One possible explanation could be that at moderately short distances from the star, the amount of O atoms is high enough to efficiently destroy SiS.

In protostellar outflows the situation is different as SiS probably forms from the silicon released during dust grain disruption. Here we aim to evaluate whether in these environments the reaction between atomic silicon and SO can be an efficient source of SiS and whether the reaction between SiS and atomic oxygen can be an important destruction process of SiS, which could potentially explain the paucity of SiS detections in outflows. Protostellar outflows lead to the generation of shocks between the ejected material and the surrounding quiescent envelope. During the shock, density and temperature are drastically enhanced and, depending on the shock strength, disruption of dust grains and dissociation of molecules can happen. After the passage of the shock, density and temperature relax and molecule formation takes place. Dedicated models have been constructed to treat this process in detail (see, e.g., Schilke et al. 1997; Gusdorf et al. 2008), although they usually focus on SiO and not SiS.

Here we have carried out a simple time-dependent chemical model to explore the influence of the SiS + O and Si + SO reactions on the chemistry of SiS in a post-shocked gas. These calculations do not aim to model a particular shocked outflow, but to evaluate how SiS abundance behaves during molecule formation after the passage of a shock. For this simple model, we consider that density and temperature remain constant with time. We adopted a density of H nuclei of 10⁸ cm⁻³, which makes molecule formation take place within a few hundreds of years. The main effect of the density is on the overall chemical timescale, which is shortened in denser gas and lengthened at lower densities. We consider three different gas kinetic temperatures: 100, 300, and 500 K. We assume that right after the shock passage, all of the molecules have been dissociated so that the initial composition of the post-shocked gas consists of neutral atoms with solar elemental abundances.

We use a large chemical network, with ∼500 gas-phase species composed of the elements H, C, O, N, S, Si, and P and linked by ∼8000 reactions, previously used to study the chemistry of IRC +10216 (Agúndez et al. 2017). The rate coefficients of most reactions involving SiS as a reactant or product in this network are mere guesses based on Si/C and S/O chemical analogies (e.g., Willacy & Cherchneff 1998).

We have checked various potentially important SiS-forming reactions involving atomic silicon as the reactant. Excluding endothermic and spin-forbidden reactions, we identify that Si + SO and Si + SO₂ are probably the most efficient ones. As discussed in the previous section, the reaction between Si and SO has two open channels, one that yields SiO, which is the most exothermic and rapid, and another one resulting in SiS, which we also find to be moderately rapid. In the case of the reaction between Si and SO₂, we assume that the channels leading to SiO and SiS are equally rapid, with a rate coefficient of 10⁻¹⁰ cm³ s⁻¹. This same value is adopted for the SiS-forming reactions SiH + S and SiH + S₂, which have been calculated to be barrierless by Rosi et al. (2018), although in our calculations these reactions are not important routes to SiS. The reason is that SiH is not efficiently formed in gas phase from Si because the reaction between Si and H₂ to yield SiH is endothermic by about 35 kJ/mol, and thus too low even at 500 K. The reaction of SiS with atomic oxygen, which is calculated to be rapid, is also included. The rate coefficients adopted for these reactions are summarized in Table 1.

In Figure 5 we show the time evolution of the abundances of SiS, SiO, and some other interesting species. We first note that SiO is formed rapidly from atomic Si reacting with O₂ and OH at any temperature. This is also found in previous models of shocked protostellar outflows (e.g., Schilke et al. 1997; Gusdorf et al. 2008). The formation of SiS, on the other hand, is found to be substantially less favored than that of SiO, both regarding the maximum abundance reached and the velocity of the formation. More concretely, SiS is formed efficiently only at temperatures higher than 100 K and its formation is not as rapid as that of SiO. At temperatures of a few hundreds of K, the formation of SiS at early postshock times occurs at the expense of atomic silicon, in a similar way as occurs for SiO. The abundance of SiS however remains 2–3 orders of magnitude below that of SiO because, on the one hand, the sulfur precursors of SiS (SO and SO₂) are not as abundant as the oxygen precursors of SiO (O₂ and OH), and on the other, because SiS is efficiently destroyed by atomic oxygen forming SiO. At late postshock times, SiS can become quite abundant, where the main formation reactions are Si + SO and Si + SO₂. We note that given the simplicity of the calculations presented here, it is not advisable to extract quantitative conclusions on abundances and timescales. Nevertheless, we can extract useful information from these calculations: (1) SiS forms with a lower abundance and later than SiO, (2) SiS is efficiently destroyed through reaction with O, and (3) the main SiS-forming reactions are Si + SO and Si + SO₂.

Zoom In Zoom Out Reset image size

We have calculated the ab initio PESs of the SiOS system to compute the rate coefficients of the reaction of destruction of SiS with atomic O and of the SiS-forming reaction Si + SO in the temperature range of 50–2500 K. We find that the reactions SiS + O and Si + SO are rapid, with rate coefficients of a few times 10⁻¹⁰ cm³ s⁻¹, almost independent of temperature. In the reaction between Si and SO, production of SiO is favored over the formation of SiS, with a 5–7 times larger rate coefficient. The reaction of SiS with O atoms can play an important role in destroying SiS in envelopes around evolved stars and definitively is an efficient destruction pathway for SiS in protostellar outflows. Through a simple chemical model we find that in this latter type of environment SiS is formed with a lower abundance than that of SiO and that the reactions of the formation are Si + SO and Si + SO₂. Dedicated observational searches for SiS and detailed shock models including SiS chemistry are needed to better understand the differential behavior of SiO and SiS in interstellar and circumstellar clouds.

The research leading to these results has received funding from the European Research Council (ERC Grant 610256: NANOCOSMOS) and from Spanish MINECO through grants FIS2014-52172-C2, FIS2017-83473-C2, and AYA2016-75066-C2-1-P. M.A. also acknowledges funding support from the Ramón y Cajal programme of Spanish MINECO (RyC-2014-16277).