© ESO, 2015

1. Introduction

To date, more than a dozen molecules containing sulphur have been detected unequivocally in interstellar gas clouds and/or circumstellar outflows: SH, SH⁺, CS, NS, SO, SiS, H₂S, SO₂, HCS⁺, OCS, C₂S, C₃S, HNCS, HSCN, H₂CS, and CH₃SH. All but one of these molecules were first detected using ground-based telescopes operating at millimeter or submillimeter wavelengths, the only exception being the mercapto radical, SH, which has its lowest rotational transition at 1383 GHz, a frequency inaccessible from ground-based observatories and also lying in a gap between Bands 5 and 6 of the Herschel Space Observatory’s HIFI instrument. The first detection of interstellar SH was obtained (Neufeld et al. 2012) in absorption towards the star-forming region W49N using the German Receiver for Astronomy at Terahertz frequencies (GREAT; Heyminck et al. 2012) instrument¹ on the Stratospheric Observatory For Infrared Astronomy (SOFIA; Young et al. 2012).

The chemistry of interstellar molecules containing sulphur is distinctive in that none of the species S, SH, S⁺, SH⁺, or H₂S⁺ can react exothermically with H₂ in a hydrogen atom abstraction reaction (i.e. a reaction of the form X + H₂ → XH + H.) This feature of the thermochemistry greatly inhibits the formation of sulphur-bearing molecules in gas at low temperature.

In this paper, we present a corrected analysis of the SH abundance presented in Neufeld et al. (2012), together with additional detections of SH obtained with SOFIA/GREAT toward four other bright continuum sources – W51, W31C, G34.3+0.1 and G29.96–0.02 – and ancillary observations of H₂S, CS, SO and H₃S⁺ obtained with the IRAM 30 m telescope at Pico Veleta toward all five sources. (The observations of H₃S⁺ led only to upper limits, not detections.) The observations and data reduction are described in Sect. 2 below, and results presented in Sect. 3. The abundances of SH, H₂S, CS, and SO – and upper limits on the H₃S⁺ abundances – derived in diffuse molecular clouds along the sight-lines to these bright continuum sources are discussed in Sect. 4 with reference to current astrochemical models.

2. Observations and data reduction

Table 1 lists the continuum sources toward which SH, H₂S, CS, SO and H₃S⁺ were observed, along with the exact positions targeted, the estimated distances to each source, the source systemic velocities relative to the local standard of rest (LSR), the dates of each SOFIA observation, and the total integration time for each transition. In Table 2, we list the transitions targeted for each species, together with their rest frequencies, telescope beam size at the transition frequency, upper and lower state energies, and spontaneous radiative decay rates. All the molecular data presented here were taken directly from the CDMS (Müller et al. 2005) line catalog, with the exception of those for the H₃S⁺ transition; for the latter transition, which is not included in the CDMS catalog, we adopted the frequency measured by Tinti et al. (2006), and obtained the radiative decay rate for an assumed dipole moment of 1.731 D (Botschwina et al. 1986).

All the SOFIA data reported here were obtained with the L1 receiver of GREAT and the AFFTS backend; the latter provides 8192 spectral channels with a spacing of 183.1 kHz. As described in Neufeld et al. (2012), the W49N SH data were acquired in the upper sideband of the L1 receiver; observations of the other four sources were performed with the SH line frequency in the lower sideband. At this frequency, the telescope beam has a diameter of ~21′′ HPBW. The observations were performed in dual beam switch mode, with a chopper frequency of 1 Hz and the reference positions located 60′′ (75′′ for W49N) on either side of the source along an east-west axis.

The raw data were calibrated to the (“forward beam brightness temperature”) scale, using an independent fit to the dry and the wet content of the atmospheric emission. Here, the assumed forward efficiency was 0.97 and the assumed main beam efficiency for the L1 band was 0.67 (except in the case of the W49N data from 2011, where the beam efficiency was 0.54). The uncertainty in the flux calibration is estimated to be ~20% (Heyminck et al. 2012). Additional data reduction was performed using CLASS². This entailed the removal of a baseline and the rebinning of the data into 10-channel bins, resulting in a 0.40 km s^-1 channel spacing for the final data product.

All the IRAM data reported here were obtained using the EMIR receivers (Carter et al. 2012) on the 30-m telescope at Pico Veleta on 2014 March 16−18, with the H₂S, CS, SO, and H₃S⁺ transitions observed respectively using the E1, E0, E0, and E3 receivers. Since two receivers can be used simultaneously (and upper and lower sideband spectra obtained separately and simultaneously for two orthogonal polarizations), we targeted the H₂S transition in all observations together with either (1) the H₃S⁺ transition when the weather conditions were excellent (precipitable water vapor <1 mm) or (2) the CS and SO transitions. The FTS backend was deployed for all observations, providing a 200 kHz channel spacing over a 4.05 GHz bandwidth. In addition, the VESPA autocorrelator was used to target the CS transition toward W49N, W31C, G34.3+0.1 and G29.96−0.02 and the SO transition toward W49N and W31C; the VESPA backend was used to obtain a 80 kHz channel spacing over a bandwidth of 160 MHz. All the data were acquired in double beam switching mode, using the wobbling secondary with a throw of ± 45 arcsec. Observations of Mars and Mercury were used for periodic checks of the telescope focus.

