11institutetext:  Laboratory for Astrophysics, Leiden Observatory, Leiden
University, PO Box 9513, 2300-RA, Leiden, The Netherlands 22institutetext:
Instituto de Física Fundamental (IFF-CSIC), C.S.I.C., Serrano 123, 28006
Madrid, Spain 33institutetext: University of Firat, Department of Physics,
23169 Elazig, Turkey 44institutetext: Institute of Physics, Faculty of
Physics, Astronomy and Informatics, Nicolaus Copernicus University in Torun,
Grudziadzka 5, 87-100 Torun, Poland 55institutetext: Université Paris-Saclay,
CEA, AIM, Dèpartement d’Astrophysique (DAp), F-91191 Gif-sur-Yvette, France
66institutetext: Observatorio Astronómico Nacional (IGN), c/ Alfonso XII 3,
28014 Madrid, Spain. 77institutetext: Laboratoire d’astrophysique de Bordeaux,
Univ. Bordeaux, CNRS, B18N, allée Geoffroy Saint-Hilaire, 33615 Pessac, France
88institutetext: Institut des Sciences Moléculaires (ISM), CNRS, Univ.
Bordeaux, 351 cours de la Libération, F-33400, Talence, France
99institutetext: Sorbonne Université, Observatoire de Paris, Univerité PSL,
CNRS, LERMA, 92190 Meudon, France 1010institutetext: Centre for Astrochemical
Studies, Max-Planck-Institute for Extraterrestrial Physics,Giessenbachstrasse
1, 85748, Garching, Germany 1111institutetext: Institut de Planétologie et
d’Astrophysique de Grenoble (IPAG), Université Grenoble Alpes, CNRS, F-38000
Grenoble, France 1212institutetext: Institut de Radioastronomie Millimétrique
(IRAM), 300 Rue de la Piscine, F-38406 Saint-Martin d’Hères, France

# Gas phase Elemental abundances in Molecular cloudS (GEMS)

VIII. Unlocking the CS chemistry: the CH + S$\rightarrow$ CS + H and C2 \+
S$\rightarrow$ CS + C reactions

Carlos M. R. Rocha 11 Octavio Roncero 22 Niyazi Bulut 33 Piotr Zuchowski 44
David Navarro-Almaida 55 Asunción Fuente 66 Valentine Wakelam 77 Jean-
Christophe Loison 88 Evelyne Roueff 99 Javier R. Goicoechea 22 Gisela
Esplugues 66 Leire Beitia-Antero 66 Paola Caselli 1010 Valerio Lattanzi 1010
Jaime Pineda 1010 Romane Le Gal 11111212 Marina Rodríguez-Baras 66 Pablo
Riviere-Marichalar 66

###### Abstract

Context. Carbon monosulphide (CS) is among a few sulphur-bearing species that
has been widely observed in all environments, including the most extreme ones
such as diffuse clouds. Moreover, it has been widely used as a tracer of the
gas density in the interstellar medium in our Galaxy and external galaxies.
Therefore, the full understanding of its chemistry in all environments is of
paramount importance for the study of the interstellar matter.

Aims. Our group is revising the rates of the main formation and destruction
mechanisms of CS. In particular, we focus on those which involve open-shell
species for which the classical capture model might not be accurate enough. In
this paper, we revise the rates of reactions CH + S $\rightarrow$ CS + H and
C2 \+ S $\rightarrow$ CS + C. These reactions are important CS formation
routes in some environments such as dark and diffuse warm gas.

Methods. We performed ab initio calculations to characterize the main features
of all the electronic states correlating to the open shell reactants. For CH+S
we have calculated the full potential energy surfaces (PES) for the lowest
doublet states and the reaction rate constant with a quasi-classical method.
For C2+S, the reaction can only take place through the three lower triplet
states, which all present deep insertion wells. A detailed study of the long-
range interactions for these triplet states allowed to apply a statistic
adiabatic method to determine the rate constants.

Results. Our detailed theoretical study of the CH + S $\rightarrow$ CS + H
reaction shows that its rate is nearly independent on the temperature in a
range of 10$-$ 500 K with an almost constant value of 5.5 $\times$ $10^{-11}$
cm3 s-1 at temperatures above 100 K. This is a factor $\sim$ 2$-$3 lower than
the value obtained with the capture model. The rate of the reaction C2 \+ S
$\rightarrow$ CS + C does depend on the temperature taking values close to 2.0
$\times$ $10^{-10}$ cm3 s-1 at low temperatures and increasing to $\sim$ 5.0
$\times$ $10^{-10}$ cm3 s-1 for temperatures higher than 200 K. In this case,
our detailed modeling taking into account the electronic and spin states
provides a rate higher than the one currently used by factor of $\sim$2.

Conclusions.  These reactions were selected for involving open-shell species
with many degenerate electronic states, and, unexpectedly, the results
obtained in the present detailed calculations provide values which differ a
factor of $\sim$2$-$3 from the simpler classical capture method. We have
updated the sulphur network with these new rates and compare our results in
the prototypical case of TMC1 (CP). We find a reasonable agreement between
model predictions and observations with a sulphur depletion factor of 20
relative to the sulphur cosmic abundance. However, it is not possible to fit
all sulphur-bearing molecules better than a factor of 10 at the same chemical
time.

###### Key Words.:

Astrochemistry – ISM: abundances – ISM: molecules

##  1 Introduction

