High-resolution SOFIA/EXES Spectroscopy of SO₂ Gas in the Massive Young Stellar Object MonR2 IRS3: Implications for the Sulfur Budget

The models assume SO₂ and ¹³CO are present in a uniform slab perpendicular to and covering the line of sight. Given a column density (N, cm⁻²), an excitation temperature (T_ex, K), an intrinsic line width (or Doppler parameter, b, km s⁻¹, related to the FWHM by $\mathrm{FWHM}=2b\sqrt{\mathrm{ln}2}$), and a Gaussian instrumental line spread function, the LTE model generates an absorption spectrum.

The molecular line parameters (Einstein A coefficient, partition function table, and quantum numbers) were retrieved from the HITRAN database. These parameters were used to calculate the population in each energy level for a gas with temperature T_ex and total column N following the standard Boltzmann equation. Subsequently, line equivalent widths were calculated and then converted to line profiles of width b at infinite resolution. Each line has a Voigt profile that is then convolved with a Gaussian profile to the instrumental resolution and used to compute a continuum normalized intensity.

The LTE model was fit across the whole spectrum simultaneously and the best fit was found by using a χ² value as our maximum likelihood estimator. To obtain the uncertainties, we used a Markov Chain Monte Carlo (MCMC) sampling to determine the posterior distributions of our model parameters. From these posterior distributions we obtained the 68% credibility interval. To do the sampling we used the open source Python package emcee (Foreman-Mackey et al. 2013). Figure 1 shows selected subsets of our best fits plotted over the data. The Doppler shift of the target is not a parameter we fit for; instead, we adopted a value based off the V_LSR of 10 km s⁻¹ measured in previous submillimeter studies of MonR2 IRS3 (van der Tak et al. 2003) and in good agreement with our SO₂ and ¹³CO observations (Figure 2).

Zoom In Zoom Out Reset image size

Concurring with previous studies (Giannakopoulou et al. 1997; Smith et al. 2016), we fit two temperature components to the observations for both SO₂ and ¹³CO. The priors common to all of our fits were the restrictions that T_ex > 1, N > 0.0, and b > 0.0. We also folded an uncertainty on spectral resolution into our determination of b, the prior for this was a Gaussian distribution. For EXES the peak probability occurred at R = 55,000 with a standard deviation of 1100, and for NIRSPEC the peak probability occurred at R = 25,000 and was uncertain to 10% at the 3σ level.

Fitting for the warm SO₂ alone, we found a T_ex of 234 ± 15 K, an N of ${4.95}_{-0.30}^{+0.29}\times {10}^{16}\,{\mathrm{cm}}^{-2}$, and an upper limit of b < 3.20 km s⁻¹. The cold component was much less well determined, and so we only quote upper limits at the 95% confidence level: T_ex < 88 K and N < 1.3 × 10¹⁴ cm⁻², assuming the same upper limit on the Doppler parameter.

For the ¹³CO we first measured the Doppler parameter by stacking the absorption features in our spectrum and measuring the line width assuming a Gaussian line profile and R of 25,000 (this approach is not possible for the crowded SO₂ spectrum). This yielded b values of 7.4 ± 2.2 and 5.3 ± 2.1 km s⁻¹ for the warm and cold ¹³CO components, respectively. These values were then used as additional priors. The other key difference from our SO₂ analysis was that we fit the warm ¹³CO to the high-J level transitions, where it is the only contributor, before fitting a cold component on top of the now determined warm component in the low-J level transitions. A single temperature component was unable to produce deep absorption features in both low-J level transitions and high-J level transitions. This yielded a best fit for the warm ¹³CO component of T_ex = 240 ± 25 K and N = 1.1 ± 0.2 × 10¹⁷ cm⁻², and a cold component with T_ex = 10 ± 7 K and $N={3.7}_{-1.0}^{+0.6}\times {10}^{16}\,{\mathrm{cm}}^{-2}$. These values are consistent with those found by a curve of growth analysis (Smith et al. 2016, R. L. Smith et al. 2018, in preparation).

