Quantum dynamical study of rotational excitations in SiS (X 1Σ+) molecule by H collisions.

The knowledge of accurate rate coefficients for collisional excitation of molecules by the abundant chemical species like He, H2 and H is important in modeling the conditions of interstellar medium. In the present paper, we computed the inelastic rotational cross sections and the corresponding rate coefficients of SiS molecule in its ground vibrational state in collisions with atomic hydrogen, H. We computed ab initio two-dimensional (rigid-rotor) potential energy surface for the H + SiS system at high accuracy for this purpose. The excitation cross sections are performed by numerically exact close-coupling quantum mechanical formalism up to collision energy of 2000 cm-1. The corresponding rate coefficients are obtained in the temperature range 5-300 K.

Introduction

Silicon sulfide (SiS) molecule was detected in the circumstellar envelope of carbon star IRC +10216 by its microwave emissions corresponding to J = 5 → 4 and J = 6 → 5 rotational transitions [1]. Grasshoff et al [2] detected another set of rotational transitions, J = 1 → 0 of SiS in IRC +10216 in ν = 0 as well as ν = 1 vibrational manifolds and J = 2 → 1 rotational line in ground vibrational state. Olofsson et al. [3] reported the detection of J = 5 → 4 transition of SiS and estimated its abundance to be 5.9 × 10−3 with respect to CO and 2.8 × 10−6 with respect to H2. A large number of rotational transitions including J = 16 → 15, J = 15 → 14, J = 14 → 13, J = 13 → 12, J = 12 → 11, J = 7 → 6, J = 6 → 5, J = 5 → 4 and J = 1 → 0 in IRC +10216. J = 7 → 6 and J = 6 → 5 transitions were also detected in other carbon-rich stars CIT 6, CRL 2688 and IRC +20370 [4]. SiS is supposed to be formed in thermochemical equilibrium conditions and is found to be more abundant in carbon-rich stars than in the oxygen-rich regions [5, 6, 7, 8, 9, 10, 11, 12, 13]. Although SiS is detected predominantly in the molecular envelopes of evolved stars it has also been found in interstellar medium (ISM). Dickinson et al. [14] and Ziurys [15, 16] reported the J = 6 → 5 and J = 5 → 4 rotational transitions in the Orion-KL and Sgr B2 which are star-forming regions. The column density is estimated to be 4 × 1015 cm−2 with a fractional abundance of about 3.9 × 10−9 with respect to H2. The maser action of J = 1 → 0 in SiS has also been detected in IRC +10216 [17]. Rosi et al. [18] have recently proposed two reactions by which SiS is thought to be formed in the ISM. Both the reactions considered, SiH + S and SiH + S2 are found to be barrierless and presumably form SiS in the ISM. Zanchet et al. [5] have studied both the formation and destruction pathways for SiS in space. The SiS forming reactions studied were Si + SiO and Si + SO2 while the SiS decomposition reaction considered was SiS + O.

Molecular and spectroscopic properties of SiS have been determined from various experimental studies [19, 20, 21, 22, 23]. Theoretically, SiS molecule has been studied thoroughly and the spectroscopic parameters like bond distance, dissociation energy, dipole moment and vibrational frequency were determined [24, 25, 26, 27, 28, 29, 30, 31]. Upadhyay et al. [32] computed a comprehensive list of rotation-vibrational lines and partition functions for 12 isotopologues of SiS employing the existing potential energy curve (PEC) [33] and a newly computed dipole moment curve.

Rotationally inelastic cross sections for excitation and de-excitation among first 26 rotational states of SiS by He collisions up to collision energies 1500 cm−1 have been computed in exact close-coupling method by Vincent et al. [34]. The corresponding rate coefficients in the temperature range 5–200 K have also been reported. Cross sections for rotational excitations among the first 51 rotational levels by para-H2 collisions have been carried out using coupled-states approximation by Lique et al. [35] for collision energies up to 2500 cm−1. The authors also tested the validity of the assumption that the rate coefficients for collisions with para-H2 can be obtained by multiplying the rate coefficients of He collisions by appropriate reduced mass correction. A comparison of the rate coefficients for He and para-H2 collisions with SiS found that the assumption is a reasonable first estimate. Lique and Kłos [36, 37] computed cross sections for rotational excitations of SiS with both para-H2 and ortho-H2 for collision energies up to 2500 cm−1 amongst the first 41 rotational levels of SiS and the corresponding state -to-state rate coefficients in the temperature range 5–300 K.

