11institutetext: Instituto de Física Fundamental (CSIC). Calle Serrano
121-123, 28006, Madrid, Spain. 11email<EMAIL_ADDRESS>22institutetext: Facultad de Ciencias. Universidad Autónoma de Madrid, 28049
Madrid, Spain. 33institutetext: Institut de Radioastronomie Millimétrique
(IRAM), Grenoble, France. 44institutetext: LERMA, Observatoire de Paris, PSL
Research University, CNRS, Sorbonne Universités, 92190 Meudon, France.
55institutetext: Observatorio Astronómico Nacional (OAN), Alfonso XII, 3,
28014 Madrid, Spain. 66institutetext: Max-Planck-Institut für Radioastronomie,
Auf dem Hügel 69, 53121 Bonn, Germany. 77institutetext: OASU/LAB-UMR5804,
CNRS, Université Bordeaux, 33615 Pessac, France. 88institutetext: European
Southern Observatory, Alonso de Cordova 3107, Vitacura, Santiago, Chile.

# Bottlenecks to interstellar sulfur chemistry

Sulfur-bearing hydrides in UV-illuminated gas and grains

J. R. Goicoechea 11 A. Aguado 22 S. Cuadrado 11 O. Roncero 11 J. Pety 33 E.
Bron 44
A. Fuente 55 D. Riquelme 66 E. Chapillon 3377 C. Herrera 33 C. A. Duran 6688

(Received 25 October 2020 / Accepted 23 December 2020)

Hydride molecules lie at the base of interstellar chemistry, but the synthesis
of sulfuretted hydrides is poorly understood and their abundances often
crudely constrained. Motivated by new observations of the Orion Bar
photodissociation region (PDR) – 1′′ resolution ALMA images of SH+; IRAM 30m
detections of bright H${}_{2}^{32}$S, H${}_{2}^{34}$S, and H${}_{2}^{33}$S
lines; H3S+ (upper limits); and SOFIA/GREAT observations of SH (upper limits)
– we perform a systematic study of the chemistry of sulfur-bearing hydrides.
We self-consistently determine their column densities using coupled
excitation, radiative transfer as well as chemical formation and destruction
models. We revise some of the key gas-phase reactions that lead to their
chemical synthesis. This includes ab initio quantum calculations of the
vibrational-state-dependent reactions $\rm
SH^{+}+H_{2}({\it{v}})\rightleftarrows H_{2}S^{+}+H$ and $\rm
S\,+\,H_{2}\,({\it{v}})\rightleftarrows SH\,+\,H$. We find that reactions of
UV-pumped H2($v$ $\geq$ 2) molecules with S+ ions explain the presence of SH+
in a high thermal-pressure gas component, $P_{\rm th}/k$ $\approx$ 108 cm-3 K,
close to the H2 dissociation front (at $A_{V}$ $<$ 2 mag). These PDR layers
are characterized by no or very little depletion of elemental sulfur from the
gas. However, subsequent hydrogen abstraction reactions of SH+, H2S+, and S
atoms with vibrationally excited H2, fail to form enough H2S+, H3S+, and SH to
ultimately explain the observed H2S column density ($\sim$2.5$\times$1014
cm-2, with an ortho-to-para ratio of 2.9 $\pm$ 0.3; consistent with the high-
temperature statistical value). To overcome these bottlenecks, we build PDR
models that include a simple network of grain surface reactions leading to the
formation of solid H2S (s-H2S). The higher adsorption binding energies of S
and SH suggested by recent studies imply that S atoms adsorb on grains (and
form s-H2S) at warmer dust temperatures ($T_{d}$ $<$ 50 K) and closer to the
UV-illuminated edges of molecular clouds. We show that everywhere s-H2S
mantles form(ed), gas-phase H2S emission lines will be detectable.
Photodesorption and, to a lesser extent, chemical desorption, produce roughly
the same H2S column density (a few 1014 cm-2) and abundance peak (a few 10-8)
nearly independently of $n_{\rm H}$ and $G_{0}$. This agrees with the observed
H2S column density in the Orion Bar as well as at the edges of dark clouds
without invoking substantial depletion of elemental sulfur abundances.

###### Key Words.:

Astrochemistry — line: identification — ISM: clouds — (ISM:) photon-dominated
region (PDR) — ISM: clouds

##  1 Introduction

Hydride molecules play a pivotal role in interstellar chemistry (e.g., Gerin
et al., 2016), being among the first molecules to form in diffuse interstellar
clouds and at the UV-illuminated edges of dense star-forming clouds, so-called
photodissociation regions (PDRs; Hollenbach & Tielens, 1997). Sulfur is on the
top ten list of most abundant cosmic elements and it is particularly relevant
for astrochemistry and star-formation studies. Its low ionization potential
(10.4 eV) makes the photoionization of S atoms a dominant source of electrons
in molecular gas at intermediate visual extinctions $A_{V}$ $\simeq$ 2 - 4 mag
(Sternberg & Dalgarno, 1995; Goicoechea et al., 2009; Fuente et al., 2016).

The sulfur abundance, [S/H], in diffuse clouds (e.g., Howk et al., 2006) is
very close to the [S/H] measured in the solar photosphere
(${\rm[S/H]}_{\odot}$ $\simeq$ 1.4$\times$10-5; Asplund et al., 2009). Still,
the observed abundances of S-bearing molecules in diffuse and translucent
molecular clouds ($n_{\rm H}$ $\simeq$ $10^{2}$ $-$ $10^{3}$ cm-3) make up a
very small fraction, $<$ 1 $\%$, of the sulfur nuclei (mostly locked as S+;
Tieftrunk et al., 1994; Turner, 1996; Lucas & Liszt, 2002; Neufeld et al.,
2015). In colder dark clouds and dense cores shielded from stellar UV
radiation, most sulfur is expected in molecular form. However, the result of
adding the abundances of all detected gas-phase S-bearing molecules is
typically a factor of $\sim$102-103 lower than [S/H]⊙ (e.g., Fuente et al.,
2019). Hence, it is historically assumed that sulfur species deplete on grain
mantles at cold temperatures and high densities (e.g., Graedel et al., 1982;
Millar & Herbst, 1990; Agúndez & Wakelam, 2013). However, recent chemical
models predict that the major sulfur reservoir in dark clouds can be either
gas-phase neutral S atoms (Vidal et al., 2017; Navarro-Almaida et al., 2020)
or organo-sulfur species trapped on grains (Laas & Caselli, 2019).
Unfortunately, it is difficult to overcome this dichotomy from an
observational perspective. In particular, no ice carrier of an abundant sulfur
reservoir other than solid OCS (hereafter s-OCS, with an abundance of
$\sim$10-8 with respect to H nuclei; Palumbo et al., 1997) has been
convincingly identified. Considering the large abundances of water ice (s-H2O)
grain mantles in dense molecular clouds and cold protostellar envelopes (see
reviews by van Dishoeck, 2004; Gibb et al., 2004; Dartois, 2005), one may also
expect hydrogen sulfide (s-H2S) to be the dominant sulfur reservoir. Indeed,
s-H2S is the most abundant S-bearing ice in comets such as
67P/Churyumov–Gerasimenko (Calmonte et al., 2016). However, only upper limits
to the s-H2S abundance of $\lesssim$1 % relative to water ice have so far been
estimated toward a few interstellar sightlines (e.g., Smith, 1991; Jiménez-
Escobar & Muñoz Caro, 2011). These values imply a maximum s-H2S ice abundance
of several 10-6 with respect to H nuclei. Still, this upper limit could be
higher if s-H2S ices are well mixed with s-H2O and s-CO ices (Brittain et al.,
2020).

The bright rims of molecular clouds illuminated by nearby massive stars are
intermediate environments between diffuse and cold dark clouds. Such
environments host the transition from ionized S+ to neutral atomic S, as well
as the gradual formation of S-bearing molecules (Sternberg & Dalgarno, 1995).
In one prototypical low-illumination PDR, the edge of the Horsehead nebula,
Goicoechea et al. (2006) inferred very modest gas-phase sulfur depletions. In
addition, the detection of narrow sulfur radio recombination lines in dark
clouds (implying the presence of S+; Pankonin & Walmsley, 1978) is an argument
against large sulfur depletions in the mildly illuminated surfaces of these
clouds. The presence of new S-bearing molecules such as S2H, the first (and so
far only) doubly sulfuretted species detected in a PDR (Fuente et al., 2017),
suggests that the chemical pathways leading to the synthesis of sulfuretted
species are not well constrained; and that the list of S-bearing molecules is
likely not complete.

Figure 1: Overview of the Orion Bar. The (0′′, 0′′) position corresponds to
$\mathrm{\alpha_{2000}=05^{h}\,35^{m}\,20.1^{s}\,}$;
$\mathrm{\delta_{2000}=-\,05^{\circ}25^{\prime}07.0^{\prime\prime}}$. Left
panel: Integrated line intensity maps in the 13CO $J$ = 3-2 (color scale) and
SO 89-78 emission (gray contours; from 6 to 23.5 K km s-1 in steps of 2.5 K km
s-1) obtained with the IRAM 30 m telescope at 8′′ resolution. The white dotted
contours delineate the position of the H2 dissociation front as traced by the
infrared H2 $v$ = 1–0 $S$(1) line (from 1.5 to 4.0 $\times$ 10-4 erg s-1 cm-2
sr-1 in steps of 0.5 $\times$ 10-4 erg s-1 cm-2 sr-1; from Walmsley et al.,
2000). The black-dashed rectangle shows the smaller FoV imaged with ALMA (Fig.
3). The DF position has been observed with SOFIA, IRAM 30 m, and Herschel.
Cyan circles represent the $\sim$15′′ beam at 168 GHz. Right panel: H2S lines
lines detected toward three positions of the Orion Bar.

Interstellar sulfur chemistry is unusual compared to that of other elements in
that none of the simplest species, X=S, S+, SH, SH+, or H2S+, react
exothermically with H2 ($v$ = 0) in the initiation reactions X + H2
$\rightarrow$ XH + H (so-called hydrogen abstraction reactions). Hence, one
would expect a slow sulfur chemistry and very low abundances of SH+
(sulfanylium) and SH (mercapto) radicals in cold interstellar gas. However,
H2S (Lucas & Liszt, 2002), SH+ (Menten et al., 2011; Godard et al., 2012), and
SH (Neufeld et al., 2012, 2015) have been detected in low-density diffuse
clouds ($n_{\rm H}\lesssim 100$ cm-3) through absorption measurements of their
ground-state rotational lines111SH was first reported by IR spectroscopy
toward the cirumstellar envelope around the evolved star R Andromedae
(Yamamura et al., 2000).. In UV-illuminated gas, most sulfur atoms are
ionized, but the very high endothermicity of reaction

${\rm S^{+}\,(^{4}{\it{S}})+H_{2}\,(^{1}\Sigma^{+},\nu=0)\rightleftarrows
SH^{+}\,(^{3}\Sigma^{-})+H\,(^{2}{\it{S}})}$ (1)

($E\,/\,k$ = 9860 K, e.g., Zanchet et al., 2013a, 2019) prevents this reaction
from being efficient unless the gas is heated to very high temperatures. In
diffuse molecular clouds (on average at $T_{\rm k}$ $\sim$ 100 K), the
formation of SH+ and SH only seems possible in the context of local regions of
overheated gas subjected to magnetized shocks (Pineau des Forets et al., 1986)
or in dissipative vortices of the interstellar turbulent cascade (Godard et
al., 2012, 2014). In these tiny pockets ($\sim$100 AU in size), the gas would
attain the hot temperatures ($T_{\rm k}\,$$\simeq$ 1000 K) and/or ion-neutral
drift needed to overcome the endothermicities of the above hydrogen
abstraction reactions (see, e.g., Neufeld et al., 2015).

Dense PDRs ($n_{\rm H}$ $\simeq 10^{3}-10^{6}$ cm-3) offer a complementary
environment to study the first steps of sulfur chemistry. Because of their
higher densities and more quiescent gas, fast shocks or turbulence dissipation
do not contribute to the gas heating. Instead, the molecular gas is heated to
$T_{\rm k}\lesssim$ 500 K by mechanisms that depend on the flux of far-UV
photons (FUV; $E$ $<$13.6 eV). A different perspective of the H2 ($v$)
reactivity emerges because certain endoergic reactions become exoergic and
fast when a significant fraction of the H2 reagents are radiatively pumped to
vibrationally excited states $v\geq 1$ (Stecher & Williams, 1972; Freeman &
Williams, 1982; Tielens & Hollenbach, 1985; Sternberg & Dalgarno, 1995). In
this case, state-specific reaction rates for H2 ($\nu,J$) are needed to make
realistic predictions of the abundance of the product XH (Agúndez et al.,
2010; Zanchet et al., 2013b; Faure et al., 2017). The presence of abundant
FUV-pumped H2 ($v$ $\geq$ 1) triggers a nonthermal “hot” chemistry. Indeed,
CH+ and SH+ emission lines have been detected in the Orion Bar PDR (Nagy et
al., 2013; Goicoechea et al., 2017) where H2 lines up to $v$ = 10 have been
detected as well (Kaplan et al., 2017).

In this study we present a systematic (observational and modeling) study of
the chemistry of S-bearing hydrides in FUV-illuminated gas. We try to answer
the question of whether gas-phase reactions of S atoms and SH+ molecules with
vibrationally excited H2 can ultimately explain the presence of abundant H2S,
or if grain surface chemistry has to be invoked.

The paper is organized as follows. In Sects. 2 and 3 we report on new
observations of H${}_{2}^{32}$S, H${}_{2}^{34}$S, H${}_{2}^{33}$S, SH+, SH,
and H3S+ emission lines toward the Orion Bar. In Sect. 4 we study their
excitation and derive their column densities. In Sect. 6 we discuss their
abundances in the context of updated PDR models, with emphasis on the role of
hydrogen abstraction reactions

${\rm SH^{+}\,(^{3}\Sigma^{-})+H_{2}\,(^{1}\Sigma^{+})\rightleftarrows
H_{2}S^{+}\,(^{2}A^{\prime})+H\,(^{2}{\it{S}})},$ (2)
$\rm{H_{2}S^{+}\,(^{2}{\it{A^{\prime}}})+H_{2}\,(^{1}\Sigma^{+})\rightleftarrows
H_{3}S^{+}\,(X^{1}{\it{A_{\rm 1}}})+H\,(^{2}{\it{S}})},$ (3)
$\rm{S\,({\it{{}^{3}P}})+H_{2}\,(^{1}\Sigma^{+})\rightleftarrows
SH\,(X^{2}\Pi)+H\,({\it{{}^{2}S}})},$ (4)

photoreactions, and grain surface chemistry. In Sect. 5 we summarize the ab
initio quantum calculations we carried out to determine the state-dependent
rates of reactions (2) and (4). Details of these calculations are given in
Appendices A and B.

##  2 Observations of S-bearing hydrides

###  2.1 The Orion Bar

At an adopted distance of $\sim$414 pc, the Orion Bar is an interface of the
Orion molecular cloud and the Huygens H ii region that surrounds the Trapezium
cluster (Genzel & Stutzki, 1989; O’Dell, 2001; Bally, 2008; Goicoechea et al.,
2019, 2020; Pabst et al., 2019, 2020). The Orion Bar is a prototypical
strongly illuminated dense PDR. The impinging flux of stellar FUV photons
($G_{0}$) is a few 104 times the mean interstellar radiation field (Habing,
1968). The Bar is seen nearly edge-on with respect to the FUV illuminating
sources, mainly $\theta^{1}$ Ori C, the most massive star in the Trapezium.
This favorable orientation allows observers to spatially resolve the H+-to-H
transition (the ionization front or IF; see, e.g., Walmsley et al., 2000;
Pellegrini et al., 2009) from the H-to-H2 transition (the dissociation front
or DF; see, e.g., Allers et al., 2005; van der Werf et al., 1996, 2013;
Wyrowski et al., 1997; Cuadrado et al., 2019). It also allows one to study the
stratification of different molecular species as a function of cloud depth
(i.e., as the flux of FUV photons is attenuated; see, e.g., Tielens et al.,
1993; van der Wiel et al., 2009; Habart et al., 2010; Goicoechea et al., 2016;
Parikka et al., 2017; Andree-Labsch et al., 2017).

Regarding sulfur222Sulfur has four stable isotopes, in decreasing order of
abundance: 32S ($I_{\rm N}$ = 0), 34S ($I_{\rm N}$ = 0), 33S ($I_{\rm N}$ =
3/2), and 36S ($I_{\rm N}$ = 0), where $I_{\rm N}$ is the nuclear spin. The
most abundant isotope is here simply referred to as S., several studies
previously reported the detection of S-bearing molecules in the Orion Bar.
These include CS, C34S, SO, SO2, and H2S (Hogerheijde et al., 1995; Jansen et
al., 1995), SO+ (Fuente et al., 2003), C33S, HCS+, H2CS, and NS (Leurini et
al., 2006), and SH+ (Nagy et al., 2013). These detections refer to modest
angular resolution pointed observations using single-dish telescopes. Higher-
angular-resolution interferometric imaging of SH+, SO, and SO+ (Goicoechea et
al., 2017) was possible thanks to the Atacama Compact Array (ACA).

###  2.2 Observations of H2S isotopologues and H3S+

We observed the Orion Bar with the IRAM 30 m telescope at Pico Veleta (Spain).
We used the EMIR receivers in combination with the Fast Fourier Transform
Spectrometer (FTS) backends at 200 kHz resolution ($\sim$0.4 km s-1, $\sim$0.3
km s-1, and $\sim$0.2 km s-1 at $\sim$168 GHz, $\sim$217 GHz, and $\sim$293
GHz, respectively). These observations are part of a complete line survey
covering the frequency range 80 $-$ 360 GHz (Cuadrado et al., 2015, 2016,
2017, 2019) and include deep integrations at 168 GHz toward three positions of
the PDR located at a distance of 14′′, 40′′, and 65′′ from the IF (see Fig.
1). Their offsets with respect to the IF position at
$\mathrm{\alpha_{2000}=05^{h}\,35^{m}\,20.1^{s}\,}$,
$\mathrm{\delta_{2000}=-\,05^{\circ}25^{\prime}07.0^{\prime\prime}}$ are
(+10′′, -10′′), (+30′′, -30′′’), and (+35′′, -55′′). The first position is the
DF.