Astrochemistry has become a necessary tool for understanding the interstellar
medium of our Galaxy and external galaxies. Nowadays, we are aware of the
existence of nearly 300 molecules in the interstellar and circumstellar
medium, as well as of around 70 molecules in external galaxies (for a complete
list, see the Cologne Database for Molecular
Spectroscopy111https://cdms.astro.uni-koeln.de/). Although the sulphur cosmic
elemental abundance is only ten times lower than that of carbon
(S/H$\approx$1.5 $\times$ 10-5), only 33 out of the currently detected
interstellar molecules contain sulphur atoms. This apparent lack of chemical
diversity in astrophysical sulphur-bearing molecules is the consequence of a
greater problem in astrochemistry: there is an unexpected paucity of sulphur-
bearing species in dense molecular clouds and star-forming regions. In such
dense regions, the sum of the observed gas-phase abundances of sulphur-bearing
species (the most abundant are SO, SO2, H2S, CS, HCS+, H2CS, C2S, C3S, and NS
) constitutes only $<$1% of the expected amount (Agúndez & Wakelam 2013;
Vastel et al. 2018; Rivière-Marichalar et al. 2019; Hily-Blant et al. 2022) .
One could think that most of the sulphur is locked on the icy grain mantles,
but a similar trend is encountered within the solid phase, where s-OCS
($"$s-$"$ indicates that the molecule is in the solid phase) and s-SO2 are the
only sulphur-bearing species detected thus far (Palumbo et al. 1995, 1997;
Boogert et al. 1997; Ferrante et al. 2008), and only upper limits to the s-H2S
abundance have been derived (Jiménez-Escobar & Muñoz Caro 2011). Recent
observations with James Webb Space Telescope (JWST) did not detect s-H2S
either (McClure et al. 2023). According to these data, the abundances of the
observed icy species account for $<$ 5% of the total expected sulphur
abundance. This means that 94% of the sulphur is missing in our counting. It
has been suggested that this missing sulphur may be locked in hitherto
undetected reservoirs in gas and icy grain mantles, or as refractory material
Shingledecker et al. (2020). In particular, laboratory experiments and
theoretical work shows that sulphur allotropes, such as S8, could be an
important refractory reservoir (Jiménez-Escobar et al. 2014; Shingledecker et
al. 2020; Cazaux et al. 2022).

Gas phase Elemental abundances in Molecular CloudS (GEMS) is an IRAM 30m Large
Program designed to estimate the S, C, N, and O depletions and the gas
ionization fraction as a function of visual extinction in a selected set of
prototypical star-forming filaments in low-mass (Taurus), intermediate-mass
(Perseus), and high-mass (Orion) star forming regions (Fuente et al. 2019;
Navarro-Almaida et al. 2020; Bulut et al. 2021; Rodríguez-Baras et al. 2021;
Esplugues et al. 2022; Spezzano et al. 2022; Fuente et al. 2023). Determining
sulphur depletion is probably the most challenging goal of this project. The
direct observation of the potential main sulphur reservoirs (s-H2S, s-OCS,
gas-phase atomic S) remains difficult even in the JWST era. Therefore, sulphur
elemental abundance needs to be estimated by comparing the observed abundances
of rarer species such as CS, SO, HCS+, H2S, SO2, and H2CS, with the
predictions of complex gas-grain chemical models. All this causes that the
sulphur chemistry in cold dark clouds remains as a puzzling problem. The
development of accurate and complete sulphur chemical networks is therefore a
requisite to disentangle the sulphur elemental abundance. Within the context
of GEMS project, we have carried out a large theoretical effort to improve the
accuracy of key reaction rates of the sulphur chemical network, with a special
interest in those associated with the formation and destruction paths of SO
and CS which are the gas-phase sulphur bearing species observed in more
different environments. We estimated the rates of the reactions S + O2
$\rightarrow$ SO + O (Fuente et al. 2016) and SO + OH $\rightarrow$ SO2 \+ H
(Fuente et al. 2019) at the low temperatures prevailing in dark clouds. These
reactions drives the SO chemistry in these cold environments. Bulut et al.
(2021) estimated the rate constant of the CS + O $\rightarrow$ SO + O
reaction, that had been proposed as efficient CS destruction mechanisms in
molecular clouds. In this work we shall study two formation reactions of CS,
which are thought to be important in regions with low ionization fraction:
CH(${}^{2}\Pi$)+S(3P) and C2(${}^{1}\Sigma_{g}^{+}$)+ S(3P). In the two cases,
the two reactants are radicals presenting several degenerate or quasi-
degenerate electronic states: for CH + S there are 36 degenerate states and
for C2+S the first excited C2(${}^{3}\Pi_{u}$) states are only 0.089 eV above
the C2(${}^{1}\Sigma_{g}^{+}$) ground state . This makes the experimental
determination of their rates difficult, because of the low densities in which
two radical species are obtained, and because of the possibility of self-
reactions. Thus, the most accurate theoretical determination of the reaction
rate constants is desirable to proper bound the abundance of CS in chemical
models.

The reactions rates currently available for these two reactions were obtained
with a classical capture method (Vidal et al. 2017). Dealing with open-shell
systems, there are several degenerate electronic states for the reactants, not
all of them leading to the desired product, and the crossing among them
originate barriers. All these effects are ignored in the classical capture
method. Precise calculations like the ones in this article are needed to
improve the accuracy of our chemical networks and to estimate the
uncertainties associated with the less demanding classical capture methods.

##  2 CH+S reaction

The reaction

$\displaystyle\textrm{CH}(^{2}\Pi)+\textrm{S}(^{3}\textrm{P})$
$\displaystyle\rightarrow$
$\displaystyle\textrm{CS}(X^{1}\Sigma^{+},a^{3}\Pi)+\textrm{H}\quad\quad(a)$
$\displaystyle\rightarrow$
$\displaystyle\textrm{SH}(X^{2}\Pi)+\textrm{C}(^{3}\textrm{P})\quad\quad(b)$

presents several rearrangement channels (CS and SH products) with several
electronic states in each case. There are 36 degenerate states (neglecting
spin-orbit couplings) in the CH(${}^{2}\Pi$)+S(${}^{3}P$) entrance channel.
The same number of degenerate states are in the SH($X^{2}\Pi$) + C$(^{3}P)$
rearrangement channel, which is only $\approx$ 0.08 eV below
CH(${}^{2}\Pi$)+S(${}^{3}P$). The reaction SH($X^{2}\Pi$) +
C(${}^{3}P$)$\rightarrow$ CS(X${}^{1}\Sigma^{+},a^{3}\Pi$) + H has already
been studied theoretically for the ground ${}^{2}A^{\prime}$ state (Stoecklin
et al. 1988, 1990a, 1990b; Song et al. 2016; Voronin 2004) and for the excited
${}^{2}A^{\prime\prime}$ state (Zhang et al. 2018).

The 36 degenerate electronic states (neglecting spin-orbit couplings) involved
in Eq. (2.a) consist of 6 spin states times 6 orbital states. The spin states
are quartet and doublets, while the orbital states are splitted in A’ and A”
states, $i.e.$, symmetric or antisymmetric with respect to the inversion
through the plane of the molecule. Of all these 36 states, only a doublet
(i.e. 2 states) correlates to the CS($X^{1}\Sigma^{+}$) states, which
correspond to the ground adiabatic ${}^{2}A^{\prime}$ states. Reaction (2.a)
towards the $X^{1}\Sigma^{+}$ is exothermic by $\approx$ 3.5 eV. However, the
excited CS($a^{3}\Pi$) is about 3.42 eV above CS($X^{1}\Sigma^{+}$), and the
reaction is nearly thermoneutral.

###  2.1 Ab initio calculations

To describe the electronic correlation along the reaction, we use here the
internally contracted multi-reference configuration interaction (ic-MRCI)
method (Werner & Knowles 1988a, b) including the Davidson correction
(hereafter called MRCI+Q) (Davidson 1975) and the calculations are performed
with the MOLPRO suite of programs (Werner et al. 2012). Three electronic
states are calculated for each symmetry, 3 ${}^{2}A^{\prime}$ and 3
${}^{2}A^{\prime\prime}$, and the same for the quartet states.

In these calculations, the molecular orbitals are optimized using a state-
averaged complete active space self-consistent field (SA-CASSCF) method, with
an active space of 10 orbitals (7 and 3 of $a^{\prime}$ and $a^{\prime\prime}$
symmetry, respectively). Five ${}^{2,4}A^{\prime}$ and four
${}^{2,4}A^{\prime\prime}$ electronic states are calculated and simultaneously
optimized at CASSCF level, using a dynamical weighting factor of 10. In all
these calculations the aug-cc-pVTZ (aVTZ) basis set is used (Dunning & Jr.
1989). For the ic-MRCI calculations, 6 orbitals are kept doubly occupied,
giving rise to $\approx$ $5\times 10^{6}$ $(382\times 10^{6})$ contracted
(uncontracted) configurations. Checks made with the aug-cc-pVQZ basis gave
nearly parallel results along some of the minimum energy paths (MEPs) shown
below. For this reason, we kept the aVTZ basis set to build the PES.
Nevertheless, to better describe the long range part, we shall use a AV5Z
basis, as discussed below.

In Fig. 1 the first 5 A’ and 4 A” states for doublet (bottom) and quartet
(top) multiplicities are shown, calculated at CASSCF level. With no spin-orbit
couplings, the six electronic states for each multiplicity are degenerate for
long distances between S(${}^{3}P$) and CH($X^{2}\Pi$), and this asymptote is
taken as the origin of energy. When they approach, they cross with the excited
states correlating to S(${}^{3}P$)+CH($a^{4}\Sigma^{-}$). After the crossing,
there are three (two) curves for the doublet (quartet) states that become
negative. For the doublets, one of these states correlates to
CS($X^{1}\Sigma^{+}$) and two, nearly degenerate, to CS($a^{3}\Pi_{r}$)
states.

Figure 1:  CASSCF energies along the reaction coordinate defined as RCS-RCH
distance difference, obtained for an angle $\gamma_{SCH}$= 179.9o for the
doublet (bottom panel) and quartet (top panel) states, for the ${\rm CH}+{\rm
S}\rightarrow{\rm CS}+{\rm H}$ reaction. Five A’ and 4 A” electronic states
are considered

This degeneracy is recovered in the MRCI+Q calculations shown in Fig. 2. The
two degenerate states, 1${}^{2}A^{\prime}$ and 1${}^{2}A^{\prime\prime}$ (and
1${}^{4}A^{\prime}$ and 1${}^{4}A^{\prime\prime}$ for quartets), correspond to
the HCS(${}^{2}\Pi$) radical (and the excited HCS(${}^{4}\Pi$)) state, which
experience Renner-Teller effects, as studied by Senekowitsch et al. (1990). In
the case of doublets, the $X^{2}A^{\prime}$ state (of ${}^{2}\Pi$ character at
collinear geometry) crosses with $2^{2}A^{\prime}$ (of ${}^{2}\Sigma$
character at collinear geometry), which correlate with the
CS(X${}^{1}\Sigma^{+}$) state of the products.

The crossing at collinear geometry is a conical intersection, between $\Sigma$
and $\Pi$ states. As long as the system bents, there are couplings and the
crossing is avoided. As a consequence, the ground $X^{2}A^{\prime}$ state
correlates to the CS(X${}^{1}\Sigma^{+}$) products, $i.e.$ only these 2
doublets among the 36 states correlating to the CH(${}^{2}\Pi$)+S(${}^{3}P$)
reactants. Therefore, only the ground electronic state is needed to describe
the reaction 2. The energy difference between the CH(${}^{2}\Pi$)+S(${}^{3}P$)
and CS(X${}^{1}\Sigma^{+}$) + H, at their corresponding equilibrium geometries
is $D_{e}$ = 3.52 eV, and $D_{0}$ = 3.62 eV when including zero-point energy
(using the fit described below). Since there is no direct measurements for
this reaction, the experimental exothermicity is estimated from the
dissociation energies D0=D${}^{CH}_{0}$ \- D${}^{CS}_{0}$ = 3.89 eV, about a
10% higher than the present results using the values reported by Huber &
Herzberg (1979).

Figure 2:  MRCI+Q energies along the reaction coordinate defined as RCS-RCH
distance difference, obtained for an angle $\gamma_{SCH}$= 179.9o for the
doublet (bottom panel) and quartet (top panel) states, for the ${\rm CH}+{\rm
S}\rightarrow{\rm CS}+{\rm H}$ reaction. Three A’ and three A” electronic
states are considered in each case.

###  2.2 Long range interaction

MRCI calculations are not size consistent and the Davidson correction (+Q) is
not adequate to describe close lying electronic states. Therefore MRCI+Q
method introduces inaccuracies in the long range region, which need to be
described accurately to obtain good rate constants at low temperatures.
Coupled Cluster methods are size consistent and are typically considered as a
holy grail to describe long range interactions, but in the presence of two
open shell reactants, is not expected to yield good accuracy either.

As an alternative, here we use the CASSCF method with a larger aV5Z basis set
(Dunning & Jr. 1989). Several points have been calculated in Jacobi
coordinates, for CH($r_{CH})$ = 1.1199 Å and R= 10-50 Å where $R$ is the
distance of CH center-of-mass and the S atom, as a function of the angle
$\theta$, with $\cos\theta={\bf R}\cdot{\bf r}_{CH}/Rr_{CH}$. At long
distances, the system behaves as dipole-quadrupole interaction, with the
dipole of CH interacting with the quadrupole of a P sulphur atom, whose
analytical form is (Zeimen et al. 2003)

$\displaystyle V_{LR}(R,\theta)=M_{QD}(\theta){Q_{A}d_{B}\over
R^{4}}+M_{QQ}(\theta){Q_{A}Q_{B}\over R^{5}},$ (2)

with $M_{QD}(\theta)$ and $M_{QQ}(\theta)$ being 3$\times$3 matrices,
depending on Legendre polynomials (Zeimen et al. 2003). The eigenvalues of
this matrix properly describe the angular dependence of the adiabatic states.
Thus, only 2 effective parameters are needed to fit the ab initio points,
$Q_{A}d_{B}$ = 0.293 hartree Å-4 and $Q_{A}Q_{B}$ = 0.092 hartree Å-5. Two
families of three states are considered, separately,
$1^{2}A^{\prime\prime},1^{2}A^{\prime},2^{2}A^{\prime\prime}$ and
$2^{2}A^{\prime},3^{2}A^{\prime},3^{2}A^{\prime\prime}$, corresponding to the
two $\Pi$ states of CH interacting with the 3 $P$ states of sulphur. The
excellent agreement between calculated and fitted analytical expressions, as
shown in Fig. 3, demonstrates the adequacy of the analytical fit in the
asymptotic region.

Figure 3:  Points are the CASSCF(aV5Z) energies (in cm-1) as a function of the
Jacobi angle $\theta$ (defined in the text) for $r_{CH}$= 1.1199Å and $R$= 10,
15, 20 and 30Å, as indicated in each panel, for the $1^{2}A^{\prime\prime}$,
$1^{2}A^{\prime}$ and $2^{2}A^{\prime\prime}$ states. Lines correspond to the
long range analytical fit of Eq. (2), taken from (Zeimen et al. 2003).

###  2.3 Analytical fit of the ground state

The analytical representation of the potential energy surface of the ground
1${}^{2}A^{\prime}$ electronic state is described by two terms

$\displaystyle
V(r_{CH},r_{CS},r_{SH})=E_{g}^{FF}(r_{CH},r_{CS},r_{SH})+V^{3B}(r_{CH},r_{CS},r_{SH}),$
(3)

where $V^{3B}$ is the three-body term added to the zero order description
provided by the lowest root, $E_{g}^{FF}(r_{CH},r_{CS},r_{SH})$, of the
3$\times$3 reactive force-field matrix defined as (Zanchet et al. 2018;
Roncero et al. 2018; Goicoechea et al. 2021)

$\displaystyle
H^{FF}=\left({\scriptsize\begin{array}[]{ccc}V_{CH}+W^{1}_{CS}+W^{1}_{SH}+V_{LR}&V_{12}&V_{13}\\\
V_{12}&V_{CS}+W^{2}_{CH}+W^{2}_{SH}&V_{23}\\\
V_{13}&V_{23}&V_{SH}+W^{3}_{CH}+W^{3}_{CS}\end{array}}\right).$ (7)

The diagonal terms describe each of the rearrangement channels, in which
$V_{CH}(r_{CH})$, $V_{CS}(r_{CS})$ and $V_{SH}(r_{SH})$, are fitted using the
diatomic terms of Aguado & Paniagua (1992). In the reactants channel 1, the
long-range term given by Eq. 2 is included. $W_{AB}$ are Morse potentials
whose parameters are determined to describe each channel independently.
Finally, the non-diagonal terms $V_{ij}$ essentially are build as Gaussian
functions $exp(-\alpha(H^{FF}_{ii}-H^{FF}_{jj})^{2})$ depending on the energy
difference between the diagonal terms of the corresponding force-field matrix.
The parameter $\alpha$ is determined to fit the transition between
rearrangements.

The three-body term, $V^{3B}$ in Eq. (3), is described with the method of
Aguado & Paniagua (1992) using a modification of the program GFIT3C (Aguado et
al. 1998). In the fit, about 7500 ab initio points, calculated at MRCI+Q
level, have been included. These points are mainly composed by a grid of 22
points for $0.6\leq r_{CH}\leq$ 8 Å, 27 points for 1$\leq r_{CS}\leq$ 8 Å, and
19 points in the $\theta_{HCS}$ angle, in interval of 10o. The points have
been weighted, with a weight of 1, for those with energy up to 1 eV, from the
entrance CH+S channel at 50 Å, and with a Gaussian function for higher
energies, and considering a minimum weight of 10-4. The final fit uses
polynomials up to order 10, with an overall root-mean-square error of 0.07 eV.

Figure 4: Contour plots of the analytical potential energy surfaces fitted to
ab initio points, corresponding to the reactant channel for $R_{CH}$=1.1199
Å(bottom panel) and to the products channel for $R_{CS}$=1.568 Å, . Energies
are in eV, and for the products the energies are shifted 3.526 eV, so that its
zero of energy corresponds to the CS($r_{eq}$) at an infinite distance from H.
The contours correspond to 0 and $\pm$ 1 meV, to show the dependence of the
potential at long distances. The units of the color box is in eV.

The main features of the potential fit is shown in Fig. 4, where the contour
plots for the reactants (bottom panel, for $R_{CH}$=1.1199 Å) and products
(top panel, for $R_{CS}$=1.568 Å, and shifted by 3.526 eV, the exoergicity of
the reaction) are shown. The CH+S reactant channel (bottom panel) is
attractive for $R_{CS}<$ 7 Å, leading to the products channel at
$R_{CS}\approx$ 1.5 Å, with an energy of -5.5 eV. These energies correspond to
the CS-H well in products channel, which is about 2 eV below the CS + H
products, shown in the top panel of Fig. 4, shifted 3.526 eV to show the
details. At $R_{CH}\approx$ 2Å, there is a barrier which arises from the curve
crossing discussed above.

###  2.4 Quasi-Classical versus Quantum wave packet dynamics in the
$1^{2}A^{\prime}$ state

To check the validity of the quasi-classical method, we first compare quantum
and quasi-classical trajectory (QCT) calculations for total angular momentum
$J$=0. The quantum wave packet (WP) calculations are performed with the
MADWAVE3 code (Zanchet et al. 2009) and the parameters used are listed in
Table 1. The WP method is considered numerically exact, but it is very
demanding computationally. The QCT calculations are performed with the MDwQT
code (Sanz-Sanz et al. 2015; Zanchet et al. 2016; Ocaña et al. 2017). Initial
conditions are sampled with the usual Monte Carlo method (Karplus et al.
1965). In this first set of calculations CH is in its ground vibrational (v)
and rotational state, (v,j)=(0,0), and the initial internuclear distance and
velocity distributions are obtained with the adiabatic switching method
(Grozdanov & Solov’ev 1982; Qu & Bowman 2016; Nagy & Lendvay 2017). An impact
parameter $b=0$ is considered, that corresponds to $J$=0. The initial distance
between sulfur and the CH center of mass is set to 50 Bohr, and the
trajectories are stopped when any internuclear distance is longer than 60
Bohr. For each energy $N_{tot}=$ 104 trajectories are run to calculate the
reaction probability as $P_{R}(E)=N_{r}/N_{tot}$, where $N_{r}$ are the number
of reactive trajectories.

The WP and QCT reactions probabilities are compared in Fig. 5, and the two
show methods show a very similar behavior: a reaction probability slightly
larger than 0.9 at energies below 0.01 eV, decreasing as a function of
collision energy, down to a probability lower than 0.2 at 1 eV. In the two
cases there are oscillations, which do not match perfectly but that show
similar envelopes. Since these oscillations are expected to wash out when
considering the partial wave summation over total angular momentum, $J$, we
consider this agreement as satisfactory. This leads us to conclude, that QCT
method is sufficiently accurate to determine the reaction rate constant, as
described below.

Figure 5: CH+S $\rightarrow$ CS + H reaction probabilities versus collision energy for $J=0$ in the 1 ${}^{2}A^{\prime}$ electronic state using quantum wave packet (WP) and quasi-classical trajectory (QCT) methods.  Table 1:  Parameters used in the wave packet calculations in reactant Jacobi coordinates: $r_{min}\leq r\leq r_{max}$ is the CH internuclear distance, $R_{min}\leq R\leq R_{max}$ is the distance between CH center-of-mass and the sulphur atom, $0\leq\gamma\leq\pi$ is the angle between ${\@vec{r}}$ and ${\@vec{R}}$ vectors. The initial wave packet is described in $R$ by a Gaussian centered at $R=R_{0}$, and at a translational energy of $E=E_{0}$, and width $\Delta E$. The total reaction probability is obtained by analyzing the total flux at $r=r_{\infty}$.  $r_{min}$, $r_{max}=$ 0.1, 17 Å |  $N_{r}$=420
---|---
$r_{abs}$= 5 Å |
$R_{min}$, $R_{max}=$ 0.001, 17.5Å |  $N_{R}$=840
$R_{abs}$= 10.5 Å |
$N_{\gamma}$ = 168 | in $[0,\pi]$
$R_{0}$ = 10 Å |  $E_{0},\Delta E$= 0.35,0.2 eV
$r_{\infty}$ = 4 Å |

The reaction rate constant for CH($v$=0 and 1) in the 1${}^{2}A^{\prime}$
electronic state is evaluated according to

$K_{v}^{1^{2}A^{\prime}}(T)=\sqrt{\frac{8k_{B}T}{\pi\mu}}\pi
b_{max}^{2}(T)P_{r}(T),$ (8)

where $b_{max}(T)$ and $P_{r}(T)$ are the maximum impact parameter and
reaction probability at constant temperature, respectively. In this case,
about 105 trajectories are run for each temperature, the same for translation
and rotation degrees of freedom, fixing the vibrational state of CH to $v$=0
or 1.

###  2.5 Thermal rate

Considering that only the double degenerate 1${}^{2}A^{\prime}$ electronic
state reacts to form CS($X^{1}\Sigma^{+}$), the electronic partition function
has to be considered. Neglecting the spin-orbit splitting, the electronic
partition function would be 2/36, $i.e.$, the thermal rate constant is about
1/18 the reaction rate, $K^{1^{2}A^{\prime}}$, associated to the
1${}^{2}A^{\prime}$ state. Including the spin-orbit splitting of
S(${}^{3}P_{J_{S}}$) sulphur atom, and assuming that only the lowest two spin-
orbit react (having an individual rate constant equal to
$K^{1^{2}A^{\prime}}$), the vibrational selected thermal rate constant is
given by

$\displaystyle K_{v}(T)={2K_{v}^{1^{2}A^{\prime}}(T)\over
4\left[5+3exp(-539.83/T)+exp(-825.34/T)\right]},$ (9)

and it is shown in Fig. 6.

Figure 6: Vibrational-selected rate constant for the CH(X${}^{2}\Pi$,v=0)+
S(${}^{3}P$) $\rightarrow$ CS($X^{1}\Sigma$) + H reaction, obtained here
according to Eq. (9).

In this figure the red line corresponds to the rate constant (1.4 $\times$
$10^{-10}$ cm3s-1) obtained from the KIDA data base, as obtained with a
classical capture model (Vidal et al. 2017), using analytical formulas
(Georgievskii & Klippenstein 2005; Woon & Herbst 2009) with dipole moments and
polarizabilities taken from the literature or calculated using Density
Functional theory.

The $K_{v=0}^{1^{2}A^{\prime}}(T)/18$ rate constant at 10 K is about 4
$\times$ $10^{-11}$ cm3/s, increasing to a nearly constant value of $\approx$
5.5 $\times$ $10^{-11}$ cm3/s at temperatures above 100 K. This is a factor
between 2 and 3 lower than the value obtained with the capture model (Vidal et
al. 2017).

When including the spin-orbit splitting, the rate $K_{v=0}(T)$ is larger at 10
K, simply because the populations of the excited sulphur spin-orbit states,
$J_{S}$= 1 and 0, are negligible. As temperature increases, their populations
increase, leading to a reduction of the rate constant $K_{v=0}(T)$, which
decreases tending to $K^{1^{2}A^{\prime}}(T)/18$ at high temperature. The rate
with spin-orbit splittings, $K_{v=0}(T)$, is only a factor $\approx$ 1/2
smaller than that of KIDA at 10 K, and is considered the most accurate
obtained here. The $K_{v=0}(T)$ rate constant has been fitted and the
parameters are listed in Table 2.

In Fig. 6, $K_{v=1}^{1^{2}A^{\prime}}(T)/18$ is also shown. Its contribution
at low temperatures is small, because the vibration energy of CH(v=1) is about
0.337 eV higher than v=0, $i.e.$ 3916 K. Therefore, $K_{v=0}(T)$ is a good
approximation to the thermal rate constant, which includes the spin-orbit
splitting.

##  3 C2+S reaction

Figure 7: CASSCF/aV$T$Z optimized reaction paths for the conversion between
C2+S and CS+C which occur on several potential energy surfaces of CCS. The
reaction coordinate is defined as $R_{\mathrm{CS}}\\!-\\!R_{\mathrm{CC}}$,
while the valence angle, $\gamma_{\mathrm{CCS}}$, was held fixed at
$179.9$\mathrm{\SIUnitSymbolDegree}$$ in the geometry optimizations. Five
$A^{\prime\prime}$ and four $A^{\prime}$ electronic states of CCS were
considered for each multiplicity: triplet (bottom panel), singlet (middle
panel) and quintet (top panel). Figure 8: Optimized reaction paths for
$\mathrm{C_{2}}(X^{1}\Sigma_{g}^{+})\\!+\\!\mathrm{S}(^{3}P)\rightarrow\mathrm{CS}(X^{1}\Sigma^{+})\\!+\\!\mathrm{C}(^{3}P)$
as obtained from MRCI+Q/aV$Q$Z//CASSCF/aV$Q$Z calculations. The reaction
coordinate is defined as $R_{\mathrm{CS}}\\!-\\!R_{\mathrm{CC}}$. In the
bottom panel, the valence angle, $\gamma_{\mathrm{CCS}}$, was held fixed at
$179.9$\mathrm{\SIUnitSymbolDegree}$$ in the geometry optimizations, while, in
the top panel, it has been freely optimized along with $R_{\mathrm{CS}}$ or
$R_{\mathrm{CC}}$ bond distances. Two ${}^{3}A^{\prime\prime}$ and one
${}^{3}A^{\prime}$ triplet electronic states of CCS were considered in each
case.

The reaction

$\mathrm{C_{2}}(X^{1}\Sigma_{g}^{+})+\mathrm{S}(^{3}P)\rightarrow\mathrm{CS}(X^{1}\Sigma^{+})+\mathrm{C}(^{3}P)$
(10)

is exothermic by $\sim\\!1\,$\mathrm{eV}$$ (Visser et al. 2019; Reddy et al.
2003). Due to the open shell nature of the involved atoms [$\mathrm{S}(^{3}P)$
and $\mathrm{C}(^{3}P)$], reactants collision and product formation can take
place adiabatically on three triplet CCS potential energy surfaces
($2^{3}A^{\prime\prime}\\!+\\!1^{3}A^{\prime}$); see,  _e.g._ , Figures 7 and
8. To the best of our knowledge, no dedicated theoretical studies are yet
available on this reaction. Vidal et al. (2017) reported a theoretical upper
limit for its rate coefficient
($k\\!\sim\\!$2\text{\times}{10}^{-10}\text{\,}$\,\rm
cm^{3}\,molecule^{-1}\,s^{-1}$ at $10\,$\mathrm{K}$$) using classical capture
rate theory. From an experimental viewpoint, presently available techniques
for the production of reactant dicarbon often generate a mixture of both
$\mathrm{C_{2}}(X^{1}\Sigma_{g}^{+})$ and $\mathrm{C_{2}}(a^{3}\Pi_{u})$ (Gu
et al. 2006; Páramo et al. 2008) (recall that this electronically excited
state lies only $\sim\\!716\,\mathrm{cm^{-1}}$ above the ground $X$ form)
which, together with the expected high reactivity of C and S atoms, make the
laboratory characterization of this specific reaction (10) extremely
cumbersome.

###  3.1 Ab initio calculations

The methodology employed to obtain optimized energy paths for the
C2+S$\rightarrow$CS+C reaction closely resembles the one utilized for the CH+S
system. Preliminary PES explorations and geometry optimizations were all
performed at the SA-CASSCF level of theory, followed by single-point MRCI+Q
calculations. The CASSCF active space involves a total of 14 correlated
electrons in 12 active orbitals (9$a^{\prime}$+3$a^{\prime\prime}$). For each
multiplicity considered (triplet, singlet and quintet), five
$A^{\prime\prime}$ and four $A^{\prime}$ electronic states were simultaneously
treated in the SA-CASSCF wave functions. The aV$X$Z ($X\\!=\\!T,Q$) basis sets
of Dunning and co-workers (Dunning 1989; Kendall et al. 1992) were employed
throughout, with the calculations done with MOLPRO.

Our calculated CASSCF/aV$T$Z optimized path for reaction (10) is shown at the
bottom panel of Figure 7. As seen, the
$\mathrm{C_{2}}(X^{1}\Sigma_{g}^{+})+\mathrm{S}(^{3}P)$ reactants collision
involves only triplet CCS PESs and can happen on two ${}^{3}A^{\prime\prime}$
and one ${}^{3}A^{\prime}$ electronic states. Proceeding through the ground-
state PES of $\mathrm{CCS}(1^{3}A^{\prime\prime})$, reaction (10) does not
encounter any activation barriers for collinear atom-diatom approaches, being
exothermic by $\sim\\!0.6\,$\mathrm{eV}$$ at CASSCF/aV$T$Z level. Note that
this process occurs via the formation of a strongly-bound intermediate complex
corresponding to the linear global minimum of CCS,
$\ell$-$\mathrm{CCS}(X^{3}\Sigma^{-})$ (Saito et al. 1987); from this
structure, the ground-state $\mathrm{CS}(X^{1}\Sigma^{+})+\mathrm{C}(^{3}P)$
products can be directly accessed, without an exit barrier.

A close look at Figure 7 (bottom panel) also reveals that the excited
$2^{3}A^{\prime\prime}$ and $1^{3}A^{\prime}$ electronic states are degenerate
along $C_{\infty v}$ atom-diatom collisions. These PESs form the Renner-Teller
components of the strongly-bound $\ell$-$\mathrm{CCS}(A^{3}\Pi)$ complex,
showing a conical intersection with $\ell$-$\mathrm{CCS}(X^{3}\Sigma^{-})$ at
$R_{\mathrm{CS}}\\!-\\!R_{\mathrm{CC}}\\!\approx\\!+0.5\,{\AA}$ (Riaplov et
al. 2003; Tarroni et al. 2007). Differently from the ground
$1^{3}A^{\prime\prime}$ state, the conversion from reactants to products as
proceeding adiabatically through the $2^{3}A^{\prime\prime}$ and
$1^{3}A^{\prime}$ PESs entails a large activation barrier
($\approx\\!1\,$\mathrm{eV}$$ at CASSCF/aV$T$Z level) which is located at
$R_{\mathrm{CS}}\\!-\\!R_{\mathrm{CC}}\\!\approx\\!-0.5\,{\AA}$; see Figure 7.
As shown, this region of the nuclear configuration space is extremely
congested by the existence of several low-lying excited triplet states
correlating with $\mathrm{C_{2}}(a^{3}\Pi_{u})+\mathrm{S}(^{3}P)$. Note that
the $\mathrm{C_{2}}(a^{3}\Pi_{u})+\mathrm{S}(^{3}P)$ reactants can approach
each other in six triplet ($3^{3}A^{\prime\prime}\\!+\\!3^{3}A^{\prime}$), six
singlet ($3^{1}A^{\prime\prime}\\!+\\!3^{1}A^{\prime}$) and six quintet
($3^{5}A^{\prime\prime}\\!+\\!3^{5}A^{\prime}$) electronic states. For
completeness, their corresponding optimized reaction paths towards CS+C
formation are also plotted in Figure 7; see bottom, middle and top panels
therein. Accordingly, when proceeding adiabatically, the reactions involving
$\mathrm{C_{2}}(a^{3}\Pi_{u})+\mathrm{S}(^{3}P)$ are all endothermic, leading
ultimately to excited-state CS+C products. So, they are expected to be highly
inefficient at the low temperature regimes here envisaged (unless non-
adiabatic transitions play a role) and, for this reason, will not be
considered further in this work.

Keeping now our focus on the target
$\mathrm{C_{2}}(X^{1}\Sigma_{g}^{+})+\mathrm{S}(^{3}P)\rightarrow\mathrm{CS}(X^{1}\Sigma^{+})+\mathrm{C}(^{3}P)$
process [Reaction (10)] and to better estimate its overall attributes, we have
performed MRCI+Q/aV$Q$Z//CASSCF/aV$Q$Z calculations along the underlying
reaction paths on both $1^{3}A^{\prime\prime}$, $2^{3}A^{\prime\prime}$ and
$1^{3}A^{\prime}$ electronic states. The results are plotted in Figure 8.
Accordingly, at this level of theory, our best estimate for the exothermicity
of reaction (10) is $1.05\,$\mathrm{eV}$$ (without zero-point energy), a value
that matches nearly perfectly the corresponding experimental estimate of
$1.04\,$\mathrm{eV}$$ (Visser et al. 2019; Reddy et al. 2003). The
stabilization energies of the $\ell$-$\mathrm{CCS}(X^{3}\Sigma^{-})$ and
$\ell$-$\mathrm{CCS}(A^{3}\Pi)$ complexes are herein predicted to be $-5.3$
and $-4.2\,$\mathrm{eV}$$, respectively, relative to the
$\mathrm{C_{2}}(X^{1}\Sigma_{g}^{+})+\mathrm{S}(^{3}P)$ reactant channel. Most
notably, Figure 8 (bottom panel) shows that, at MRCI+Q/aV$Q$Z//CASSCF/aV$Q$Z
level, the predicted activation barriers along linear $C_{\infty v}$ paths for
the $2^{3}A^{\prime\prime}$ and $1^{3}A^{\prime}$ states are largely reduced
with respect to CASSCF/aV$T$Z values (Figure 7; bottom panel), going from
$\approx\\!1$ to less than $0.1\,$\mathrm{eV}$$. Indeed, by allowing the
valence C–C–S angle ($\gamma_{\mathrm{CCS}}$) to also be freely optimized in
the MRCI+Q/aV$Q$Z//CASSCF/aV$Q$Z calculations, Figure 8 (top panel)
unequivocally pinpoints that such barriers actually become submerged, thence
lying below the corresponding
$\mathrm{C_{2}}(X^{1}\Sigma_{g}^{+})+\mathrm{S}(^{3}P)$ reactant channel; note
therein the existence of small discontinuities on the _ab initio_ curves which
are associated with abrupt changes in $\gamma_{\mathrm{CCS}}$ near the top of
these barriers. This clearly indicates that all such PESs
($1^{3}A^{\prime\prime}$, $2^{3}A^{\prime\prime}$ and $1^{3}A^{\prime}$)
contribute to the overall dynamics/kinetics of reaction (10), even at low
temperatures. In the following, we describe the methodology employed to obtain
rate coefficients for this target reaction.

###  3.2 Long range interactions

Restricting the calculations to the equilibrium distance of
C2(${}^{1}\Sigma_{g}^{+}$), allow us to consider C2(${}^{1}\Sigma_{g}^{+}$)
and C2(${}^{3}\Pi$) separately. Under this approximation, the situation
simplifies to a closed shell diatom plus an open shell atom for the two
asymptotes: C2+S and CS+C. For both cases let us adopt the Jacobi coordinate
system and expansion in terms of orthogonal functions as defined eg. in Flower
& Launay (1977); Dubernet & Hutson (1994):

$V(R,\theta,\phi,\theta_{a},\phi_{a})=\sum_{\lambda\lambda_{a}\mu}V_{\lambda\lambda_{a}\mu}(R)C_{\lambda\mu}(\theta,\phi)C_{\lambda_{a},-\mu}(\theta_{a}\phi_{a})$
(11)

where
$C_{\lambda\mu}(\theta,\phi)=\left(\frac{4\pi}{2\lambda+1}\right)^{1/2}Y_{\lambda\mu}(\theta,\phi)$
are spherical harmonics in Racah normalization. The $\theta,\phi$ angles
correspond to the orientation of diatomic molecules with respect to vector
connecting the center-of mass of the diatomic molecule with the atom in a
Jacobi frame, while $\theta_{a},\phi_{a}$ angles to the orientation of doubly
occupied $p$ orbital of S and C atom, respectively. The $R$ is the distance
between the COM of diatom and atom. For both systems in $C_{s}$ symmetry the
solution of electronic Schrödinger equation is hard, since one cannot use
single-reference methods since two solutions will always belong to the same
irreducible representation, and are nearly degenerate. The situation is
simpler in case of symmetric configurations: C∞v (linear) - for both systems,
and C2v (T-shaped) for C2+S. For these cases atom+diatom states belong to
distinct irreducible representations: for linear geometry $\Pi$ and
$\Sigma^{-}$ states, while for C2v these symmetries are $B_{1}$, $B_{2}$ and
$A_{2}$. For these particular configurations of the CCS system one can use
single-reference gold-standard CCSD(T) method for calculations of interaction
energy in each symmetry (Atahan et al. 2006; Alexander 1998; Klos et al.
2004).

Using linear and T-shape geometries we have calculated UCCSD(T) interaction
energies in the range of 6-30 $\AA$ and converted it to
$V_{\lambda\lambda_{r}\mu}(R)$ potentials using the following formula for C2+S

$\displaystyle V_{000}$ $\displaystyle=$
$\displaystyle(2V_{\Pi}+V_{\Sigma}+2V_{B_{1}}+2V_{B_{2}}+2V_{A_{2}})/9$ (12)
$\displaystyle V_{020}$ $\displaystyle=$
$\displaystyle-2(V_{\Pi}-V_{\Sigma}+V_{B_{1}}+V_{B_{2}}-2V_{A_{2}})/9$ (13)
$\displaystyle V_{200}$ $\displaystyle=$ $\displaystyle
2(2V_{\Pi}+V_{\Sigma}-2V_{B_{1}}-V_{B_{2}}-V_{A_{2}})/9$ (14) $\displaystyle
V_{220}$ $\displaystyle=$
$\displaystyle-2(2V_{\Pi}-2V_{\Sigma}-V_{B_{1}}-V_{B_{2}}+2V_{A_{2}})/9$ (15)

and CS+C

$\displaystyle V_{000}$ $\displaystyle=$
$\displaystyle(2V_{\Pi,0}+V_{\Sigma,0}+2V_{\Pi,180}+V_{\Sigma,180})/6$ (16)
$\displaystyle V_{100}$ $\displaystyle=$
$\displaystyle(2V_{\Pi,0}+V_{\Sigma,0}-(2V_{\Pi,180}+V_{\Sigma,180}))/6$ (17)
$\displaystyle V_{020}$ $\displaystyle=$
$\displaystyle(V_{\Sigma,0}-V_{\Pi,0}+V_{\Sigma,180}-V_{\Pi,180})/3$ (18)
$\displaystyle V_{120}$ $\displaystyle=$
$\displaystyle(V_{\Sigma,0}-V_{\Pi,0}-(V_{\Sigma,180}-V_{\Pi,180}))/3$ (19)

where in the latter equation $V_{\Pi,0}$/$V_{\Sigma,0}$ denote potential for
appropriate symmetry in C-S-C configuration and $V_{\Pi,180}$/$V_{\Sigma,180}$
in C-C-S alignment. The above equations can be obtained by calculating Eq. 11
for $\theta,\phi,\theta_{a},\phi_{a}$ corresponding to symmetric
configurations (Flower & Launay 1977).

Then, the $V_{\lambda\lambda_{r}\mu}(R)$ potentials were carefully fitted to
analytical form to get inverse power expansion form. It is important to
realize, that collision-induced rotations of CS molecules are driven directly
by $V_{100}$ and $V_{120}$ terms, while for C2 molecules colliding with atoms
such terms are $V_{200}$ and $V_{220}$. When the atoms are assumed to be
spherically symmetric, the terms $V_{120}$ and $V_{220}$ can be ignored. Thus,
from now on we will skip the dependence on $\theta_{a},\phi_{a}$ in the Eq.
11. Thus our model of the potential used for the statistical method includes
only isotropic term $V_{000}$ and leading anisotropies $V_{100}$ and $V_{200}$
which were fitted to analytical forms of van der Waals expansion:
$\sum_{i=0}^{3}C_{6+i}R^{-(6+i)}$ for $V_{000}$ and $V_{200}$, and
$\sum_{i=0}^{3}C_{7+i}R^{-(7+i)}$ for $V_{100}$. These potentials can be
viewed as averaged over all orientations of $P$-state atoms. Moreover, for the
leading coefficients we also performed similar calculations using open-shell
Symmetry Adapted Perturbation theory (Hapka et al. 2012) and confirmed the
values for the leading coefficients $C_{6}$ and $C_{7}$ of the reactants and
products (they agreed to within 10%).

The final analytical form of the potential used for C2+S reads (in atomic
units of distance and energy)

$\displaystyle V(R,\theta)$ $\displaystyle=$ $\displaystyle
A\,R^{-6}+B\,R^{-8}+C\,R^{-10}$ $\displaystyle+$
$\displaystyle\frac{1}{2}(3\cos^{2}(\theta)-1)(D\,R^{-6}+E\,R^{-8}+F\,R^{-10})$

with A=$-125.8$ , B=$-9.444\times 10^{3}$,C=$1.743\times 10^{5}$, D=$-13.66$ ,
E=$-2.848\times 10^{3}$ and F=$-1.929\times 10^{5}$ while for CS+C

$\displaystyle V(R,\theta)$ $\displaystyle=$ $\displaystyle
A\,R^{-6}+B\,R^{-8}+C\,R^{-10}$ $\displaystyle+$ $\displaystyle
D\,R^{-12}+{\cos(\theta)\over R^{-7}}\,(E\,+F\,R^{-2}+G\,R^{-4}+H\,R^{-6})$

with A= $-147.5$, B=$-207.3$, C=$-3.834\times 10^{6}$, D=$-2.805\times
10^{8}$, E=$-487.3$, F=$-4.309\times 10^{3}$, G=$-2.394\times 10^{7}$ and
H=$-1.997\times 10^{9}$.

Figure 9: Reactive rate constant for the C2+S $\rightarrow$ CS + S collisions obtained with the present AS method and compared to that available in the KIDA data base from Vidal et al. (2017) Table 2: Parameters used to fit the total reaction rates calculated for CH+S and C2 \+ S, according to the expression $K(T)=A(T/300)^{B}e^{-C/T}$ reaction | A (cm${}^{3}/s$) | B | C (K)
---|---|---|---
CH + S | 7.52 $\times$ 10-11 | -0.117 | 5.42
C2 \+ S | 5.08 $\times$ 10-10 | 0.121 | 2.70

###  3.3 Rate constant calculations

Due to the deep well appearing for the three lowest adiabatic states of C2\+ S
reaction, and the large masses of the three atoms involved, this system has a
high density of resonances near the thresholds. Therefore, this reaction is
expected to be governed by a statistical mechanism. In this work we shall use
the Adiabatic Statistical (AS) method (Quack & Troe 1974), using the AZTICC
code recently implemented (Gómez-Carrasco et al. 2022). Here, we shall use the
rigid rotor approach, similar to the method already used to treat several
reactions and inelastic processes (Konings et al. 2021). We consider the
experimental exoergic values of D0=D${}^{C_{2}}_{0}$ \- D${}^{CS}_{0}$ = 1.14
eV (Huber & Herzberg 1979), which already includes the vibrational zero-point
energy (ZPE) of reactants and products. This exothermicity is very close to
that shown in the ab initio calculations, of $\approx$ 1\. eV in Fig. 8, which
does not include ZPE.

In the rigid rotor approach, the diatomic distances for the diatomic molecules
are frozen to their equilibrium values, $r_{e}$= 1.2425 and 1.568 Å for C2 and
CS, respectively. Since C2 is homonuclear, the even and odd rotational
channels are not coupled, and are treated independently. The large
exothermicity requires the inclusion of vibrational levels for the CS
products, up to v=15, to describe adequately the products density of states.
To this end, the adiabatic potentials are generated for each final $v$,
independently, using the rigid rotor approximation at $r_{eq}$ with the same
long-range potential but shifting the energy using the anharmonic vibration
constants, $\omega_{e}=$ 0.15933 eV $\omega_{e}\chi_{e}=$ 0.801 meV from
Herzberg (1950).

Table 3: Model predictions using the old and new reaction rate coefficients Molecule | Age (Myr) | TMC1-old | TMC1-new | HC-old | HC-new | PDR-old | PDR-new
---|---|---|---|---|---|---|---
CS | 0.1 | 1.92$\times 10^{-8}$ | 1.61$\times 10^{-8}$ | 6.33$\times 10^{-11}$ | 6.13$\times 10^{-11}$ | 2.69$\times 10^{-9}$ | 2.88$\times 10^{-9}$
| 1.0 | 6.47$\times 10^{-8}$ | 6.47$\times 10^{-8}$ | 2.04$\times 10^{-12}$ | 1.66$\times 10^{-12}$ | 2.89$\times 10^{-9}$ | 3.06$\times 10^{-9}$
SO | 0.1 | 1.08$\times 10^{-7}$ | 4.70$\times 10^{-8}$ | 6.36$\times 10^{-7}$ | 6.28$\times 10^{-7}$ | 1.70$\times 10^{-9}$ | 1.70$\times 10^{-9}$
| 1.0 | 1.70$\times 10^{-9}$ | 1.67$\times 10^{-9}$ | 4.28$\times 10^{-7}$ | 3.50$\times 10^{-7}$ | 1.82$\times 10^{-9}$ | 1.82$\times 10^{-9}$
SO2 | 0.1 | 7.68$\times 10^{-9}$ | 1.41$\times 10^{-8}$ | 1.32$\times 10^{-7}$ | 1.56$\times 10^{-7}$ | 9.99$\times 10^{-14}$ | 1.33$\times 10^{-13}$
| 1.0 | 4.28$\times 10^{-11}$ | 2.20$\times 10^{-10}$ | 3.71$\times 10^{-7}$ | 4.49$\times 10^{-7}$ | 1.06$\times 10^{-13}$ | 1.43$\times 10^{-13}$

All product ro-vibrational levels are then considered to obtain the square of
the S-matrix at each total angular momentum $J$ (see Ref.Gómez-Carrasco et al.
(2022) for more details). The calculations are done using a body-fixed frame,
in which the z-axis is parallel to the Jacobi vector R, joining the diatomic
center-of-mass to the atom, and the three atoms are considered in the body-
fixed xz-plane. In the present calculations, a maximum rotational quantum
number $j_{max}=$ 200 and 250 have been considered for C2 and CS, and a
maximum helicity quantum number of $\Omega_{max}$= 15. The individual state-
to-state reactive cross sections are obtained performing the summation over
all $J$ in the partial wave expression up to $J_{max}=$ 200.

Finally, integrating over the translational energy according to a Boltzmann
energy distribution, and summing over all accessible states, the thermal rate
constants, $K_{\alpha}(T)$, are obtained, with $\alpha$=
$1^{3}A^{\prime\prime}$, $1^{3}A^{\prime}$and $1^{3}A^{\prime}$. The
$2^{3}A^{\prime\prime}$ and $1^{3}A^{\prime}$ state present a submerged
barrier which could reduce the reactivity, but here we consider that the three
$\alpha$ electronic states have the same reactive rate constant. For this
reason, when the spin-orbit splitting of S(${}^{3}P_{J}$) is included, the
final thermal rate constant averaging over the spin-orbit states, $K^{AS}(T)$,
is equal to $K_{\alpha}(T)$, shown in Fig. 9, whose parameters are fitted and
listed in Table 2.

The C2\+ S $\rightarrow$ CS + S available in KIDA data
base222https://kida.astrochem-tools.org/ is that of Vidal et al. (2017)
obtained by a capture method, which corresponds to a value of $K^{C}(T)$= 2
10-10 cm3/s, independent of temperature, and was calculated using a capture
method, $i.e.$ assuming that the reaction is exothermic, but without
considering the deep well described in this work. Such approach neglects the
possibility that the C2S complex formed under the statistical assumption can
exit back to C2 \+ S products. This probability is small because it is
proportional to the density of states in each channel, and this density is
much lower in the C2+S due to the exothermicity, the homonuclear symmetry and
the larger rotational constant.

##  4 Discussion

In this section we evaluate the impact of the new reaction rates in our
understanding of interstellar chemistry. This is not straightforward since the
formation and destruction routes of the different species depend on the local
physical and chemical conditions as well as the chemical time. Therefore, we
need to consider different environments to be able to have a comprehensive
view of the impact of the new reactions rates on astrochemical calculations,
and in particular in our ability to reproduce the abundance of CS.

We performed chemical calculations using Nautilus 1.1 (Ruaud et al. 2016), a
three-phase model, in which gas, grain surface and grain mantle phases, and
their interactions, are considered. We used the code upgraded as described by
Wakelam et al. (2021) with the chemical network of KIDA2 that has been
modified to account for the reaction rates estimated by Fuente et al. (2016),
Fuente et al. (2019), and this paper (reaction rates in Table 2). In Nautilus,
desorption into the gas phase is only allowed for the surface species,
considering both thermal and non-thermal mechanisms. In the regions where the
temperature of grain particles is below the sublimation temperature, non-
thermal desorption processes become important to calculate the number of
molecules in gas phase. The latter include desorption induced by cosmic rays
(Hasegawa & Herbst 1993), direct (UV field) and indirect (secondary UV field
induced by the cosmic-ray flux) photo-desorption, and reactive chemical
desorption (Garrod et al. 2007; Minissale et al. 2016). In the following
calculations, we use the prescription proposed by Minissale et al. (2016) for
ice coated grains to calculate the reactive chemical desorption. The physical
and chemical conditions associated with these three simulations are detailed
below:

Table 4: Chemical abundances in TMC 1 (CP)

Species | Observed | Model (0.1 Myr) | Model (1 Myr)
---|---|---|---
H2 | 0.5 |  |  |
CS | 6.5$\times 10^{-9}$(1) | 6.4$\times 10^{-9}$(2) | 1.6$\times 10^{-8}$ | 6.5$\times 10^{-8}$
CO | 4.8$\times 10^{-5}$(1) | 4.2$\times 10^{-5}$(2) | 1.0$\times 10^{-4}$ | 6.8$\times 10^{-6}$
HCO+ | 5.0$\times 10^{-9}$(1) | 8.4$\times 10^{-9}$(2) | 1.0$\times 10^{-8}$ | 2.0$\times 10^{-9}$
HCS+ | 1.1$\times 10^{-10}$(∗) | 7.0$\times 10^{-11}$(1) | 2.3$\times 10^{-11}$ | 1.4$\times 10^{-10}$
HCN | 4.8$\times 10^{-8}$(∗) | 3.5$\times 10^{-8}$(2) | 5.0$\times 10^{-8}$ | 2.1$\times 10^{-7}$
SO | 9.0$\times 10^{-10}$(1) | 5.3$\times 10^{-10}$(2) | 4.7$\times 10^{-8}$ | 1.7$\times 10^{-9}$
OCS | 1.2$\times 10^{-10}$(2) |  | 3.2$\times 10^{-9}$ | 4.5$\times 10^{-9}$
H2S | 8.0$\times 10^{-10}$(2) | 1.1$\times 10^{-9}$(3) | 1.4$\times 10^{-9}$ | 1.3$\times 10^{-8}$
H2CS | 1.1$\times 10^{-9}$(4) |  | 2.1$\times 10^{-9}$ | 6.7$\times 10^{-8}$
C2S | 2.8$\times 10^{-9}$(4) |  | 2.6$\times 10^{-10}$ | 9.4$\times 10^{-9}$
C3S | 3.8$\times 10^{-10}$(4) |  | 3.9$\times 10^{-10}$ | 3.3$\times 10^{-9}$
HCSCN | 8.9$\times 10^{-12}$(5) |  |  |
HCSCCH | 2.2$\times 10^{-12}$(5) |  |  |
NCS | 2.2$\times 10^{-11}$(5) |  |  |
HCCS | 1.9$\times 10^{-11}$(5) |  | 1.0$\times 10^{-11}$ | 1.1$\times 10^{-8}$
H2CCS | 2.2$\times 10^{-11}$(5) |  | 2.1$\times 10^{-11}$ | 4.7$\times 10^{-11}$
H2CCCS | 1.0$\times 10^{-11}$(5) |  |  |
C4S | 1.0$\times 10^{-12}$(5) |  | 6.2$\times 10^{-12}$ | 3.1$\times 10^{-10}$
C5S | 1.4$\times 10^{-12}$(5) |  |  |
HC3S+ | 5.5$\times 10^{-12}$(5) |  | 9.3$\times 10^{-12}$ | 6.5$\times 10^{-11}$
HSCN | 1.9$\times 10^{-12}$(6) |  | 6.3$\times 10^{-13}$ | 1.3$\times 10^{-10}$
HNCS | 1.9$\times 10^{-12}$(6) |  | 8.8$\times 10^{-13}$ | 2.0$\times 10^{-10}$
HCCCCS | 2.6$\times 10^{-12}$(7) |  |  |
HCCS+ | 3.0$\times 10^{-11}$(8) |  | 1.4$\times 10^{-11}$ | 1.3$\times 10^{-10}$
NS+ | 1.4$\times 10^{-12}$(9) |  | 1.7$\times 10^{-11}$ | 1.5$\times 10^{-12}$
NS | 4.7$\times 10^{-11}$(9) |  | 2.2$\times 10^{-10}$ | 1.0$\times 10^{-9}$
HCS | 1.5$\times 10^{-10}$(5) |  | 2.1$\times 10^{-11}$ | 3.1$\times 10^{-10}$

References: (1) Fuente et al. (2019) 2) Rodríguez-Baras et al. (2021) 3)
Navarro-Almaida et al. (2020) 4) Gratier et al. (2016) 5) Cernicharo et al.
(2021a) 6) Adande et al. (2010) 7) Fuentetaja et al. (2022) 8) Cernicharo et
al. (2021b) 9) Cernicharo et al. (2018)

