The Extraordinary Outburst in the Massive Protostellar System NGC 6334 I-MM1: Strong Increase in Mid-Infrared Continuum Emission

We used the Stratospheric Observatory for Infrared Astronomy (SOFIA; Temi et al. 2014) during Cycle 7 to image NGC 6334 I with the Faint Object infraRed CAmera for the SOFIA Telescope (FORCAST; Herter et al. 2012). Observations were performed in the asymmetric chop-and-nod imaging mode C2NC2. This field was also observed with the High-resolution Airborne Wideband Camera-plus (HAWC+; Harper et al. 2018) in Lissajous scanning mode during Guaranteed Time Observations. Further details are listed in Table 1. We used the Level 4 and Level 3 archival data products for FORCAST and HAWC+, respectively. To facilitate closer comparison with higher resolution data, we also generated a deconvolved version of the FORCAST 25 μm image with ≈2'' resolution.

ALMA observations with the 12 m array were obtained with comparable angular resolution in Bands 4, 7, 8, and 9 (2.2, 1.0, 0.76, and 0.43 mm) within 8 months during ALMA Cycle 5 (Table 1). We calibrated the data with the Cycle 7 release of the ALMA pipeline. NGC 6334 I is an extremely line-rich source (Section 1). The method described in Brogan et al. (2018c) was employed to identify line-free channels; the resulting net continuum bandwidths are listed in Table 1. Iterative self-calibration was performed on the continuum data to improve the image quality (Brogan et al. 2018a). The Briggs robust weighting parameter (smaller values down-weight shorter UV-spacings) was used to achieve similar angular-scale sensitivity between the different bands beyond that afforded by the projected UV-ranges of the data (see Table 1); these ALMA data are not sensitive to smooth emission larger than about 8''. Spectral cubes of SiO (3-2) and CS (6-5) with 2 km s⁻¹ channels were created from the Band 4 and 7 continuum-subtracted data, respectively. With the exception of the C43-6 Band 4 data (McGuire et al. 2018; Xue et al. 2019; Ligterink et al. 2020), results from these ALMA data sets have not previously been published.

We compiled the latest maser monitoring spectra of NGC 6334 I from the HartRAO 26 m telescope. Observation details are provided by MacLeod et al. (2018). For each maser transition, we extracted a time series of the intensity in the velocity channel closest to the LSR velocity of the thermal gas (−7.25 km s⁻¹), which also represents the peak methanol maser emission from MM1 (Hunter et al. 2018).

To establish accurate astrometry for the infrared images, we extracted an archival pre-outburst 06 K_s image from the High Acuity Wide field K-band Imager (HAWKI) on the Very Large Telescope (VLT) on 2012 June 9. Applying 2MASS photometric calibration, the faintest detected object is ∼21 mag. Despite the large difference in wavelength, the morphological features from IRS-I-1 and IRS-I-3 are remarkably similar to the 1998 18 μm image from the Cerro Tololo Inter-American Observatory (CTIO) 4 m telescope (Figure 1(a); De Buizer et al. 2000). We calibrated the HAWKI image astrometry to Gaia DR2 (Lindegren et al. 2018), then aligned the CTIO image to it. We then aligned the three subfields observed in 1999 at W. M. Keck Observatory (De Buizer et al. 2002) to the CTIO image. Next, we aligned the FORCAST images to the CTIO image, and finally the HAWC+ image to the 37 μm image. We also extracted the 70 μm PACS image (Tigé et al. 2017) from the Herschel archive and aligned it to the HAWC+ image. Table 2 lists the shifts applied to each image. Finally, we also extracted 3.4 and 4.6 μm images from the (NEO)WISE survey covering pre- and post-outburst (Mainzer et al. 2011), and K_s images from the VVV/VVVX survey (Minniti et al. 2010), all post-outburst.

Prior to the 2015 outburst, mid-IR images of NGC 6334 I were dominated by dust emission from the UCHII region MM3 (IRS-I-1; Figure 1(a)). Aside from MM3, faint emission at 10 and 18 μm was detected toward and slightly southwest of MM2 in the CTIO and Keck images (IRS-I-2; Figures 1(a) and (b)). The elongation of this feature coincides with high-velocity blueshifted emission in the ALMA images of SiO 3-2 and CS 6-5 (Figures 1(c) and (d)) leading toward H₂ knot C (Eislöffel et al. 2000), providing another example of mid-IR emission tracing outflow cavities in massive star-forming regions (De Buizer et al. 2017). MM1 is undetected in pre-burst mid-IR imaging (De Buizer et al. 2000), with a 4σ upper limit of 0.36 Jy at 18 μm (see Figure 2(d) contours).