We carried out these observations in the position switching mode taking a
distant reference position at ($-$600′′, 0′′). The half power beam width
(HPBW) at $\sim$168 GHz, $\sim$217 GHz, and $\sim$293 GHz is $\sim$15′′,
$\sim$11′′, and $\sim$8′′, respectively. The latest observations (those at 168
GHz) were performed in March 2020. The data were first calibrated in the
antenna temperature scale $T^{*}_{\rm A}$ and then converted to the main beam
temperature scale, $T_{\rm mb}$, using $T_{\rm mb}$ = $T^{*}_{\rm
A}/\upeta_{\rm mb}$, where $\upeta_{\rm mb}$ is the antenna efficiency
($\upeta_{\rm mb}$ = 0.74 at $\sim$168 GHz). We reduced and analyzed the data
using the GILDAS software as described in Cuadrado et al. (2015). The typical
rms noise of the spectra is $\sim$3.5, 5.3, and 7.8 mK per velocity channel at
$\sim$168 GHz, $\sim$217 GHz, and $\sim$293 GHz, respectively. Figures 1 and 2
show the detection of $o$-H2S $1_{1,0}-1_{0,1}$ (168.7 GHz), $p$-H2S
$2_{2,0}-2_{1,1}$ (216.7 GHz), and $o$-H234S $1_{1,0}-1_{0,1}$ lines (167.9
GHz) (see Table 7 for the line parameters), as well as several $o$-H233S
$1_{1,0}-1_{0,1}$ hyperfine lines (168.3 GHz).

We complemented our dataset with higher frequency H2S lines detected by the
Herschel Space Observatory (Nagy et al., 2017) toward the “CO+ peak” position
(Stoerzer et al., 1995), which is located at only $\sim$4′′ from our DF
position (i.e., within the HPBW of these observations). These observations
were carried out with the HIFI receiver (de Graauw et al., 2010) at a
spectral-resolution of 1.1 MHz (0.7 km s-1 at 500 GHz). HIFI’s HPBW range from
$\sim$42′′ to $\sim$20′′ in the 500 - 1000 GHz window (Roelfsema et al.,
2012). The list of additional hydrogen sulfide lines detected by Herschel
includes the $o$-H2S $2_{2,1}-2_{1,2}$ (505.5 GHz), $2_{1,2}-1_{0,1}$ (736.0
GHz), and $3_{0,3}-2_{1,2}$ (993.1 GHz), as well as the $p$-H2S
$2_{0,2}-1_{1,1}$ (687.3 GHz) line. We used the line intensities, in the
$T_{\rm mb}$ scale, shown in Table A.1 of Nagy et al. (2017).

In order to get a global view of the Orion Bar, we also obtained 2.5′ $\times$
2.5′ maps of the region observed by us with the IRAM 30 m telescope using the
330 GHz EMIR receiver and the FTS backend at 200 kHz spectral-resolution
($\sim$0.2 km s-1). On-the-fly (OTF) scans were obtained along and
perpendicular to the Bar. The resulting spectra were gridded to a data cube
through convolution with a Gaussian kernel providing a final resolution of
$\sim$8′′. The total integration time was $\sim$6 h. The achieved rms noise is
$\sim$1 K per resolution channel. Figure 1 shows the spatial distribution of
the 13CO $J$=3-2 (330.5 GHz) and SO 89-78 (346.5 GHz) integrated line
intensities.

###  2.3 ALMA imaging of Orion Bar edge in SH+ emission

We carried out mosaics of a small field of the Orion Bar using twenty-seven
ALMA 12 m antennas in band 7 (at $\sim$346 GHz). These unpublished
observations belong to project 2012.1.00352.S (P.I.: J. R. Goicoechea) and
consisted of a 27-pointing mosaic centered at $\alpha$(2000) = 5h35m20.6s;
$\delta$(2000) = -05o25′20′′. The total field-of-view (FoV) is
58′′$\times$52′′ (shown in Fig. 1). The two hyperfine line components of the
SH+ $N_{J}$ = $1_{0}-0_{1}$ transition were observed with correlators
providing $\sim$500 kHz resolution (0.4 km s-1) over a 937.5 MHz bandwidth.
The total observation time with the ALMA 12 m array was $\sim$2h. In order to
recover the large-scale extended emission filtered out by the interferometer,
we used deep and fully sampled single-dish maps, obtained with the total-power
(TP) antennas at 19′′ resolution, as zero- and short-spacings. Data
calibration procedures and image synthesis steps are described in Goicoechea
et al. (2016). The synthesized beam is $\sim$1′′. This is a factor of $\sim$4
better than previous interferometric SH+ observations (Goicoechea et al.,
2017). Figure 3 shows the resulting image of the SH+ $1_{0}-0_{1}$ $F$ =
1/2-3/2 hyperfine emission line at 345.944 GHz. We rotated this image 37.5o
clockwise to bring the FUV illumination in the horizontal direction. The
typical rms noise of the final cube is $\sim$ 80 mK per velocity channel and
1′′-beam. As expected from their Einstein coefficients, the other $F$ =
1/2-1/2 hyperfine line component at 345.858 GHz is a factor of $\sim$2 fainter
(see Table 8) and the resulting image has low signal-to-noise (S/N).

We complemented the SH+ dataset with the higher frequency lines observed by
HIFI (Nagy et al., 2013, 2017) at $\sim$526 GHz and $\sim$683 GHz (upper
limit). These pointed observations have HPBWs of $\sim$41′′ and $\sim$32′′
respectively, thus they do not spatially resolve the SH+ emission. To
determine their beam coupling factors ($f_{\rm b}$), we smoothed the bigger
4′′-resolution ACA + TP SH+ image shown in Goicoechea et al. (2017) to the
different HIFI’s HPBWs. We obtain $f_{\rm b}$ $\simeq$ 0.4 at $\sim$526 GHz
and $f_{\rm b}$ $\simeq$ 0.6 at $\sim$683 GHz. The corrected intensities are
computed as $W_{\rm corr}$ = $W_{\rm HIFI}$ / $f_{\rm b}$. These correction
factors are only a factor of $\lesssim$ 2 lower than simply assuming uniform
SH+ emission from a 10′′ width filament.

###  2.4 SOFIA/GREAT search for SH emission

We finally used the GREAT receiver (Heyminck et al., 2012) on board the
Stratospheric Observatory For Infrared Astronomy (SOFIA; Young et al., 2012)
to search for the lowest-energy rotational lines of SH (${}^{2}\Pi_{3/2}$ $J$
= 3/2-1/2) at 1382.910 and 1383.241 GHz (e.g., Klisch et al., 1996; Martin-
Drumel et al., 2012). These lines lie in a frequency gap that Herschel/HIFI
could not observe from space. These SOFIA observations belong to project
$07\\_0115$ (P.I.: J. R. Goicoechea). The SH lines were searched on the lower
side band of 4GREAT band 3. We employed the 4GREAT/HFA frontends and 4GFFT
spectrometers as backends. The HPBW of SOFIA at 1.3 THz is $\sim$20′′, thus
comparable with IRAM 30 m/EMIR and Herschel/HIFI observations. We also
employed the total power mode with a reference position at ($-$600′′,0′′). The
original plan was to observe during two flights in November 2019 but due to
bad weather conditions, only $\sim$70 min of observations were carried out in
a single flight.

After calibration, data reduction included: removal of a first order spectral
baseline, dropping scans with problematic receiver response, rms weighted
average of the spectral scans, and calibration to $T_{\rm mb}$ intensity scale
($\upeta_{\rm mb}$ = 0.71). The final spectrum, smoothed to a velocity-
resolution of 1 km s-1 has a rms noise of $\sim$50 mK (shown in Fig. 4). Two
emission peaks are seen at the frequencies of the $\Lambda$-doublet lines.
Unfortunately, the achieved rms is not enough to assure the unambiguous
detection of each component of the doublet. Although the stacked spectrum does
display a single line (suggesting a tentative detection) the resulting line-
width ($\Delta$v $\simeq$ 7 km s-1) is a factor of $\sim$3 broader than
expected in the Orion Bar (see Table 9). Hence, this spectrum provides
stringent upper limits to the SH column density but deeper integrations would
be needed to confirm the detection.

##  3 Observational results

Figure 2: Detection of H233S (at $\sim$168.3 GHz) toward the DF position of
the Orion Bar. Red lines indicate hyperfine components. Blue lines show
interloping lines from 13CCH. The length of each line is proportional to the
transition line strength (taken from the Cologne Database for Molecular
Spectroscopy, CDMS; Endres et al., 2016). Figure 3: ALMA 1′′-resolution images
zooming into the edge of the Orion Bar in 12CO 3-2 (left panel, Goicoechea et
al., 2016) and SH+ 10-01 $F$ = 1/2-3/2 line (middle panel, integrated line
intensity). The right panel shows the H2 $v$ = 1–0 $S$(1) line (Walmsley et
al., 2000). We rotated these images (all showing the same FoV) with respect to
Fig. 1 to bring the FUV illuminating direction in the horizontal direction
(from the right). The circle shows the DF position targeted with SOFIA in SH
(20′′ beam) and with the IRAM 30m telescope in H2S and H3S+.

###  3.1 H232S, H234S, and H233S across the PDR

Figure 1 shows an expanded view of the Orion Bar in the 13CO ($J$ = 3-2)
emission. FUV radiation from the Trapezium stars comes from the upper-right
corner of the image. The FUV radiation field is attenuated in the direction
perpendicular to the Bar. The infrared H2 $v$ = 1–0 $S$(1) line emission
(white contours) delineates the position of the H-to-H2 transition, the DF.
Many molecular species, such as SO, specifically emit from deeper inside the
PDR where the flux of FUV photons has considerably decreased. In contrast,
H2S, and even its isotopologue H${}_{2}^{34}$S, show bright 11,0-10,1 line
emission toward the DF (right panels in Fig. 1; see also Jansen et al., 1995).
Rotationally excited H2S lines have been also detected toward this position
(Nagy et al., 2017), implying the presence of warm H2S close to the irradiated
cloud surface (i.e., at relatively low extinctions). The presence of
moderately large H2S column densities in the PDR is also demonstrated by the
unexpected detection of the rare isotopologue H${}_{2}^{33}$S toward the DF
(at the correct LSR velocity of the PDR: vLSR $\simeq$ 10.5 km s-1). Figure 2
shows the H${}_{2}^{33}$S 11,0-10,1 line and its hyperfine splittings
(produced by the 33S nuclear spin). To our knowledge, H${}_{2}^{33}$S lines
had only been reported toward the hot cores in Sgr B2 and Orion KL before
(Crockett et al., 2014).

The observed $o$-H2S/$o$-H${}_{2}^{34}$S 11,0-10,1 line intensity ratio toward
the DF is 15 $\pm$ 2, below the solar isotopic ratio of 32S/34S $=$ 23 (e.g.,
Anders & Grevesse, 1989). The observed ratio thus implies optically thick
$o$-H2S line emission at $\sim$168 GHz. However, the observed
$o$-H${}_{2}^{34}$S/$o$-H${}_{2}^{33}$S 11,0-10,1 intensity ratio is 6 $\pm$
1, thus compatible with the solar isotopic ratio (34S/33S $=$ 5.5) and with
H${}_{2}^{34}$S and H${}_{2}^{33}$S optically thin emission.

Figure 4: Search for the SH 2$\Pi_{3/2}$ $J$=5/2-3/2 doublet (at $\sim$1383
GHz) toward the DF position with SOFIA/GREAT. Vertical magenta lines indicate
the position of hyperfine splittings taken from CDMS.

###  3.2 SH+ emission from the PDR edge

Figure 3 zooms into a small field of the Bar edge. The ALMA image of the CO
$J$ = 3-2 line peak temperature was first presented by Goicoechea et al.
(2016). Because the CO $J$ = 3-2 emission is nearly thermalized and optically
thick from the DF to the molecular cloud interior, the line peak temperature
scale ($T_{\rm peak}$) is a good proxy of the gas temperature ($T_{\rm k}$
$\simeq$ $T_{\rm ex}$ $\simeq$ $T_{\rm peak}$). The CO image implies small
temperature variations around $T_{\rm k}$ $\simeq$ 200 K. The middle panel in
Fig. 3 shows the ALMA image of the SH+ $N_{J}$ = 10-01 $F$ = 1/2-3/2 hyperfine
line at 345.944 GHz. Compared to CO, the SH+ emission follows the edge of the
molecular PDR, akin to a filament of $\sim$10′′ width (for the spatial
distribution of other molecular ions, see, Goicoechea et al., 2017). The SH+
emission shows localized small-scale emission peaks (density or column density
enhancements) that match, or are very close to, the vibrationally excited H2
($v$ = 1-0) emission (Fig. 3). We note that while some H2 ($v$ = 1-0) emission
peaks likely coincide with gas density enhancements (e.g., Burton et al.,
1990), the region also shows extended emission from FUV-pumped H2 ($v$ = 2-1)
(van der Werf et al., 1996) that does not necessarily coincide with the H2
($v$ = 1-0) emission peaks.

###  3.3 Search for SH, H3S+, and H2S $\nu_{2}$ = 1 emission

We used SOFIA/GREAT to search for SH 2$\Pi_{3/2}$ $J$=5/2-3/2 lines toward the
DF (Fig. 4). This would have been the first time that interstellar SH
rotational lines were seen in emission. Unfortunately, the achieved rms of the
observation does not allow a definitive confirmation of these lines, so here
we will only discuss upper limits to the SH column density. The red, green,
and blue curves in Fig. 4 show radiative transfer models for $n_{\rm H}$ = 106
cm-3, $T_{\rm k}$ = 200 K, and different SH column densities (see Sect. 4 for
more details).

Our IRAM 30 m observations toward the DF neither resulted in a detection of
H3S+, a key gas-phase precursor of H2S. The $\sim$293.4 GHz spectrum around
the targeted H3S+ $1_{0}$-$0_{0}$ line is shown in Fig. 5. Again, the achieved
low rms allows us to provide a sensitive upper limit to the H3S+ column
density. This results in $N$(H3S+) =(5.5-7.5)$\times$1010 cm-2 (5$\sigma$)
assuming an excitation temperature range $T_{\rm ex}$ = 10-30 K and extended
emission. Given the bright H2S emission close to the edge of the Orion Bar,
and because H2S formation at the DF might be driven by very exoergic
processes, we also searched for the 11,0-10,1 line of vibrationally excited
H2S (in the bending mode $\nu_{2}$). The frequency of this line lies at
$\sim$181.4 GHz (Azzam et al., 2013), thus at the end of our 2 mm-band
observations of the DF (rms $\simeq$ 16 mK). However, we do not detect this
line either.

Figure 5: Search for H3S+ toward the Orion Bar with the IRAM 30 m telescope.
The blue curve shows the expected position of the line.

##  4 Coupled nonlocal excitation and chemistry

In this section we study the rotational excitation of the observed S-bearing
hydrides333Readers interested only in the chemistry of these species and in
depth-dependent PDR models could directly jump to Section 6.. We determine the
SH+, SH (upper limit), and H2S column densities in the Orion Bar, and the
“average” gas physical conditions in the sense that we search for the
combination of single $T_{\rm k}$, $n_{\rm H}$, and $N$ that better reproduces
the observed line intensities (so-called “single-slab” approach). In Sect. 6
we expand these excitation models to multi-slab calculations that take into
account the expected steep gradients in a PDR.

In the ISM, rotationally excited levels are typically populated by inelastic
collisions. However, the lifetime of very reactive molecules can be so short
that the details of their formation and destruction need to be taken into
account when determining how these levels are actually populated (Black,
1998). Reactive collisions (collisions that lead to a reaction and thus to
molecule destruction) influence the excitation of these species when their
timescales become comparable to those of nonreactive collisions. The lifetime
of reactive molecular ions observed in PDRs (e.g., Fuente et al., 2003; Nagy
et al., 2013; van der Tak et al., 2013; Goicoechea et al., 2017, 2019) can be
so short that they do not get thermalized by nonreactive collisions or by
absorption of the background radiation field (Black, 1998). In these cases, a
proper treatment of the molecule excitation requires including chemical
formation and destruction rates in the statistical equilibrium equations
(d$n_{i}$ / d$t$ = 0) that determine the level populations:

$\displaystyle\sum_{j>i}n_{j}\,A_{ji}+\sum_{j\neq
i}n_{i}\left(B_{ji}\,\bar{J}_{ji}+C_{ji}\right)+F_{i}=$ (5)
$\displaystyle=n_{i}\left(\sum_{j<i}A_{ij}+\sum_{j\neq
i}\left(B_{ij}\,\bar{J}_{ij}+C_{ij}\right)\,+\,D_{i}\right),$ (6)

where $n_{i}$ $[\rm cm^{-3}]$ is the population of rotational level $i$,
$A_{ij}$ and $B_{ij}$ are the Einstein coefficients for spontaneous and
induced emission, $C_{ij}$ $[\rm s^{-1}]$ is the rate of inelastic
collisions444We use the following inelastic collision rate coefficients
$\gamma_{ij}$:
$\bullet$ SH+– $e^{-}$, including hyperfine splittings (Hamilton et al.,
2018).
$\bullet$ SH+– $o$-H2 and $p$-H2, including hyperfine splittings (Dagdigian,
2019).
$\bullet$ SH+– H, including hyperfine splittings (Lique et al., 2020).
$\bullet$ $o$-H2S and $p$-H2S with $o$-H2 and $p$-H2 (Dagdigian, 2020).
$\bullet$ SH– He, including fine-structure splittings (Kłos et al., 2009).
($C_{ij}$ = $\sum_{k}\gamma_{ij,\,k}\,n_{k}$, where $\gamma_{ij,\,k}(T)$ $[\rm
cm^{3}s^{-1}]$ are the collisional rate coefficients and $k$ stands for H2, H,
and $e^{-}$), and $\bar{J}_{ij}$ is the mean intensity of the total radiation
field over the line profile. In these equations, $n_{i}\,D_{i}$ is the
destruction rate per unit volume of the molecule in level $i$, and $F_{i}$ its
formation rate per unit volume (both in $\rm cm^{-3}s^{-1}$). When state-to-
state formation rates are not available, and assuming that the destruction
rate is the same in every level ($D_{i}$ = $D$), one can use the total
destruction rate $D\,[{\rm s^{-1}}]$ ($=\,\sum_{k}n_{k}\,k_{k}(T)$ \+
photodestruction rate, where $k_{k}$ $[\rm cm^{3}s^{-1}]$ is the state-
averaged rate of the two-body chemical reaction with species $k$) and consider
that the level populations of the nascent molecule follow a Boltzmann
distribution at an effective formation temperature $T_{\rm form}$:

$F_{i}=F\,g_{i}\,e^{-E_{i}/kT_{\rm form}}\,/\,Q(T_{\rm form}).$ (7)

In this formalism, $F$ $[\rm cm^{-3}\,s^{-1}]$ is the state-averaged formation
rate per unit volume, $g_{i}$ the degeneracy of level $i$, and $Q(T_{\rm
form})$ is the partition function at $T_{\rm form}$ (van der Tak et al.,
2007).

This “formation pumping” formalism has been previously implemented in large
velocity gradient codes to treat, for example, the local excitation of the
very reactive ion CH+ (Nagy et al., 2013; Godard & Cernicharo, 2013; Zanchet
et al., 2013b; Faure et al., 2017). However, interstellar clouds are
inhomogeneous and gas velocity gradients are typically modest at small spatial
scales. This means that line photons can be absorbed and reemitted several
times before leaving the cloud. Here we implemented this formalism in a Monte
Carlo code that explicitly models the nonlocal behavior of the excitation and
radiative transfer problem (see Appendix of Goicoechea et al., 2006).

Although radiative pumping by dust continuum photons does not generally
dominate in PDRs, for completeness we also included radiative excitation by a
modified blackbody at a dust temperature of $\sim$50 K and a dust opacity
$\tau_{\lambda}$ = 0.03 (150/$\lambda[\upmu{\rm m}]$)1.6 (which reproduces the
observed intensity and wavelength dependence of the dust emission in the Bar;
Arab et al., 2012). The molecular gas fraction, $f$(H2) = 2$n$(H2)/$n_{\rm
H}$, is set to 2/3, where $n_{\rm H}$ = $n$(H) + 2$n$(H2) is the total density
of H nuclei. This choice is appropriate for the dissociation front and implies
$n$(H2) = $n$(H). As most electrons in the DF come from the ionization of
carbon atoms, the electron density $n_{e}$ is set to $n_{e}$ $\simeq$ $n$(C+)
= 1.4$\times$10-4 $n_{\rm H}$ (e.g., Cuadrado et al., 2019). For the inelastic
collisions with $o$-H2 and $p$-H2, we assumed that the H2 ortho-to-para (OTP)
ratio is thermalized to the gas temperature.

###  4.1 SH+ excitation and column density

We start by assuming that the main destruction pathway of SH+ are reactions
with H atoms and recombinations with electrons (see Sect. 6.1). Hence, the SH+
destruction rate is $D$ $\simeq$ $n_{e}$ $k_{\rm e}$($T$) + $n$(H) $k_{\rm
H}$($T$) (see Table 1 for the relevant chemical destruction rates). For
$T_{\rm k}$ = $T_{\rm e}$ = 200 K and $n_{\rm H}$ = 106 cm-3 (e.g., Goicoechea
et al., 2016) this implies $D$ $\simeq$ 10-4 s-1 (i.e., the lifetime of an SH+
molecule in the Bar is less than 3 h). At these temperatures and densities,
$D$ is about ten times smaller than the rate of radiative and inelastic
collisional transitions that depopulate the lowest-energy rotational levels of
SH+. Hence, formation pumping does not significantly alter the excitation of
the observed SH+ lines, but it does influence the population of higher-energy
levels. Formation pumping effects have been readily seen in CH+ because this
species is more reactive555CH+ is more reactive than SH+ because CH+ does
react with H2($v$=0) exothermically producing CH${}_{2}^{+}$ at $k$ =
1.2$\times$10-9 cm3 s-1 (Anicich, 2003) and also because reaction of CH+ with
H is faster, $k$ = 7.5$\times$10-10 cm3 s-1. and its rotationally excited
levels lie at higher-energy (i.e., their inelastic collision pumping rates are
slower, e.g., Zanchet et al., 2013b)

Figure 6: Non-LTE excitation models of SH+. The horizontal lines mark the
observed line intensities in the Orion Bar. Dotted curves are for a standard
model ($F=D=0$). Continuous curves are for a model that includes chemical
destruction by H atoms and $e^{-}$ (model $F,D$). Dashed lines are for a model
in which destruction rates are multiplied by ten (model $F,D$ $\times$10). The
vertical black line marks the best model.

Figure 6 shows results of several models: without formation pumping (dotted
curves for model “$F=D$ = 0”), adding formation pumping with SH+ destruction
by H and $e^{-}$ (continuous curves for model “$F,D$”), and using a factor of
ten higher SH+ destruction rates (simulating a dominant role of SH+
photodissociation or destruction by reactions with vibrationally excited H2;
dashed curves for model “$F,D$ $\times$10”). Since the formation of SH+ is
driven by reaction (1) when H2 molecules are in $v$ $\geq$ 2, here we adopted
$T_{\rm form}$ $\simeq$ $E$($v$ = 2, $J$ = 0) / $k$ $-$ 9860 K $\approx$ 2000
K. Because these are constant column density $N$(SH+) excitation and radiative
transfer models, we used a normalized formation rate $F$ = $\sum F_{i}$ that
assumes steady-state SH+ abundances consistent with the varying gas density in
each model. That is, $F$ = $\sum F_{i}$ = $x$(SH+) $n_{\rm H}$ $D$ $[\rm
cm^{-3}s^{-1}]$, where $x$ refers to the abundance with respect to H nuclei.

The detected SH+ rotational lines connect the fine-structure levels $N_{J}$ =
10-01 (345 GHz) and 12-01 (526 GHz). Upper limits also exist for the 11-01
(683 GHz) lines. SH+ critical densities ($n_{\rm cr}$ = $A_{ij}$ /
$\gamma_{ij}$) for inelastic collisions with H or H2 are of the same order and
equal to several 106 cm-3. As for many molecular ions (e.g., Desrousseaux et
al., 2021), SH+–H2 (and SH+–H) inelastic collisional rate coefficients4 are
large ($\gamma_{ij}$ $\gtrsim$ 10-10 cm3 s-1). Thus, collisions with H (at low
$A_{V}$) and H2 (at higher $A_{V}$) generally dominate over collisions with
electrons ($\gamma_{ij}$ of a few 10-7 cm3 s-1). At low densities (meaning
$n_{\rm H}$ $<$ $n_{\rm cr}$) formation pumping increases the population of
the higher-energy levels (and their $T_{\rm ex}$), but there are only minor
effects in the low-energy submillimeter lines. At high densities, $n_{\rm H}$
$>$ 107 cm-3, formation pumping with $T_{\rm form}$ = 2000 K produces lower
intensities in these lines because the lowest-energy levels ($E_{\rm u}/k$ $<$
$T_{\rm k}$ $<$ $T_{\rm form}$) are less populated.

The best fit to the observed lines in model F, D is for $N$(SH+) $\simeq$
1.1$\times$1013 cm-3, $n_{\rm H}$ $\simeq$ 3$\times$105 cm-3, and $T_{\rm k}$
$\simeq$ 200 K. This is shown by the vertical dotted line in Fig. 6. This
model is consistent with the upper limit intensity of the 683 GHz line (Nagy
et al., 2013). In this comparison, and following the morphology of the SH+
emission revealed by ALMA (Fig. 3), we corrected the line intensities of the
SH+ lines detected by Herschel/HIFI with the beam coupling factors discussed
in Sec. 2.3, The observed 12-01/10-01 line ratio ($R$ =
$W$(526.048)/$W$(345.944) $\simeq$ 2) is sensitive to the gas density. In
these models, $R$ is 1.1 for $n_{\rm H}$ = 105 cm-3 and 3.0 for $n_{\rm H}$
=106 cm-3. We note that $n_{\rm H}$ could be lower if SH+
formation/destruction rates were faster, as in the $F,D$ $\times$10 model.
This could happen if SH+ photodissociation or destruction reactions with
H2($v$ $\geq$2) were faster than reactions of SH+ with H atoms or with
electrons. In Sec. 6 we show that this is not the case.

###  4.2 SH excitation and column density

SH is a ${}^{2}\Pi$ open-shell radical with fine-structure,
$\Lambda$-doubling, and hyperfine splittings (e.g., Martin-Drumel et al.,
2012). However, the frequency separation of the SH ${}^{2}\Pi_{3/2}$ $J$ =
5/2-3/2 hyperfine components is too small to be spectrally resolved in
observations of the Orion Bar (see Fig. 4).The available rate coefficients for
inelastic collisions of SH with helium atoms do not resolve the hyperfine
splittings. Hence, we first determined line frequencies, level degeneracies,
and Einstein coefficients of an SH molecule without hyperfine structure. To do
this, we took the complete set of hyperfine levels tabulated in CDMS. Lacking
specific inelastic collision rate coefficients, we scaled the available SH– He
rates of Kłos et al. (2009) by the square root of the reduced mass ratios and
estimated the SH– H and SH– H2 collisional rates.

The scaled rate coefficients are about an order of magnitude smaller than
those of SH+. However, the chemical destruction rate of SH at the PDR edge
(reactions with H, photodissociation, and photoionization, see Sect. 6.1) is
also slower (we take the rates of SH–H reactive collisions from Zanchet et
al., 2019). We determine $D$ $\simeq$ 3$\times$10-6 s-1 for $n_{\rm H}$ = 106
cm-3, $T_{\rm k}$ = 200 K, and $A_{V}$ $\simeq$ 0.7 mag. Models in Fig. 7
include these chemical rates for $T_{\rm form}$ =$T_{\rm k}$ (a lower limit to
the unknown formation temperature). Formation pumping enhances the intensity
of the ${}^{2}\Pi_{3/2}$ $J$ = 5/2-3/2 ground-state lines by a few percent
only.

Figure 7: Non-LTE excitation models of SH emission lines targeted with
SOFIA/GREAT. Horizontal dashed lines refer to observational limits, assuming
extended emission (lower intensities) and for a 10′′ width emission filament
at the PDR surface (higher intensities).

To estimate the SH column density in the Orion Bar we compare with the upper
limit intensities of the SH lines targeted by SOFIA. If SH and SH+ arise from
roughly the same gas at similar physical conditions ($n_{\rm H}$ $\simeq$ 106
cm-3 and $T_{k}$ $\simeq$ 200 K) the best model column density is for $N$(SH)
$\leq$ (0.6-1.6)$\times$1014 cm-2. If densities were lower, around $n_{\rm H}$
$\simeq$ 105 cm-3, the upper limit $N$(SH) column densities will be a factor
ten higher.

Figure 8: Non-LTE excitation models for $o$-H2S and $p$-H2S. Thin horizontal
lines show the observed intensities assuming either extended emission (lower
limit) or emission that fills the 15′′ beam at 168.7 GHz. The vertical line
marks the best model, resulting in an OTP ratio of 2.9 $\pm$ 0.3.

###  4.3 H2S excitation and column density

H2S has a $X^{2}A$ ground electronic state and two nuclear spin symmetries
that we treat separately, $o$-H2S and $p$-H2S. Previous studies of the H2S
line excitation have used collisional rates coefficients scaled from those of
the H2O – H2 system. Dagdigian (2020) recently carried out specific
calculations of the cross sections of $o$-H2S and $p$-H2S inelastic collisions
with $o$–H2 and $p$-H2 at different temperatures. The behavior of the new and
the scaled rates is different and it depends on the H2 OTP ratio (e.g., on gas
temperature) because the collisional cross sections are different for
$o$-H2–H2S and $p$-H2–H2S systems. At the warm temperatures of the PDR,
collisions with $o$-H2 dominate, resulting in rate coefficients for the
$\sim$168 GHz $o$-H2S line that are a factor up to $\sim$2.5 smaller than
those scaled from H2O–H2.

H2S is not a reactive molecule. At the edge of the PDR its destruction is
driven by photodissociation. We determine that the radiative and collisional
pumping rates are typically a factor of $\sim$100 higher than $D$ $\approx$
2$\times$10-6 s-1 (for $n_{\rm H}$ = 106 cm-3, $T_{\rm k}$ = 200 K, $G_{0}$
$\simeq$104, and $A_{V}$ $\simeq$ 0.7 mag). Figure 8 shows non-LTE $o$-H2S and
$p$-H2S excitation and radiative transfer models. As H2S may have its
abundance peak deeper inside the PDR and display more extended emission than
SH+ (e.g., Sternberg & Dalgarno, 1995), we show results for $T_{\rm k}$ = 200
and 100 K. When comparing with the observed line intensities, we considered
either emission that fills all beams, or a correction that assumes that the
H2S emission only fills the 15′′ beam of the IRAM 30m telescope at 168 GHz.
The vertical dotted lines in Fig. 8 show the best model, $N$(H2S) =
$N$($o$-H2S)+$N$($p$-H2S) = 2.5$\times$1014 cm-2, with an OTP ratio of 2.9
$\pm$ 0.3, thus consistent with the high-temperature statistical ratio of 3/1
(see discussion at the end of Sect. 6.4). Models with lower densities, $n_{\rm
H}$ $\simeq$ 105 cm-3, show worse agreement, and would translate into even
higher $N$(H2S) of $\gtrsim$ 1015 cm-2. In either case, these calculations
imply large columns of warm H2S toward the PDR. They result in a limit to the
SH to H2S column density ratio of $\leq$ 0.2-0.6. This upper limit is already
lower than the $N$(SH)/$N$(H2S) = 1.1-3.0 ratios observed in diffuse clouds
(Neufeld et al., 2015). This difference suggests an enhanced H2S formation
mechanism in FUV-illuminated dense gas.

Figure 9: Minimum energy paths for reactions (1), (2), and (3). Points
correspond to RCCSD(T)-F12a calculations and lines to fits (Appendix A). The
reaction coordinate, $s$, is defined independently for each path. The
geometries of each species at $s$=0 are different.

##  5 New results on sulfur-hydride reactions

In this section we summarize the ab initio quantum calculations we carried out
to determine the vibrationally-state-dependent rates of gas-phase reactions of
H2($v$ $>$ 0) with several S-bearing species. We recall that all hydrogen
abstraction reactions,
${\rm S}^{+}\xrightarrow[(1)]{{{+{\rm H_{2}}}}}{\rm
SH}^{+}\xrightarrow[(2)]{{+{\rm H_{2}}}}{\rm
H_{2}S^{+}}\xrightarrow[(3)]{{+{\rm H_{2}}}}{\rm H_{3}S^{+},}\hskip
28.45274pt{\rm S}\xrightarrow[(4)]{{{+{\rm H_{2}}}}}{\rm SH,}\vspace{0.1cm}$
are very endoergic for H2 ($v$ = 0), with endothermicities in Kelvin units
that are significantly higher than $T_{\rm k}$ even in PDRs. This is markedly
different to O+ chemistry, for which all hydrogen abstraction reactions
leading to H3O+ are exothermic and fast (Gerin et al., 2010; Neufeld et al.,
2010; Hollenbach et al., 2012).

The endothermicity of reactions involving HnS+ ions decreases as the number of
hydrogen atoms increases. The potential energy surfaces (PES) of these
reactions possess shallow wells at the entrance and products channels (shown
in Fig. 9). In addition, these PESs show saddle points between the energy
walls of reactants and products whose heights increase with the number of H
atoms. For reaction (2), the saddle point has an energy of 0.6 eV
($\simeq$7,000 K) and is slightly below the energy of the products. However,
for reaction (3), the saddle point is above the energy of the products and is
a reaction barrier. These saddle points act as a bottleneck in the gas-phase
hydrogenation of S+.

If one considers the state dependent reactivity of vibrationally excited H2,
the formation of SH+ through reaction (1) becomes exoergic666If one considers
H2 rovibrational levels, reaction (1) becomes exoergic for $v$ = 0, $J$ $\geq$
11 and for $v$ = 1, $J$ $\geq$ 7 (Zanchet et al., 2019). when $v$ $\geq$ 2
(Zanchet et al., 2019). The detection of bright H2S emission in the Orion Bar
(Figs. 1 and 4) might suggest that subsequent hydrogen abstraction reactions
with H2 ($v$ $\geq$ 2) proceed as well. Motivated by these findings, and
before carrying out any PDR model, we studied reaction (2) and the reverse
process in detail. This required to build a full dimensional quantum PES of
the H3S+ (X1A1) system (see Appendix A).

In addition, we studied reaction (4) (and its reverse) through quantum
calculations. Details of these ab initio calculations and of the resulting
reactive cross sections are given in Appendix B. Table 1 summarizes the
updated reaction rate coefficients that we will include later in our PDR
models.

Table 1:  Relevant rate coefficients from a fit of the Arrhenius-like form
$k\,(T)$ = $\alpha\,(T/300\,{\rm K})^{\beta}\,{\rm exp}(-\gamma/T)$ to the
calculated reaction rates.

