Sulfur-bearing Molecules in Orion KL

The ALMA archive data were obtained from ALMA science verification⁸ (SV) data in band 6 (project ID 2011.0.00009.SV). The observations were performed on 2012 January 20, with baselines ranging from 17 to 265 m. The phase-tracking center was at R.A. = 05^h35^m1435 and decl. = −05°22'350. The ALMA spectral coverage is from 213.719 to 246.619 GHz, with a spectral resolution of 0.488 MHz (corresponding to ∼0.7 km s⁻¹ at 230 GHz). Bandpass and flux calibration were performed with Callisto. Quasar J0607-085 was observed for phase calibration. Continuum subtraction was performed, and the images were reduced with the MIRIAD (Sault et al. 1995) software. The images were deconvolved with natural weighting using the Clean algorithm. The synthesized beam size is 164 × 120 at 230 GHz. The 1σ root mean square (rms) noise level of the continuum and lines are ∼10 mJy beam⁻¹ and 30 mJy beam⁻¹ per channel, respectively. Figure 1 shows the 1.3 mm continuum map from the ALMA-only data. Adopting the same nomenclature as that given by Wu et al. (2014), we label the brightest condensation as the hot core and the resolved condensations as mm2-mm7.

Zoom In Zoom Out Reset image size

The ALMA-only observations filtered out ∼30%–50% of the extended emission in general compared to the 1.3 mm line emission detected with the IRAM-30 m single-dish observations (see the observational details in Feng et al. 2015). Therefore, we convert the IRAM-30 m data into the MIRIAD data format and use the following procedure to combine the ALMA and the IRAM-30 m data.

First, with the task UVMODEL, model visibility data from the single-dish spectral line are generated in the UV plane. The ALMA visibilities (red dot) and the IRAM-30 m visibilities (black dot) in the amplitude-UV-distance plane are shown in Figure 2 as an example. Then, using the Clean algorithm with natural weighting, the ALMA and IRAM-30 m visibilities are combined into deconvolved images containing large-scale emission (see Figure 3 for an example of the complementation of missing flux). The combined data cube has a synthesized beam of 186 × 153 (P.A. = −19°) in the upper sideband and 203 × 176 (P.A. = −126) in the lower sideband. The rms noise level of spectral lines is ∼60 mJy beam⁻¹ per km s⁻¹.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The broad bandwidth of the combined data set covers multiple transitions from a large number of S-bearing species, allowing us to perform an unbiased S-bearing study by excluding excitation effects. Adopting the eXtended CASA Line Analysis Software Suite (XCLASS, Möller et al. 2017), we are able to identify all lines of a particular species in our data simultaneously, based on molecular databases such as the Cologne Database for Molecular Spectroscopy (CDMS; Müller et al. 2005) or the database from from the Jet Propulsion Laboratory (JPL; Pickett et al. 1998).

The spectra extracted from the peak position toward each substructure are shown in Figure 4. The hot core contains the largest number of line detections, showing more intensive emission than the remainder of Orion KL. Therefore, using XCLASS, we identified 74 lines from 20 S-bearing isotopologues toward the hot core (listed in Table 1), including 13 lines of SO₂, 14 lines of ³⁴SO₂, 8 lines of ³³SO₂, 9 lines of OS¹⁷O, 9 lines of OS¹⁸O, 4 lines of SO, 1 line of ³⁴SO, 1 line of ³³SO, 1 line of S¹⁸O, 1 line of H₂S, 1 line of ${{\rm{H}}}_{2}^{34}{\rm{S}}$, 1 line of ${{\rm{H}}}_{2}^{33}{\rm{S}}$, 2 lines of OCS, 1 line of OC³³S, 2 lines of O¹³CS, 1 line of ¹⁸OCS, 1 line of ¹³CS, 1 line of H₂CS, 2 lines of H₂C³⁴S, and 1 line of ${{\rm{H}}}_{2}{{\rm{C}}}^{33}{\rm{S}}$.

For each species, assuming that all lines of that species are under local thermodynamic equilibrium (LTE), we can use the Modeling and Analysis Generic Interface for eXternal numerical codes (MAGIX; Möller et al. 2013) package to perform the fitting process. Therein, a synthetic spectrum is modeled from an isothermal object in one dimension by taking the optical depth, line blending, source size, velocity, and linewidth into account.