The CLASS software package was used to coadd all the data obtained for each source, with a weighting inversely proportional to the square of the radiometric noise in each scan. The final coadded spectra show the antenna temperature, averaged over both polarizations, as a function of the frequency in the LSR frame.

3. Results

In Table 3, we report the continuum brightness temperatures measured at each frequency, together with the rms noise in the coadded spectra. The continuum spectra of H₂S, SO and CS are of a particularly high signal-to-noise ratio, providing exquisite sensitivity to weak spectral features. In Figs. 1 through 5, we show the spectra obtained toward each source. The left panels are scaled to show the full range of observed antenna temperatures. The H₂S, CS, and SO spectra show strong emission at the systemic velocity of each source, together with multiple narrower absorption features. The SH spectra show only absorption. The H₃S⁺ spectra, which show multiple interloper emission lines (i.e. emission lines other than the targeted line), show no evidence for H₃S⁺ emission or absorption. The right panels show the same data with the vertical scaling adjusted to reveal weaker absorption features.

In each panel in Figs. 1−5, the red parts of each histogram denote spectral regions that were assumed to be devoid of foreground absorption. These regions were subsequently fit with a flat continuum and one or more Gaussian emission features, to provide a background spectrum of the source itself, T_bg(v_LSR). The observed spectra were then divided by that fitted background spectrum (including emission lines), to obtain a transmission spectrum, t = T_A/T_bg, and an optical depth spectrum, τ(v_LSR) = −ln t(v_LSR). In the case of the SH transition, which is split by lambda doubling and hyperfine splitting, the velocity scale is computed for the F = 3⁻ ← 2⁺ transition. As is clear from Figs. 1−5, the absorption pattern is therefore repeated with a velocity shift of +71.8 km s^-1, corresponding to the magnitude of the lambda doubling, and is broadened by the hyperfine structure.

Figures A.1 through A.5 (black histogram) show the optical depth as a function of LSR velocity for the SH, H₂S, CS, and SO transitions. The vertical axis shows the total molecular column density per velocity interval. The latter is computed for an excitation temperature of 2.73 K, the temperature of the cosmic microwave background. The validity of this assumption is justified by excitation calculations described in Appendix B. In the case of H₂S, an ortho-to-para ratio of 3 was assumed. Blue curves show multigaussian fits to the column density per velocity interval. In the case of SH, the conversion to column density is complicated by the presence of hyperfine structure and lambda doubling. Here, the blue curve represents the fitted column density per unit velocity interval, while the red curve represents the corresponding optical depth (after convolution with the hyperfine structure and lambda doubling). Thus, in this case, it is the red curve that should be compared with the optical depths shown in the black histogram. The blue symbols above the spectra indicate the centroids and widths (full-widths-at-half-maximum) for each fitted Gaussian. In our fits, the linewidths were constrained to be identical for all four molecules, and the centroids allowed to vary by at most 0.2 km s^-1.

In Table 4, we present the foreground molecular column densities inferred for selected velocity intervals toward each source, along with column density ratios. The selected velocity windows are marked by the red vertical lines in Figs. A.1−A.5, and the numbers appearing above each window give the molecular column densities in units of 10¹² cm^-2. The average abundances for each velocity range appear in Table 5, computed relative to molecular hydrogen (Cols. 4−7) and relative to hydrogen nuclei (Cols. 8−11). For this purpose, we adopted atomic hydrogen column densities, N(H), and hydrogen nucleus column densities, N_H = N(H) + 2 N(H₂), obtained with the assumptions given in the table footnote. The values given in Table 4 for the W51 v_LSR = 62−75 km s^-1 window must be regarded with caution. Material absorbing in this velocity interval may lie in close proximity to W51, where radiative excitation is important: except for the case of SH, the derived column densities, computed for the conditions in foreground diffuse clouds, may be unreliable. For this reason, the W51 v_LSR = 62−75 km s^-1 has not been included in Table 5.

Our observations did not lead to any detection of H₃S⁺. The line emission spectrum is quite rich in this spectral region, and the presence of interloper emission lines (i.e. lines from species other than H₃S⁺) reduces our sensitivity to broad H₃S⁺ features.In particular, an emission feature lying ~7 km s^-1 blueward of the H₃S⁺ line – most likely the 3₂₁−4₁₄ transition of methanol at 293.464203 GHz – greatly limits our sensitivity to H₃S⁺ emission from the background sources. However, good upper limits can be placed on the H₃S⁺ optical depths associated with certain narrow H₂S absorption components. The best limits on the optical depth ratio, τ(H₃S⁺) /τ(H₂S), are obtained for the 40 km s^-1 absorption component in the foreground of W49N and the 30 km s^-1 absorption component in the foreground of W31C. The 3σ upper limits on τ(H₃S⁺) /τ(H₂S) obtained for these components are 1.06 × 10^-2 and 1.02 × 10^-2 respectively. Figure A.6 shows the H₃S⁺ (green) and H₂S (red, scaled by a factor 0.1) optical depth spectra obtained for these two components, together with the 3σ upper limit on H₃S⁺ (blue). These limits correspond to a column density ratio N(H₃S⁺) /N(H₂S) ~ 0.028(T_rot/ 100 K). Here, T_rot is the rotational temperature for the metastable (J = K) states of (the symmetric top molecule) H₃S⁺, a negligible population being assumed for all other states within each K-ladder.