Reaction | $\alpha$ | $\beta$ | $\gamma$
---|---|---|---
| (cm3 s-1) |  | (K)
SH+ \+ H2 ($v$=1) $\rightarrow$ H2S+ \+ H | 4.97e-11 | 0 | 1973.4 a
SH+ \+ H2 ($v$=2) $\rightarrow$ H2S+ \+ H | 5.31e-10 | -0.17 | 0 a
SH+ \+ H2 ($v$=3) $\rightarrow$ H2S+ \+ H | 9.40e-10 | -0.16 | 0 a
SH+ \+ H $\rightarrow$ S+ \+ H2 | 1.86e-10 | -0.41 | 27.3 b
SH+ \+ $e^{-}$ $\rightarrow$ S + H | 2.00e-07 | -0.50 | c
H2S+ \+ H $\rightarrow$ SH+ \+ H2 | 6.15e-10 | -0.34 | 0 a
S + H2 ($v$=2) $\rightarrow$ SH + H |  $\sim$8.6e-13 |  $\sim$2.3 |  $\sim$2500 a
S + H2 ($v$=3) $\rightarrow$ SH + H |  $\sim$1.7e-12 |  $\sim$2.0 |  $\sim$1500 a
SH + H $\rightarrow$ S + H2 | 5.7e-13 | 2.48 | 1600a,†
| 7.7e-14 | 0.39 |  $-$1.3a,†
S+ \+ H2 ($v$=2) $\rightarrow$ SH+ \+ H | 2.88e-10 | -0.15 | 42.9 b
S+ \+ H2 ($v$=3) $\rightarrow$ SH+ \+ H | 9.03e-10 | -0.11 | 26.2 b
S+ \+ H2 ($v$=4) $\rightarrow$ SH+ \+ H | 1.30e-09 | -0.04 | 40.8 b
S+ \+ H2 ($v$=5) $\rightarrow$ SH+ \+ H | 1.21e-09 | 0.09 | 34.5 b
777(a) This work. (b) From Zanchet et al. (2019). (c) From Prasad & Huntress
(1980). †Total rate is the sum of the two expressions.

The H2S+ formation rate through reaction (2) with H2 ($v$ = 0) is very slow.
For H2 ($v$ = 1), the rate constant increases at $\approx$ 500 K,
corresponding to the opening of the H2S+ \+ H threshold. For H2 ($v$ = 2) and
H2 ($v$ = 3), the reaction rate is much faster, close to the Langevin limit
(see Appendix A.2). However, our estimated vibrational-state specific rates
for SH formation through reaction (4) (S + H2) are considerably smaller than
for reactions (1) and (2), and show an energy barrier even for H2 ($v$ = 2)
and H2 ($v$ = 3). We anticipate that this reaction is not a relevant formation
route for SH.

In FUV-illuminated environments, collisions with H atoms are very important
because they compete with electron recombinations in destroying molecular
ions, and also they contribute to their excitation. An important result of our
calculations is that the destruction rate of H2S+ (SH+) in reactions with H
atoms are a factor of $\geq$ 3.5 ($\geq$ 1.7) faster (at $T_{\rm k}$ $\leq$
200 K) than those previously used in astrochemical models (Millar et al.,
1986). Conversely, we find that destruction of SH in reactions with H atoms
(Appendix B) is slower than previously assumed.

##  6 PDR models of S-bearing hydrides

We now investigate the chemistry of S-bearing hydrides and the effect of the
new reaction rates in PDR models adapted to the Orion Bar conditions. In this
analysis we used version 1.5.4. of the Meudon PDR code (Le Petit et al., 2006;
Bron et al., 2014). Following our previous studies, we model the Orion Bar as
a stationary PDR at constant thermal-pressure (i.e., with density and
temperature gradients). When compared to time-dependent hydrodynamic PDR
models (e.g., Hosokawa & Inutsuka, 2006; Bron et al., 2018; Kirsanova & Wiebe,
2019), stationary isobaric models seem a good description of the most exposed
and compressed gas layers of the PDR, from $A_{V}$ $\approx$ 0.5 to $\approx$
5 mag (Goicoechea et al., 2016; Joblin et al., 2018).

In our models, the FUV radiation field incident at the PDR edge is $G_{0}$ =
2$\times$104 (e.g., Marconi et al., 1998). We adopted an extinction to color-
index ratio, $R_{V}$ = $A_{V}$/$E_{B-V}$, of 5.5 (Joblin et al., 2018),
consistent with the flatter extinction curve observed in Orion (Lee, 1968;
Cardelli et al., 1989). This choice implies slightly more penetration of FUV
radiation into the cloud (e.g., Goicoechea & Le Bourlot, 2007). The main input
parameters and elemental abundances of these PDR models are summarized in
Table 2. Figure 10 shows the resulting H2, H, and electron density profiles,
as well as the $T_{\rm k}$ and $T_{\rm d}$ gradients.

Our chemical network is that of the Meudon code updated with the new reaction
rates listed in Table 1. This network includes updated photoreaction rates
from Heays et al. (2017). To increase the accuracy of our abundance
predictions, we included the explicit integration of wavelength-dependent SH,
SH+, and H2S photodissociation cross sections ($\sigma_{\rm diss}$), as well
as SH and H2S photoionization cross sections ($\sigma_{\rm ion}$). These cross
sections are shown in Fig 23 of the Appendix. The integration is performed
over the specific FUV radiation field at each position of the PDR. In
particular, we took $\sigma_{\rm ion}$(SH) from Hrodmarsson et al. (2019) and
$\sigma_{\rm diss}$(H2S) from Zhou et al. (2020), both determined in
laboratory experiments. Figure 11 summarizes the relevant chemical network
that leads to the formation of S-bearing hydrides and that we discuss in the
following sections.

Table 2: Main parameters used in the PDR models of the Orion Bar.

Model parameter | Value | Note
---|---|---
FUV illumination, $G_{0}$ | 2$\times$104 Habing | $(a)$
Total depth $A_{\rm V}$ | 10 mag |
Thermal pressure $P_{\rm th}/k$ | 2$\times$108 cm-3K |
Density $n_{\rm H}$ = $n$(H) + 2$n$(H2) |  $n_{\rm H}$ = $P_{\rm th}\,/\,kT_{\rm k}$ | Varying
Cosmic Ray $\zeta_{\rm CR}$ | 10-16 H2 s-1 | $(b)$
$R_{\rm V}$ = $A_{\rm V}$/$E_{\rm B-V}$ | 5.5 | Orionc
$M_{\rm gas}/M_{\rm dust}$ | 100 | Local ISM
Abundance O / H | 3.2$\times$10-4 |
Abundance C / H | 1.4$\times$10-4 | Oriond
Abundance S / H | 1.4$\times$10-5 | Solare
888aMarconi et al. (1998). bIndriolo et al. (2015). cCardelli et al. (1989).
dSofia et al. (2004). eAsplund et al. (2009).

Figure 10: Structure of an isobaric PDR representing the most FUV-irradiated
gas layers of the Orion Bar (see Table 2 for the adopted parameters). This
plot shows the H2, H, and electron density profiles (left axis scale), and the
gas and dust temperatures (right axis scale). Figure 11: Main gas and grain
reactions leading to the formation of sulfur hydrides. Red arrows represent
endoergic reactions (endothermicity given in units of K). Dashed arrows are
uncertain radiative associations (see Sect. A.3), $\gamma$ stands for a FUV
photon, and “s-” for solid.

###  6.1 Pure gas-phase PDR model results

Figure 12 shows results of the “new gas-phase” model using the reaction rates
in Table 1. The continuous curves display the predicted fractional abundance
profiles as a function of cloud depth in magnitudes of visual extinction
($A_{V}$). The dashed curves are for a model that uses the standard thermal
rates previously adopted in the literature (see, e.g., Neufeld et al., 2015).
As noted by Zanchet et al. (2013a, 2019), the inclusion of H2 ($v$ $\geq$ 2)
state-dependent quantum rates for reaction (1) enhances the formation of SH+
in a narrow layer at the edge of the PDR ($A_{V}$ $\simeq$ 0 to 2 mag). This
agrees with the morphology of the SH+ emission revealed by ALMA images (Fig.
3). For H2 ($v$ $=$ 2), the reaction rate enhancement with respect to the
thermal rate $\Delta k$ = $k_{2}(T)/k_{0}(T)$ (see discussion by Agúndez et
al., 2010) is about 4$\times$108 at $T_{\rm k}$ = 500 K (Millar et al., 1986).
Indeed, when the fractional abundance of H2 ($v$ = 2) with respect to H2 ($v$
= 0), defined as $f_{\,2}$ = $n$ (H2 $v$ = 2)/$n$ (H2 $v$ = 0), exceeds a few
times 10-9, meaning $\Delta k\cdot f_{\,2}>1$, reaction (1) with H2 ($v$
$\geq$ 2) dominates SH+ formation. This reaction enhancement takes place only
at the edge of the PDR, where FUV-pumped H2 ($v$ $\geq$ 2) molecules are
abundant enough (gray dashed curves in Fig. 12) and drive the formation of
SH+. The resulting SH+ column density increases by an order of magnitude
compared to models that use the thermal rate.

In this isobaric model, the SH+ abundance peak occurs at $A_{V}$ $\simeq$ 0.7
mag, where the gas density has increased from $n_{\rm H}$ $\simeq$
6$\times$104 cm-3 at the PDR edge (the IF) to $\sim$5$\times$105 cm-3 (at the
DF). At this point, SH+ destruction is dominated by recombination with
electrons and by reactive collisions with H atoms. This implies $D$(SH+)
$[$s-1$]$ $\sim$ $n_{e}$ $k_{e}$ $\simeq$ $n_{\rm H}$ $k_{\rm H}$ $\gg$ $n$(H2
$v$$\geq$ 2) $k_{2}$, as we assumed in the single-slab SH+ excitation models
(Sec. 4.1). Therefore, only a small fraction of SH+ molecules further react
with H2 ($v$ $\geq$ 2) to form H2S+. The resulting low H2S+ abundances limit
the formation of abundant SH from dissociative recombinations of H2S+ (recall
that we estimated that reaction S + H2 ($v$ $\geq$2) $\rightarrow$ SH + H is
very slow). The SH abundance peak is shifted deeper inside the cloud, at about
$A_{V}$ $\simeq$ 1.8 mag, where SH forms by dissociative recombination of H2S+
and it is destroyed by FUV photons and reactions with H atoms. In these gas-
phase models the H2S abundance peaks even deeper inside the PDR, at AV
$\simeq$ 5 mag, where it forms by recombinations of H2S+ and H3S+ with
electrons as well as by charge exchange S + H2S+. However, the new rate of
reaction H2S+ \+ H is higher than assumed in the past, so the new models
predict lower H2S+ abundances at intermediate PDR depths (thus, less H3S+ and
H2S; see Fig. 12).

Table 3: Column density predictions from different PDR models (up to $A_{V}$ =
10 mag) and estimated values from observations (single-slab approach).

| log $N$ (cm-2) |  |  |  |
---|---|---|---|---|---
Type of PDR modela | SH+ | SH | H2S | H2S+ | H3S+
Standard gas-phase | 11.0a–12.2b | 11.4a–12.5b | 11.3a–12.4b | 9.9a–11.1b | 7.8a–9.0b
New gas-phase (Table 1) | 12.1a–13.2b | 11.4a–12.5b | 10.6a–11.7b | 9.9a–11.0b | 7.7a–8.9b
Gas-grain (low $E_{\rm b}$, $\epsilon$=1%) | 12.0a–13.2b | 13.2a–14.4b | 12.9a–14.1b | 9.6a–10.7b | 10.1a–11.2b
Gas-grain (high $E_{\rm b}$, $\epsilon$=1%) | 12.0a–13.1b | 13.6a–14.8b | 13.7b–14.8b | 9.9b–11.0b | 10.8b–12.0b
Estimated from observations |  $\sim$13.1 |  $<$13.8 |  $\sim$14.4 | – |  $<$ 10.7

999aColumn densities for a face-on PDR. bEdge-on PDR with a tilt angle
$\alpha$ = 4o, leading to the maximum expected geometrical enhancement.

Figure 12: Pure gas-phase PDR models of the Orion Bar. Continuous curves show
fractional abundances as a function of cloud depth, in logarithm scale to
better display the irradiated edge of the PDR, using the new reaction rates
listed in Table 1. The gray dotted curve shows $f_{2}$, the fraction of H2
that is in vibrationally excited levels $v$ $\geq$ 2 (right axis scale).
Dashed curves are for a model using standard reaction rates.

The SH column density predicted by the new gas-phase model is below the upper
limit determined from SOFIA. However, the predicted H2S column density is much
lower than the value we derive from observations (Table 3) and the predicted
H2S line intensities are too faint (see Sect. 6.4).

Because the cross sections of the different H2S photodissociation channels
have different wavelength dependences (Zhou et al., 2020), the H2S and SH
abundances between $A_{V}$ $\approx$ 2 and 6 mag are sensitive to the specific
shape of the FUV radiation field (determined by line blanketing, dust
absorption, and grain scattering; e.g., Goicoechea & Le Bourlot, 2007). Still,
we checked that using steeper extinction curves does not increase H2S column
density any closer to the observed levels. This disagreement between the
observationally inferred $N$(H2S) column density and the predictions of gas-
phase PDR models is even worse101010Older gas-phase PDR models previously
predicted low H2S column densities (Jansen et al., 1995; Sternberg & Dalgarno,
1995). if one considers the uncertain rates of radiative association reactions
S+ \+ H2 $\rightarrow$ H2S+ \+ $h\nu$ and SH+ \+ H2 $\rightarrow$ H3S+ \+
$h\nu$ included in the new gas-phase model. For the latter reaction, the main
problem is that the electronic states of the reactants do not correlate with
the ${}^{1}A_{1}$ ground electronic state of the activated complex H3S+∗
(denoted by $*$). Instead, H3S+∗ forms in an excited triplet state
(${}^{3}A$). Herbst et al. (1989) proposed that a spin-flip followed by a
radiative association can occur in interstellar conditions and form
H3S+∗($X^{1}A_{1}$) (Millar & Herbst, 1990). In Appendix A.3, we give
arguments against this mechanism. For similar reasons, Prasad & Huntress
(1982) avoided to include the S+ \+ H2 radiative association in their models.
Removing these reactions in pure gas-phase models drastically decreases the
H2S+ and H3S+ abundances, and thus those of SH and H2S (by a factor of
$\sim$100 in these models). The alternative H2S+ formation route through
reaction SH+ \+ H2($v$ = 2) is only efficient at the PDR surface ($A_{V}$ $<$
1 mag). This is due to the large H2($v$ = 2) fractional abundances, $f_{2}$
$>$ 10-6 at $T_{k}$ $>$ 500 K, required to enhance the H2S+ production.
Therefore, and contrary to S+ destruction, reaction of SH+ with H2 is not the
dominant destruction pathway for SH+. Only deeper inside the PDR, reactions of
S with H${}_{3}^{+}$ produce small abundances of SH+ and H2S+, but the
hydrogenation of HnS+ ions is not efficient and limits the gas-phase
production H2S.

###  6.2 Grain surface formation of solid H2S

Similarly to the formation of water ice (s-H2O) on grains (e.g., Hollenbach et
al., 2009, 2012), the formation of H2S may be dominated by grain surface
reactions followed by desorption back to the gas (e.g., Charnley, 1997).
Indeed, water vapor is relatively abundant in the Bar ($N$(H2O) $\approx$ 1015
cm-2; Choi et al., 2014; Putaud et al., 2019) and large-scale maps show that
the H2O abundance peaks close to cloud surfaces (Melnick et al., 2020).

To investigate the s-H2S formation on grains, we updated the chemical model by
allowing S atoms to deplete onto grains as the gas temperature drops inside
the molecular cloud (for the basic grain chemistry formalism, see, Hollenbach
et al., 2009). The timescale of this process ($\tau_{\rm gr,\,S}$) goes as
$x({\rm S})^{-1}$ $n_{\rm H}^{-1}\,T_{\rm k}^{-1/2}$, where $x$(S) is the
abundance of neutral sulfur atoms with respect to H nuclei. In a PDR, the
abundance of H atoms is typically higher than that of S atoms111111We only
consider the depletion of neutral S atoms. S+ ions are expected to be more
abundant than S atoms at the edge of the Orion Bar ($A_{V}$ $\lesssim$ 2 mag)
where $T_{\rm k}$ and $T_{\rm d}$ are too high, and the FUV radiation field
too strong, to allow the formation of abundant grain mantles. and H atoms
stick on grains more frequently than S atoms unless $x$(H) $<$
$x$(S)$\cdot$0.18. An adsorbed H atom (s-H) is weakly bound, mobile, and can
diffuse throughout the grain surface until it finds an adsorbed S atom (s-S).
If the timescale for a grain to be hit by a H atom ($\tau_{\rm gr,\,H}$) is
shorter that the timescale for a s-S atom to photodesorb ($\tau_{\rm
photdes,\,S}$) or sublimate ($\tau_{\rm subl,\,S}$) then reaction of s-H with
s-S will proceed and form a s-SH radical roughly upon “collision” and without
energy barriers (e.g., Tielens & Hagen, 1982; Tielens, 2010). Likewise, if
$\tau_{\rm gr,\,H}$ $<$ $\tau_{\rm photdes,\,SH}$ and $\tau_{\rm gr,\,H}$ $<$
$\tau_{\rm subl,\,SH}$, a newly adsorbed s-H atom can diffuse, find a grain
site with an s-SH radical and react without barriers to form s-H2S. In these
surface processes, a significant amount of S is ultimately transferred to
s-H2S (e.g., Vidal et al., 2017), which can subsequently desorb: thermally, by
FUV photons, or by cosmic rays. In addition, laboratory experiments show that
the excess energy of certain exothermic surface reactions can promote the
direct desorption of the product (Minissale et al., 2016). In particular,
reaction s-H + s-SH directly desorbs H2S with a maximum efficiency of $\sim$
60 % (as observed in experiments, Oba et al., 2018). Due to the high flux of
FUV photons in PDRs, chemical desorption may not always compete with
photodesorption. However, it can be a dominant process inside molecular clouds
Garrod et al. (2007); Esplugues et al. (2016); Vidal et al. (2017); Navarro-
Almaida et al. (2020).

The photodesorption timescale of an ice mantle is proportional to $Y^{-1}$
$G_{0}^{-1}$ exp (+$b$ $A_{V}$), where $Y$ is the photodesorption yield (the
number of desorbed atoms or molecules per incident photon) and $b$ is a dust-
related FUV field absorption factor. The timescale for mantle sublimation
(thermal desorption) goes as $\nu_{\rm ice}^{-1}$ exp (+$E_{\rm b}$ /
$k\,T_{\rm d}$), where $\nu_{\rm ice}$ is the characteristic vibrational
frequency of the solid lattice, $T_{\rm d}$ is the dust grain temperature, and
$E_{\rm b}/k$ is the adsorption binding energy of the species (in K). Binding
energies play a crucial role in model predictions because they determine the
freezing temperatures and sublimation timescales. Table 4 lists the $E_{\rm
b}/k$ and $Y$ values considered here.