The abundance of molecular species in ISM can be determined by assuming local thermodynamic equilibrium (LTE) from molecular line observations in the absence of collisional rate coefficients [38]. It has been established [39, 40] that the LTE model fails in interpreting the molecular lines. If the collisional rate coefficients with the most abundant species in ISM H2, He and atomic hydrogen are available then non-LTE calculations can be used to accurately determine the column densities of molecular species. To the best of our knowledge there are no studies for collisions of SiS with atomic hydrogen. We believe that the rate coefficient information of collisional excitation of SiS with H will be useful in determining the molecular abundance of SiS in ISM. In the present study the SiS scattering by H atom is considered on a newly computed ab initio rigid-rotor PES. The paper is organized as follows. In section 2, the details of ab initio computations and PES calculations are discussed. Scattering calculations are described in section 3 and a brief summary and conclusions are given in section 4.

Methodology

Potential energy surface

The structural parameters of SiS are computed first using Dunning's aug-cc-pVQZ [41] basis set at MRCI level of accuracy and are compared with available literature results in Table 1 and they are found to be in good agreement. The HSiS and HSSi structures are optimized at AVQZ/MRCI level of accuracy. Both the isomers are found to have bent structures. The optimized parameters are given in Table 2. The HSSi isomer is found to be higher in energy than the HSiS isomer by 0.171 eV. Before computing the rigid rotor potential energy surface of H–SiS, we have computed the potential energy curves for collinear and off-collinear approach of H towards SiS.

Table 1

Spectroscopic parameters of SiS (X1Σ+).

	re (a0)	De (eV)	ωe (cm−1)	Be (cm−1)
Present	3.6762	6.24	734.79	0.2985
Chattopadhyaya et al[30]	3.6987	5.85	733	-
Shi et al[31]	3.6615	6.39	747.07	0.3010
Exp. [42]	3.6464	6.40	749.64	0.3035

Table 2

Optimized equilibrium geometries of HSiS and HSSi (X2A′).

	R(HSi)/R(HS) (bohr)	R(SiS) (bohr)	θ(H–Si–S)/θ(H–S–Si) (degree)
HSiS	2.8338	3.7086	120.4
HSSi	2.5385	3.9934	101.0

In Fig. 1, the one-dimensional cuts of potential energy surfaces of 2Σ+ and 2Π electronic states for collinear approach of H towards Si end of SiS (γ = 0°), 12A'/12A″ electronic states for off-collinear (γ = 60°) approach of H are presented. The 2Σ+ and 12A′ electronic states asymptotically correlate to H + SiS (X 1Σ+) while the 2Π, 22A′ and 12A″ asymptotically correlate to H + SiS (A 3Π). The transition states (TS) to the following reactive processes have also been optimized:

Fig. 1

1D-cuts of potential energy surfaces for collinear (γ = 0°) and off-collinear (γ = 60°) approach of H towards SiS. The asymptotic correlation of all the surfaces is also shown.

The saddle points for R1 and R2 have been optimized at AVQZ/MRCI level of accuracy. The TS for R1 has coordinates, R(Si–S) = 3.963 a0, R(Si–H) = 3.172 a0 and θ(H–S–Si) = 47.7° whereas for TS of R2 the coordinates are R(Si–S) = 3.963 a0, R(S–H) = 2.973 a0 and θ(H–S–Si) = 52.1°. The barrier for both the reactions R1 and R2 is found to be 3138 cm−1 which is much higher than the collision energy employed in the present study. It can be safely assumed that in the collision energy range considered here, only inelastic scattering occurs and reactive collisions can be ruled out. The validity of the rigid-rotor model in the present study can be discussed in the light of the computed barriers to the reactive processes. Lique [43] established the validity of the rigid-rotor model for H and H2 collisions with interstellar species if the reactions involved proceeds through a barrier and for low collision energies. Hence, the two-dimensional rigid-rotor PES of ground state of H–SiS collision system (2Σ+/2A′) correlating asymptotically to H + SiS (X 1Σ+) has been computed which was employed in the present dynamics study.