4. Discussion

4.1. Measured column densities, abundances, and correlations among the molecules observed

All the molecular abundances in Table 5 lie far below the solar abundance of sulphur, (n_S/n_H)_⊙ = 1.3 × 10^-5 (Asplund et al. 2009), indicating that the four molecules together account for at most 0.14% of the sulphur nuclei within any of the velocity ranges we defined. Moreover, sulphur is notable in showing no measurable depletion in diffuse atomic clouds (Sofia et al. 1994), and only moderate levels of sulphur depletion (by a factor ≲ 4) have been inferred even in the dense (n_H ~ 10⁵ cm^-3) gas within the Horsehead PDR (Goicoechea et al. 2006), suggesting that the molecules that we have observed contain a similarly small fraction of the gas-phase sulphur nuclei.

Although the chemical pathways leading to the formation of SH, H₂S, CS and SO are clearly inefficient in converting sulphur atoms and ions into sulphur-bearing molecules, they evidently operate in a rather uniform manner, with the relative abundances among the sulphur-bearing molecules showing a remarkable constancy. In Figs. 6, we show (top panels) the column densities of SH, CS and SO for each velocity interval in Table 3, as a function of N(H₂S). The dependence is nearly linear in every case, with the column density ratios (bottom panels in Figs. 6 and Cols. 6−9 in Table 4) spanning a fairly narrow range; the N(SH) /N(H₂S), N(CS) /N(H₂S), and N(SO) /N(H₂S) ratios lie in the ranges 1.12−2.96, 0.32−0.61, and 0.08−0.30 respectively. Furthermore, these column density ratios exhibit no obvious trend with N(H₂S) as the latter varies over a factor ~13. In Figs. 7, we present the column densities (left panel) and abundances (right panel) of SH, H₂S, CS and SO relative to H₂, as a function of N(H₂). Similar results are presented in Figs. 8, but with abundances given relative to H nuclei as a function of N_H. Once again, there is no strong trend in abundance with either N(H₂) or N_H.

4.2. Principal component analysis

We have used Principal Component Analysis (PCA; e.g. Jolliffe 2002) to explore differences and similarities among the absorption spectra observed toward W31C and W49N. These sources are suited particularly well to this analysis because their absorption spectra exhibit a complex structure involving multiple foreground velocity components, detected at high signal-to-noise ratio, and unassociated with the source itself. In our analysis, we made use of the optical depths, as a function of LSR velocity, for the SH, H₂S, CS and SO transitions observed here, together with analogous results for seven additional transitions of interest: the HI 21 cm transition (Winkel et al., in prep.); the CH N=1, J=3 / 2−1 / 2 lambda doublet (Gerin et al. 2010a); the H₂O⁺1₁₁−0₀₀ transition (Gerin et al. 2010b; Neufeld et al. 2010); the ¹³CH⁺J = 1−0 transition (Godard et al. 2012); the SH⁺J_N = 1₂−0₁ transition; the ArH⁺J = 1−0 transition (Schilke et al. 2014), and the CN N = 1−0 transition (Godard et al. 2010). In the case of those transitions (the SH, CH, SH⁺, ¹³CH⁺, and CN transitions) that show hyperfine structure, we used the hyperfine-deconvolved spectra (i.e. the optical depths inferred for the strongest hyperfine component).

These eleven “optical depth spectra” were rebinned on a common velocity scale covering LSR velocities in the range 10 to 50 km s^-1 for W31C and 35 to 75 km s^-1 for W49N. They are shown by the black histograms in Fig. 9 (for W31C) and Fig. C.1 for W49N. Prior to implementing the PCA, we followed the standard procedure of scaling and shifting each optical depth spectrum so that the variances were unity and the means were zero. This ensures that each spectrum is given equal weight. The PCA method represents each optical depth spectrum as a linear combination of eleven orthogonal (i.e. uncorrelated) components, with the optical depth for the ith transition given by (1)where the S_j(v) are the 11 principal component optical depth spectra shown in Fig. 10 (for W31C) and Fig. C.2 (for W49N), and a_i, b_i, and C_ij are coefficients. The S_j(v) all have a mean of zero and a variance of unity, and the C_ij obey the normalization condition .