Figure 13: Representative timescales relevant to the formation of s-H2S and
s-H2O as well as their freeze-out depths. Upper panel: The continuous black
curve is the timescale for a grain to be hit by an H atom. Once in the grain
surface, the H atom diffuses and can react with an adsorbed S atom to form
s-SH. The dashed magenta curves show the timescale for thermal desorption of
an s-S atom ($E_{\rm b}/k$ (S) = 1100 K left curve, and 2600 K right curve)
and of an s-O atom (blue curve; $E_{\rm b}/k$ (O) = 1800 K). The gray dotted
curve is the photodesorption timescale of s-S. At $G_{0}$ values where the
continuous line is below the dashed and dotted lines, s-O and s-S atoms remain
on grain surfaces sufficiently long to combine with an adsorbed H atom and
form s-OH and s-SH (and then s-H2O and s-H2S). These timescales are for
$n_{\rm H}$ = 105 cm-3 and $n$(H) = 100 cm-3. Bottom panel: Freeze-out depth
at which most O and S are incorporated as s-H2O and s-H2S (assuming no
chemical desorption and $T_{\rm k}$ = $T_{\rm d}$).

Representative timescales of the basic grain processes described above are
summarized in the upper panel of Fig. 13. In this plot, $T_{d}$ is a
characteristic dust temperature inside the PDR, $T_{d}$ = (3$\cdot$104 \+
2$\cdot$103 $G_{0}^{1.2}$)0.2, taken from Hollenbach et al. (2009). In the
upper panel, the continuous black curve is the timescale for a grain to be hit
by an H atom ($\tau_{\rm gr,\,H}$). The dashed magenta curves show the
timescale for thermal desorption of an s-S atom ($\tau_{\rm subl,\,S}$) (left
curve for $E_{\rm b}/k$ (S) = 1100 K and right curve for $E_{\rm b}/k$ (S) =
2600 K), and the same for an s-O atom (blue curve). The gray dotted curve is
the timescale for s-S atom photodesorption ($\tau_{\rm photodes,\,S}$) at
$A_{V}$ = 5 mag. At $G_{0}$ strengths where the continuous line is below the
dashed and dotted lines, an adsorbed s-S atom remains on the grain surface
sufficiently long to react with a diffusing s-H atom, form s-SH, and
ultimately s-H2S.

Figure 13 shows that, if one takes $E_{\rm b}/k$ (S) = 1100 K (the most common
value in the literature; Hasegawa & Herbst, 1993), the formation of s-H2S is
possible inside clouds illuminated by modest FUV fields, when grains are
sufficiently cold ($T_{\rm d}$ $<$ 22 K). However, recent calculations of s-S
atoms adsorbed on water ice surfaces suggest higher binding energies ($\sim$
2600 K; Wakelam et al., 2017). This would imply that S atoms freeze at higher
$T_{\rm d}$ ($\lesssim$ 50 K) and that s-H2S mantles form in more strongly
illuminated PDRs (the observed $T_{\rm d}$ at the edge of the Bar is $\simeq$
50 K and decreases to $\simeq$35 K behind the PDR; see, Arab et al., 2012).

The freeze-out depth for sulfur in a PDR, the $A_{V}$ at which most sulfur is
incorporated as S-bearing solids (s-H2S in our simple model) can be estimated
by equating $\tau_{\rm gr,\,S}$ and $\tau_{\rm photdes,\,H_{2}S}$. This
implicitly assumes the H2S chemical desorption does not dominate in FUV-
irradiated regions, which is in line with the particularly large FUV
absorption cross section of s-H2S measured in laboratory experiments (Cruz-
Diaz et al., 2014). With these assumptions, the lower panel of Fig. 13 shows
the predicted s-H2S and s-H2O freeze-out depths. Owing to the lower abundance
and higher atomic mass of sulfur atoms (i.e., grains are hit slower by S atoms
than by O atoms), the H2S freeze-out depth appears slightly deeper than that
of water ice. For the FUV-illumination conditions in the Bar, the freeze-out
depth of sulfur is expected at $A_{V}$ $\gtrsim$ 6 mag. This implies that
photodesorption of s-H2S can produce enhanced abundances of gaseous H2S at
$A_{V}$ $<$ 6 mag.

Table 4: Adopted binding energies and photodesorption yields.

Species | $E_{\rm b}/k$ | Yield
---|---|---
| (K) | (FUV photon)-1
S | 1100 a/2600 b | 10-4
SH | 1500 a/2700 b | 10-4
H2S | 2700 b,c | 1.2$\times$10-3 g (as H2S)
CO | 1300 d | 3$\times$10-3 h
O | 1800 e | 10-4 h
O2 | 1200 d | 10-3 h
OH | 4600 a | 10-3 h
H2O | 4800 f | 10-3 h (as H2O)
|  | 2$\times$10-3 h (as OH)
121212aHasegawa & Herbst (1993). bWakelam et al. (2017). cCollings et al.
(2004). dMinissale et al. (2016). eHe et al. (2015).fSandford & Allamandola
(1988). gFuente et al. (2017) hSee, Hollenbach et al. (2009).

Figure 14: Gas-grain PDR models leading to the formation of s-H2S (shown as
black curves). Continuous colored curves show gas-phase fractional abundances
as a function of depth into the cloud. $\epsilon$ refers to the efficiency of
the chemical desorption reaction s-H + s-H2S $\rightarrow$ SH + H2 (see text).
Left panel: Gas-grain high $E_{\rm b}$ model (high adsorption binding energies
for S and SH, see Table 4). Right panel: Low $E_{\rm b}$ model.

FUV-irradiation and thermal desorption of H2S ice mantles have been studied in
the laboratory (e.g., Cruz-Diaz et al., 2014; Jiménez-Escobar & Muñoz Caro,
2011). These experiments show that pure s-H2S ices thermally desorb around 82
K, and at higher temperatures for H2S–H2O ice mixtures. These experiments
determine a photodesorption yield of $Y_{\rm H_{2}S}$ $\sim\,$1.2$\times$10-3
molecules per FUV photon (see also Fuente et al., 2017). Regarding surface
grain chemistry, experiments show that reaction s-H + s-SH $\rightarrow$ s-H2S
is exothermic (Oba et al., 2018), whereas reaction s-H + s-H2S, although it
has an activation energy barrier of $\sim$1500 K, it may directly desorb
gaseous SH. Finally, reaction s-SH + s-SH $\rightarrow$ s-H2S2 may trigger the
formation of doubly sulfuretted species, but it requires mobile s-SH radicals
(e.g., Jiménez-Escobar & Muñoz Caro, 2011; Fuente et al., 2017). Here we will
only consider surface reactions with mobile s-H.

###  6.3 Gas-grain PDR model results

Here we show PDR model results in which we add a simple network of gas-grain
reactions for a small number of S-bearing (S, SH, and H2S) and O-bearing (O,
OH, H2O, O2, and CO) species. These species can adsorb on grains as
temperatures drop, photodesorb by FUV photons (stellar and secondary), desorb
by direct impact of cosmic-rays, or sublimate at a given PDR depth (depending
on $T_{\rm d}$ and on their $E_{\rm b}$). Grain size distributions ($n_{\rm
gr}$ $\propto$ $a^{-3.5}$, where $a$ is the grain radius) and gas-grain
reactions are treated within the Meudon code formalism (see, Le Petit et al.,
2006; Goicoechea & Le Bourlot, 2007; Le Bourlot et al., 2012; Bron et al.,
2014). As grain surface chemistry reactions we include s-H + s-X $\rightarrow$
s-XH and s-H + s-XH $\rightarrow$ s-H2X, where s-X refers to s-S and s-O. In
addition, we add the direct chemical desorption reaction s-H + s-SH
$\rightarrow$ H2S with an efficiency of 50 $\%$ per reactive event, and also
tested different efficiencies ($\epsilon$) for the chemical desorption process
s-H + s-H2S $\rightarrow$ SH + H2.

In our models we compute the relevant gas-grain timescales and atomic
abundances at every depth $A_{V}$ of the PDR. If the timescale for a grain to
be struck by an H atom ($\tau_{\rm gr,\,H}$) is shorter than the timescales to
sublimate or to photodesorb an s-X atom or a s-XH molecule; and if H atoms
stick on grains more frequently than X atoms, we simply assume these surface
reactions proceed instantaneously. At large $A_{V}$, larger than the freeze-
out depth, this grain chemistry builds abundant s-H2O and s-H2S ice mantles.

Figure 14 shows results of two types of gas-grain models. The only difference
between them is the adopted adsorption binding energies for s-S and s-SH. Left
panel is for a “high $E_{\rm b}$” model and right panel is for a “low $E_{\rm
b}$” model (see Table 4). We note that these models do not include the gas-
phase radiative association reactions S+ \+ H2 $\rightarrow$ H2S+ \+ $h\nu$
and SH+ \+ H2 $\rightarrow$ H3S+ \+ $h\nu$; although their effect is smaller
than in pure gas-phase models.

The chemistry of the most exposed PDR surface layers ($A_{V}$ $\lesssim$ 2
mag) is the same to that of the gas-phase models discussed in Sect. 6.1.
Photodesorption keeps dust grains free of ice mantles, and fast gas-phase ion-
neutral reactions, photoreactions, and reactions with FUV-pumped H2 drive the
chemistry. The resulting SH+ abundance profile is nearly identical and there
is no need to invoke depletion of elemental sulfur from the gas-phase to
explain the observed SH+ emission (see Fig. 15). Beyond these first PDR
irradiated layers, the chemistry does change because the formation of s-H2S on
grains and subsequent desorption alters the chemistry of the other S-bearing
hydrides.

Figure 15: Line intensity predictions for different isobaric PDR models.
Calculations were carried out in a multi-slab Monte Carlo code (Sect. 4) that
uses the output of the PDR model. Blue stars show the line intensities
observed toward the Bar (corrected by beam dilution). Left panel: SH+ emission
models for PDRs of different $P_{\rm th}$ values and $\alpha$ = 5o. Right
panel: SH and H2S (adopting an OTP ratio of 3) emission from: high $E_{\rm b}$
(magenta squares), low $E_{\rm b}$ (gray triangles), and gas-phase (cyan
circles) PDR models, all with $P_{\rm th}$ / $k$ = 2$\times$108 K cm-3. Upper
limit intensity predictions are for a PDR with an inclination angle of
$\alpha$ = 5o with respect to a edge-on geometry. Lower limit intensities
refer to a face-on PDR model.

In model high $E_{\rm b}$, S atoms start to freeze out closer to the PDR edge
($T_{\rm d}$ $<$ 50 K). Because of the increasing densities and decreasing
temperatures, the s-H2S abundance with respect to H nuclei reaches $\sim$10-6
at $A_{V}$ $\simeq$ 4 mag. In model low $E_{\rm b}$, this level of s-H2S
abundance is only reached beyond an $A_{V}$ of 7 mag. At lower $A_{V}$, the
formation of s-H2S on bare grains and subsequent photodesorption produces more
H2S than pure-gas phase models independently of whether H2S chemical
desorption is included or not. In these intermediate PDR layers, at $A_{V}$
$\simeq$ 2-7 mag for the strong irradiation conditions in the Bar, the flux of
FUV photons drives much of the chemistry, desorbing grain mantles, preventing
complete freeze out, and dissociating the gas-phase products.

There are two H2S abundance peaks at $A_{V}$ $\simeq$ 4 and 7 mag. The H2S
abundance in these “photodesorption peaks” depends on the amount of s-H2S
mantles formed on grains and on the balance between s-H2S photodesorption and
H2S photodissociation (which now becomes the major source of SH). The enhanced
H2S abundance modifies the chemistry of H2S+ and H3S+ as well: H2S
photoionization (with a threshold at $\sim$10.4 eV) becomes the dominant
source of H2S+ at $A_{V}$ $\simeq$ 4 mag because the H2 ($v$ $\geq$2)
abundance is too low to make reaction (2) competitive. Besides, reactions of
H2S with abundant molecular ions such as HCO+, H${}_{3}^{+}$, and H3O+
dominate the H3S+ production.

Our gas-grain models predict that other S-bearing molecules, such as SO2 and
SO, can be the major sulfur reservoirs at these intermediate PDR depths.
However, their abundances strongly depend on those of O2 and OH through
reactions S + O2 $\rightarrow$ SO + O and SO + OH $\rightarrow$ SO2 \+ H (see
e.g., Sternberg & Dalgarno, 1995; Fuente et al., 2016, 2019). These reactions
link the chemistry of S- and O-bearing neutral molecules (Prasad & Huntress,
1982) and are an important sink of S atoms at $A_{V}$ $\gtrsim$ 5 mag.
However, while large column densities of OH have been detected in the Orion
Bar ($\gtrsim$ 1015 cm-2; Goicoechea et al., 2011), O2 remains undetected
despite deep searches (Melnick et al., 2012). Furthermore, the inferred upper
limit $N$(O2) columns are below the expectations of PDR models (Hollenbach et
al., 2009). This discrepancy likely implies that these gas-grain models miss
details of the grain surface chemistry leading to O2 (for other environments
and modeling approaches see, e.g., Ioppolo et al., 2008; Taquet et al., 2016).
Here we will not discuss SO2, SO, or O2 further.

At large cloud depths, $A_{V}$ $\gtrsim$ 8 mag, the FUV flux is largely
attenuated, temperatures drop, the chemistry becomes slower, and other
chemical processes dominate. The H2S abundance is controlled by the chemical
desorption reaction s-H + s-SH $\rightarrow$ H2S. This process keeps a floor
of detectable H2S abundances ($>$10-9) in regions shielded from stellar FUV
radiation. In addition, and although not energetically favorable, the chemical
desorption s-H + s-H2S $\rightarrow$ SH + H2 enhances the SH production at
large $A_{V}$ (the enhancement depends on the desorption efficiency
$\epsilon$), which in turn boosts the abundances of other S-bearing species,
including that of neutral S atoms.

The H2S abundances predicted by the high $E_{\rm b}$ model reproduce the H2S
line intensities observed in the Bar (Sect. 6.4). In this model s-H2S becomes
the main sulfur reservoir. However, we stress that here we do not consider the
formation of more complex S-bearing ices such as s-OCS, s-H2S2, s-Sn, s-SO2 or
s-HSO (Jiménez-Escobar & Muñoz Caro, 2011; Vidal et al., 2017; Laas & Caselli,
2019). Together with our steady-state solution of the chemistry, this implies
that our predictions are not precise deep inside the PDR. However, we recall
that our observations refer to the edge of the Bar, so it is not plausible
that the model conditions at $A_{V}$ $\gtrsim$ 8 mag represent the line of
sight we observe.

Model low $E_{\rm b}$ produces less H2S in the PDR layers below $A_{V}$
$\lesssim$ 8 mag because S atoms do not freeze until the dust temperature
drops deep inside the PDR. Even beyond these layers, thermal desorption of s-S
maintains higher abundances of S atoms at large depths. Indeed, model low
$E_{\rm b}$ predicts that the major sulfur reservoir deep inside the cloud are
gas-phase S atoms. This agrees with recent chemical models of cold dark clouds
(Vidal et al., 2017; Navarro-Almaida et al., 2020).

Figure 16: Constant density gas-grain PDR models using the high $E_{\rm b}$
chemical network and undepleted sulfur elemental abundances. Left panel:
Effects of changing the FUV radiation field. Right panel: Effects of varying
the gas density.

###  6.4 Line intensity comparison and H2S ortho-to-para ratio

We now specifically compare the SH+, SH, and H2S line intensities implied by
the different PDR models with the intensities observed toward the DF position
of the Bar. We used the output of the PDR models – $T_{\rm k}$, $T_{\rm d}$,
$n$(H2), $n$(H), $n_{e}$, $n$(SH+), $n$(SH), and $n$(H2S) profiles from
$A_{V}$ = 0 to 10 mag – as input for a multi-slab Monte Carlo model of their
line excitation, including formation pumping (formalism presented in Sect. 4)
and radiative transfer. As the Orion Bar is not a perfectly edge-on, this
comparison requires a knowledge of the tilt angle ($\alpha$) with respect to a
pure edge-on PDR. Different studies suggest $\alpha$ of $\approx$ 5o (e.g.,
Jansen et al., 1995; Melnick et al., 2012; Andree-Labsch et al., 2017). This
inclination implies an increase in line-of-sight column density, compared to a
face-on PDR, by a geometrical factor (sin $\alpha$)-1. It also means that
optically thin lines are limb-brightened.

The left panel of Fig. 15 shows SH+ line intensity predictions for isobaric
PDR models of different $P_{\rm th}$ values (leading to different $T_{\rm k}$
and $n_{\rm H}$ profiles). Since the bulk of the SH+ emission arises from the
PDR edge ($A_{V}$ $\simeq$ 0 to 2 mag) all models (gas-phase or gas-grain)
give similar results. The best fit is for $P_{\rm th}$ $\simeq$
(1–2)$\times$108 cm-3 K and $\alpha$ $\simeq$5o. These high pressures, at
least close to the DF, agree with those inferred from ALMA images of HCO+ ($J$
= 4-3) emission (Goicoechea et al., 2016), Herschel observations of high-$J$
CO lines (Joblin et al., 2018), and IRAM 30 m detections of carbon
recombination lines (Cuadrado et al., 2019).

Right panel of Fig. 15 shows SH and H2S line emission predictions for the high
$E_{\rm b}$ gas-grain model (magenta squares), low $E_{\rm b}$ gas-grain model
(gray triangles), and a pure gas-phase model (cyan circles). For each model,
the upper limit intensities refer to radiative transfer calculations with an
inclination angle $\alpha$ = 5o. The lower intensity limits refer to a face-on
PDR. Gas-phase models largely underestimate the observed H2S intensities.
Model low $E_{\rm b}$ produces higher H2S columns and brighter H2S lines, but
still below the observed levels (by up to a factor of ten). Model high $E_{\rm
b}$ provides a good agreement with observations; the two possible inclinations
bracket the observed intensities, and it should be considered as the reference
model of the Bar. It is also consistent with the observational SH upper
limits.