(*) Calculated in this work.

  * •

TMC 1: This case represents the physical and chemical conditions prevailing in
molecular cloud complexes where low-mass stars are formed. The physical
properties of these regions are typically described with moderate number
densities of atomic hydrogen nuclei n${}_{\rm H}=3\times 10^{4}$ cm-3, cold
gas and dust temperatures ${\rm T}=10$ K, a moderate visual extinction ${\rm
A}_{V}=20$ mag, a cosmic-ray H2 ionization rate of $\zeta_{{\rm
H}_{2}}=10^{-16}$ s-1 (Fuente et al. 2019, 2023), and an intensity of the far-
ultraviolet field (FUV) equal to $\chi=5$ in Draine units. We assume that
sulphur is depleted by a factor of 20 relative to cosmic abundance as derived
by Fuente et al. (2019) and Fuente et al. (2023) in Taurus and Perseus.

  * •

Hot core: This case represents the physical and chemical conditions in the
warm interior of young protostars where the gas and dust temperature is $>$100
K and the icy grains mantles are sublimated. Typical physical conditions in
these regions are: nH = 3 $\times$ 106 cm-3, Tk= 200 K, AV= 20 mag,
$\zeta_{\rm H_{2}}$= 1.3 $\times$ 10-17 s-1, $\chi$=1 in Draine units. Sulphur
depletion in hot cores is not well established. We adopted [S/H]=8$\times$10-7
since this value is commonly used to model massive hot cores (Gerner et al.
2014).

  * •