A key feature of PCA is that the terms within the sum in Eq. (1) above make successively smaller contributions to the right-hand-side, with the first term accounting for most of the variance in the optical depth spectra, the second term the next largest fraction of the variance, and so on. In this particular case, the first two principal components are found to account for 72% of the total variance, and the red histograms in Figs. 9 and B.1 show that a fairly good approximation to most of the observed spectra is obtained by retaining just the first two principal components³. The only notable exception is ArH⁺, for which the first two terms in equation (1) yield a poor fit to the observed spectra; for that case, the blue histograms in Figs. 9 and C.1 show the sum of the first three terms.

In Fig. 11, we present, for each optical depth spectrum, i, the first three coefficients, C_i1, C_i2, and C_i3 in the sum in Eq. (1). The upper panel shows C_i2 versus C_i1, while the lower panel shows C_i3 versus C_i2. Each transition occupies a specific position within these plots, which provides a valuable graphical representation of the similarities and differences between the spectra. Color-coded squares show the coefficients for each transition, and lines are drawn from the origin to each point. The results plotted in Fig. 11 were obtained from a simultaneous analysis of the W31C and W49N data, but similar results were obtained by analyzing each source individually. Moreover, the relative positions of the SH, H₂S, CS, SO, CH, CN and HI transitions in the upper panel do not depend strongly upon whether the molecular ion transitions are included in the analysis. Several features are noteworthy:

• First, with the notable exception of ArH⁺ transition, most transitions lie fairly close to the unit circle (black circle) in the upper panel, indicating that the first two principal components contain most of the information about their spectra. (If the coefficients for the third and higher principal components were zero, then all points would lie exactly on the unit circle.) ArH⁺ is unusual in that its spectrum is largely accounted for by the third principal component. The fact that it is poorly correlated with most of the other molecules may reflect its unusual chemistry (Schilke et al. 2014): it is believed to be a tracer of gas with an extremely small (≲ 10^-3) molecular fraction.

• Second, the transitions of the four neutral sulphur-bearing molecules lie very close to each other in the plot, are significantly separated from the HI, CH, and all the molecular ion transitions, and are relatively close to the CN transition. In the limit where the points lie very close to the unit circle, the correlation coefficient between any pair of transitions is equal to cosine of the angle between the relevant vectors (i.e. the lines drawn from the origin to each point). Indeed, the linear correlation coefficients among the neutral sulphur-bearing molecules range from 0.75 to 0.96 (with those among CS, H₂S, and SH lying in the range 0.90 to 0.96). By contrast, the correlation coefficients between CH and the neutral sulphur-bearing molecules are smaller, in the range 0.52 and 0.67, and those between HI and the neutral sulphur-bearing molecules are smaller yet: 0.32 to 0.39.

• Third, the molecular ions – other than ArH⁺ – have transitions that typically lie close together and close to HI in the upper panel. In particular, the SH⁺ molecular ion shows a distribution closer to that of CH⁺ and H₂O⁺ than to that of the neutral sulphur-bearing molecules observed in the present study. The molecular ions do, however, occupy a range of positions in the lower panel (i.e. they show a range of coefficients for the third principal component).

• Fourth, the vector corresponding to CH in Fig. 11 lies between the vectors for HI and those for CN and the neutral sulphur-bearing molecules. CH is believed to be a good tracer for H₂, both in gas that is mainly atomic and in gas that is mainly molecular, while CN is thought to be abundant only in regions of large molecular fraction, so this behavior may indicate that the neutral sulphur-bearing species are preferentially found in regions of large molecular fraction and are less prevalent in partially-molecular regions where CH is presentbut CN is absent. We note, in this context, that the relative positions of CH⁺, CH, and CN in the upper panel of Fig. 11 are consistent with previous UV absorption line observations (Lambert et al. 1990; Pan et al. 2004, 2005), which showed that CH line profiles could be seen as the sum of a broad component similar to the CH⁺ profile, and a narrow one, similar to that of the CN line.

4.3. Comparison with previous observations of CS, H₂S and SO

Sulphur-bearing molecules have previously been studied within diffuse molecular clouds, both through observations of line emission (e.g. Drdla et al. 1989; Heithausen et al. 1998; Turner 1996a,b) and absorption (Tieftrunk et al. 1994; Lucas & Liszt 2002; Destree et al. 2009). Absorption-line observations can provide the most reliable estimates of molecular column densities, because – as noted in Appendix B – the latter often depend only weakly upon the assumed physical conditions in the absorbing material or upon the adopted rate coefficients for collisional excitation. Accordingly, we focus here on previous absorption-line observations, which provide the most direct quantitative comparison with results obtained in the present study.

Millimeter-wavelength observations of absorption by S-bearing molecules in diffuse clouds have been reported by Lucas & Liszt (2002), who detected CS, H₂, and/or SO toward six extragalactic radio sources. The typical N(CS) /N(H₂) ratio inferred from these observations was 4 × 10^-9, a value lying close to the middle of the range of N(CS) /N(H₂) obtained in the present study (Table 5). Absorption line observations of CS in diffuse clouds are also possible at ultraviolet wavelengths through observations of the C–X (0, 0) band near 1401 Å, which has been detected along 13 sight-lines toward hot stars (Destree et al. 2009). The observed absorption line strengths are broadly consistent with the CS abundances inferred by Lucas & Liszt (2002) and in the present study, but a detailed comparison is limited by a lack of knowledge of the relevant f-values.

