Formation of interstellar SH$^+$ from vibrationally excited H$_2$: Quantum study of S$^+$ + H$_2$ $\rightleftarrows$ SH$^+$ + H reactions and inelastic collisions

TL;DR

This study uses quantum mechanical calculations to determine new rate constants for the formation and destruction of interstellar SH$^+$, revealing higher rates that better match recent astronomical observations.

Contribution

The paper introduces two new potential energy surfaces and provides updated rate constants for SH$^+$ chemistry, improving astrochemical models of dense photodissociation regions.

Findings

New rate constants are 2-6 times larger than previous estimates.

Updated models produce four times larger SH$^+$ column densities.

Results align better with ALMA observations of the Orion bar.

Abstract

The rate constants for the formation, destruction, and collisional excitation of SH$^+$ are calculated from quantum mechanical approaches using two new SH$_2^+$ potential energy surfaces (PESs) of $^4A''$ and $^2A''$ electronic symmetry. The PESs were developed to describe all adiabatic states correlating to the SH$^+$ ($^3\Sigma^-$) + H($^2S$) channel. The formation of SH$^+$ through the S$^+$ + H$_2$ reaction is endothermic by $\approx$ 9860 K, and requires at least two vibrational quanta on the H$_2$ molecule to yield significant reactivity. Quasi-classical calculations of the total formation rate constant for H$_2$($v=2$) are in very good agreement with the quantum results above 100K. Further quasi-classical calculations are then performed for $v=3$, 4, and 5 to cover all vibrationally excited H$_2$ levels significantly populated in dense photodissociation regions (PDR). The new calculated formation and destruction rate constants are two to six times larger than the previous ones and have been introduced in the Meudon PDR code to simulate the physical and illuminating conditions in the Orion bar prototypical PDR. New astrochemical models based on the new molecular data produce four times larger SH$^+$ column densities, in agreement with those inferred from recent ALMA observations of the Orion bar.

Tables

Table 1. Table 1: Diatomic constants of SH + determined theoretically (top) and experimentally.

re(Å)	De(eV)	D0(eV)	Reference
1.365	3.72	3.56	This work
1.365	3.67	3.52	Zanchet et al. (2013a)
1.365	3.68		Song et al. (2018)
1.365	3.52	3.36	Stancil et al. (2000)
1.354	3.84	3.68	McMillan et al. (2016)
1.364	3.65 $\pm$ 0.13		Huber & Herzberg (1979)
1.364		3.48	Dunlavey et al. (1979)
1.363	3.70	3.54	Rostas et al. (1984)

Table 2. Table 2: Parameters used in the abc calculations.

jtot	0–100	Total angular momentum quantum number $J$.
jmax	24	Maximum rotational quantum number of any channel.
kmax	6	Helicity truncation parameter.
rmax	15.9	Maximum hyper-radius ${\rho }_{\text{max}}$ (in Å).
mtr	600	Number of log derivative propagation sectors.
emax	3.5	Maximum internal energy in any channel (in eV).
ipar	-1, 1	Triatomic parity eigenvalue $P$.
jpar	-1, 1	Diatomic parity eigenvalue $p$.

Table 3. Table 3: Parameters used in the wave packet calculations for the doublet PES

Reactant scattering coordinate range:	$R_{min}$=0.001Å;	$R_{max}$=40.0Å
Number of grids points in $R$:	600
Diatomic coordinate range:	$r_{min}$=0.001Å	$r_{max}$=32.0Å
Number of grid points in $r$:	480
Number of angular basis functions:	120
Number of projection of total angular momentum:	${\Omega }_{max}$=29
Center of initial wave packet:	19.0Å
Initial translational kinetic energy/eV:	0.495
Distance for flux determination:	$r$=12.0Å
Number of Chebychev iterations:	201000 for $J=0$	185000 for $J=10-60$	5000 for $J>60$

Table 4. Table 4: Parameters of the fit of the formation and destruction rate constants in the quartet electronic state

Reaction	$\alpha$	$\beta$	$\gamma$
	(cm3 s-1)		(K)