Our observations and models provide a (line-of-sight)
$N$($o$-H2S)/$N$($p$-H2S) OTP ratio of 2.9 $\pm$ 0.3, consistent with the
(gas-phase) high-temperature statistical equilibrium value. However, the cold
“nuclear-spin-temperatures” ($T_{\rm spin}$ $\ll$ $T_{\rm k}$; see definition
in eq. 16) implied by the low water vapor OTP ratios observed in some sources
($<$ 2.5) have been associated with the temperature of the ice mantles where
H2O molecules might have formed (i.e., $T_{\rm spin}$ $\simeq$ $T_{\rm d}$;
Mumma et al., 1987; Lis et al., 2013). In the case of H2S, our derived OTP
ratio toward the DF position implies any $T_{\rm spin}$ above 30 $\pm$10 K
(see Fig. 24). Hence, this temperature might be also compatible with s-H2S
formation131313Crockett et al. (2014) inferred $N$($o$-H2S)/$N$($p$-H2S) = 2.5
$\pm$ 0.8 in the hot core of Orion KL using LTE rotational digrams. However,
they favored an OTP ratio of 1.7 $\pm$ 0.8 based on the column density ratio
of selected pairs of rotational levels with similar energies. This latter OTP
ratio implies $T_{\rm spin}$(H2S) $\simeq$ 12 K (Fig. 24), perhaps related to
much colder dust grains than in PDRs or to colder gas conditions just before
the hot core phase; so that reactive collisions did not have time to establish
the statistical equilibrium value. We note that the observed OTP ratios of
H2CO, H2CS, and H2CCO in the Bar are also $\sim$3 (Cuadrado et al., 2017). in
warm grains if $T_{\rm spin}$ $\simeq$ $T_{\rm d}$ upon formation is preserved
in the gas-phase after photodesorption (e.g., Guzmán et al., 2013).
Interestingly, the H2O OTP ratio derived from observations of the Orion Bar is
2.8 $\pm$ 0.1 (Putaud et al., 2019) and implies $T_{\rm spin}$(H2O) = 35 $\pm$
2 K. This value is compatible with $T_{\rm spin}$(H2S) and might reflect the
similar $T_{\rm d}$ of the PDR layers where most s-H2O and s-H2S form and
photodesorb. Nevertheless, laboratory experiments have challenged this $T_{\rm
spin}$ $\simeq$ $T_{\rm d}$ association, at least for s-H2O: cold water ice
surfaces, at 10 K, photodesorb H2O molecules with an OTP ratio of $\sim$3
(Hama et al., 2016). Follow up observations of $p$-H2S lines across the Bar
will allow us to study possible variations of the OTP ratio as $G_{0}$
diminishes and grains get colder.

###  6.5 Generalization to different G0 and nH conditions

In this section we generalize our results to a broader range of gas densities
and FUV illumination conditions (i.e., to clouds with different $G_{0}$ /
$n_{\rm H}$ ratios). We run several PDR models using the high $E_{\rm b}$ gas-
grain chemistry. The main difference compared to the Orion Bar models is that
here we model constant density clouds with standard interstellar grain
properties ($R_{V}$ = 3.1). Figure 16 (left panel) shows models of clouds with
constant $n_{\rm H}$ = 104 cm-3 and varying FUV radiation fields, while Fig.
16 (right panel) show models of constant FUV illumination ($G_{0}$ = 100) and
varying gas densities141414In these models we consider undepleted [S/H]
abundances and only the chemical desorption s-H + s-SH $\rightarrow$ H2S (with
a 50 % efficiency).. The main result of this study is the similar gas-phase
H2S column density (a few 1014 cm-2 up to $A_{V}$ = 10) and H2S abundance peak
(a few 10-8 close to the FUV-irradiated cloud edge) predicted by these models
nearly irrespective of $G_{0}$ and $n_{\rm H}$. A similar conclusion was
reached previously for water vapor in FUV-illuminated clouds (Hollenbach et
al., 2009, 2012). Increasing $G_{0}$ shifts the position of the H2S abundance
peak to larger $A_{V}$ until the rate of S atoms sticking on grains balances
the H2S photodissociation rate (the dominant H2S destruction mechanism except
in shielded gas; see also Fig.13). Since s-H2S photodesorption and H2S
photodissociation rates depend on $G_{0}$, the peak H2S abundance in the PDR
is roughly the same independently of $G_{0}$. On the other hand, the formation
rate of s-H2S mantles depends on the product $n$(S) $n_{\rm gr}$ $\propto$
$n_{\rm H}^{2}$, whereas the H2S photodesorption rate depends on $n_{\rm gr}$
$\propto$ $n_{\rm H}$. Hence, the H2S abundance peak moves toward the cloud
surface for denser PDRs (like the Orion Bar). The exact abundance value
depends on the adopted grain-size distribution and on the H2S photodesorption
yield (which is well constrained by experiments; see, Cruz-Diaz et al., 2014;
Fuente et al., 2017).

The role of chemical desorption increases and can dominate beyond the
photodesorption peak as the flux of stellar FUV photons is attenuated. Here we
do not carry out an exhaustive study of this mechanism, which is hard to model
in full detail because its efficiency decreases considerably with the
properties of grain surfaces (bare vs. icy; see e.g., Minissale & Dulieu,
2014). In our models, and depending on $\zeta_{\rm CR}$, photodesorption by
secondary FUV photons can also be important in cloud interiors. These
processes limit the conversion of most of the sulfur reservoir into S-bearing
ices and increase the abundance of other gas-phase species deep inside clouds,
notably S atoms and H2S molecules.

The H2S abundance in shielded gas depends on the destruction rate by gas-phase
reactions different than photodissociation, in particular H2S reactions with
H${}_{3}^{+}$. The H${}_{3}^{+}$ abundance increases with $\zeta_{\rm CR}$ and
decreases with the electron density. Figure 16 (right) shows models of
constant $G_{0}$ and constant $\zeta_{\rm CR}$ in which the H2S abundance at
large depths increases with decreasing density (more penetration of FUV
photons, more ionization, more electrons, less H${}_{3}^{+}$). The lowest gas
density model, $n_{\rm H}$ = 103 cm-3, shows the highest H2S abundance at
large $A_{V}$. Because S freeze-out is less efficient at low densities, the
low-density model shows higher gas-phase S abundances at large depths, making
atomic S a dominant gas-phase sulfur reservoir. Unfortunately, direct
observation of atomic S in cold gas is complicated, which makes it difficult
to benchmark this prediction.

In warm PDRs, in addition to S radio recombination lines (e.g., Smirnov et
al., 1995), the ${}^{3}P$ fine-structure lines of atomic sulfur, the [S i] 25,
56 $\upmu$m lines, can be interesting diagnostics of gas physical conditions
and of [S/H] abundances. Unfortunately, the low sensitivity of previous
infrared telescopes was not sufficient to detect the [S i] 25 $\upmu$m line
($\Delta E_{12}$ = 570 K) in the Orion Bar (Rosenthal et al., 2000); although
it is detected in protostellar outflows (e.g., Neufeld et al., 2009;
Goicoechea et al., 2012). Moreover, the ${}^{3}P_{2}$-${}^{1}D_{2}$ forbidden
line of atomic sulfur at 1.082 $\upmu$m can be an interesting tracer of the
ionization and dissociation fronts in PDRs. Some of these lines will be
accesible to high-angular-resolution and high sensitivity observations with
JWST.

####  6.5.1 The origin of H2S emission in other environments

Irrespective of $n_{\rm H}$ and $G_{0}$, grain surface formation of s-H2S and
photodesorption back to the gas-phase lead to H2S column densities of a few
1014 cm-2 in PDRs. This is in agrement with the observed column in the Bar
($G_{0}$ $\approx$ 104) as well as at the mildly illuminated rims of TMC-1 and
Barnard 1b clouds ($G_{0}$ $\approx$ 10; Navarro-Almaida et al., 2020). The
inferred H2S abundance in the shielded interior of these dark clouds ($A_{V}$
$>$ 10 mag) drops to a few 10-9, but the species clearly does not disappear
from the gas ($N$(H2S) of a few 1013 cm-2; Navarro-Almaida et al., 2020).
Interestingly, neither in the Bar the H2S line emission at $\sim$168 GHz
decreases much behind the PDR (Fig. 1) even if the flux of FUV photons is
largely attenuated compared to the irradiated PDR edge.

Despite oxygen is $\sim$25 times more abundant than sulfur, the H2O to H2S
column density ratio in the Orion Bar PDR is only about $\sim$ 5\. This
similarity must also reflect the higher abundances of CO compared to CS.
Furthermore, the H2S column density in cold cores is strikingly similar to
that of water vapor (Caselli et al., 2010, 2012). This coincidence points to a
more efficient desorption mechanism of s-H2S compared to s-H2O in gas shielded
from stellar FUV photons. Navarro-Almaida et al. (2020) argues that chemical
desorption is able to reproduce the observed H2S abundance floor if the
efficiency of this process diminishes as ice grain mantles get thicker inside
cold dense cores.

Turning back to warmer star-forming environments, our predicted H2S abundance
in FUV-illuminated gas is comparable to that observed toward many hot cores
($\sim$10-9-10-8; van der Tak et al., 2003; Herpin et al., 2009). In these
massive protostellar environments, thermal desorption of icy mantles, suddenly
heated to $T_{\rm d}$ $\gtrsim$ 100 K by the luminosity of the embedded
massive protostar, drives the H2S production. Early in their evolution, young
hot cores ($\lesssim$ 104 yr) can show even higher abundances of recently
desorbed H2S (before further chemical processing takes place in the gas-phase;
e.g., Charnley, 1997; Hatchell et al., 1998; Jiménez-Serra et al., 2012;
Esplugues et al., 2014). Indeed, Crockett et al. (2014) reports a gas-phase
H2S abundance of several 10-6 toward the hot core in Orion KL. This high value
likely reflects the minimum s-H2S abundance locked as s-H2S mantles just
before thermal desorption. In addition, the H2S abundance in the Orion Bar is
only slightly lower than that inferred in protostellar outflows (several
10-8). In these regions, fast shocks erode and sputter the grain mantles,
releasing a large fraction of their molecular content and activating a high-
temperature gas-phase chemistry that quickly reprocesses the gas (e.g.,
Holdship et al., 2019). All in all, it seems reasonable to conclude that
everywhere s-H2S grain mantles form, or already formed in a previous
evolutionary stage, emission lines from gas-phase H2S will be detectable.

In terms of its detectability with single-dish telescopes, H2S rotational
lines are bright in hot cores ($T_{\rm peak,\,168\,GHz}$ $\simeq$ 30 K in
Orion KL but $\simeq$ 1-3 K toward most hot cores; Tercero et al., 2010; van
der Tak et al., 2003; Herpin et al., 2009), in strongly irradiated PDRs
($\simeq$ 6 K, this work), and in lower-illumination PDRs such as the
Horsehead ($\simeq$ 1 K; Rivière-Marichalar et al., 2019). The H2S emission is
fainter toward cold dark clouds ($\simeq$ 0.2 K in TMC-1; Navarro-Almaida et
al., 2020) and protostellar outflows ($\simeq$ 0.6 K in L1157; Holdship et
al., 2019). These line intensity differences are mostly produced by different
gas physical conditions and not by enormous changes of the H2S abundance.

Finally, H2S is also detected outside the Milky Way (firstly by Heikkilä et
al., 1999). Lacking enough spatial-resolution it is more difficult to
determine the origin of the extragalactic H2S emission. The derived abundances
in starburst galaxies such as NGC 253 ($\sim$10-9; Martín et al., 2006) might
be interpreted as arising from a collection of spatially unresolved hot cores
(Martín et al., 2011). However, hot cores have low filling factors at star-
forming cloud scales. Our study suggests that much of this emission can arise
from (the most common) extended molecular gas illuminated by stellar FUV
radiation (e.g., Goicoechea et al., 2019).

##  7 Summary and conclusions

We carried out a self-consistent observational and modeling study of the
chemistry of S-bearing hydrides in FUV-illuminated gas. We obtained the
following results:
– ALMA images of the Orion Bar show that SH+ is confined to narrow gas layers
of the PDR edge, close to the H2 dissociation front. Pointed observations
carried out with the IRAM 30m telescope show bright H${}_{2}^{32}$S,
H${}_{2}^{34}$S, H${}_{2}^{33}$S emission toward the PDR (but no H3S+, a key
gas precursor of H2S) as well as behind the Bar, where the flux of FUV photons
is largely attenuated. SOFIA observations provide tight limits to the SH
emission.
– The SH+ line emission arises from a high-pressure gas component, $P_{\rm
th}$ $\simeq$ (1–2)$\times$108 cm-3 K, where SH+ ions are destroyed by
reactive collisions with H atoms and electrons (as most HnS+ ions do). We
derive $N$(SH+) $\simeq$ 1013 cm-2 and an abundance peak of several
$\sim$10-9. H2S shows larger column densities toward the PDR, $N$(H2S) =
$N$($o$-H2S) + $N$($p$-H2S) $\simeq$ 2.5$\times$1014 cm-2. Our tentative
detection of SH translates into an upper limit column density ratio
$N$(SH)/$N$(H2S) of $<$ 0.2-0.6, already lower than the ratio of 1.1-3.0
observed in low-density diffuse molecular clouds (Neufeld et al., 2015). This
implies an enhanced H2S production mechanism in FUV-illuminated dense gas.
– All gas-phase reactions X + H2($v$=0) $\rightarrow$ XH + H (with X = S+, S,
SH+, or H2S+) are highly endoergic. While reaction of FUV-pumped H2($v$ $\geq$
2) molecules with S+ ions becomes exoergic and explains the observed levels of
SH+, further reactions of H2($v$ $\geq$ 2) with SH+ or with neutral S atoms,
both reactions studied here through ab initio quantum calculations, do not
form enough H2S+ or H3S+ to ultimately produce abundant H2S. In particular,
pure gas-phase models underestimate the H2S column density observed in the
Orion Bar by more than two orders of magnitude. This implies that these models
miss the main H2S formation route. The disagreement is even worse as we favor,
after considering the potential energy surfaces of the H2S+∗ and H3S+∗
complexes, that the radiative associations S+ \+ H2 $\rightarrow$ H2S+ \+
$h\nu$ and SH+ \+ H2 $\rightarrow$ H3S+ \+ $h\nu$ may actually not occur or
possess slower rates than considered in the literature.
– To overcome these bottlenecks, we built PDR models that include a simple
network of gas-grain and grain surface reactions. The higher binding energies
of S and SH suggested by recent studies imply that bare grains start to grow
s-H2S mantles not far from the illuminated edges of molecular clouds. Indeed,
the observed $N$(H2S) in the Orion Bar can only be explained by the freeze-out
of S atoms, grain surface formation of s-H2S mantles, and subsequent
photodesorption back to the gas phase. The inferred H2S OTP ratio of 2.9 $\pm$
0.3 (equivalent to $T_{\rm spin}$ $\geq$ 30 K) is compatible with the high-
temperature statistical ratio as well as with warm grain surface formation if
$T_{\rm spin}$ $\simeq$ $T_{\rm d}$ and if $T_{\rm spin}$ is preserved in the
gas-phase after desorption.
– Comparing observations with chemical and excitation models, we conclude that
the SH+-emitting layers at the edge of the Orion Bar ($A_{V}$ $<$ 2 mag) are
charaterized by no or very little depletion of sulfur from the gas-phase. At
intermediate PDR depths ($A_{V}$ $<$ 8 mag) the observed H2S column densities
do not require depletion of elemental (cosmic) sulfur abundances either.
– We conclude that everywhere s-H2S grain mantles form (or formed) gas-phase
H2S will be present in detectable amounts. Independently of $n_{\rm H}$ and
$G_{0}$, FUV-illuminated clouds produce roughly the same H2S column density (a
few 1014 cm-2) and H2S peak abundances (a few 10-8). This agrees with the H2S
column densities derived in the Orion Bar and at the edges of mildly
illuminated clouds. Deep inside molecular clouds ($A_{V}$ $>$ 8 mag), H2S
still forms by direct chemical desorption and photodesorption by secondary FUV
photons. These processes alter the abundances of other S-bearing species and
makes difficult to predict the dominant sulfur reservoir in cloud interiors.

In this study we focused on S-bearing hydrides. Still, many subtle details
remain to be fully understood: radiative associations, electron
recombinations, and formation of multiply sulfuretted molecules. For example,
the low-temperature (Tk $<$ 1000 K) rates of the radiative and dielectronic
recombination of S+ used in PDR models may still be not accurate enough
(Badnell, 1991). In addition, the main ice-mantle sulfur reservoirs are not
fully constrained observationally. Thus, some of the narrative may be subject
to speculation. Similarly, reactions of S+ with abundant organic molecules
desorbed from grains (such as s-H2CO, not considered in our study) may
contribute to enhance the H2S+ abundance through gas-phase reactions (e.g., S+
\+ H2CO $\rightarrow$ H2S+ \+ CO; Prasad & Huntress, 1982). Future
observations of the abundance and freeze out depths of the key ice carriers
with JWST will clearly help in these fronts.

###### Acknowledgements.

We warmly thank Prof. György Lendvay for interesting discussions and for
sharing the codes related to their S(${}^{3}P$) + H2(${}^{1}\Sigma_{g}^{+}$,v)
PES. We thank Paul Dagdigian, François Lique, and Alexandre Faure for sharing
their H2S–H2, SH+–H, and SH+–e- inelastic collisional rate coefficients and
for interesting discussions in Grenoble and Salamanca. We thank Helgi
Hrodmarsson for sending his experimental SH photoionization cross section in
tabulated format. We finally thank our referee, John H. Black, for encouraging
and insightful suggestions. This paper makes use of the ALMA data
ADS/JAO.ALMA#2012.1.00352.S. ALMA is a partnership of ESO (representing its
member states), NSF (USA), and NINS (Japan), together with NRC (Canada), and
NSC and ASIAA (Taiwan), in cooperation with the Republic of Chile. The Joint
ALMA Observatory is operated by ESO, AUI/NRAO, and NAOJ. It also includes IRAM
30 m telescope observations. IRAM is supported by INSU/CNRS (France), MPG
(Germany), and IGN (Spain). We thank the staff at the IRAM 30m telescope and
the work of the USRA and NASA staff of the Armstrong Flight Research Center in
Palmdale and of the Ames Research Center in Mountain View (California), and
the Deutsches SOFIA Institut. We thank the Spanish MICIU for funding support
under grants AYA2016-75066-C2-2-P, AYA2017-85111-P, FIS2017-83473-C2
PID2019-106110GB-I00, and PID2019-106235GB-I00 and the French-Spanish
collaborative project PICS (PIC2017FR). We finally acknowledge computing time
at Finisterrae (CESGA) under RES grant ACCT-2019-3-0004.