Lucas & Liszt (2002) reported N(CS) /N(SO) ratios in the range 1.7 ± 0.8, which overlaps the range (2−6) implied by Table 4, and N(CS) /N(H₂S) ratios in the range 6 ± 1. The H₂S column densities derived by Lucas & Liszt (2002), however, appear to be based on a column-density-to-optical-depth ratio (their Table 2) that is a factor of ~12 too small; applying this correction brings their N(CS) /N(H₂S) ratio down to 0.48 ± 0.08, in excellent agreement with the range (0.32−0.61) reported here in Table 4. In addition, Lucas & Liszt also detected the 85.348 GHz J = 2−1 transition of HCS⁺ in absorption toward one of the six sources, 3C 111, with an implied N(CS) /N(HCS⁺) ratio of 13.3 ± 0.1. (The HCS⁺J = 2−1 frequency was not covered in the present study.)

4.4. Comparison with astrochemical models

4.4.1. Models without “warm chemistry”

Standard models of photodissociation regions (PDRs) are known to greatly underpredict the abundances of sulfur-bearing molecules (e.g. Tieftrunk et al. 1994). For the densities (n_H ~ few × 100 cm^-3) and extinctions (A_V ~ few × 0.1 mag) typical of the diffuse foreground clouds we observed (e.g. Gerin et al. 2015), we find the abundances to be underpredicted by a factor of ~100, even in models that include gas-grain chemistry (Vasyunin & Herbst 2011) and the photodesorption and chemical desorption⁴ of molecules from grain surfaces.

4.4.2. Models for turbulent dissipation regions

As discussed by Neufeld et al. (2012), the discrepancy described in the previous subsection may provide evidence for a warm chemistry driven by turbulent dissipation in magnetized shocks (e.g. Pineau des Forêts et al. 1986) or in vortices (TDR model, Godard et al. 2009) where the gas temperature and/or ion-neutral drift velocity are sufficient to activate highly endothermic formation pathways. Such a scenario has indeed proven successful (Godard et al. 2014) in reproducing the large abundances of CH⁺ and SH⁺ observed in the diffuse interstellar gas assuming that turbulence dissipates in regions with density n_H ≤ 50 cm^-3, ion-neutral drift velocity u_in ~ 3 km s^-1, and at an average rate , a value in agreement with measurements of the kinetic energy transfer rate at large scales (Hennebelle & Falgarone 2012).

In Fig. 12, we show the main production and destruction routes for 19 sulfur- and carbon-bearing molecules predicted by the PDR, turbulent dissipation region (TDR) and shock models for a gas with density n_H = 50 cm^-3 and shielding A_V = 0.4 mag. For each species, only the processes that together contribute to more than 70% of the total destruction and formation rates are displayed. Solid blue arrows indicate reactions that are only important in PDRs while solid red arrows indicate reactions that are only important in TDRs or shocks. Green and dashed red arrows denote additional processes discussed in Sect. 4.3.4 below. Activation energy barriers are indicated by temperatures (ΔE/k) marked next to the relevant reaction arrow. The chemical network we used here is an updated version of that used in the Meudon PDR code (Le Petit et al. 2006⁵).Table D.1 lists the rates adopted for all the processes shown in Fig. 12, together with those for two additional reactions for which updated rate coefficients have been used. Here, many of the rates adopted for endothermic reactions were obtained from those for the exothermic reverse reactions using the principle of detailed balance. This procedure may introduce inaccuracies, because it is based upon the assumption that the reactants are in local thermodynamic equilibrium with the kinetic temperature. At the low densities typically of relevance in astrophysical applications, the degree of vibrational excitation – in particular – may be highly subthermal. Thus, the reaction rates for endothermic reactions with vibrationally-cold reactants may be significantly overestimated if vibrational excitation is substantially more efficient than translational energy in overcoming the energy barrier. Unfortunately, detailed scattering calculations with a quantum chemistry potential are needed to determine which reactions are most affected; these are not typically available in the literature.

TDR model predictions are shown in Fig. 13, where the predicted column density ratios are plotted againstN(SH) /N_H. Grey shaded regions indicate the range of observed values listed in Tables 4 and 5 (or, in the case of N(CS) /N(HCS⁺), reported by Lucas & Liszt 2002). On the one hand,the observed N(CS) /N(SH) and N(SH) /N_H ratios can be accounted for by models with an ion-neutral drift velocity of at least 6 km s^-1, somewhat larger than that invoked to explain the CH⁺/SH⁺ column density ratio. This result suggests that turbulent energy dissipates over a distribution of structures with different strengths rather than a single type of event. On the other hand, however, the TDR model predictions systematically underestimate the column densities of both H₂S and SO by roughly an order of magnitude unless A_V is unrealistically large (A_V ≥ 1)⁶.Moreover, the CS/HCS⁺ ratios reported by Lucas & Liszt (2002) for the sight-line toward 3C111 are larger than the model predictions except for the highest assumed density (n_H = 300 cm^-3).