Photon-dominated regions (PDRs) are those environments where the far-
ultraviolet photons emitted by hot stars are determining the physical and
chemical conditions of gas and dust. PDRs can be found on the surfaces of
protoplanetary disks and molecular clouds, globules, planetary nebulae, and
starburst galaxies. As representative of the physical and chemical conditions
in the PDRs associated with massive star forming regions, we select: nH = 5
$\times$ 105 cm-3, Tk= 100 K, AV= 4 mag, $\zeta_{H_{2}}$= 10-16 s-1,
$\chi$=104 in Draine units. The amount of sulphur in gas phase in PDRs is
still an open question. Based on observations of sulphur recombination lines
in the Orion Bar, Goicoechea & Cuadrado (2021) obtained that the abundance of
sulphur should be close to the solar value in this prototypical PDR. We are
aware that sulphur recombination lines are arising from PDR layers close to
the S+/S transition, at $A_{v}\sim$ 4 mag (Simon et al. 1997; Goicoechea et
al. 2006)) and sulphur depletion might be higher towards more shielded regions
from where the emission of most sulphur-bearing molecules comes. In spite of
this, we adopt [S/H]=1.5$\times$10-5 for our calculations.

Figure 10: Comparison between Nautilus predictions obtained using our updated
chemical network and the abundances observed towards TMC 1 (CP) (see Table 4).
Grey regions indicate the observed values. We assumed an uncertainty of a
factor of 2 in the measured abundances. The physical and chemical conditions
used in our calculations are: nH = 3 $\times$ 104 cm-3, Tk= 10 K, AV= 20 mag,
$\zeta_{H_{2}}$= 10-16 s-1, $G_{0}$=1 in Draine units and
[S/H]=7.5$\times$10-7 (sulphur depleted by a factor of $\sim$20).