Zoom In Zoom Out Reset image size

Figures 2(a) and (b) show ALMA continuum images of NGC 6334 I, highlighting the four primary protocluster members (MM1–MM4). The UCHII region MM3 is prominent at 2.2 mm due to its free–free emission, but becomes less apparent at 1.0 mm. The spectral index (α) image (S_ν ∝ ν^α ) formed from these images (Figure 2(c)) demonstrates that MM1, MM2, and MM4 are dominated by dust emission, with the most optically thick sources (α approaching +2) being MM4 and locations within MM1. In contrast, MM3 shows α ∼−0.1 typical of optically thin free–free emission.

As of mid-2019, MM1 is the brightest source at 25 and 37 μm, outshining MM3 (Figures 2(d) and (e); this is also true in the mid-2018 native 53 μm HAWC+ image, not shown). This remarkable result confirms a large increase in luminosity, first inferred from millimeter data (Hunter et al. 2017). Separated by only 3'', MM1 and MM2 are difficult to distinguish in the native resolution FORCAST images. However, the deconvolved 25 μm image (Figure 2(f)) demonstrates that both objects contribute at this wavelength. Interestingly, the 25 μm emission peaks not at the center of MM1, where the column density is likely highest, but rather ∼1''southwest toward a cluster of 6.7 GHz masers that were not present prior to the outburst (Hunter et al. 2018). The 25 μm emission also extends northward following an arc of masers along the N–S outflow (McGuire et al. 2018; Brogan et al. 2018c; Chibueze et al. 2021) and westward, including the MM2 area, supporting the conclusion of De Buizer et al. (2012) that these masers trace mid-IR-bright outflow cavities surrounding massive protostars. The coincidence of all the methanol masers with 25 μm emission also reflects the expectation of pumping models that require T_dust ≳ 120 K (Cragg et al. 2005).

The millimeter-wavelength flux ratio of MM1:MM2 increased from r_preburst = 2.1 in the SMA ¹² 0.87 mm data to r_burst = 10.5 in the ALMA 1.0 mm data. However, because MM1 and MM2 cannot be clearly separated in the SOFIA images, the most accurate spectral energy distribution (SED) we can produce is for the combination of the two sources. Thus, we used the ALMA images to define a polygonal aperture that generally follows a low intensity contour of elliptical shape covering MM1 and MM2 but avoiding the bright northern edge of MM3 (Figures 2(a) and (b)). Primary beam correction was applied to all ALMA images prior to measuring flux densities. The ALMA flux density uncertainties are the estimated rms noise in the aperture (rms per beam $\times \sqrt{\mathrm{beams}\,\mathrm{in}\,\mathrm{aperture}}$) added in quadrature with the absolute flux uncertainty (10% for Bands 4 and 7; 15% for Bands 8 and 9).

Considering the coarser resolution of the mid-IR images, we needed to model and remove the UCHII region emission to avoid contamination. To construct a model, we masked the non-UCHII emission in the CTIO 18 μm image, then smoothed it to the resolution of each of the longer wavelength images (Figure 2(g)–(i)). Next, each smoothed 18 μm image was iteratively scaled upward to minimize the residuals when subtracting it from the corresponding longer wavelength image. For the 25, 37, 53, and 70 μm images, we found optimal scale factors of 3.33 (α₁₈₋₂₅ = 3.53), 7.0 (α₁₈₋₃₇ = 2.7), 12 (α₁₈₋₅₃ = 2.3), and 2.5, respectively. These values are comparable to intermediate values of α₁₉₋₃₇ computed for massive protostar models by Liu et al. (2020), and the 37:20 μm and 70:20 μm flux ratios of the cometary H ii region M17UC1 (5.7 and 2.25: Lim et al. 2020).