4.4.3. Shock models

A similar discrepancy is found when the predictions of standard shock models are compared with the data. Using the model described by Lesaffre et al. (2013), we have conducted a parameter study in which the column densities of sulphur-bearing molecules behind a shock were computed as a function of the shock velocity, v_s, the preshock magnetic field, B₀, the preshock density, n_H, and the UV shielding. The latter is parameterized by the visual extinction, A_V, to the cloud surface, which we assumed to be irradiated by the typical Galactic UV radiation field (Draine 1978). As long as v_s exceeds a value sufficient to drive the endothermic production reactions of relevance (about 10 km s^-1 for n_H = 100 cm^-3 and A_V = 0.1), the SH/H₂S ratio is virtually independent of v_s and B₀. As in the TDR case, the predicted H₂S and SO column densities relativeto those of SH are underestimated by a factor of about 10 for clouds with n_H = 100 cm^-3 and small shielding A_V ≤ 1. Similarly, the CS/HCS⁺ ratio is also underestimated except for high density values. Unlike TDR models, however, shock models underestimate the CS/SH ratio by a factor of about 5 for A_V ≤ 1.

4.4.4. Additional effects on the chemistry

As a possible explanation of the small observed N(SH) /N(H₂S) ratios relative to the predictions of standard shock models, we have investigated two effects. First, we considered the dependence of the predicted SH/H₂S column density ratio upon the branching ratio for dissociative recombination of H₃S⁺. In the standard model, we adopted a branching ratio for H₂S production of 0.17, the value measured by Kamińska et al. (2008) for the production of D₂S by dissociative recombination of D₃S⁺. If we assume instead that every recombination of H₃S⁺ leads to the production of H₂S, the SH/H₂S ratio decreases to a value of about 5, i.e. about twice the typically-observed value.

Second, we have investigated a physical effect not previously included in standard models for the chemistry of sulphur-bearing molecules in interstellar shock waves. After being produced by dissociative recombination of H₂S⁺ or H₃S⁺, SH molecules will initially carry an imprint of the velocity of the ionized parents from which they formed; this velocity drift, relative to other neutral species, might significantly enhance the production of H₂S via the endothermic neutral-neutral reaction SH(H₂,H)H₂S (for which ΔE/k = 6862 K). As shown by Flower & Pineau des Forêts (1998), who followed the velocity of CH molecules within a C-type shock, computing the drift between SH and H₂ only requires the rate coefficient for momentum transfer between SH and H₂, k_mt(SH,H₂). In principle, this rate could be obtained by subtracting the inelastic collisions rates (including both reactive and non reactive processes) from the total SH + H₂ collision rate. However, as far as we know, the collision frequency for this system has never been estimated; moreover the rates of elastic and inelastic processes depend strongly on the gas kinetic temperature. As a result, we have explored the shock model predictions over a wide range of momentum transfer rate coefficient: 10^-12 ≤ k_mt(SH,H₂) ≤ 5 × 10^-10 cm³ s^-1.

In Fig. 14, we show the drift velocity of SH relative to H₂, as a function of distance through the shock. As expected, the SH velocity lies somewhere between the ion velocity (red) and the H₂ velocity (black), with the exact behavior depending upon the value adopted for k_mt(SH,H₂).

Figure 15 shows predictions for the shock model, with “velocity memory effects” included for SH. The reaction pathway enabled by such effects is shown by the green arrow in Fig. 12. The results are shown as a function of n_H and for three values of A_V (0.1, 1, and 2). Here, the top panel shows the SH column density for a single shock, while the otherfour panels showthe CS/SH, H₂S/SH, SO/SH, and CS/HCS⁺ column density ratios, all computed for a single C-type shock with a velocity of 14 km s^-1 and a preshock magnetic field B₀ = (n_H/ cm^-3)^{1 / 2}μG. Results are shown for k_mt(SH,H₂) of 1 × 10^-10 cm³ s^-1 (dashed curves) and 5 × 10^-10 cm³ s^-1 (solid curves). The grey shaded areas in each panel indicate the range of observed values presented in Table 4. In the case of the top panel, the shaded area refers to the SH column density per magnitude of visual extinction; a comparison with the model predictions then leads to the number of shocks per magnitude of visual extinction required to reproduce the observations.For A_V ≤ 1, the inclusion of velocity memory effects for SH, with an assumed k_mt(SH,H₂) of 1 × 10^-10 cm³ s^-1 (dashed curves), significantlyincreases the predicted H₂S/SH ratio, although the observed values still lie a factor ~3 above the predictions unless A_V is unrealistically large (≳ 1).This process has however a small impact on all the other ratios displayed in Fig. 15⁷.The SO/SH and CS/SH column density ratios are similarly underpredicted. The effect of velocity memory on the TDR model predictions will be investigated in future studies but is expected to have an impact similar to that found for shock modelsin increasing the predicted H₂S/SH column density ratio.