## References

  * Aguado et al. (2010) Aguado, A., Barragan, P., Prosmiti, R., et al. 2010, J. Chem. Phys., 133, 024306
  * Aguado & Paniagua (1992) Aguado, A. & Paniagua, M. 1992, J. Chem. Phys., 96, 1265
  * Aguado et al. (2001) Aguado, A., Tablero, C., & Paniagua, M. 2001, Comput. Phys. Comm., 134, 97
  * Agúndez et al. (2010) Agúndez, M., Goicoechea, J. R., Cernicharo, J., Faure, A., & Roueff, E. 2010, ApJ, 713, 662
  * Agúndez & Wakelam (2013) Agúndez, M. & Wakelam, V. 2013, Chemical Reviews, 113, 8710
  * Allers et al. (2005) Allers, K. N., Jaffe, D. T., Lacy, J. H., Draine, B. T., & Richter, M. J. 2005, ApJ, 630, 368
  * Anders & Grevesse (1989) Anders, E. & Grevesse, N. 1989, Geochim. Cosmochim. Acta., 53, 197
  * Andree-Labsch et al. (2017) Andree-Labsch, S., Ossenkopf-Okada, V., & Röllig, M. 2017, A&A, 598, A2
  * Anicich (2003) Anicich, V. G. 2003, JPL Publication 03-19, 1-1194
  * Arab et al. (2012) Arab, H., Abergel, A., Habart, E., et al. 2012, A&A, 541, A19
  * Asplund et al. (2009) Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481
  * Azzam et al. (2013) Azzam, A. A. A., Yurchenko, S. N., Tennyson, J., Martin-Drumel, M.-A., & Pirali, O. 2013, J. Quant. Spec. Radiat. Transf., 130, 341
  * Badnell (1991) Badnell, N. R. 1991, ApJ, 379, 356
  * Bally (2008) Bally, J. 2008, Overview of the Orion Complex, ed. B. Reipurth, 459
  * Bañares et al. (2003) Bañares, L., Aoiz, F. J., Honvault, P., Bussery-Honvault, B., & Launay, J.-M. 2003, J. Chem. Phys., 118, 565
  * Bañares et al. (2004) Bañares, L., Aoiz, F. J., Honvault, P., & Launay, J.-M. 2004, J. Phys. Chem., 108, 1616
  * Black (1998) Black, J. H. 1998, Faraday Discussions, 109, 257
  * Bonnet & Rayez (1997) Bonnet, L. & Rayez, J.-C. 1997, Chem. Phys. Lett., 277, 183
  * Bonnet & Rayez (2004) Bonnet, L. & Rayez, J.-C. 2004, Chem. Phys. Lett., 397, 106
  * Brittain et al. (2020) Brittain, A., Coolbroth, K., & Boogert, A. 2020, in American Astronomical Society Meeting Abstracts, Vol. 236, American Astronomical Society Meeting Abstracts #236, 247.08
  * Bron et al. (2018) Bron, E., Agúndez, M., Goicoechea, J. R., & Cernicharo, J. 2018, ArXiv e-prints
  * Bron et al. (2014) Bron, E., Le Bourlot, J., & Le Petit, F. 2014, A&A, 569, A100
  * Buckinghan (1967) Buckinghan, A. D. 1967, Adv. Chem. Phys., 12, 107
  * Burton et al. (1990) Burton, M. G., Hollenbach, D. J., & Tielens, A. G. G. M. 1990, ApJ, 365, 620
  * Calmonte et al. (2016) Calmonte, U., Altwegg, K., Balsiger, H., et al. 2016, MNRAS, 462, S253
  * Cardelli et al. (1989) Cardelli, J. A., Clayton, G. C., & Mathis, J. S. 1989, ApJ, 345, 245
  * Caselli et al. (2012) Caselli, P., Keto, E., Bergin, E. A., et al. 2012, ApJ, 759, L37
  * Caselli et al. (2010) Caselli, P., Keto, E., Pagani, L., et al. 2010, A&A, 521, L29
  * Charnley (1997) Charnley, S. B. 1997, ApJ, 481, 396
  * Choi et al. (2014) Choi, Y., van der Tak, F. F. S., Bergin, E. A., & Plume, R. 2014, A&A, 572, L10
  * Collings et al. (2004) Collings, M. P., Anderson, M. A., Chen, R., et al. 2004, MNRAS, 354, 1133
  * Crockett et al. (2014) Crockett, N. R., Bergin, E. A., Neill, J. L., et al. 2014, ApJ, 781, 114
  * Cruz-Diaz et al. (2014) Cruz-Diaz, G. A., Muñoz Caro, G. M., Chen, Y. J., & Yih, T. S. 2014, A&A, 562, A119
  * Cuadrado et al. (2017) Cuadrado, S., Goicoechea, J. R., Cernicharo, J., et al. 2017, A&A, 603, A124
  * Cuadrado et al. (2015) Cuadrado, S., Goicoechea, J. R., Pilleri, P., et al. 2015, A&A, 575, A82
  * Cuadrado et al. (2016) Cuadrado, S., Goicoechea, J. R., Roncero, O., et al. 2016, A&A, 596, L1
  * Cuadrado et al. (2019) Cuadrado, S., Salas, P., Goicoechea, J. R., et al. 2019, A&A, 625, L3
  * Dagdigian (2019) Dagdigian, P. J. 2019, MNRAS, 487, 3427
  * Dagdigian (2019) Dagdigian, P. J. 2019, J. Chem. Phys., 150, 084308
  * Dagdigian (2020) Dagdigian, P. J. 2020, MNRAS, 494, 5239
  * Dartois (2005) Dartois, E. 2005, Space Sci. Rev., 119, 293
  * Davidson (1975) Davidson, E. R. 1975, J. Comp. Phys., 17, 87
  * de Graauw et al. (2010) de Graauw, T., Helmich, F. P., Phillips, T. G., et al. 2010, A&A, 518, L6
  * Desrousseaux et al. (2021) Desrousseaux, B., Lique, F., Goicoechea, J. R., Quintas-Sánchez, E., & Dawes, R. 2021, A&A, 645, A8
  * Endres et al. (2016) Endres, C. P., Schlemmer, S., Schilke, P., Stutzki, J., & Müller, H. S. P. 2016, Journal of Molecular Spectroscopy, 327, 95
  * Esplugues et al. (2016) Esplugues, G. B., Cazaux, S., Meijerink, R., Spaans, M., & Caselli, P. 2016, A&A, 591, A52
  * Esplugues et al. (2014) Esplugues, G. B., Viti, S., Goicoechea, J. R., & Cernicharo, J. 2014, A&A, 567, A95
  * Farah et al. (2012) Farah, K., Muller-Plathe, F., & Bohm, M. C. 2012, Chem. Phys. Chem., 13, 1127
  * Faure et al. (2017) Faure, A., Halvick, P., Stoecklin, T., et al. 2017, MNRAS, 469, 612
  * Freeman & Williams (1982) Freeman, A. & Williams, D. A. 1982, Ap&SS, 83, 417
  * Fuente et al. (2016) Fuente, A., Cernicharo, J., Roueff, E., et al. 2016, A&A, 593, A94
  * Fuente et al. (2017) Fuente, A., Goicoechea, J. R., Pety, J., et al. 2017, ApJ, 851, L49
  * Fuente et al. (2019) Fuente, A., Navarro, D. G., Caselli, P., et al. 2019, A&A, 624, A105
  * Fuente et al. (2003) Fuente, A., Rodrıguez-Franco, A., Garcıa-Burillo, S., Martın-Pintado, J., & Black, J. H. 2003, A&A, 406, 899
  * Garrod et al. (2007) Garrod, R. T., Wakelam, V., & Herbst, E. 2007, A&A, 467, 1103
  * Genzel & Stutzki (1989) Genzel, R. & Stutzki, J. 1989, ARA&A, 27, 41
  * Gerin et al. (2010) Gerin, M., de Luca, M., Black, J., et al. 2010, A&A, 518, L110
  * Gerin et al. (2016) Gerin, M., Neufeld, D. A., & Goicoechea, J. R. 2016, ARA&A, 54, 181
  * Gibb et al. (2004) Gibb, E. L., Whittet, D. C. B., Boogert, A. C. A., & Tielens, A. G. G. M. 2004, ApJS, 151, 35
  * Godard & Cernicharo (2013) Godard, B. & Cernicharo, J. 2013, A&A, 550, A8
  * Godard et al. (2012) Godard, B., Falgarone, E., Gerin, M., et al. 2012, A&A, 540, A87
  * Godard et al. (2014) Godard, B., Falgarone, E., & Pineau des Forêts, G. 2014, A&A, 570, A27
  * Goicoechea et al. (2012) Goicoechea, J. R., Cernicharo, J., Karska, A., et al. 2012, A&A, 548, A77
  * Goicoechea et al. (2017) Goicoechea, J. R., Cuadrado, S., Pety, J., et al. 2017, A&A, 601, L9
  * Goicoechea et al. (2011) Goicoechea, J. R., Joblin, C., Contursi, A., et al. 2011, A&A, 530, L16
  * Goicoechea & Le Bourlot (2007) Goicoechea, J. R. & Le Bourlot, J. 2007, A&A, 467, 1
  * Goicoechea et al. (2020) Goicoechea, J. R., Pabst, C. H. M., Kabanovic, S., et al. 2020, A&A, 639, A1
  * Goicoechea et al. (2016) Goicoechea, J. R., Pety, J., Cuadrado, S., et al. 2016, Nature, 537, 207
  * Goicoechea et al. (2009) Goicoechea, J. R., Pety, J., Gerin, M., Hily-Blant, P., & Le Bourlot, J. 2009, A&A, 498, 771
  * Goicoechea et al. (2006) Goicoechea, J. R., Pety, J., Gerin, M., et al. 2006, A&A, 456, 565
  * Goicoechea et al. (2019) Goicoechea, J. R., Santa-Maria, M. G., Bron, E., et al. 2019, A&A, 622, A91
  * Gómez-Carrasco & Roncero (2006) Gómez-Carrasco, S. & Roncero, O. 2006, J. Chem. Phys., 125, 054102
  * Graedel et al. (1982) Graedel, T. E., Langer, W. D., & Frerking, M. A. 1982, ApJS, 48, 321
  * Grozdanov & Solov’ev (1982) Grozdanov, T. P. & Solov’ev, E. A. 1982, J. Phys. B, 15, 1195
  * Guzmán et al. (2013) Guzmán, V. V., Goicoechea, J. R., Pety, J., et al. 2013, A&A, 560, A73
  * Habart et al. (2010) Habart, E., Dartois, E., Abergel, A., et al. 2010, A&A, 518, L116
  * Habing (1968) Habing, H. J. 1968, Bull. Astron. Inst. Netherlands, 19, 421
  * Hama et al. (2016) Hama, T., Kouchi, A., & Watanabe, N. 2016, Science, 351, 65
  * Hamilton et al. (2018) Hamilton, J. R., Faure, A., & Tennyson, J. 2018, MNRAS, 476, 2931
  * Hasegawa & Herbst (1993) Hasegawa, T. I. & Herbst, E. 1993, MNRAS, 261, 83
  * Hatchell et al. (1998) Hatchell, J., Thompson, M. A., Millar, T. J., & MacDonald, G. H. 1998, A&A, 338, 713
  * He et al. (2015) He, J., Shi, J., Hopkins, T., Vidali, G., & Kaufman, M. J. 2015, ApJ, 801, 120
  * Heays et al. (2017) Heays, A. N., Bosman, A. D., & van Dishoeck, E. F. 2017, A&A, 602, A105
  * Heikkilä et al. (1999) Heikkilä, A., Johansson, L. E. B., & Olofsson, H. 1999, A&A, 344, 817
  * Herbst et al. (1989) Herbst, E., DeFrees, D. J., & Koch, W. 1989, Mon. Not. R. Astyr. Soc., 237, 1057
  * Herpin et al. (2009) Herpin, F., Marseille, M., Wakelam, V., Bontemps, S., & Lis, D. C. 2009, A&A, 504, 853
  * Heyminck et al. (2012) Heyminck, S., Graf, U. U., Güsten, R., et al. 2012, A&A, 542, L1
  * Hogerheijde et al. (1995) Hogerheijde, M. R., Jansen, D. J., & van Dishoeck, E. F. 1995, A&A, 294, 792
  * Holdship et al. (2019) Holdship, J., Jimenez-Serra, I., Viti, S., et al. 2019, ApJ, 878, 64
  * Hollenbach et al. (2009) Hollenbach, D., Kaufman, M. J., Bergin, E. A., & Melnick, G. J. 2009, ApJ, 690, 1497
  * Hollenbach et al. (2012) Hollenbach, D., Kaufman, M. J., Neufeld, D., Wolfire, M., & Goicoechea, J. R. 2012, ApJ, 754, 105
  * Hollenbach & Tielens (1997) Hollenbach, D. J. & Tielens, A. G. G. M. 1997, ARA&A, 35, 179
  * Hosokawa & Inutsuka (2006) Hosokawa, T. & Inutsuka, S.-i. 2006, ApJ, 646, 240
  * Howk et al. (2006) Howk, J. C., Sembach, K. R., & Savage, B. D. 2006, ApJ, 637, 333
  * Hrodmarsson et al. (2019) Hrodmarsson, H. R., Garcia, G. A., Nahon, L., Loison, J.-C., & Gans, B. 2019, Physical Chemistry Chemical Physics (Incorporating Faraday Transactions), 21, 25907
  * Indriolo et al. (2015) Indriolo, N., Neufeld, D. A., Gerin, M., et al. 2015, ApJ, 800, 40
  * Ioppolo et al. (2008) Ioppolo, S., Cuppen, H. M., Romanzin, C., van Dishoeck, E. F., & Linnartz, H. 2008, ApJ, 686, 1474
  * Jansen et al. (1995) Jansen, D. J., Spaans, M., Hogerheijde, M. R., & van Dishoeck, E. F. 1995, A&A, 303, 541
  * Jiménez-Escobar & Muñoz Caro (2011) Jiménez-Escobar, A. & Muñoz Caro, G. M. 2011, A&A, 536, A91
  * Jiménez-Serra et al. (2012) Jiménez-Serra, I., Zhang, Q., Viti, S., Martín-Pintado, J., & de Wit, W. J. 2012, ApJ, 753, 34
  * Joblin et al. (2018) Joblin, C., Bron, E., Pinto, C., et al. 2018, A&A, 615, A129
  * Johnson (1987) Johnson, B. R. 1987, J. Chem. Phys., 86, 1445
  * Kaplan et al. (2017) Kaplan, K. F., Dinerstein, H. L., Oh, H., et al. 2017, ApJ, 838, 152
  * Karplus et al. (1965) Karplus, M., Porter, R. N., & Sharma, R. D. 1965, J. Chem. Phys., 43, 3259
  * Kirsanova & Wiebe (2019) Kirsanova, M. S. & Wiebe, D. S. 2019, MNRAS, 486, 2525
  * Klisch et al. (1996) Klisch, E., Klaus, T., Belov, S. P., et al. 1996, ApJ, 473, 1118
  * Kłos et al. (2009) Kłos, J., Lique, F., & Alexander, M. H. 2009, Chemical Physics Letters, 476, 135
  * Knizia et al. (2009) Knizia, G., Adler, T. B., & Werner, H. J. 2009, J. Chem. Phys., 130, 054104
  * Laas & Caselli (2019) Laas, J. C. & Caselli, P. 2019, A&A, 624, A108
  * Le Bourlot et al. (2012) Le Bourlot, J., Le Petit, F., Pinto, C., Roueff, E., & Roy, F. 2012, A&A, 541, A76
  * Le Petit et al. (2006) Le Petit, F., Nehmé, C., Le Bourlot, J., & Roueff, E. 2006, ApJS, 164, 506
  * Lee (1968) Lee, T. A. 1968, ApJ, 152, 913
  * Leurini et al. (2006) Leurini, S., Rolffs, R., Thorwirth, S., et al. 2006, A&A, 454, L47
  * Levine & Bernstein (1987) Levine, R. D. & Bernstein, R. B. 1987, Molecular Reaction Dynamics and Chemical Reactivity (Oxford University Press)
  * Lique et al. (2020) Lique, F., Zanchet, A., Bulut, N., Goicoechea, J. R., & Roncero, O. 2020, A&A, 638, A72
  * Lis et al. (2013) Lis, D. C., Bergin, E. A., Schilke, P., & van Dishoeck, E. F. 2013, Journal of Physical Chemistry A, 117, 9661
  * Lucas & Liszt (2002) Lucas, R. & Liszt, H. S. 2002, A&A, 384, 1054
  * Maiti et al. (2004) Maiti, B., Schatz, G. C., & Lendvay, G. 2004, Journal of Physical Chemistry A, 108, 8772
  * Marconi et al. (1998) Marconi, A., Testi, L., Natta, A., & Walmsley, C. M. 1998, A&A, 330, 696
  * Martín et al. (2011) Martín, S., Krips, M., Martín-Pintado, J., et al. 2011, A&A, 527, A36
  * Martín et al. (2006) Martín, S., Mauersberger, R., Martín-Pintado, J., Henkel, C., & García-Burillo, S. 2006, ApJS, 164, 450
  * Martin-Drumel et al. (2012) Martin-Drumel, M. A., Eliet, S., Pirali, O., et al. 2012, Chemical Physics Letters, 550, 8
  * Melnick et al. (2012) Melnick, G. J., Tolls, V., Goldsmith, P. F., et al. 2012, ApJ, 752, 26
  * Melnick et al. (2020) Melnick, G. J., Tolls, V., Snell, R. L., et al. 2020, ApJ, 892, 22
  * Menten et al. (2011) Menten, K. M., Wyrowski, F., Belloche, A., et al. 2011, A&A, 525, A77
  * Millar et al. (1986) Millar, T. J., Adams, N. G., Smith, D., Lindinger, W., & Villinger, H. 1986, MNRAS, 221, 673
  * Millar & Herbst (1990) Millar, T. J. & Herbst, E. 1990, Astron. Astrophys., 231, 466
  * Millar & Herbst (1990) Millar, T. J. & Herbst, E. 1990, A&A, 231, 466
  * Minissale & Dulieu (2014) Minissale, M. & Dulieu, F. 2014, J. Chem. Phys., 141, 014304
  * Minissale et al. (2016) Minissale, M., Dulieu, F., Cazaux, S., & Hocuk, S. 2016, A&A, 585, A24
  * Mumma et al. (1987) Mumma, M. J., Weaver, H. A., & Larson, H. P. 1987, A&A, 187, 419
  * Nagy & Lendvay (2017) Nagy, T. & Lendvay, G. 2017, J. Phys. Chem. Lett., 8, 4621
  * Nagy et al. (2017) Nagy, Z., Choi, Y., Ossenkopf-Okada, V., et al. 2017, A&A, 599, A22
  * Nagy et al. (2013) Nagy, Z., Van der Tak, F. F. S., Ossenkopf, V., et al. 2013, A&A, 550, A96
  * Navarro-Almaida et al. (2020) Navarro-Almaida, D., Le Gal, R., Fuente, A., et al. 2020, A&A, 637, A39
  * Neufeld et al. (2012) Neufeld, D. A., Falgarone, E., Gerin, M., et al. 2012, A&A, 542, L6
  * Neufeld et al. (2015) Neufeld, D. A., Godard, B., Gerin, M., et al. 2015, A&A, 577, A49
  * Neufeld et al. (2010) Neufeld, D. A., Goicoechea, J. R., Sonnentrucker, P., et al. 2010, A&A, 521, L10
  * Neufeld et al. (2009) Neufeld, D. A., Nisini, B., Giannini, T., et al. 2009, ApJ, 706, 170
  * Oba et al. (2018) Oba, Y., Tomaru, T., Lamberts, T., Kouchi, A., & Watanabe, N. 2018, Nature Astronomy, 2, 228
  * O’Dell (2001) O’Dell, C. R. 2001, ARA&A, 39, 99
  * Pabst et al. (2019) Pabst, C., Higgins, R., Goicoechea, J. R., et al. 2019, Nature, 565, 618
  * Pabst et al. (2020) Pabst, C. H. M., Goicoechea, J. R., Teyssier, D., et al. 2020, A&A, 639, A2
  * Palumbo et al. (1997) Palumbo, M. E., Geballe, T. R., & Tielens, A. G. G. M. 1997, ApJ, 479, 839
  * Pankonin & Walmsley (1978) Pankonin, V. & Walmsley, C. M. 1978, A&A, 64, 333
  * Parikka et al. (2017) Parikka, A., Habart, E., Bernard-Salas, J., et al. 2017, A&A, 599, A20
  * Pellegrini et al. (2009) Pellegrini, E. W., Baldwin, J. A., Ferland, G. J., Shaw, G., & Heathcote, S. 2009, ApJ, 693, 285
  * Peterson et al. (2008) Peterson, K. A., Adler, T. B., & Werner, H. J. 2008, J. Chem. Phys., 128, 084102
  * Pineau des Forets et al. (1986) Pineau des Forets, G., Flower, D. R., Hartquist, T. W., & Dalgarno, A. 1986, MNRAS, 220, 801
  * Prasad & Huntress (1980) Prasad, S. S. & Huntress, W. T., J. 1980, ApJS, 43, 1
  * Prasad & Huntress (1982) Prasad, S. S. & Huntress, W. T., J. 1982, ApJ, 260, 590
  * Putaud et al. (2019) Putaud, T., Michaut, X., Le Petit, F., Roueff, E., & Lis, D. C. 2019, A&A, 632, A8
  * Qu & Bowman (2016) Qu, C. & Bowman, J. M. 2016, J. Phys. Chem. A, 120, 4988
  * Rivière-Marichalar et al. (2019) Rivière-Marichalar, P., Fuente, A., Goicoechea, J. R., et al. 2019, A&A, 628, A16
  * Roelfsema et al. (2012) Roelfsema, P. R., Helmich, F. P., Teyssier, D., et al. 2012, A&A, 537, A17
  * Roncero et al. (2018) Roncero, O., Zanchet, A., & Aguado, A. 2018, Phys. Chem. Chem. Phys., 20, 25951
  * Rosenthal et al. (2000) Rosenthal, D., Bertoldi, F., & Drapatz, S. 2000, A&A, 356, 705
  * Sandford & Allamandola (1988) Sandford, S. A. & Allamandola, L. J. 1988, Icarus, 76, 201
  * Sanz-Sanz et al. (2013) Sanz-Sanz, C., Roncero, O., Paniagua, M., & Aguado, A. 2013, J. Chem. Phys., 139, 184302
  * Shiozaki & Werner (2013) Shiozaki, T. & Werner, H.-J. 2013, Mol. Phys., 111, 607, mRCI-F12
  * Smirnov et al. (1995) Smirnov, G. T., Sorochenko, R. L., & Walmsley, C. M. 1995, A&A, 300, 923
  * Smith (1991) Smith, R. G. 1991, MNRAS, 249, 172
  * Sofia et al. (2004) Sofia, U. J., Lauroesch, J. T., Meyer, D. M., & Cartledge, S. I. B. 2004, ApJ, 605, 272
  * Stecher & Williams (1972) Stecher, T. P. & Williams, D. A. 1972, ApJ, 177, L141
  * Sternberg & Dalgarno (1995) Sternberg, A. & Dalgarno, A. 1995, ApJS, 99, 565
  * Stoerzer et al. (1995) Stoerzer, H., Stutzki, J., & Sternberg, A. 1995, A&A, 296, L9
  * Stowe et al. (1990) Stowe, G. F., Schultz, R. H., Wright, C. A., & Armentrout, P. B. 1990, Int. J. Mass Spectrom. Ion Proc., 100, 377
  * Tablero et al. (2001) Tablero, C., Aguado, A., & Paniagua, M. 2001, Comput. Phys. Comm., 140, 412
  * Taquet et al. (2016) Taquet, V., Furuya, K., Walsh, C., & van Dishoeck, E. F. 2016, MNRAS, 462, S99
  * Tercero et al. (2010) Tercero, B., Cernicharo, J., Pardo, J. R., & Goicoechea, J. R. 2010, A&A, 517, A96
  * Tieftrunk et al. (1994) Tieftrunk, A., Pineau des Forets, G., Schilke, P., & Walmsley, C. M. 1994, A&A, 289, 579
  * Tielens (2010) Tielens, A. G. G. M. 2010, The Physics and Chemistry of the Interstellar Medium
  * Tielens & Hagen (1982) Tielens, A. G. G. M. & Hagen, W. 1982, A&A, 114, 245
  * Tielens & Hollenbach (1985) Tielens, A. G. G. M. & Hollenbach, D. 1985, ApJ, 291, 722
  * Tielens et al. (1993) Tielens, A. G. G. M., Meixner, M. M., van der Werf, P. P., et al. 1993, Science, 262, 86
  * Turner (1996) Turner, B. E. 1996, ApJ, 468, 694
  * van der Tak et al. (2007) van der Tak, F. F. S., Black, J. H., Schöier, F. L., Jansen, D. J., & van Dishoeck, E. F. 2007, A&A, 468, 627
  * van der Tak et al. (2003) van der Tak, F. F. S., Boonman, A. M. S., Braakman, R., & van Dishoeck, E. F. 2003, A&A, 412, 133
  * van der Tak et al. (2013) van der Tak, F. F. S., Nagy, Z., Ossenkopf, V., et al. 2013, A&A, 560, A95
  * van der Werf et al. (2013) van der Werf, P. P., Goss, W. M., & O’Dell, C. R. 2013, ApJ, 762, 101
  * van der Werf et al. (1996) van der Werf, P. P., Stutzki, J., Sternberg, A., & Krabbe, A. 1996, A&A, 313, 633
  * van der Wiel et al. (2009) van der Wiel, M. H. D., van der Tak, F. F. S., Ossenkopf, V., et al. 2009, A&A, 498, 161
  * van Dishoeck (2004) van Dishoeck, E. F. 2004, ARA&A, 42, 119
  * Velilla et al. (2008) Velilla, L., Lepetit, B., Aguado, A., Beswick, J., & Paniagua, M. 2008, J. Chem. Phys., 129, 084307
  * Vidal et al. (2017) Vidal, T. H. G., Loison, J.-C., Jaziri, A. Y., et al. 2017, MNRAS, 469, 435
  * Wakelam et al. (2017) Wakelam, V., Loison, J. C., Mereau, R., & Ruaud, M. 2017, Molecular Astrophysics, 6, 22
  * Walmsley et al. (2000) Walmsley, C. M., Natta, A., Oliva, E., & Testi, L. 2000, A&A, 364, 301
  * 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
  * Wyrowski et al. (1997) Wyrowski, F., Schilke, P., Hofner, P., & Walmsley, C. M. 1997, ApJ, 487, L171
  * Yamamura et al. (2000) Yamamura, I., Kawaguchi, K., & Ridgway, S. T. 2000, ApJ, 528, L33
  * Young et al. (2012) Young, E. T., Becklin, E. E., Marcum, P. M., et al. 2012, ApJ, 749, L17
  * Zanchet et al. (2013a) Zanchet, A., Agúndez, M., Herrero, V. J., Aguado, A., & Roncero, O. 2013a, AJ, 146, 125
  * Zanchet et al. (2018) Zanchet, A., del Mazo, P., Aguado, A., et al. 2018, PCCP, 20, 5415
  * Zanchet et al. (2013b) Zanchet, A., Godard, B., Bulut, N., et al. 2013b, ApJ, 766, 80
  * Zanchet et al. (2019) Zanchet, A., Lique, F., Roncero, O., Goicoechea, J. R., & Bulut, N. 2019, A&A, 626, A103
  * Zanchet et al. (2009) Zanchet, A., Roncero, O., González-Lezana, T., et al. 2009, Journal of Physical Chemistry A, 113, 14488
  * Zhou et al. (2020) Zhou, J., Zhao, Y., Hansen, C. S., et al. 2020, Nature Communications, 11, 1547