Formation
S+ + H2($v$=2) $\to$ SH+ + H (*)	2.88 10-10	-0.15	42.93
S+ + H2($v$=3) $\to$ SH+ + H	9.03 10-10	-0.11	26.15
S+ + H2($v$=4) $\to$ SH+ + H	12.96 10-10	-0.04	40.80
S+ + H2($v$=5) $\to$ SH+ + H	12.09 10-10	0.09	34.51
Destruction
SH+($v$=0) + H $\to$ S+ + H2 (*)	1.86 10-10	-0.41	27.38

Equations

$${\mathrm{M}}^{+}+{\mathrm{H}}_2\to {\mathrm{MH}}^{+}+\mathrm{H},$$

$$S{H}^{+}{(}^3{\Sigma }^{-})+H{(}^{2}S)$$

$$V_{ABC}=\sum _A{V}_A^{(1)}+\sum _{AB}{V}_{AB}^{(2)}(r_{AB})+V_{ABC}^{(3)}(r_{AB},r_{AC},r_{BC}),$$

$${\rho }_{AB}(r_{AB})=r_{AB}{e}^{-{\beta }_{AB}^{(2)}{r}_{AB}},AB=S{H}^{+},S{H}^{+},HH,$$

$$V_{ABC}^{(3)}(r_{AB},r_{AC},r_{BC})=\sum _{ijk}^K{d}_{ijk}{\rho }_{AB}^i{\rho }_{AC}^j{\rho }_{BC}^k,$$

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

Full text

1 Introduction

Molecular hydrogen is by far the most abundant molecule in the Universe and triggers the chemistry in the interstellar medium (ISM) through reactions with the most abundant atoms and ions. Owing to their high reactivity, the formation of molecular hydride cations is an important first step toward the synthesis of more complex interstellar molecules.

In UV-illuminated environments such as the low-density diffuse clouds or the edges of dense molecular clouds close to massive stars, the so-called photodissociation regions (PDRs, Hollenbach & Tielens 1997), the formation of hydride cations can start through the following reaction.

[TABLE]

where M+ is an ion typically formed by photo-ionization of element M, such as carbon or sulfur, with an ionization potential below 13.6 eV.

The launch of the Herschel satellite in 2009 opened a new frequency window in the far-IR/submillimeter wavelengths allowing for the detection of molecular hydrides in the ISM (for a review, see Gerin et al. 2016). The Herschel satellite allowed for the detection of rotational line emission from CH+ or SH+ in dense interstellar PDRs (Nagy et al. 2013; Naylor et al. 2010; Pilleri et al. 2014; Joblin et al. 2018), in the irradiated walls of protostellar outflows (Falgarone et al. 2010; Benz et al. 2010, 2016), and in the circumstellar envelopes around hot planetary nebulae (da Silva Santos et al. 2018).

For CH+ and SH+ ions, reaction (1) is very endothermic ($\Delta E/k=4300$ and 9860 K, respectively) and the above interstellar detections may seem surprising. In UV-irradiated environments however, H2 can be radiatively pumped to vibrationally excited levels (e.g., Black & Dalgarno 1976). These vibrational states have internal energies high enough to overcome reaction endothermicities (Sternberg & Dalgarno 1995). Indeed, laboratory experiments showed that for CH+, reaction (1) becomes exothermic and fast if H2 ($v\geq 1$) (Hierl et al. 1997). State-to-state rate constants for this reaction have have been calculated from quantum calculations (Zanchet et al. 2013b) and PDR and excitation models using these data predict CH+ abundances and rotational line intensities close to the observed ones (e.g., Agúndez et al. 2010; Godard & Cernicharo 2012; Nagy et al. 2013; Faure et al. 2017; Joblin et al. 2018). In addition, mapping observations of the Orion molecular cloud have revealed spatial correlation between the [C ii] 158 $\mu$m emission (from ionized carbon C+), the IR H2 ($\geq$1) emission, and the CH+ ($j=1-0$) rotational emission (Goicoechea et al. 2019). This work observationally demonstrates that the interstellar CH+ emission is widespread in UV-irradiated dense gas, and that its main formation route is reaction (1).