Since the measured line width for the warm ¹³CO gas (7.4 ± 2.2 km s⁻¹) is broader than that of the SO₂ gas (<3.20 km s⁻¹), we are possibly including additional gas not associated with the reservoir of SO₂ we observed. Assuming a typical ¹²CO/¹³CO = 80 and H₂/¹²CO = 5000 (Lacy et al. 1994), we derive a lower limit on the abundance of SO₂ relative to N_H (=N(H) + 2N(H₂)) of (5.6 ± 0.5) × 10⁻⁷ for the warm SO₂. Comparing this lower limit to the cosmic sulfur abundance (S/H = 1.3 × 10⁻⁵; Asplund et al. 2009) shows that this SO₂ gas accounts for >4% of the sulfur budget. We place an upper limit on the cold SO₂ abundance of 4.4 × 10⁻⁹ by using a cold ¹²CO gas column of N = 80 × 3.7 × 10¹⁶ = 3.0 × 10¹⁸ cm⁻². Frozen CO contributes little. Using data from Gibb et al. (2004) we derived a ¹²CO ice column upper limit of 0.5 × 10¹⁷ cm⁻².

We also calculate an abundance relative to H₂O, allowing for direct comparison with cometary results. Boonman et al. (2003) reported a column density ${N}_{{{\rm{H}}}_{2}{\rm{O}}}=5\pm 2\times {10}^{17}\,{\mathrm{cm}}^{-2}$ with ${250}_{-100}^{+200}\,{\rm{K}}$. This yields a warm abundance, SO₂/H₂O, of >10 ± 3%. Lacking measurements of the foreground H₂O, we could not do the same for our cold SO₂ measurements. All abundances are summarized in Table 2.

Table 2.  Abundances of Sulfur-bearing Molecules toward MonR2 IRS3 and Comet 67/P

Species	Hot Corea		Foreground Gasb	Foreground Ice		Comet 67/P's Coma
	XH	${X}_{{{\rm{H}}}_{2}{\rm{O}}}$	XH	XHc	${X}_{{{\rm{H}}}_{2}{\rm{O}}}$ d	${X}_{{{\rm{H}}}_{2}{\rm{O}}}$
	10−7	%	10−7	10−7	%	%
SO2	$\gt 5.6\pm 0.5$	10 ± 3	$\lt 0.044$	$\lt 5.7$ (1)	$\lt 0.6$ (1)	0.127 ± 0.003 (4)
${{\rm{H}}}_{2}{\rm{S}}$	⋯	⋯	⋯	$\lt 2.8$ (2)	$\lt 1.1$ (2)	1.10 ± 0.05 (4)
$\mathrm{OCS}$	⋯	⋯	⋯	<0.18 (3)	<0.07 (3)	0.041 ± 0.001 (4)
S2	⋯	⋯	⋯	⋯	⋯	0.197 ± 0.003 (4)

Notes. Sources: (1) derived from data in Gibb et al. (2004), (2) Smith (1991), (3) Palumbo et al. (1997), (4) Calmonte et al. (2016).

^aCalculated from an SO₂ column of 4.95 × 10¹⁶ cm⁻², and a derived hydrogen column of 8.8 × 10²² cm⁻² or an H₂O column of 5 × 10¹⁷ cm⁻² (Section 3.3). ^bCalculated from an SO₂ column of 1.3 × 10¹⁴ cm⁻², and a hydrogen column of 3.0 × 10²² cm⁻² (Section 3.3). ^cRelative to a hydrogen column of 3.0 × 10²² cm⁻² derived from our cold ¹³CO gas column (Section 3.3). ^dRelative to an ice column ${N}_{{{\rm{H}}}_{2}{\rm{O}}}$ = 1.9 × 10¹⁸ cm⁻² (Gibb et al. 2004).

Download table as:  ASCIITypeset image

Millimeter-wave observations of MonR2 IRS3 have measured the SO₂ column density at beam sizes of ∼15'', probing the cool envelope (van der Tak et al. 2003). The column density of 1.5 ± 0.3 × 10¹⁴ cm⁻² is consistent with our cold SO₂ component's upper bound of 1.3 × 10¹⁴ cm⁻². Additionally, our measured warm SO₂ component is consistent with infrared measurements by Keane et al. (2001), who reported ${T}_{{ex}}={225}_{-70}^{+50}$ and N = 4.0 ± 0.8 × 10¹⁶ cm⁻² using an adopted b of 3 km s⁻¹.