Table 3 shows the fractional abundances of CS, SO, and SO2 predicted using the
”old” and ”new” chemical network and the physical conditions above described.
In regions where the ionization fraction is large, CS is essentially produced
from the electronic dissociative recombination of HCS+, where HCS+ is formed
by reactions of S+ and CH (Sternberg & Dalgarno 1995; Lucas & Liszt 2002).
Vidal et al. (2017) showed that the production of HCS+ with S+ \+ CH
$\rightarrow$ CS+ \+ H, and CS+ \+ H2 $\rightarrow$ HCS+ \+ H, is the most
efficient HCS+ formation route at the cloud surface. In more shielded regions
where the sulphur is mainly in neutral atomic form, CS is also produced by
neutral-neutral reactions, with significant contributions of the reactions
studied in this paper. For dark cloud conditions, the reaction C2 \+ S forms
CS more efficiently than CH + S. At low temperatures, the calculated C2 \+ S
rate is slightly higher than previous value (see Fig. 9). However, its
possible effect is canceled by the lower value of the new CH + S reaction rate
(see Fig. 6). As a result, the impact of our new rates on the CS abundance for
the TMC1 case is negligible. Not surprisingly, the major impact of the new
reactions rates calculated on the CS abundance is observed for the hot core
case with variations of $\sim$ 20% due to the significantly higher C2 \+ S
reaction rate at temperatures $>$100 K (see Fig. 9). The S + O2 $\rightarrow$
SO + O and SO + OH $\rightarrow$ SO2 \+ H reactions rates published by Fuente
et al. (2016) and Fuente et al. (2019) produce the maximum variations in the
case TMC1, with variations of the SO and SO2 abundances of a factor of $>$2,
which demonstrates the need of performing this type of calculations.