Formation of SH+ through reaction (1) is much more endothermic than that of CH+ and reaction (1) only becomes exothermic when H2 ($v\geq 2$) (e.g., Zanchet et al. 2013a). This is true considering H2 vibrational levels alone. However, we note that taking into account H2 rovibrational levels, reaction (1) becomes exothermic for v=0, j$\geq$11 and v=1, j$\geq$7. These rotationally excited levels within a given vibrational state can be populated, but their observed column densities are much lower than those of the low-j levels of vibrationally excited H2 (e.g., Habart et al. (2011); Kaplan et al. (2017)). . Still, SH+ has been detected in low-density diffuse clouds, for example through absorption lines measurements (Menten et al. 2011; Godard et al. 2012) and in denser and more strongly UV-irradiated environments through emission lines (Benz et al. 2010; Nagy et al. 2013). Contrary to CH+, SH+ can be observed from ground-based millimeter-wave telescopes (Müller et al. 2014; Halfen & Ziurys 2015). Then, high-angular interferometric images of the Orion Bar PDR obtained with ALMA have revealed that the SH+ emission nicely delineates the H2 dissociation front at the edge of the PDR (Goicoechea et al. 2017), the same gas layers where H2 molecules are UV-pumped to highly excited vibrationally levels (Kaplan et al. 2017). These results further confirm the need for accurate state-to-state reaction rates for reaction (1).

In addition, the starting hydrogen abstraction reactions involving H2 and S, S+, and SH+ are all highly endothermic (Neufeld et al. 2015) and this implies that SH+ is destroyed in reactions with atomic hydrogen and electron recombinations. Both H and $e^{-}$ are abundant in UV-irradiated gas, and therefore SH+ is expected to be relatively reactive and a short-lived hydride in PDRs.

Therefore, it is very relevant to study the first steps of sulfur chemistry. In particular, Goicoechea et al. (2017) concluded that the abundances of SH+ inferred from the ALMA images of the Orion Bar were $\sim$ 3-30 times higher than those predicted by a state-of-the-art PDR model using specific rate constants, computed by some of us for the S+ + H2 ($v$) reaction (with $v$ from 0 to 4) using the ground-quartet-state potential energy surface (PES) and a quasi-classical trajectories method (Zanchet et al. 2013a).

In this work, we focus on a ab initio quantum study of the relevant rate constants for the formation, destruction, and excitation of SH+. By studying the reaction dynamics using a quasi-classical trajectory (QCT) method on a new PES, Zanchet et al. (2013a) demonstrate that the S+(${}^{4}S$) + H2(v$\geq$ 2) collisions was an efficient way to produce SH+. The PES was also used to perform a quantum study (Zanchet et al. 2016) that showed that at low collisional energies some resonances appear, associated to roaming, and that the quantum reaction cross section for H${}_{2}(v=2,j=0)$ was considerably larger than the classical ones. These differences disappeared for higher rotational excitation of H2. The same PES was also used to estimate the destruction rate constants of SH+ colliding with H. However, this rate of destruction can be considered as not very accurate since the PES was designed essentially to study the formation.

More recently, a new PES has been calculated using a larger basis set in the ab initio calculations (Song et al. 2018). From this PES, these latter authors derived slightly larger QCT cross sections than those previously reported. These latter works show the necessity for high-level  ab initio calculations to correctly describe the interactions between the three atoms and to take in account quantum effects in the reaction dynamics in order to get accurate estimations of the rate constants.

In order to describe the reaction

[TABLE]

not only is the quartet PES needed, but also a doublet PES, since the spins of both SH+(${}^{2}\Sigma^{-}$) and H(${}^{2}S$) can give rise to both doublet and quartet states which are asymptotically degenerate, as shown in Fig. 1. In this work the objective is to present a new set of PESs for the quartet and doublet states using a larger electronic basis set than that used by Zanchet et al. (2013a) and Song et al. (2018).It should be noted that during the preparation of this work, a new PES for the doublet electronic state was also published (Zhang et al. 2018) and that QCT calculations were performed for the formation from the excited state of S${}^{+}(^{2}D)$. Reactions of H2 molecules with S+ ions in electronic excited states are expected to be negligible inside molecular clouds and only relevant at the PDR/HII interface layers where a small fraction of H2 might exist.

The paper is organized as follows. First the construction of the new PESs for the two electronic states (quartet and doublet) is described in Sect. 2. Sections 3 and 4 are devoted to the reaction dynamics on the quartet and doublet states, respectively. Finally, in Sect. 5, the astrophysical implications of the new results obtained for both the formation and destruction of SH+ are discussed.