The PES computations were performed in the Jacobi scattering coordinates (R, r, γ), where r is the interatomic distance of SiS, R is the distance of H from the center of mass of SiS and angle between R and r is given by γ = cos−1(R.r). The approach of H towards Si end of SiS is taken to be γ = 0°. Calculations were carried out in the C2υ symmetry for collinear geometry and in the Cs symmetry for the off-collinear geometries. The ground state surface was computed for the 2Σ+ for the collinear approach of H and 2A′ for off-collinear approaches of H. Computations were carried out at internally contracted multi-reference configuration interaction (icMRCI) [44, 45, 46] with Dunning's aug-cc-pVQZ basis set using MOLPRO 2010.1 [47] suite of programs. At the SCF level, the basis set produced 214 molecular orbitals (MOs) from contracted Gaussian atomic orbitals and they are listed as [85a1, 51b1, 51b2, 27a2] and [136a′, 78a''] in the C2υ and the Cs point groups respectively. The [6a1, 2b1, 2b2] MOs in the C2υ and [8a′, 2a''] MOs in the Cs were treated as core orbitals. The ground state electronic configuration for the C2υ is [13a1, 3b1] with the thirteenth a1 orbital singly occupied to give the 2Σ+ state accounting for 31 total electrons. The electronic configuration is [13a′, 3a''] with the thirteenth a' orbital singly occupied in the Cs point group calculation. In the complete active space self-consistent field (CASSCF) calculations [6-11a1, 3–4b1, 3–4b2] for C2υ and [9-15a′, 3-4a''] for the Cs are the active orbitals respectively. The wavefunction in a typical calculation consisted of 1568 configuration state functions (CSFs) with 2640 Slater determinants. The configuration interaction (CI) calculations were performed with a reference space of 696 configurations; N, N-1 and N-2 internal configurations are 2304, 2907 and 3139. The threshold value of the CSFs selection was kept at 3.2 × 10−7 a.u. The energy values reported are Davidson corrected [48]. The size-inconsistency inherent to MRCI method could be corrected to some extent by including Davidson correction (MRCISD + Q) [49]. The PES was obtained on the following grid points: γ = 0°–180° (10°); R = 0.8–2.2 (0.2), 2.3–3.5 (0.05), 3.6–6.0 (0.1), 6.2–10.0 (0.2), 11.0–20.0 (0.5), 21.0–30.0 (1.0) and 35.0–80.0 (5.0). R is in atomic units and the numbers in the parentheses indicate the increment value. The value of r is frozen at experimental equilibrium bond length value of 3.6459 a0.

The basis set superposition error (BSSE) is accounted for at all geometries in computing interaction potential of H–SiS , by employing the counterpoise correction method of Boys and Bernandi [50]:where E (R, r, γ) is the total energy of HSiS triatomic system, EH and ESiS are the energies of the subsystems H and SiS, respectively, computed in the full basis set of the system.

The minimum of the rigid rotor PES is found at R = 4.10 a0 and γ = 40°. The well depth of the PES is found to be 1.366 eV. The ab initio points of the PES are fitted using cubic spline interpolation method on a finer grid. The multipolar expansion coefficients are computed using the splined interaction potential.

The rigid-rotor interaction potential has been fitted to the following analytic expression:where V’s are expansion coefficients with λ varying from 0 to 18 (number of γ values) and P’s are Legendre polynomials. The interaction anisotropy can be examined in terms of V’s which are plotted as a function of R and shown in Fig. 2. The V1 term shows the deepest attractive interaction potential well followed by V2, V0, and V3. V4, V5 and V6 are repulsive.

Fig. 2

Legendre radial expansion coefficients of interaction potential.

Scattering calculations