##  Appendix A H2S+ formation and destruction

In this Appendix we give details about how we calculated the H2 vibrational-
state-dependent rates of reaction (2) and of the reverse reaction, the
destruction of H2S+ (${}^{2}A^{\prime}$) by reactive collisons with H
(${}^{2}S$) atoms (summarized in Fig. 17).

We first built a full dimensional potential energy surface (PES) of the
triplet H3S+ (${}^{3}A$) system by fitting more than 150,000 ab initio points,
including the long range interactions in the reactants and products channels.
The main topological features of the PES are summarized in the minimum energy
path between reactants and products (see middle panel of Fig. 9). These ab
initio points were calculated with an explicitly correlated restricted coupled
cluster including a single, double, and (perturbatively) triple excitations
(RCCSD(T)-F12a) method (Knizia et al. 2009). The analytical fit has a overall
rms error of $\simeq$ 0.01 eV (Fig. 18). Appendix A.1 provides more details.

Reaction (2) is endothermic by 0.672 eV, and the PES of the triplet state
shows two shallow wells in the H2 \+ SH+ entrance channel (named
${}^{3}W_{1a}$ and ${}^{3}W_{1b}$, with a depth of $\simeq$ 0.118 eV) and
another one near the H + H2S+ products (named ${}^{3}W_{2}$, with a depth of
0.08 eV). Between the reactants and products wells there is a saddle point,
with an energy of 0.601 eV. This saddle point, slightly below the products,
has a geometry similar to ${}^{3}W_{2}$ in which the H–H distance is strongly
elongated compared to that of H2. These features are also present in the
maximum multiplicity PES of reactions H2 \+ S${}^{+}(^{4}S)$ and H2 \+
H2S+(${}^{2}A$) (see Fig. 9). We determine the state-dependent rates of
reaction (2) and of the reverse reaction using a quasi-classical trajectory
(QCT) method on our ground triplet PES. We provide more details on how the
reactive cross sections for fixed collision energies were calculated in
Appendix. A.2.

The formation rate of H2S+ from H2 ($v$ = 0) is very slow. For H2 ($v$ = 1),
the rate constant significantly increases at $\approx$ 500 K, corresponding
with the opening of the H2S+ \+ H threshold. At this point, it is important to
consider the zero-point energy (ZPE) of the products (see next section for
details). For H2 ($v$ = 2) and H2 ($v$ = 3), reaction rates are faster, close
to the Langevin limit. Finally, the H2S+ destruction rate constant is very
similar to that of its formation from H2 ($v$ = 2). In Appendix A.3 we provide
more information about the destruction of HnS+ ions through radiative
association and spin flip mechanisms.

Figure 17: Calculated rate constants as a function of temperature (for
translation and rotation) for SH+ ($v$ = 0, $j$ = 0) + H2 ($v$ = 1, 2, 3, $j$
= 0) and H2S+ ($v$ = 0, $j$ = 0) + H reactions (lavender) using ZPE corrected
QCT method. Dotted curves are fits of the form $k(T)$ = $\alpha$ ($T$/300)β
exp$(-\gamma/T)$. Rate coefficients are listed in Table 1. Figure 18: Rms
error as a function of total energy, showing the number of ab initio points
used to evaluate the error in the PES calculation. Arrows indicate selected
critical points in the PES and provide an estimate of the error in each
region. TS means transition state.

###  A.1 Ab initio calculations and PES

Dagdigian (2019) presented a PES for the SH+-H2 system that includes
4-dimensions and is based on RCCSD(T)-F12a ab initio calculations. This PES
was used to study SH+–H2 inelastic collisions using a rigid rotor approach in
which the two diatomic molecules are kept fixed at their equilibrium
distances. However, in order to study the reactivity of the collision, the two
diatomic distances have to be included to account for the breaking and
formation of new bonds.

Reaction (2) corresponds to a triplet state H3S+ (${}^{3}A$). The H2S+
(${}^{2}A^{\prime}$) + H (${}^{2}S$) products can form a triplet and a singlet
state. The triplet state can lead to the destruction of H2S+ through reaction
with H atoms. The singlet state, however, produces very excited states of the
reactants. Thus, it only leads to inelastic collisions but not not to the
destruction of H2S${}^{+}(^{2}A^{\prime})$. In consequence, here we only
consider the ground triplet electronic state of the system. In addition, the
H${}_{3}^{+}$ \+ S (${}^{3}P$) channel is about 2.4 eV above the H2 \+ SH+
asymptote, and will not be included in the present study.

In order to study the regions where several electronic states intersect, we
performed a explicitly correlated internally contracted multireference
configuration interaction (ic-MRCI-F12) calculation (Shiozaki & Werner 2013;
Werner & Knowles 1988a, b) including the Davidson correction (icMRCI-F12+Q;
Davidson 1975). The ic-MRCI-F12 calculations were carried out using state-
averaged complete active space self-consistent field (SA-CASSCF) orbitals with
all the CAS configurations as the reference configuration state functions. We
used a triple zeta correlation consistent basis set for explicitly correlated
wave functions (cc-pVTZ-F12; Peterson et al. 2008). In order to avoid orbital
flipping between core and valence orbitals. SA-CASSCF calculations with three
lowest triplet states were carried out including the core and valence orbitals
as active space (18 electrons in 11 orbitals). For the ic-MRCI-F12
calculation, the core orbitals was kept doubly occupied, resulting in about
$2.5\times 10^{6}$ $(9\times 10^{7})$ contracted (uncontracted)
configurations. All ab initio calculations were performed with MOLPRO (Werner
et al. 2012).

Our ic-MRCI-F12 calculations show that the crossings with electronic excited
states are 2 eV above the energy of the reactants. The energy interval below 2
eV is enough to study reaction 2. In these low-energy regions, RCCSD(T)-F12a
calculations were also performed. They are in good agreement with the ic-
MRCI-F12 results and the t1 diagnostic is always below 0.03. This allows us to
conclude that for energies below 2 eV, the RCCSD(T)-F12a method performs well,
presents a simple convergence, and being size consistent, is well adapted to
the present case. This method is the same one employed in the inelastic
collision calculations by Dagdigian (2019).

We performed extensive RCCSD(T)-F12a calculations in all accessible regions to
properly describe the six-dimensional phase space. 150000 ab initio points
were fitted to a multidimensional analytic function, that generates the six-
dimensional PES represented as

$\displaystyle H=H^{diab}+H^{MB}$ (8)

(Aguado et al. 2010; Sanz-Sanz et al. 2013; Zanchet et al. 2018; Roncero et
al. 2018), where $H^{diab}$ is an electronic diabatic matrix in which each
diagonal matrix element describes a rearrangement channel – six in this case,
three equivalent for SH+ \+ H2 channels, and three equivalent for H2S+ \+ H
fragments (we omitted the H${}_{3}^{+}$ \+ S channel) – as an extension of the
reactive force field approach (Farah et al. 2012). In each diagonal term, the
molecular fragments (SH+, H2 and H2S+) are described by 2 or 3 body fits
(Aguado & Paniagua 1992), and the interaction among them is described by a sum
of atom-atom terms plus the long range interaction. The non diagonal terms of
$H^{diab}$ are described as previously (Zanchet et al. 2018; Roncero et al.
2018) and the parameters are fitted to approximately describe the saddle
points along the minimum energy path in the right geometry.

In the reactants channel, the leading long range interaction
$\mbox{SH}^{+}(X^{3}\Sigma^{-})+\mbox{H}_{2}(X^{1}\Sigma_{g}^{+})$ corresponds
to charge-quadrupole and charge-induced dipole interactions (Buckinghan 1967):

$\displaystyle
V_{\mbox{charge}}(\mathbf{r}_{HH},\mathbf{R})=\Theta_{2}(r_{HH})P_{2}(\cos\theta_{2})R^{-3}$