2 Ab initio calculations and analytical fit

The ${}^{4}A^{\prime\prime}$ and ${}^{2}A^{\prime\prime}$ electronic states of the SH${}_{2}^{+}$ system both correlate with the SH+ (${}^{3}\Sigma^{-}$) + H(${}^{2}S$) asymptote where they are degenerate. In the other asymptotic channel, each state correlates to a different electronic state of the S+ cation, the quartet state connects to the ground state S+(${}^{4}S$) + H2(${}^{1}\Sigma^{+}$) while the doublet connects to the excited state S+(${}^{2}D$) + H2(${}^{1}\Sigma^{+}$) reactants. The ${}^{4}S$ and ${}^{2}D$ states of S+ are separated by $\approx$ 1.8 eV as shown in Fig. 1. These states are only coupled by the spin-orbit coupling which is neglected in this work. In the S+(${}^{2}D$) asymptote, five electronic states are degenerate but only one correlates to the ground state of SH+ while the four others correlate to excited states of SH+ products, and hence their PESs are not considered in this work. However, to get an accurate and homogeneous description of the doublet and quartet PESs in all the configuration space and in particular in the asymptotic channels, the five doublet states as well as the quartet state were taken into account in the ab initio calculations. The state-average complete active space (SA-CASSCF) method (Werner & Knowles 1985) was employed to calculate the first ${}^{4}A^{\prime\prime}$ together with the two first ${}^{2}A^{\prime}$ and the three first ${}^{2}A^{\prime\prime}$ electronic states. The active space considered consists in seven electrons distributed in six orbitals (5 - 9a’ and 2a”) in order to include all valence orbitals of sulfur and the 1s orbitals from both hydrogen atoms. The obtained state-average orbitals and multireference configurations were then used to calculate both the lowest ${}^{4}A^{\prime\prime}$ and the lowest ${}^{2}A^{\prime\prime}$ state energies with the internally contracted multireference configuration interaction method (ic-MRCI) including simple and double excitation (Werner et al. 1988) and Davidson correction (Davidson 1975). The 1s orbital of sulfur was kept frozen. For both sulfur and hydrogen atoms, the aug-cc-pV5Z basis set (aV5Z) was used, including spdfgh and spdfg basis functions and all calculations were done using the MOLPRO suite of programs (MOLPRO is a package of ab initio programs designed by H. -J. Werner and P. J. Knowles and with contributions from version 2012).

Ab initio calculations were performed over 3800 geometries of the SH${}_{2}^{+}$ system. To sample the geometries, three sets of coordinates were used. Jacobi coordinates associated to the S+ + H2 channel were used for a good description of this asymptotic channel while internal coordinates H-H-S+ and H-S+-H were employed to sample the channel associated to SH++H. These icMRCI+Q energies for the electronic states ${}^{4}A^{\prime\prime}$ and ${}^{2}A^{\prime\prime}$ were then fitted separately using the GFIT3C procedure (Aguado & Paniagua 1992; Aguado et al. 1993, 1998), in which a global PES is represented by a many-body expansion:

[TABLE]

where $V_{A}^{(1)}$ represents the energy of the atom A ($A$ = S+, H, H) in its corresponding electronic state, $V_{AB}^{(2)}$ the diatomic terms ($AB$ = SH+, SH+, HH) in the corresponding electronic state, and $V_{ABC}^{(3)}$ the three-body term($ABC$ = SHH+).

The diatomic terms are written as a sum of short- and long-range contributions. The short-range potential is defined as a shielded Coulomb potential, whereas the long-range term is a linear combination of modified Rydberg functions (Rydberg 1931) defined as :

[TABLE]

with $\beta_{AB}^{(2)}>0$. The root-mean-square (rms) error of the fitted SH+ potential which is common for both PESs since they share the same SH++H asymptote is $\approx 0.045$ kcal/mol. The rms of the fitted H2 potential is $\approx 0.029$ kcal/mol for the quartet state where H2 is the pure diatomic and $\approx 0.346$ kcal/mol for the doublet state where the H2 implicitly considers the avoided crossing arising from the crossing of the two SH+ electronic states correlating to S+(${}^{4}S$) and S+(${}^{2}D$).

The three-body term is expressed as an expansion:

[TABLE]