Finally, as a possible explanation of the underprediction of SO column densities by TDR and shock models, we have considered a formation pathway of SO not yet included in standard chemical networks and denoted by dashed red arrows in Fig. 12: H₂S⁺ + O → HSO⁺ + H followed by the dissociative recombination of HSO⁺. So far, standard chemical networks available online (e.g. the UMIST database; McElroy et al. 2013) only include the channels H₂S⁺ + O → SH⁺ + OH and H₂S⁺ + O → SO⁺ + H₂ with identical rate coefficients of 3.1 × 10^-10 cm³ s^-1 (Prasad & Huntress 1980), despite the fact that the production of HSO⁺ is also a highly exothermic process. Implementing this channel in the shock model with an assumed rate coefficient of 6.2 × 10^-10 cm³ s^-1 enhances the predicted column density of SO by a factor of 10. Such a scenario has yet to be confirmed by theoretical or experimental estimates of the reactivity of H₂S⁺ with oxygen atoms and a detailed analysis of the possible products.

4.4.5. Line widths and velocity shifts

In the case of TDRs, an interpretation of the observed line profiles requires a convolution of the kinematic signatures of large scale turbulence with those of active regions of dissipation – taking into account their orientations, strengths and distribution along the line of sight – and of the relaxation stages that necessarily follow the dissipative bursts. In TDRs, SH⁺, CH⁺, and H₃S⁺ are formed in the active phase (Godard et al. 2014), while SH, H₂S, SO, and CS are mainly produced by the relaxation phase. The linewidths of the former species are therefore expected to be wider than those of the latter; this effect is indeed apparent in the optical depth spectra (Figs. 9 and C.1) and shows up in the PCA results (Fig. 11). In C-type shocks, the chemical production may induce a velocity shift between the line profiles of species. Flower & Pineau des Forêts (1998) estimate the shift between CH and CH⁺ induced by a single shock of 14 km s^-1 propagating along the line of sight to be roughly 2 km s^-1. The same could be expected between SH and the neutral fluid, yet no significant line shifts are apparent in the observed data. We note, however, that taking a random distribution of shocks with different orientations will drastically reduce the line shifts. A detailed calculation of the line profiles in the framework of both the TDR and shock models remains to be done.

5. Summary

• 1.

  With the use of the GREAT instrument on SOFIA, we observed the SH 1383 GHz lambda doublet towards the star-forming regions W31C, G29.96–0.02, G34.3+0.1, W49N and W51. We detected foreground absorption toward all five sources.

• 2.

  Using the IRAM 30 m telescope at Pico Veleta, we observed the H₂S 1₁₀−1₀₁ (169 GHz), CS J = 2−1 (98 GHz), SO 3₂−2₁ (99 GHz) and H₃S⁺1₀−0₀ (293 GHz) transitions toward the same five sources. The H₃S⁺1₀−0₀ transition was not detected in any source, but, for the other three transitions, we detected emission from the background source, and absorption by foreground clouds, toward every source.

• 3.

  In nine foreground absorption components detected towards these sources, the inferred column densities of the four detected molecules showed relatively constant ratios, with N(SH) /N(H₂S) in the range 1.1−3.0, N(CS) /N(H₂S) in the range 0.32−0.61, and N(SO) /N(H₂S) in the range 0.08−0.30. The observed SH/H₂ ratios – in the range 5 to 26 × 10^-9 – indicate that SH (and other sulphur-bearing molecules) account for ≪ 1% of the gas-phase sulphur nuclei.

• 4.

  We have used Principal Component Analysis to investigate similarities and differences in the distributions of SH, H₂S, CS, SO, and seven other species present in foreground diffuse clouds along the sight-lines to W31C and W49N. The neutral sulphur-bearing molecules we observed are very well correlated amongst themselves and with CN, moderately well correlated with CH (a surrogate tracer for H₂), and poorly correlated with atomic hydrogen and the molecular ions SH⁺, CH⁺, H₂O⁺, and ArH⁺. A natural, but not unique, interpretation of these correlations is that SH, H₂S, CS, SO are preferentially found in material of large molecular fraction.

• 5.

  The observed abundances of neutral sulphur-bearing molecules greatly exceed those predicted by standard models of cold diffuse molecular clouds, providing further evidence for the enhancement of endothermic reaction rates by elevated temperatures or ion-neutral drift.

• 6.

  Standard models for C-type shocks or turbulent dissipation regions (TDRs) can successfully account for the observed abundances of SH and CS in diffuse clouds, but they fail to reproduce the column densities of H₂S and SO by a factor of about 10. More elaborate shock models, which include “velocity memory effects” for SH, reduce this discrepancy to a factor ~3. Moreover, a possible reaction pathway not included in previous chemical networks, H₂S⁺(O,H)HSO⁺(e,H)SO, can potentially explain the larger-than-expected SO abundances.