Using MAGIX, we fit the synthetic spectra of all 20 isotopologues toward each substructure. The input parameters of a particular molecule in MAGIX include the source size (in arcsecond), the rotational temperature T_rot (K), the molecular total column density N_tot (cm⁻²), the FWHM linewidth Δυ (km s⁻¹), and the central velocity υ_lsr (km s⁻¹). These input parameters are assumed to be the same for different transitions as initial guesses. By minimizing χ² in the given parameter space, MAGIX yields optimized results as output. For the isotopologues of the same species, we assume that all lines have the same centroid velocity. Moreover, we treat cases of line blending in which two or more possible lines contribute more than 5% of the observed intensity. An optimal fit is obtained by using three fitting algorithms, genetic, Levenberg–Marquardt, and errorestim-ins. Figure 4 shows an example of the synthetic spectra produced by MAGIX fitting to the lines in the frequency range of 215,380–217,218 MHz toward individual substructures.

Zoom In Zoom Out Reset image size

View extended figure (10 panels)

Analyzing the fitting results, we obtain the following results:

• 1.  
  In general, all the S-bearing species identified in our data set are well fit, exhibiting small optical depths (Figure 1). One exception is SO₂ toward the hot core, mm2, mm3, and mm4, where the fits for some low-J lines indicate optical depths greater than 1. These lines may originate from different temperature components (i.e., from both the central protostellar objects and the outer envelope). Photons from the inner hot component will be absorbed by molecules from the colder envelope, which will lead to overfitting of the low-J lines. (See the similar results given in Ahmadi et al. 2018, Appendix B). In such cases, we only use the high-J lines for the fitting process. The exclusion of the optically thick (low-J) lines will yield more accurate rotational temperatures.
• 2.  
  The LTE condition seems to be a reasonable assumption for the substructures we study here. Given that a successful fitting of a particular species with MAGIX requests at least three confirmed transitions and that SO₂ is the only species for which we detected more than three unblended lines toward all substructures (see the example of its unblended lines toward the hot core in Figure 3), we can use its fitting results to validate the assumptions in our source environment when using MAGIX. These SO₂ lines in our data set cover an E_u range of 19.03–1126.34 K. The rotation temperature T_rot of SO₂ from the fitting results (Table 2) is shown in Figure 5. We found that SO₂ exhibits the highest T_rot toward the hot core region (∼176 K) and the lowest T_rot toward the mm7 region (∼68 K). The T_rot of SO₂ is ∼106 K toward mm2, mm3, mm4, and mm5, without significant variations. This result is consistent with Feng et al. 2015. At such high temperature, the critical densities of these lines are 10⁴ ∼ 10⁶ cm⁻³, which is less than the number density of the substructures (10⁷ ∼ 10⁸ cm⁻³). For the rest of the S-bearing species (given in Table 1), we assume that their excitation temperatures toward individual substructures are the same as those of SO₂. Thus, the LTE condition is also a good approximation for estimating their column densities.
• 3.  
  The S-bearing lines in our data set seem to show different line profiles. Figure 6 shows line profiles of representative lines for each species. Figure 7 gives the statistical FWHM linewidth of these species toward individual substructures. Most lines from the carbon-free S-bearing species (including SO₂, ³⁴SO₂, ³⁴SO, and H₂S) exhibit a single peak toward the hot core, with a central velocity of ∼7 km s⁻¹. Those lines have large line wings toward mm2, mm3, mm4, and mm5. Due to overlap of different velocity components along the line-of-sight, the line profiles of H₂S show multiple peaks toward mm2, mm3, mm4, and mm5. The FWHM linewidth toward the hot core is 7–14 km s⁻¹ and toward mm7 is 3–8 km s⁻¹ (The FWHM linewidth of ³⁴SO₂ is 15 km s⁻¹, due to its weak emission), which is significantly narrower than the values for the rest of the substructures (12 ∼ 26 km s⁻¹). Lines from the carbon–sulfur compounds (including OCS, ¹³CS, and H₂CS) exhibit a single velocity component toward each substructure, with a central velocity in the range of 7 ∼ 9 km s⁻¹. However, they exhibit a narrower FWHM line width toward mm2 to mm7 (2 ∼ 8 km s⁻¹) than that of the hot core (7 ∼ 11 km s⁻¹). Since Orion KL is a complex region with multiple outflows, bringing layers of time-dependent shocks (Zapata et al. 2011; Crockett et al. 2014; Bally et al. 2017), the varying line widths of different S-bearing groups may be the result of components with different chemical ages.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Integrating the intensity of the representative lines for each species in the velocity range of 0–16 km s⁻¹, we present the spatial distribution maps of six species in Figure 8. The carbon-free S-bearing species and the carbon–sulfur compounds exhibit different spatial extents. The extended emissions of the carbon-free species cover the region from northeast mm4 down to the south mm2. The extended emissions of the carbon–sulfur compounds do not cover the northern mm4 but instead shift to the southern mm7, exhibiting a "heart-shaped" morphology. This result indicates that these two groups of S-bearing species are chemically different, tracing different gas.