where $\rho_{AB}=r_{AB}\exp{-\beta_{AB}r_{AB}}$ are modified Rydberg functions (Rydberg 1931; Aguado & Paniagua 1992). For SHH+, there are only two nonlinear parameters, $\beta_{SH^{+}}$ and $\beta_{HH}$, and additional constraints in the linear parameters $d_{ijk}$ to ensure symmetry of the PES with respect to the permutation of the two H atoms  (Aguado & Paniagua 1992; Aguado et al. 1993, 1998). The linear parameters $d_{ijk},(i+j+k)$ and the two nonlinear parameters $\beta_{SH^{+}}$ and $\beta_{HH}$, are determined by fitting the approximately 3000 calculated ab initio energies after the substraction of the one- and two-body contributions. In the present case, the order $L$ is 10 for both states, giving an overall rms error of 0.430 kcal/mol and 0.392 kcal/mol for the ${}^{4}A^{\prime\prime}$ and ${}^{2}A^{\prime\prime}$ states, respectively.

The two PESs exhibit completely different topographies. The ${}^{4}A^{\prime\prime}$ state does not present any minimum out of the van der Waals wells in the asymptotic channels and does not present any barrier to reaction. The SH++H $\rightarrow$ S++H2 reaction is exothermic on this surface and reactive collisions are likely to occur in competition with the inelastic collisions. On the other hand, the ${}^{2}A^{\prime\prime}$ state presents a deep insertion HSH well and does not present any barrier either. For this state, in contrast with the previous case, the SH++H $\rightarrow$ S++H2 is endothermic and only inelastic collisions can occur (pure or involving H exchange). The main features are summarized in the minimum energy path shown in Fig. 1

The present PES for the quartet state, calculated with the aV5Z basis set, is very similar to that of Zanchet et al. (2013a), calculated with the smaller aVQZ basis set. Two relevant differences can be appreciated.

The first one is the De well depth of the SH+ diatomic, which is deeper by 5 meV for the present, larger aV5Z basis set as compared to that obtained with the AVQZ basis set (Zanchet et al. 2013a), as can be seen in Fig. 2. In Table 1 the equilibrium distance and dissociation energies, De and D0, of SH+ obtained in previous works are compared with the present results. All the theoretical values of D0 (McMillan et al. 2016; Stancil et al. 2000; Song et al. 2018; Zanchet et al. 2013a) are within the experimental uncertainty (Huber & Herzberg 1979; Dunlavey et al. 1979; Rostas et al. 1984). Assuming that the accuracy of theoretical calculations improves with the size of the basis set, we conclude that the SH+ diatomic considered in this work is slightly more accurate than the one used previously (Zanchet et al. 2013a). This sensitivity to the basis set highlights the difficulty in treating the electronic structure of third raw atoms, which remain to be challenging calculations. On the other hand, the constants of the H2 diatomic are similar using the AVQZ or the AV6Z basis set, leading to the same De. As a consequence, the increase of the SH+ dissociation energy yields a reduction in the endothermicity of the reaction by about 5meV in our present PES, but this value may still be overestimated considering that McMillan et al. (2016) find a De that is 12 meV deeper using a aV6Z basis set in a study centered on the SH+ diatomic.

The second difference is not related to the size of the basis set, but to the sampling of ab initio points on the PES. Indeed, in the PES of Zanchet et al. (2013a), a poor density of ab initio points was included in the fit of the SH++H channel. As a consequence, the barrier present in the collinear configuration H-SH+ was not sampled and is not present in the former PES while it is correctly described by the new potential. As we see in the following sections, these differences have significant implications for the formation and destruction of SH+. The barrier considerably reduces the acceptance cone of the SH+ + H reaction as it only allows destruction at low temperatures if the H atom collides on the hydrogen side of SH+. To put it in another way, the absence of the barrier in the potential of Zanchet et al. (2013a) artificially increases the probability of destruction allowing collisions on the S side to react when they are not supposed to do so.

3 Collision dynamics

3.1 Dynamics on the quartet state