• 7.

  Further theoretical studies will be needed to determine the importance of velocity memory effects in TDRs and any kinematic signatures that might be expected.

---

Acknowledgments

Based on observations made with the NASA/DLR Stratospheric Observatory for Infrared Astronomy. SOFIA Science Mission Operations are conducted jointly by the Universities Space Research Association, Inc., under NASA contract NAS2-97001, and the Deutsches SOFIA Institut under DLR contract 50 OK 0901. This research was supported by USRA through a grant for Basic Science Program 01-0039. E.F., M.G. and B.G. acknowledge support from the French CNRS/INSU Programme PCMI (Physique et Chimie du Milieu Interstellaire). E.H. thanks NASA for its support of his research in the chemistry of the diffuse interstellar medium. We thank B. Winkel for communicating HI 21 cm absorption data prior to publication, and V. Fish for providing data from Fish et al. (2003) in digital format. We thank H. Liszt for helpful discussions, and the anonymous referee for several valuable suggestions that improved the paper. We thank H. Müller for helpful discussions about the laboratory spectroscopy for SH, and for providing the line frequencies for SH – obtained from his simultaneous fit to the available data – through CDMS. We gratefully acknowledge the outstanding support provided by the SOFIA Operations Team, the GREAT Instrument Team, and the staff of the IRAM 30-m telescope.

References

Online material

Appendix A: Absorption features and molecular column densities

Figures A.1−A.5 show the column density per velocity interval for each detected molecule in each observed source. Figure A.6 presents upper limits on H₃S⁺.

Appendix B: Excitation of sulphur-bearing molecules in diffuse molecular clouds

In order to determine the equilibrium level populations for CS and SO in diffuse molecular clouds, we have solved the equations of statistical equilibrium as a function of temperature and density. Our excitation model includes radiative pumping by the cosmic microwave background (CMB), at temperature T_CMB = 2.73 K, spontaneous radiative decay, and collisional excitation by molecular hydrogen. We adopted the collisional rate coefficients recommended in the LAMDA database (Schöier et al. 2005)⁸, which were obtained by a simple scaling of the values computed for the excitation of CS (Lique et al. 2006a) and SO (Lique et al. 2006b) by He.

In Fig. B.1, we present results for the ratio of total molecular column density to optical depth, N/τ, for the CS J = 2−1 (upper panel) and SO J_N = 3₂−2₁ (lower panel) transitions that we observed. Results are presented as a function of H₂ density for four different temperatures. The left panels show the results on a log-log plot, while the right panels show the same data on a log-linear plot for a narrower range of n(H₂). In the limit of low density, the level populations are in equilibrium with the CMB, at a temperature of 2.73 K, while in the limit of high density the level populations are in local thermodynamic equilibrium (LTE) at the gas temperature, T. In the latter limit, N/τ is proportional to T² provided kT ≫ the transition energy, ΔE; this result arises because the difference between the absorption and stimulated emission rates is proportional to where N_u and N_l are the column densities in the lower and upper states, and Z is the partition function, proportional to T for diatomic molecules. At intermediate densities and temperatures ≳ 100 K, the populations are inverted.

Two caveats affecting the results shown in B.1 are our neglect of radiative trapping and of electron impact excitation in determining the level populations. While the SO transition is clearly optically-thin in all the absorption components we observed, the optical depth of the CS transition can reach values of order unity at line center; thus, for CS, radiative trapping can lead to a modest reduction in the critical density at which departures from the low-density behavior become significant. For polar molecules like CS and SO, electron impact excitation can be of comparable importance to excitation by neutral species in regions where carbon is fully-ionized. If the CS and SO absorption is occurring in regions where C is ionized, a further modest reduction in the critical density could result.

Notwithstanding these caveats, for the temperatures (T ≤ 80 K) and densities (≤ few × 10² cm^-3) typical of the foreground diffuse molecular clouds observed in this study (e.g. Gerin et al. 2015), the N/τ ratio is very comfortably in the low density limit for both CS and SO, justifying the assumption underlying our computation of column densities in Sect. 3 above. In the case of H₂S, the effects of collisional excitation upon the N/τ ratio are expected to be even smaller than those for CS and SO, both because the spontaneous radiative rate for the H₂S transition is larger than those for CS and SO, and because the H₂S absorption takes place out of thelowest rotational state of the ortho spin-species.

Appendix C: Principal component analysis for W49N

Figures C.1 and C.2 show the optical depth spectra for eleven transitions observed towards W49N, and the principal components determined from a simultaneous analysis of the W31C and W49N spectra.

Appendix D: Reaction network for sulphur-bearing molecules

The chemical network for sulphur-containing molecules adopted in the Meudon PDR code version 1.4.4 – a slightly-updated version of which is used here for TDR and shock modeling – includes 432 reactions. In Table D.1, we present an abbreviated list containing rates for those reactions and protoprocesses that play a dominant role in the formation and destruction of sulfur bearing species (Fig. 12) and/or whose rates have been updated for use in our study.

All Tables

All Figures