TMC 1 (CP) is the astrophysical object for which the highest number of
sulphur-bearing species have been detected so far, with more than ten complex
sulphur-bearing molecules detected for the first time in the last 3 years
(Cernicharo et al. 2021a; Fuentetaja et al. 2022, Table 4). The large number
of atoms in these new species ($>$5) shows that a rich and complex organo-
sulphur chemistry is going on in this dark cloud(Laas & Caselli 2019).
Although these large molecules carry a small percentage of the sulphur budget,
their detection is useful to test the predictive power of our chemical
network. Table 4 shows a compilation of the observed abundances towards TMC 1
(CP). In order to perform the most uniform and reliable comparison, the
abundances have been re-calculated assuming NH=3.6$\times$1022 cm-2 (Fuente et
al. 2019). The abundances of HCS+ and H13CN were already estimated by
Rodríguez-Baras et al. (2021). Here, they have been re-calculated using the
most recent collisional coefficients reported by Denis-Alpizar et al. (2022)
and Navarro-Almaida et al. (2023). Following the same methodology explained in
Fuente et al. (2019) and Rodríguez-Baras et al. (2021), we assumed
X(HCN)/X(H13CN)=60 to estimate the HCN abundance.

We performed chemical calculations using Nautilus 1.1 and the chemical network
modified to account for the reaction rates presented in this paper, Fuente et
al. (2019), and Fuente et al. (2016), and the physical conditions derived by
Fuente et al. (2023).These physical conditions are the same as in the TMC1
case in Table 3. Only 17 of the total number of sulphur-bearing species
detected towards TMC 1 (CP) are included in our chemical network. Fig. 10
shows the chemical predictions for all the sulphur-bearing species that have
been observed in this proto-typical source and are included in our chemical
network. In addition, we show the CO, HCO+, and HCN abundances because these
molecules are considered good tracers of the gas ionization degree and C/O
ratio (Fuente et al. 2019).