The reaction dynamics on the quartet state has been studied from a time-independent treatment based on hyperspherical coordinates. The calculations were performed with the ABC code (Skouteris et al. 2000) using the parameters listed in Table  2. The formation of SH+ through the S$(^{4}S)^{+}$ + H${}_{2}(v,j)$ $\rightarrow$ SH${}^{+}(^{3}\Sigma^{-})$ + H reaction presents some similar trends to that obtained previously (Zanchet et al. 2013a, 2016), but significant quantitative difference are found. The total reaction cross sections are shown in Fig. 2 and compared with quantum wave packet calculations of Zanchet et al. (2016) using the previous quartet state PES of Zanchet et al. (2013a). The present formation cross sections for collisions of S$(^{4}S)^{+}$ + H2($v=2,j=0$) is larger than that previously obtained, while the reactivity for the S$(^{4}S)^{+}$ + H2($v=2,j=1$) reaction is the reverse. The cross section for collisions with H2($v=2,j=1$) remains however considerably larger than that for H2($v=2,j=0$); that is, more than a factor 4 at 100 meV. This ortho/para-H2 differentiation in the reactivity is also valid for $j$ = 2 and 3, but disappears for higher $j$. At very low collisional energy, the cross sections on the present PES are always larger. When performing the full-dimensional dynamic calculation it is difficult to attribute this increase to a particular detail of the PES. We tentatively attribute this increase to the decrease of the endothermicity in the present PES as compared to the previous one.

In any case, the overall behavior is rather similar, and in both PESs, the resonances associated to roaming mediate the reactivity at low collisional energies (Zanchet et al. 2016).

The destruction cross sections, for the SH${}^{+}(^{3}\Sigma^{-},v,j)$ + H $\rightarrow$ S$(^{4}S)$ + H2 reaction, are shown in Fig. 3. At low translational energy ($E<0.1$ eV), the cross sections decrease with increasing rotation, that is, the rotational excitation inhibits the reaction. However, for $E>0.2$ eV, all the cross sections exhibit an almost constant value of 3-4 $\AA^{2}$. The order of magnitude of the formation and destruction cross sections is very similar for high energies and nearly independent of rotational excitation of the reactants.

The inelastic and exchange cross sections for the H + SH+(v=0,j=0) are shown in Fig. 4. The exchange cross sections are more than two orders of magnitude smaller than the inelastic ones, what is explained by the fact that the collision proceeds, in this case directly, without forming a collision complex, making the H-exchange very improbable. The inelastic cross sections decrease with increasing $\Delta j$ as usual.

3.2 Dynamics on the doublet state: inelastic and exchange

In order to fully describe the inelastic collision, the dynamics on the doublet state must be considered, as it accounts for one third (or 2/6) in the electronic partition function. This state correlates with the excited sulfur cation, S${}^{+}(^{2}D)$, which is nearly 2 eV higher than the S${}^{+}(^{4}S)$ and $\approx$ 1 eV above the SH+ (see Fig. 1). Therefore, at low to intermediate collision energies, only the inelastic and exchange process are possible in this electronic state, according to Eq. (1.b) if the Born-Oppenheimer approximation is assumed. However, the presence of the deep well of $\approx$ 4eV from the SH+ + H asymptote makes it difficult to apply the ABC code to study this collision. Instead, here we used a quantum wavepacket MADWAVE3 code (Gómez-Carrasco & Roncero 2006; Zanchet et al. 2009) as it was applied to study the OH++H dynamics in the doublet state (Bulut et al. 2015) with a very similar deep well. The parameters used in these calculations are listed in Table 3. The state-to-state reaction probabilities have been calculated for $J=0,10,20,30,...,100$ and for intermediate values of $J$; they are interpolated using the $J$-shifting approximation as done previously (Zanchet et al. 2013b; Gómez-Carrasco et al. 2014; Bulut et al. 2015). After summing over all partial waves, the state-to-state cross sections were obtained.

The calculated cross sections for the SH+(v=0,j=0) + H in the doublet state in the inelastic and exchange channels are shown in Fig. 5. The exchange cross section is more than one order of magnitude lower than the inelastic one. The doublet state shows a deep well and as a consequence, the calculated reaction probabilities at all the partial waves show the presence of many resonances. This necessitates propagation of the wave packet for a high number of iterations. In spite of the presence of a dense manifold of resonance, the reaction is not statistical. The same situation was obtained for OH++H dynamics in the doublet state (Bulut et al. 2015), since in the deep well the SH+ distance does not vary with respect to that of bare SH+, and the light colliding H atom is not able produce a complete randomization of the available energy between the two SH bonds, which would lead to a probability of 0.5 for each rearrangement channel.