The exact close coupling (CC) formalism, due to Arthurs & Dalgarno [51], was used in the present study. Full close-coupling calculations have been carried in the collision energy range of 0–2000 cm−1 for rotational excitations. Time-independent coupled scattering equations have been solved to compute cross-sections as implemented in the MOLSCAT code [52]. For high collision energies a parallel version of the code, PMP MOLSCAT [53] has been used. The numerical convergence of the cross sections is ensured in terms of number of rotational levels included in the basis as well as in terms of total angular momentum. Since the interaction potential well is deep, the number of rotational levels of SiS that get populated will be large for a given collision energy. So we included a large number of closed channels in the rotational basis to obtain numerically converged values for cross section. The number of closed channels varied between 15 and 30 for collision energies below 1500 cm−1. For collision energies above 1500 cm−1 we included 10 closed channels which we found to be sufficient to obtain convergence in cross sections for transitions between lowest seven rotational energy levels of SiS. The convergence is tested by including 20 closed channels in the basis. The cross sections are shown in Table 3. It can be seen that inclusion of 10 closed channels is sufficient to obtain the numerically converged cross sections. Maximum value of rotational quantum number taken is j = 91 at Ec.m = 2000 cm−1. Also, convergence of cross-sections is achieved through sufficient number of partial waves, for instance, at Ec.m = 2000 cm−1, Jtot (total angular momentum) is kept at 310. The close-coupled radial equations were numerically integrated using the Log Derivative propagator of Manolopoulos [54]. The following input parameters are taken in the calculation: rotational constant of SiS, B = 0.3035279 cm−1, reduced mass of the system (H–SiS), μ = 0.99138 a.m.u. with values of R and R as 0.8 a0 and 80.0 a0, respectively with ΔR = 0.1 a0. The energy spacing is also carefully spanned to capture the resonances at low collision energies. Below 50 cm−1 the energy spacing is 0.2 cm−1, between 50 cm−1 to 100 cm−1 the spacing is 1 cm−1, between 100 cm−1 to 500 cm−1 the spacing is 5 cm−1, between 500 cm−1 to 1000 cm−1 the spacing is 50 cm−1 and between 1000 cm−1 to 2000 cm−1 the spacing is 100 cm−1.

Table 3

Test of convergence of cross sections as a function of number of closed channels (10 and 20) included in the rotational basis for collision energies from 1500 cm−1 and above.

σ (Å2)	Ec.m = 1500 cm−1		Ec.m = 1800 cm−1		Ec.m = 2000 cm−1
	10	20	10	20	10	20
σ0→1	113.35	113.32	106.71	106.54	102.99	102.91
σ1→2	145.26	145.34	139.04	139.47	135.44	135.12
σ0→2	233.26	233.26	228.11	227.90	223.30	222.97
σ2→3	146.21	146.29	141.15	141.53	138.05	138.22
σ3→4	136.24	136.45	130.78	131.01	127.52	127.24

Results and discussion

Integral cross sections

The integral cross sections for rotational excitations j → j' have been computed using the following formula.where the total angular momentum Jtot = l + j includes orbital angular momentum of HSiS and the rotational angular momentum of the diatom. represents the wavevector where Ec.m is center-of-mass collision energy and is the S-matrix.

The close-coupling integral cross sections as a function of collision energy are shown in Fig. 3 for transitions from j = 0 to j = 1, 2, 3, 4 and 5. All the cross sections decrease with increase in collision energy. Of the five transitions shown, j = 0 to j = 4 transition has the lowest cross section. The conspicuous feature of all the cross section curves is the oscillatory structure for all the transitions up to around 100 cm−1. This structure is attributed to the resonances (shape and Feshbach) appearing as a result of the quasibound states. The shape resonances appear due to the tunneling through the centrifugal barrier [34, 55, 56, 57], whereas the Feshbach resonances appear as a result of temporary trapping of H in the attractive interaction potential well forming a quasibound states of H–SiS complex. Similar quasibound states formation is observed in He–SiS [34], ortho-H2-SiS and para-H2-SiS [36, 37]. It is proposed that both the shape and Feshbach type resonances occur in the same energy region for SiS as the rotational energy levels in SiS are very closely spaced because of its small rotational constant [34].