We find a reasonable agreement (within a factor of 10) between model and
observations for CO, HCO+, HCN, CS, HCS+, H2S, SO, C2S, C3S, C4S, H2CS, HC2S+,
HC3S+, H2CCS, NS+, HNCS, and HSCN for times between 0.1 Myr and 1 Myr.
However, as already commented by Bulut et al. (2021) and Wakelam et al.
(2021), the chemical time at which we find the best solutions depends on the
considered species. The most important restrictions respect to the chemical
time comes for CO and HCO+ whose abundances rapidly decrease for times later
than 0.4 Myr. On the opposite side, we find that the abundances of HCS+, SO,
and C2S are better reproduced for times $>$ 1 Myr. Only OCS, C4S, and NS
cannot be fitted with chemical times between 0.1 Myr and 1 Myr. We would like
to recall that the chemical time is not the same as the dynamical time, since
various physical phenomena such as turbulent motions that carry molecules to
the cloud surfaces or shocks can reset the chemical age of the gas. sulphur-
bearing species are very sensitive to the chemical time with several species
whose abundances vary in several orders of magnitude from 0.1 Myr to 1 Myr.
Therefore the chemical time is a critical parameter to fit them (see Fig.10).
In our 0D chemical calculations, the physical conditions remain fixed, which
is far from the real case of a collapsing and fragmenting cloud. In
forthcoming papers, we will explore the influence that the cloud dynamical
evolution would have on the sulphur chemistry.

##  5 Summary and conclusions

The rate constants for the formation of CS($X^{1}\Sigma^{+}$), CH$(^{2}\Pi)$
\+ S(${}^{3}P$) and C2($X^{1}\Sigma_{g}^{+}$)+ S(${}^{3}P$) have been obtained
in this work, and are tabulated in Table 2. These two reactions involve open
shell reactants, and therefore present several degenerate, or nearly
degenerate, electronic states. The role of each initial electronic state in
the formation of CS has been analyzed in detail.

For CH$(^{2}\Pi)$ \+ S(${}^{3}P$) it is found that only the
1${}^{2}A^{\prime}$ can contribute to the CS($X^{1}\Sigma^{+}$) formation
through an exothermic barrierless mechanism, $i.e.$ 2 states of the doublet
among the 36 degenerate electronic states correlating to the CH$(^{2}\Pi)$ \+
S(${}^{3}P$) asymptote. When the spin-orbit splitting is taken into account,
at the low temperatures of 10 K the electronic partition function becomes 2/9,
tending to 2/36 at high temperature. Surprisingly, the rate constant obtained
with a capture model (Vidal et al. 2017) is only a factor of two higher than
the present result at 10 K, and this difference increases very little with
increasing temperature up to 500K.

For the C2($X^{1}\Sigma_{g}^{+}$)+ S(${}^{3}P$) reaction the three triply
degenerate states connect to the CS($X^{1}\Sigma^{+}$) products. It is found
that the three states present a deep insertion well, with depths between 5.5
and 4 eV. The ground electronic state proceeds with no barrier, while the two
excited states have a barrier in the products channel, which becomes submerged
for bent configuration. The presence of the deep insertion well justifies the
use of an adiabatic statistical method to calculate the reactive rate
constant, which in turn is very similar to that obtained with a capture model
(Vidal et al. 2017) at 10 K. The difference increases with temperature, and
the present results become 2.5 times larger than the constant value obtained 2
$\times$10-10 cm3/s at 500 K.

The present results corroborate those obtained with classical capture models
at low temperatures, differing by a factor of 2.5 at most. It should be noted,
however, that this is due to different reasons for the two reactions studied
here. For open shell reactants, as those treated here, for which experiments
are difficult, it is required to do a detailed analysis of the reactivity of
all initial degenerate electronic states of the reactants before generalizing
these findings.

The new rates have been implemented in a chemical network to compare with the
observations of sulphur-bearing species towards TMC 1(CP). Model predictions
are in reasonable agreement with the observation for most the sulphur-bearing
species, except for OCS and NS, which cannot be fitted with our model.
However, it is not possible to fit all of them with an unique chemical time,
which suggests that dynamical effects are important for sulphur chemistry.

###### Acknowledgements.

The research leading to these results has received funding from MICIN (Spain)
under grant PID2021-122549NB-C21 and from the European Union’s Horizon 2020
research and innovation program under the Marie Sklodowska-Curie grant
agreement No 894321. N.B is grateful for support from the Polish National
Agency for Academic Exchange (NAWA) Grant and also acknowledges TUBITAK’s
2219-Program by scholarship no. 1059B192200348. P.S.Z. is grateful to National
Science Centre of Poland for funding the project No 2019/34/E/ST4/00407. AF
and PRM are grateful to Spanish MICIN for funding under grant
PID2019-106235GB-I00. JRG thanks the Spanish MCINN for funding support under
grant PID2019-106110GB-I00. JEP was supported by the Max-Planck Society. DNA
acknowledges funding support from Fundación Ramón Areces through its
international postdoc grant program. RLG would like to thank the ”Physique
Chimie du Milieu Interstellaire” (PCMI) programs of CNRS/INSU for their
financial supports.