3.3 Formation and destruction rate constants

Quantum rate constants are obtained from the corresponding cross sections, described above for several reactive/inelastic collisions, by numerically integrating over collision energy with a Boltzmann distribution at each temperature and summing over the rotational state of the reactants with a weight determined from a Boltzmann distribution. These rate constants are specific for each vibrational state of either H2 or SH+, for formation and destruction processes, respectively.

Quasi-classical rate constants were also calculated using the method of Karplus et al. (1965) as implemented in the miQCT code (Zanchet et al. 2013a; Dorta-Urra et al. 2015; Zanchet et al. 2016). For each temperature, five hundred thousands trajectories are run changing the initial conditions, consistent with a Boltzmann distribution of translation and rotation energy at a given temperature $T$. The initial distance between reactants is set to 21 Å, with a maximum impact parameter of 11.6 Å. The trajectories are stopped when any distance becomes greater than 24 Å. For the destruction reaction, SH+ is considered to be in $v=0$, and for the formation reaction, H2 is in $v=2$, 3, 4, and 5. A Boltzmann distribution for rotational excitation is considered in the QCT calculations of reaction rate for all initial vibrational excitations.

The calculated (quantum and QCT) reaction rate constants for formation and destruction processes are shown in Fig. 6, and compared to those of Zanchet et al. (2013a), calculated using a QCT method. The quantum and QCT calculations shown in the two lower panels are in very good agreement for temperatures higher than 100-200 K. This indicates that the QCT results for the rate constants are relatively accurate and justify the use of this less time-consuming method to calculate the formation rates for H2($v=3$, 4 and 5). Below 100 K, quantum effects are more important, and the quantum and classical results diverge.

The QCT results obtained here are larger than those reported from the previous quartet PES of Zanchet et al. (2013a). This difference is attributed to the slight decrease in the endothermicity. The factor depends on the initial vibrational state, but we can conclude that the present formation rate constants are larger than those previously reported by a factor of between three and six. The present QCT results are also similar to those reported by Song et al. (2018) using a different PES for the quartet.

The destruction rate constants of this work are also larger than the previous ones reported in Zanchet et al. (2013a) by a factor that strongly depends on the temperature: about five at 1000 K and two at 100K. The reason for this change is also attributed to the higher accuracy of the present PES.

Since the formation and destruction rate constants increase significantly with respect to those reported previously (Zanchet et al. 2013a), the predicted abundance of SH+ in interstellar conditions should be influenced by these new data. The formation and destruction quantum rate constants have been fitted to the usual expression

[TABLE]

and the $\alpha$, $\beta$ and $\gamma$ parameters obtained are listed in Table 4.

4 Updated predictions of the SH+ abundance in the Orion Bar PDR

We use the new S+ + H2 $\rightleftarrows$ SH+ + H rate constants to update the estimated abundance of SH+ in a dense PDR like the Orion Bar (Goicoechea et al. 2016, and references therein). The external and most UV-irradiated layers of dense PDRs have moderate column densities of vibrationally excited H2 that depend on the gas density and flux of UV photons. These H2($v\geq 1$) molecules greatly affect the formation of hydrides (e.g., see models of Sternberg & Dalgarno 1995; Agúndez et al. 2010). SH+ has been detected in the Orion Bar with Herschel (Nagy et al. 2013) and also from the ground (Müller et al. 2014). Zanchet et al. (2013a) have previously shown PDR models that emphasized the role of H${}_{2}\,(v\geq 2$) in the formation of SH+. Indeed, ALMA images of the Orion Bar (Goicoechea et al. 2017) have shown that the SH+ emission comes from a narrow gas layer that delineates the irradiated edge of the PDR. The SH+ column densities estimated from observations however, revealed values that are higher than model predictions by a factor of between 3 and 30 (depending on cloud geometry considerations). The latter predictions were made using the Meudon PDR model (e.g., Le Petit et al. 2006; Goicoechea & Le Bourlot 2007) adapted to the physical and chemical structure of the Orion Bar.