Zoom In Zoom Out Reset image size

Moreover, when checking the line spatial distribution maps channel by channel, we also note that a ring-like structure appears on the line maps of carbon-free S-bearing species (SO₂, ³⁴SO, H₂S) in the velocity range of 10–15 km s⁻¹ (Figure 9). This ring, centered at R.A. = 05^h35^m14235 and decl. = −05°22'327, has a radius of ∼5'' in the plane of the sky. The dust emission peaks for the hot core, mm2, mm3, and mm5 are at the edge of the ring, suggesting that the carbon-free S-bearing species may be excited by shocks from the OMC1 explosion 500 yr ago (Plambeck & Wright 2016; Bally et al. 2017).

Zoom In Zoom Out Reset image size

Using XCLASSMapfit (with the same assumptions and algorithms as MAGIX), we fit the T_rot and column density maps by modeling the synthetic spectra toward all pixels. From the map fittings to SO₂ lines (Figure 10), we note a significant temperature and column density gradient from the hottest (100 ∼180 K) and densest (7.8 × 10¹⁷ cm⁻²) center of the hot core, through the warm (100 ∼ 120 K) mm2-mm5, and to the most distant and coolest (∼60 K, 8 × 10¹⁵ cm⁻²) mm7. This gradient may be the result of radiative pumping from being externally heated (De Buizer et al. 2012; Orozco-Aguilera et al. 2017) or shocks (Zapata et al. 2011; Wright & Plambeck 2017).

Zoom In Zoom Out Reset image size

A precise measurement of the molecular column density is key to study the chemical variations in different species. Under the LTE assumption, the column density of SO₂ is obtained by fitting multiple optically thin lines using MAGIX (as stated in Section 3.2). Column densities of the other S-bearing species are estimated by using the SO₂ rotation temperature as the gas temperature. Therefore, to test if these assumptions are reasonable, we derive the relative abundance ratios between the ³²S-/³⁴S-/³³S-isotopologues as well as the ¹²C-/¹³C-isotopologues (Table 3). We compare the values with the canonic ³²S/³⁴S, ³⁴S/³³S, and ¹²C/¹³C isotopic ratios, respectively (Anders & Grevesse 1989; Langer & Penzias 1990, 1993; Chin et al. 1996; Lucas & Liszt 1998; Persson et al. 2007; Tercero et al. 2010).

We find that the average relative abundance ratio of ³²SO₂/³⁴SO₂ toward Orion KL is 20 ± 4, which is consistent with previous observations toward Orion KL (e..g, 23 ± 7 from ³²SO₂/³⁴SO₂, Persson et al. 2007; 20 ± 6 from OC³²S/OC³⁴S, Tercero et al. 2010). The result is consistent with the solar value of 23 (Anders & Grevesse 1989), the local diffuse ISM value of 19 ± 8 from C³²S/C³⁴S from absorption observations (Lucas & Liszt 1998), and the galactic average value of 24 ± 5 from C³²S/C³⁴S (Chin et al. 1996). The average relative abundance ratio of H₂C³²S/H₂C³⁴S is 53 ± 7 in our study, which is slightly greater than the value of ³²S/³⁴S from ³²SO₂/³⁴SO₂ and a previous study by Tercero et al. (2010). The ratio of H₂C³²S/H₂C³⁴S toward the hot core is 19 ± 8, which is consistent with previous results. We notice that the linewidth of H₂C³⁴S is less than 2 km s⁻¹ (2 channels in the spectra) toward mm2, mm3, mm4, and mm5, which means that the column density of H₂C³⁴S may be underestimated. Thus, the ³²S/³⁴S ratio derived from ³²SO₂/³⁴SO₂ (20 ± 4) should be more reasonable.

The average abundance ratio of ³⁴S/³³S derived from ³²SO/³⁴SO is 6 ± 1 toward Orion KL, which is consistent with the same molecular pair of Esplugues et al. (2013; 6 ± 3), Tercero et al. (2010; ∼5), and Persson et al. (2007; 4.9). Our result is consistent with the solar value of 5.6 (Anders & Grevesse 1989) and the galactic average value of 6 ± 1 from C³⁴S/C³³S (Chin et al. 1996). Moreover, the derived average ¹²C/¹³C ratio (38 ± 9 from O¹²CS/O¹³CS) is consistent with the previous work of Tercero et al. (2010; 45 ± 20 from O¹²CS/O¹³CS) and Persson et al. (2007; 57 ± 14 from 12CH₃OH/¹³CH₃OH). Langer & Penzias (1990, 1993) found a galactic gradient in the ¹²C/¹³C isotopic ratio, which increases from 30 in the inner galaxy (5 kpc) to 70 at 12 kpc. The ¹²C/¹³C ratio they derived from ¹²C¹⁸O/¹³C¹⁸O is 63 ± 6 toward Orion A, which is greater than our result.