## References

  * Adande et al. (2010) Adande, G. R., Halfen, D. T., Ziurys, L. M., Quan, D., & Herbst, E. 2010, ApJ, 725, 561
  * Aguado & Paniagua (1992) Aguado, A. & Paniagua, M. 1992, J. Chem. Phys., 96, 1265
  * Aguado et al. (1998) Aguado, A., Tablero, C., & Paniagua, M. 1998, Comput. Phys. Commun., 108, 259
  * Agúndez & Wakelam (2013) Agúndez, M. & Wakelam, V. 2013, Chemical Reviews, 113, 8710
  * Alexander (1998) Alexander, M. H. 1998, The Journal of chemical physics, 108, 4467
  * Atahan et al. (2006) Atahan, S., Kłos, J., Zuchowski, P. S., & Alexander, M. H. 2006, Physical Chemistry Chemical Physics, 8, 4420
  * Boogert et al. (1997) Boogert, A. C. A., Schutte, W. A., Helmich, F. P., Tielens, A. G. G. M., & Wooden, D. H. 1997, A&A, 317, 929
  * Bulut et al. (2021) Bulut, N., Roncero, O., Aguado, A., et al. 2021, Astron. AstroPhys., 646, A5
  * Cazaux et al. (2022) Cazaux, S., Carrascosa, H., Muñoz Caro, G. M., et al. 2022, A&A, 657, A100
  * Cernicharo et al. (2021a) Cernicharo, J., Cabezas, C., Endo, Y., et al. 2021a, A&A, 650, L14
  * Cernicharo et al. (2021b) Cernicharo, J., Cabezas, C., Endo, Y., et al. 2021b, A&A, 646, L3
  * Cernicharo et al. (2018) Cernicharo, J., Lefloch, B., Agúndez, M., et al. 2018, ApJ, 853, L22
  * Davidson (1975) Davidson, E. R. 1975, J. Comp. Phys., 17, 87
  * Denis-Alpizar et al. (2022) Denis-Alpizar, O., Quintas-Sánchez, E., & Dawes, R. 2022, MNRAS, 512, 5546
  * Dubernet & Hutson (1994) Dubernet, M. L. & Hutson, J. 1994, J. Chem. Phys., 101, 1939
  * Dunning (1989) Dunning, T. H. 1989, J. Chem. Phys., 90, 1007
  * Dunning & Jr. (1989) Dunning, T. H. & Jr. 1989, J. Chem. Phys., 90, 1007
  * Esplugues et al. (2022) Esplugues, G., Fuente, A., Navarro-Almaida, D., et al. 2022, A&A, 662, A52
  * Ferrante et al. (2008) Ferrante, R. F., Moore, M. H., Spiliotis, M. M., & Hudson, R. L. 2008, ApJ, 684, 1210
  * Flower & Launay (1977) Flower, D. R. & Launay, J. M. 1977, Journal of Physics B: Atomic and Molecular Physics, 10, 3673
  * Fuente et al. (2016) Fuente, A., Cernicharo, J., Roueff, E., et al. 2016, A&A, 593, A94
  * Fuente et al. (2019) Fuente, A., Navarro, D. G., Caselli, P., et al. 2019, A&A, 624, A105
  * Fuente et al. (2023) Fuente, A., Rivière-Marichalar, P., Beitia-Antero, L., et al. 2023, A&A, 670, A114
  * Fuentetaja et al. (2022) Fuentetaja, R., Agúndez, M., Cabezas, C., et al. 2022, A&A, 667, L4
  * Garrod et al. (2007) Garrod, R. T., Wakelam, V., & Herbst, E. 2007, A&A, 467, 1103
  * Georgievskii & Klippenstein (2005) Georgievskii, Y. & Klippenstein, S. J. 2005, J.Chem. Phys., 122, 194103
  * Gerner et al. (2014) Gerner, T., Beuther, H., Semenov, D., et al. 2014, A&A, 563, A97
  * Goicoechea et al. (2021) Goicoechea, J. R., Aguado, A., Cuadrado, S., et al. 2021, AA, 647, A10
  * Goicoechea & Cuadrado (2021) Goicoechea, J. R. & Cuadrado, S. 2021, A&A, 647, L7
  * Goicoechea et al. (2006) Goicoechea, J. R., Pety, J., Gerin, M., et al. 2006, A&A, 456, 565
  * Gómez-Carrasco et al. (2022) Gómez-Carrasco, S., Félix-González, D., Aguado, A., & Roncero, O. 2022, J. Chem. Phys., 157, 084301
  * Gratier et al. (2016) Gratier, P., Majumdar, L., Ohishi, M., et al. 2016, ApJS, 225, 25
  * Grozdanov & Solov’ev (1982) Grozdanov, T. P. & Solov’ev, E. A. 1982, J. Phys. B, 15, 1195
  * Gu et al. (2006) Gu, X., Guo, Y., Zhang, F., Mebel, A. M., & Kaiser, R. I. 2006, Faraday Discuss., 245
  * Hapka et al. (2012) Hapka, M., Żuchowski, P. S., Szczȩśniak, M. M., & Chałasiński, G. 2012, J. Chem. Phys., 137, 164104
  * Hasegawa & Herbst (1993) Hasegawa, T. I. & Herbst, E. 1993, MNRAS, 261, 83
  * Herzberg (1950) Herzberg, G. 1950, Molecular spectra and molecular structure. I. Spectra of diatomic molecules (van Nostrand Reinhold Co. (New York))
  * Hily-Blant et al. (2022) Hily-Blant, P., des Forêts, G. P., Faure, A., & Lique, F. 2022, A&A, 658, A168
  * Huber & Herzberg (1979) Huber, K. P. & Herzberg, G. 1979, Molecular Spectra and Molecular Structure. Vol IV. Constants of Diatomic Molecules (Van Nostrand, Toronto)
  * Jiménez-Escobar & Muñoz Caro (2011) Jiménez-Escobar, A. & Muñoz Caro, G. M. 2011, A&A, 536, A91
  * Jiménez-Escobar et al. (2014) Jiménez-Escobar, A., Muñoz Caro, G. M., & Chen, Y. J. 2014, MNRAS, 443, 343
  * Karplus et al. (1965) Karplus, M., Porter, R. N., & Sharma, R. D. 1965, J. Chem. Phys., 43, 3259
  * Kendall et al. (1992) Kendall, R. A., Dunning, T. H., & Harrison, R. J. 1992, J. Chem. Phys., 96, 6796
  * Klos et al. (2004) Klos, J., Szczesniak, M. M., & Chalasinski *, G. 2004, International Reviews in Physical Chemistry, 23, 541
  * Konings et al. (2021) Konings, M., Desrousseaux, B., Lique, F., & Loreau, J. 2021, J. Chem. Phys., 155, 104302
  * Laas & Caselli (2019) Laas, J. C. & Caselli, P. 2019, A&A, 624, A108
  * Lucas & Liszt (2002) Lucas, R. & Liszt, H. S. 2002, A&A, 384, 1054
  * McClure et al. (2023) McClure, M. K., Rocha, W. R. M., Pontoppidan, K. M., et al. 2023, Nature Astronomy, 7, 431
  * Minissale et al. (2016) Minissale, M., Dulieu, F., Cazaux, S., & Hocuk, S. 2016, A&A, 585, A24
  * Nagy & Lendvay (2017) Nagy, T. & Lendvay, G. 2017, J. Phys. Chem. Lett., 8, 4621
  * Navarro-Almaida et al. (2023) Navarro-Almaida, D., Bop, C. T., Lique, F., et al. 2023, A&A, 670, A110
  * Navarro-Almaida et al. (2020) Navarro-Almaida, D., Le Gal, R., Fuente, A., et al. 2020, A&A, 637, A39
  * Ocaña et al. (2017) Ocaña, A. J., Jiménez, E., Ballesteros, B., et al. 2017, AstroPhys. J., 850, 28
  * Palumbo et al. (1997) Palumbo, M. E., Geballe, T. R., & Tielens, A. G. G. M. 1997, ApJ, 479, 839
  * Palumbo et al. (1995) Palumbo, M. E., Tielens, A. G. G. M., & Tokunaga, A. T. 1995, ApJ, 449, 674
  * Páramo et al. (2008) Páramo, A., Canosa, A., Le Picard, S. D., & Sims, I. R. 2008, J. Phys. Chem. A, 112, 9591
  * Qu & Bowman (2016) Qu, C. & Bowman, J. M. 2016, J. Phys. Chem. A, 120, 4988
  * Quack & Troe (1974) Quack, M. & Troe, J. 1974, Ber. Bunsenges. Phys. Chem, 78, 240
  * Reddy et al. (2003) Reddy, R. R., Nazeer Ahammed, Y., Rama Gopal, K., & Baba Basha, D. 2003, Astrophys. Space Sci., 286, 419
  * Riaplov et al. (2003) Riaplov, E., Wyss, M., Maier, J. P., et al. 2003, J. Mol. Spectrosc., 222, 15
  * Rivière-Marichalar et al. (2019) Rivière-Marichalar, P., Fuente, A., Goicoechea, J. R., et al. 2019, A&A, 628, A16
  * Rodríguez-Baras et al. (2021) Rodríguez-Baras, M., Fuente, A., Riviére-Marichalar, P., et al. 2021, A&A, 648, A120
  * Roncero et al. (2018) Roncero, O., Zanchet, A., & Aguado, A. 2018, Phys. Chem. Chem. Phys., 20, 25951
  * Ruaud et al. (2016) Ruaud, M., Wakelam, V., & Hersant, F. 2016, MNRAS, 459, 3756
  * Saito et al. (1987) Saito, S., Kawaguchi, K., Yamamoto, S., et al. 1987, ApJL, 317, L115
  * Sanz-Sanz et al. (2015) Sanz-Sanz, C., Aguado, A., Roncero, O., & Naumkin, F. 2015, J. Chem. Phys., 143, 234303
  * Senekowitsch et al. (1990) Senekowitsch, J., Carter, S., Rosmus, P., & Werner, H.-J. 1990, Chem. Phys., 147, 281
  * Shingledecker et al. (2020) Shingledecker, C. N., Lambers, T., Laas, J. C., et al. 2020, AstroPhys. J., 888, 52
  * Simon et al. (1997) Simon, R., Stutzki, J., Sternberg, A., & Winnewisser, G. 1997, A&A, 327, L9
  * Song et al. (2016) Song, Y.-Z., Zhang, L.-L., Gao, S.-B., & Meng, Q.-T. 2016, Scientific Rep., 6, 37734
  * Spezzano et al. (2022) Spezzano, S., Fuente, A., Caselli, P., et al. 2022, A&A, 657, A10
  * Sternberg & Dalgarno (1995) Sternberg, A. & Dalgarno, A. 1995, ApJS, 99, 565
  * Stoecklin et al. (1988) Stoecklin, T., Halvick, P., & Rayez, J. C. 1988, J. Mol. Struct. (Theochem), 163, 267
  * Stoecklin et al. (1990a) Stoecklin, T., Rayez, J. C., & Duguay, B. 1990a, Chem. Phys., 148, 381
  * Stoecklin et al. (1990b) Stoecklin, T., Rayez, J. C., & Duguay, B. 1990b, Chem. Phys., 148, 399
  * Tarroni et al. (2007) Tarroni, R., Carter, S., & Handy, N. C. 2007, Mol. Phys., 105, 1129
  * Vastel et al. (2018) Vastel, C., Quénard, D., Le Gal, R., et al. 2018, MNRAS, 478, 5514
  * Vidal et al. (2017) Vidal, T. H. G., Loison, J.-C., Jaziri, A. Y., et al. 2017, Mon. Not. R. Astron. Soc., 469, 435
  * Visser et al. (2019) Visser, B., Beck, M., Bornhauser, P., et al. 2019, Mol. Phys., 117, 1645
  * Voronin (2004) Voronin, A. I. 2004, Chem. Phys., 297, 49
  * Wakelam et al. (2021) Wakelam, V., Dartois, E., Chabot, M., et al. 2021, A&A, 652, A63
  * Werner & Knowles (1988a) Werner, H. J. & Knowles, P. J. 1988a, J. Chem. Phys., 89, 5803
  * Werner & Knowles (1988b) Werner, H. J. & Knowles, P. J. 1988b, Chem. Phys. Lett., 145, 514
  * Werner et al. (2012) Werner, H.-J., Knowles, P. J., Knizia, G., Manby, F. R., & Schütz, M. 2012, WIREs Comput Mol Sci, 2, 242
  * Woon & Herbst (2009) Woon, D. E. & Herbst, E. 2009, ApJ. Sup. Series, 185, 273
  * Zanchet et al. (2018) Zanchet, A., del Mazo, P., Aguado, A., et al. 2018, PCCP, 20, 5415
  * Zanchet et al. (2016) Zanchet, A., Roncero, O., & Bulut, N. 2016, Phys. Chem. Chem. Phys., 18, 11391
  * Zanchet et al. (2009) Zanchet, A., Roncero, O., González-Lezana, T., et al. 2009, J. Phys. Chem. A, 113, 14488
  * Zeimen et al. (2003) Zeimen, W. B., Klos, J., Groenenboom, G. C., & van der Avoird, A. 2003, J. Chem. Phys., 118, 7340
  * Zhang et al. (2018) Zhang, L. L., Song, Y. Z., Gao, S. B., & Meng, Q. T. 2018, J. Phys. Chem. A, 122, 4390