The Orion Bar is illuminated by a far-UV ($<$13.6 eV) radiation field $G_{0}$ of a few times 104, where $G_{0}=1$ is equal to 1.6$\times$103 erg cm*-2* s*-1*, the far-UV flux in the solar neighborhood, integrated from $\sim$912 Å to $\sim$2400 Å  (Habing 1968). Here we used version 1.5.2 of the Meudon PDR code to model an isobaric PDR, with a constant gas thermal pressure of $P_{\rm th}=n_{\rm H}$  $T$ = 2$\cdot$108 cm*-3* K determined from previous observations (see e.g., Goicoechea et al. 2016; Joblin et al. 2018) and introduce the new SH+ formation and destruction rate constants. We adopt an impinging FUV field of $G_{0}=2\cdot 10^{4}$, an undepleted sulfur abundance of 1.4$\times$10*-5* with respect to H nuclei (Asplund et al. 2005), and dust grain properties appropriate to the flatter extinction curve observed toward Orion (i.e., an extinction-to-color-index ratio $R_{V}=A_{V}/E_{B-V}=5.5$, e.g., Cardelli et al. 1989). Figure 7 shows fractional abundances of SH+, H2 (total), H2 ($v=2$), S+, and S, and atomic hydrogen as a function of depth into the cloud (in magnitudes of visual extinction $A_{V}$). The figure also shows the gas temperature (in grey), and $f_{\rm H_{2}}$, the fraction of H2 that is in vibrationally excited levels $v\geq 2$ with respect to the ground $f_{\rm H_{2}}$ = $n\,(v\geq 2)$/$n\,(v=0)$. The SH+ abundance profiles obtained from the older S+ + H2 $\rightleftarrows$ SH+ + H rate constants of Zanchet et al. (2013a) (dashed curve) are also shown.

As observationally shown by Goicoechea et al. (2017), SH+ is only abundant in the irradiated edge of the PDR, with the abundance peak located at $A_{V}\approx 1$. The models that use the new formation and destruction rate constants produce SH+ column densities ($N$(SH+)$=$(0.2-2.7)$\cdot$1013 cm*-2*, depending on the assumed inclination of the PDR, see Goicoechea et al. 2017) that are higher than models using previous data (Nagy et al. 2013; Zanchet et al. 2013a; Goicoechea et al. 2017) by a factor of approximately four. This enhancement in the production of SH+ reconciles the SH+ column densities inferred from observations with the PDR model predictions. These results further empathize the role of UV-pumped H2 molecules in the formation of simple molecules at the edge of strongly irradiated dense PDRs.

5 Conclusions

New PESs for the ${}^{4}A^{\prime\prime}$ and ${}^{2}A^{\prime\prime}$ electronic states have been calculated for the S+ + H2 $\rightleftarrows$ SH+ + H reactions. Quantum state-to-state reactive and inelastic cross sections have been calculated for S+ + H2(v=2) and SH+(v=0) + H cases, and the corresponding rate constants were calculated and compared to QCT calculations. Good agreement was found. This justifies the use of the QCT method for the calculation of the formation rate constants for higher vibrational states, S+ + H2($v=3$, 4 and 5 ). These new rate constants are larger than those previously calculated (Zanchet et al. 2013a) by roughly a factor of between two and six.

We used the new reaction rates in the Meudon PDR code and simulated the UV-illuminating and physical conditions in the Orion Bar PDR. The new models yield column densities that are a factor of four higher than those obtained with the previous formation and destruction rate constants (Zanchet et al. 2013a). The new fractional abundances of SH+ are within the uncertainties of the SH+ column densities inferred from observations (Goicoechea et al. 2017).

The spin-rotation couplings of SH+(${}^{3}\Sigma^{-}$) have not been included, and therefore the present state-to-state rates cannot be directly used to model the fine structure observed for SH+. This is being done in a separate study using a recoupling technique (Faure & Lique 2012), in which it is assumed that the electronic spin is a spectator during the collision.

Acknowledgements.

We thank E. Bron for his help with the new version of the Meudon PDR code. We acknowledge the French-Spanish collaborative project PICS (Ref. PIC2017FR7). The research leading to these results has received funding from MICIU under grants No. FIS2017-83473-C2 and AYA2017-85111-P. FL acknowledges financial support from the Institut Universitaire de France. NB acknowledges the computing facilities by TUBITAK-TRUBA.