From the above test, we believe that the assumptions of LTE in XCLASS and gas temperature in this work are reasonable. Therefore, the chemical variations between different species can be directly indicated by comparing their relative abundances at the same position.

In this work, we estimate the relative abundance (column density ratio) of S-bearing species with respect to H₂. Instead of using continuum emission and the gas-to-dust ratio, we use C¹⁸O to estimate the H₂ column density for the following reasons: (1) The emission lines at 1 mm are so rich toward Orion KL that it is difficult to define the "line-free" part for the continuum. Moreover, we do not have the bolometric data to compensate for the missing flux of the ALMA-only continuum observations. Therefore, we do not know how significantly the dust continuum we are currently using is overestimated or underestimated. (2) The extended emission from the ALMA C¹⁸O (2-1) data is complemented using IRAM-30 m data, and our MAGIX fitting results indicate that this line is optically thin (τ ≤ 0.2). At a kinetic distance of ∼400 pc, the H₂ column density can be estimated by converting from C¹⁸O as [C¹⁸O/H₂] ∼ 2 × 10⁻⁷ toward Orion KL (Frerking et al. 1982; Plume et al. 2012; Crockett et al. 2014; Giannetti et al. 2014).

The column density maps of the S-bearing species are derived, and their relative abundance ratios with respect to H₂ are shown in Figure 11. In general, the carbon-free S-bearing species exhibit larger relative abundance with respect to H₂ to the north, i.e., from the hot core and mm4 to mm5 and mm6, with low abundances toward mm2 and mm7. In contrast, the carbon–sulfur compounds seem to have large abundance not only toward the continuum peaks but also toward mm2 and mm7. Given that the lines we used to derive the molecular column densities are observed simultaneously and cover a large range of E_U/k, we believe this chemical segregation is neither an excitation effect nor a sensitivity bias but rather the result of chemical differentiation.

Zoom In Zoom Out Reset image size

Inspecting possible correlations between both parameters pixel by pixel (Figure 12), we note that the abundances of OCS and ¹³CS decrease as the gas temperature increases from 40 to 180 K, and that of H₂CS does not seem to be sensitive to this gas temperature range. As for the carbon-free S-bearing species, their abundances exhibit a drastic increase as the gas temperature changes from 40 to 90 K, increasing by an order of magnitude. Similar observational results confirmed that the abundance of SO₂ was enhanced by more than two orders of magnitude from cold outer envelopes (T < 100 K) to hot inner envelopes (T > 100 K) (van der Tak et al. 2003). However, molecular abundances of these species exhibit no variation when the gas temperature increases from 90 to 180 K, which may be the result of depletion mechanisms (Wakelam et al. 2011).

Zoom In Zoom Out Reset image size

It has long been proposed that relative abundance ratios of different S-bearing species have the possibility to trace the chemical history of star-forming regions (e.g., Charnley 1997; Hatchell et al. 1998; Wakelam et al. 2004, 2011). In Orion KL, the abundances of two groups of S-bearing species seem to be dependent on the gas temperature (Section 4.2). Therefore, we investigate the possible correlation between the gas temperature and the relative abundance ratios of ³⁴SO/³⁴SO₂, OCS/SO₂, and H₂S/SO₂ (Figure 13).

Zoom In Zoom Out Reset image size

The ³⁴SO/³⁴SO₂ ratio decreases from 1.5 to 0.7 as the temperature increases from 60 ∼ 180 K, which is consistent with the results of Esplugues et al. (2014). It is known that SO can easily convert to SO₂ as the temperature increases (Charnley 1997). In particular, shocks can enhance the abundances of SO and SO₂ by two orders of magnitude (e.g., Pineau des Forets et al. 1993; Bachiller & Pérez Gutiérrez 1997; Wright & Plambeck 2017). As shown in Section 3.3, both SO and SO₂ may enhance by shock events in Orion KL.

The OCS/SO₂ ratio decreases from 1.0 to 0.3 when the temperature increases from 40 to 100 K and stays constant at ∼0.2 from 100 to 180 K, consistent with the model prediction by Wakelam et al. 2011 (shown as their Figure 7). SO₂ may be produced from species such as SO or evaporated to the gas phase at >100 K, whereas OCS may be destructed in such a temperature regime.

The H₂S/SO₂ ratio is ∼0.3, exhibiting no obvious variations in the temperature range of 50–180 K. This result is consistent with that of Esplugues et al. (2014), which indicates a chemical age of ∼5 × 10⁴ yr for the hot core region.