Shocks have previously been observed to lead to substantial enhancements in SO₂ abundances (Pineau des Forets et al. 1993; Bachiller & Pérez Gutiérrez 1997; Podio et al. 2015). The extreme temperatures enable a variety of gas-phase reactions of pre-shock gas or sublimated ices that enhance the formation of SO₂. Indeed, shock chemistry models successfully replicated the measured abundance of gas-phase SO₂ toward the Orion Plateau, which shows very broad lines indicative of shocks generated by Orion IRc 2 outflows (b ≳ 12–15 km s⁻¹; Blake et al. 1987). Also, shocks of tens of km s⁻¹ can shatter or sputter dust grains (May et al. 2000), possibly leading to the release of more sulfur and subsequent SO₂ formation. However, the SO₂ lines toward MonR2 IRS3 are substantially narrower (b < 3.20 km s⁻¹; Figure 2) than those in the Orion Plateau. Our results are therefore more consistent with release from the ices due to radiative heating (or perhaps mild shocks) rather than from refractory grains in strong shocks. Also, the SO₂/H abundance derived for MonR2 IRS3 (0.6 × 10⁻⁶) is somewhat higher than that measured in Orion IRc 2 (0.2 × 10⁻⁶; Blake et al. 1987) and consistent with the range reported for the outflow target HH 212 ((0.4–1.2) × 10⁻⁶; Podio et al. 2015). The formation of SO₂ in hot cores is thus at least as efficient as in shocks and, importantly, sputtering of sulfur off of grains does not seem to be a required process.

This process, for the first time observed in a massive hot core, might also be important in lower-mass YSOs. The strong enhancement of SO molecules observed in a protoplanetary disk (Booth et al. 2018) might thus relate to ice sublimation rather than grain sputtering in shocks.

There still remains the question of which molecular species lead to the efficient formation of gas-phase SO₂. The large mismatch between ice-phase detections of SO₂ (Boogert et al. 1997; Zasowski et al. 2009) and the warm component we measured indicate that SO₂ cannot be sublimating directly from the ice. Chemical models generally rely on large abundances of sublimated H₂S ice for SO₂ formation (Charnley 1997; Woods et al. 2015). Indeed, measurements by the ROSINA experiment have shown that H₂S is the largest contributor to the sulfur budget in the comet 67P/Churyumov–Gerasimenko (Calmonte et al. 2016). In stark contrast, H₂S ice measurements toward dense clouds and protostellar envelopes yielded upper limits (Smith 1991; Jiménez-Escobar & Muñoz Caro 2011) a factor of two below our gas-phase SO₂ abundance (Table 2), and we thus hesitate to conclude this is the primary source of sulfur that is driven into SO₂ molecules. Instead, we consider the release of sulfur allotropes (e.g., S₂, S₃, S₄) that have been shown to be the second most abundant sulfur carrier in the ices of comet 67P (Calmonte et al. 2016; Drozdovskaya et al. 2018). These allotropes can be broken apart by helium atoms, allowing for gas-phase reactions between OH, O₂, and S that lead to the formation of SO₂ (McElroy et al. 2013).

Chemical models suggest that at temperatures ≳230 K the oxygen is preferentially driven into H₂O (Charnley 1997; Doty et al. 2002; van Dishoeck et al. 2011). The measured SO₂ temperature is comparable to this threshold. This suggests the gas-phase formation of SO₂ in MonR2 IRS3 is suppressed in favor of the production of H₂O, though our results are not inconsistent with sub-230 K temperatures. Alternatively, the SO₂ we observe formed before further heating of the gas, or at larger, cooler radii. High H₂O abundances were measured along this line of sight, though there remains the possibility that these observations probed a warmer region closer to the hot core than the SO₂ gas we observed. This cannot be excluded following observations by Boonman et al. (2003), who measured a temperature of ${250}_{-100}^{+200}\,{\rm{K}}$. A high spectral resolution analysis by N. Indriolo et al. (2018, in preparation) is expected to shed more light on this possibility.

The basis for the LTE assumption is that collisional excitations at least match the rate of radiative de-excitations. The minimum density at which this occurs is called the critical density, and values for SO₂ are of order 10⁶ cm⁻³ to 10⁷ cm⁻³ (Wakelam et al. 2004; Williams & Viti 2014). Using a modeled temperature profile for a hot core from van der Tak et al. (1999) and the radial density profile of MonR2 IRS3 reported by van der Tak et al. (2000), we find that SO₂ with a temperature of 230 K is expected to reside at a radius of ∼500 au with a particle density of n = 6 × 10⁶ cm⁻³, comparable to the SO₂ critical densities. Also, the critical density of ¹³CO is over an order of magnitude smaller than that of SO₂, and thus ¹³CO is most likely in LTE. The similarity of the excitation temperatures of these molecules thus further indicates that the SO₂ transitions are thermalized.