SOFIA and Herschel values in Table 2 are measured from "UCHII removed" images; to account for the larger 37, 53, and 70 μm beams, we expanded the photometric aperture radially from its centroid by the relative increase in beam radius compared to 25 μm (015, 12, and 135, respectively, see expanded apertures in Figures 2(h) and (i)). Removing the UCHII region reduces the measured flux densities by 9.2%, 4.8%, 12%, and 13% at 25, 37, 53, and 70 μm, respectively. The photometric uncertainty due to astrometric error was estimated by shifting the aperture by 1 pixel in each of the four cardinal directions, repeating the photometry, then computing the rms of these four possible differences with respect to the nominal flux density: it is 5.1%, 5.5%, 1.9%, and 4.5% at 25, 37, 53, and 70 μm, respectively. The total flux density uncertainties—12% at 25 and 37 μm, 16% at 53 μm, and 20% at 70 μm—include the estimated noise in the aperture, photometric uncertainty due to astrometric error, UCHII model subtraction uncertainty (half the subtracted contribution), and flux calibration uncertainty (10% for FORCAST, 15% for HAWC+, 5% for PACS, Balog et al. 2014; Chuss et al. 2019; Lim et al. 2020), added in quadrature. The total uncertainty at 70 μm includes an additional term for the background subtraction.

No emission was detected toward MM1 or MM2 in the HAWKI image, with upper limits of 2.65 and 5.3 μJy, respectively. In the (NEO)WISE images, MM1 is within the saturated PSF core of MM3 in both bands; subtracting pre-outburst from post-outburst images did not reveal any residual toward MM1. Stacking all the K_s band VVV/VVVX images yields an upper limit of 4.2 μJy near both positions.

In order to assess the stability of the continuum emission from NGC 6334 I over the span of the post-outburst data, we use the maser light curves from HartRAO. Among the four species shown in Figure 3(a), the Class II methanol masers (6.7 and 12.2 GHz) are pumped by mid-IR photons (Cragg et al. 2005) and may be a good proxy for detecting changes in this radiation. Prior to the 2015 January outburst, the methanol masers arose only from MM2 and MM3 (Green et al. 2015). The 12.2 GHz flare lasted only a few months, but the 6.7 GHz flare persists. During the 1.5 yr interval spanning our ALMA and SOFIA observations (Table 1), the 6.7 GHz maser flux density remained constant over 44 epochs with a mean and standard deviation of 716 ± 45 Jy. Computed using pyMannKendall (Hussain & Mahmud 2019), the Mann–Kendall normalized test statistic is 0.23 and the Theil–Sen slope estimator is 0.16, both consistent with no trend.

Zoom In Zoom Out Reset image size

Combining observations spanning 18 months into an SED can be problematic for a variable source. However, massive protostars exhibiting near-IR variability (Kumar et al. 2016) show similar timescale fluctuations as the 6.7 GHz masers surrounding deeply embedded massive protostars (Goedhart et al. 2014). Also, in a sample of 6.7 GHz masers with radio continuum counterparts observed on arcsecond scales, the maser luminosity correlates significantly with the inferred ionizing photon rate (Szymczak et al. 2018). Along with theoretical expectations, these observational findings suggest that variations in the radiation field will often manifest in maser variations. Thus, the observed constancy of the methanol maser flux density from NGC 6334 I-MM1 during the interval spanning our mid-IR/millimeter observations gives us confidence in constructing the separate SEDs in Figure 3.

We fit these SEDs to protostellar radiative transfer (RT) models (Robitaille 2017) using the methodology of Towner et al. (2019). The best-fit outburst model (Figure 3b) includes a power-law envelope with an outflow cavity and a central protostar with bolometric luminosity L_burst = 49,000 ± 8000 L_⊙, where the fractional uncertainty in L_burst has been set to that of the flux density measurement nearest the SED peak (53 μm). The best-fit pre-outburst model (Figure 3c) comes from the same model family, with nearly identical outflow inclination angle but half the opening angle and has L_preburst = 4300 ± 900 L_⊙ (taking the uncertainty from the 70 μm measurement). MM1 and MM2 both harbor massive protostars, with hot core molecular spectra of comparable complexity and temperature (T_exc = 150–200 K, El-Abd et al. 2019). Therefore, we use their pre-outburst millimeter flux density ratio, r_preburst (Section 3.2), to estimate their pre-outburst luminosity ratio. Since the brightness of MM2 has remained unchanged (Hunter et al. 2017), we compute MM1's pre-outburst luminosity as L_MM1preburst = 2900 ± 600 L_⊙ and its outburst luminosity as L_MM1burst = 47,600 ± 7800 L_⊙ (assuming 20% uncertainty on r_preburst). Thus, the outburst ratio is 16.3 ± 4.4, consistent with a magnitude 3 outburst (Meyer et al. 2021).