Fig. 3

Computed exact close-coupling inelastic rotational integral cross sections for H + SiS collisions as a function of collision energy from j = 0 level.

Rate coefficients

The computed rotational cross-sections are used to calculate rate-coefficient as function of temperature:where k is the Boltzmann constant, μ is the reduced mass of the system and Ec.m is center-of-mass collision energy. The computed rate coefficients are plotted as a function of temperature in the range 5–300 K and shown in Fig. 4 for Δj = 1 and in Fig. 5 for Δj = 2 transitions.

Fig. 4

Computed rate coefficients as a function of temperature for Δj = 1 transitions.

Fig. 5

Computed rate coefficients as a function of temperature for Δj = 2 transitions.

The rate coefficients for Δj = 1 and Δj = 2 transitions show similar trend in which they increase with increasing temperature up to 300 K. But the rate coefficient corresponding to j = 0 to j = 2 transition increases up to around 150 K and then starts decreasing. The profile of rate coefficient for j = 0 to j = 2 transition could be due to the almost flat cross section between the collision energies 200 cm−1 and 1000 cm−1, whereas the cross sections for other transitions decrease monotonically with collision energy. The rate coefficients for j = 1 to j = 2 is the leading one in Δj = 1 transitions followed by j = 2 to j = 3 and j = 0 to j = 1.

The rate coefficients are of the order of ∼10−9 cm3 s−1 for the collision of H with SiS which are two orders of magnitude greater than the values reported in He – SiS [34] and in ortho/para-H2 – SiS [36, 37] systems in the same temperature range. The rate coefficient for ortho-H2-SiS is j = 1 → 0 is the largest and is an order of magnitude lower than that that for H–SiS. The rate coefficients are large for low j transitions in case of H–SiS collisions and progressively decrease for transitions amongst higher j values, whereas in H2–SiS and He–SiS collisions rate coefficients for transitions between higher rotational levels dominate the rate coefficients of transitions amongst lower j levels. In He–SiS and para-H2-SiS collisions a propensity for Δj = 2 is found while in ortho-H2-SiS a propensity for Δj = 1 is observed. In H–SiS collisions no such propensity has been observed.

Conclusions

The hydrogenated isomers of SiS, HSiS and HSSi have been optimized and the structural parameters are reported. To the best of our knowledge there are no previous literature reports for the two isomers. Both the isomers have bent geometries and have 2A′ as their ground state symmetry. The HSiS isomer is found to be the global minimum (in the two-dimensional rigid rotor context). The HSSi isomer is higher in energy than the HSiS isomer by 5.71 kcal/mol. Ab initio two-dimensional (rigid-rotor) PES for H - SiS collision system of the ground electronic state 2Σ+/2A′ asymptotically correlating to H + SiS (X 1Σ+) has been computed at MRCI/aug-cc-pVQZ level of accuracy for 19 angles between 0° and 180°. The potential anisotropy has been examined in terms of radial multipolar expansion.

Rotational inelastic cross sections are computed by exact close-coupling method for collision energies up to 2000 cm−1. The cross sections are examined for transitions among first 7 rotational levels of SiS molecule. Rate coefficients are computed in the temperature range 5–300 K. The rate coefficients are found to be two orders of magnitude greater than those reported for He - SiS and ortho/para-H2 – SiS systems in the same temperature range. Thus even if atomic hydrogen is not the dominant collision partner, the rotational population of SiS may be significantly affected by collisions with H. We believe these rate coefficients are useful for accurately determining the molecular abundance of SiS where the LTE conditions are not fulfilled in modeling the conditions of the interstellar medium accurately.

Declarations

Author contribution statement

Bhargava Anusuri: Conceived and designed the experiments; Performed the experiments; Analyzed and interpreted the data; Contributed reagents, materials, analysis tools or data; Wrote the paper.

Funding statement

This work was supported by a fellowship from IIT Madras and UGC, New Delhi.

Competing interest statement

The authors declare no conflict of interest.

Additional information

No additional information is available for this paper.