© ESO 2021

1 Introduction

The collapse of dense molecular cloud cores gives rise to the birth of stars, a process which was thought to proceed in a smooth and continuous fashion. However, the first piece of evidence for the unsteady growth of forming stars emerged by recognizing that the outburst of FU Orionis, which was thought to be a rare phase of early stellar evolution (Herbig 1966), was instead an episode of enhanced disk accretion (Hartmann & Kenyon 1985). Since then, it has been realized that episodic accretion is an intrinsic feature of forming young stars (Hartmann & Kenyon 1996; Audard et al. 2014). While this knowledge had been established exclusively from observations of low-mass stars, which become optically visible while still accreting, it was unknown until recently whether high-mass stars (M_⋆≳8 M_⊙) show the same behavior during their formation. Their scarcity and fast formation timescales imply that they are still deeply embedded in their parental core while reaching the main sequence. This suggests that similar outbursts during high-mass star formation, if present at all, might be difficult to detect.

A candidate young massive eruptive variable, V723 Car, was identified by Tapia et al. (2015). The object brightened in the K band by 3.7 mag between 1993 and 2003. Since the outburst was found a posteriori, no information on the possible accretion luminosity is available. The post-burst luminosities range from 2.5 × 10³ L_⊙ (Povich et al. 2011) to ~ 4 × 10³ L_⊙ (Tapia et al. 2015) which correspond to a mass of 8–9 M_⊙ (Tout et al. 1996) for a zero-age main sequence (ZAMS) object. In the context of the forthcoming discussion, it has to be emphasized that this source is not associated with masers. So, while the V723 Car event might be considered to be an accretion burst from a massive young stellar object (MYSO), the lack of observational coverage before and at the time of its incidence precludes major conclusions with regard to high-mass star formation.

Recently the situation concerning MYSO accretion bursts abruptly changed following the discoveries of the almost coincident events from the MYSOs S255IR-NIRS3 (Stecklum et al. 2016; Caratti o Garatti et al. 2017; Liu et al. 2018) and NGC6334I-MM1 (Hunter et al. 2017, 2018). The luminosity increase, seen at infrared (IR) and (sub)mm wavelengths for S255IR-NIRS3 and in the (sub)mm regime for NGC6334I-MM1, provided direct evidence for enhanced accretion rates. Most notably, these outbursts were accompanied by flares of Class II methanol masers (hereafter methanol masers; Fujisawa et al. 2015; Szymczak et al. 2018a; MacLeod et al. 2018). This confirmed a radiative pumping mechanism of this kind of methanol masers (Menten 1991a; Sobolev et al. 1997; Cragg et al. 2005), which is consistent with variability studies of the maser emission (Szymczak et al. 2018b; Durjasz et al. 2019). Since methanol masers trace the very early stages of massive star formation (for example Breen et al. 2013), maser flares might be taken as a proxy for accretion variability of the protostellar host. Keeping this in mind, the international maser community established the Maser Monitoring Organization (M2O)¹ to coordinate single-dish monitoring of masers and interferometric follow-up measurements.

G358.93-0.03 (hereafter G358, RA: 17^h43^m10.^s02, Dec: −29°51″45.′′8, J2000) is a hitherto little explored massive star forming site as evident from just eight SIMBAD (Wenger et al. 2000) entries until 2018. The kinematic distance amounts to 6.75 kpc (Brogan et al. 2019) which implies a galactocentric one of 1.6 kpc. It is consistent with Gaia distances of visible stars in the G358 region of ≲ 5 kpc which impose a lower limit to the distance of the molecular cloud hosting G358 (Burns et al. 2020). In mid-January 2019, flaring of the 6.7 GHz CH₃ OH maser line (Menten 1991b) in G358 was announced (Sugiyama et al. 2019). Thus, for the first time, M2O orchestrated an extensive observing campaign which became extremely successful.

Thanks to the immediate response, an unprecedented wealth of masering lines, including numerous new transitions, could be observed (Breen et al. 2019; Brogan et al. 2019; MacLeod et al. 2019) and new maser species were discovered (Chen et al. 2020a,b). Interferometric imaging in the (sub)millimeter spectral range using the Atacama Large Millimeter/submillimeter Array (ALMA) and the Submillimeter Array (SMA) dissected the star forming region and pinpointed the MYSO which hosts the flaring masers (Brogan et al. 2019). In all likelihood, the brightest continuum source MM1, which turned out to be a hot molecular core, experienced an accretion burst. For the first time, a spectacular confirmation of the event was achieved by high-resolution, multi-epoch observations of the 6.7-GHz methanol maser emission which revealed outward maser spot propagation, tracing the spread of the thermal radiation emanating from the burst (Burns et al. 2020). However, without evidence for a significant rise in (sub)millimeter dust continuum emission from MM1 (Brogan et al. 2019), the energetics of the burst remained an open issue. Therefore, we aimed for observations to identify the IR counterpart of MM1 and to verify its luminosity increase, thus independently confirming the third MYSO accretion burst witnessed so far.

At present, the Stratospheric Observatory for Infrared Astronomy (SOFIA, Erickson & Davidson 1993; Young et al. 2012) is the only facility which offers the capability to trace the far-infrared (FIR) flux increase caused by such an event. Consequently, an attempt wasmade for observing G358 with SOFIA which turned out to be successful. The observational results in context with supplementary data, the analysis and interpretation are the subjects of the present paper and will be outlined in the following.

The paper is organized as follows. At the beginning the observational foundation will be explained. The next part deals with deriving constraints on the spectral energy distribution (SED) of the bursting source. The estimation of the luminosity increase, the central quantity for assessing the accretion burst, is performed in two steps. First, a simplified treatment using graybody functions is applied. Then, a more thorough analysis is performed, utilizing dust continuum radiative transfer. The discussion section concludes the paper in which results of the present investigation are put in context with regard to previous observational and theoretical findings.

We note that the term “luminosity”, if not declared otherwise, refers to the bolometric luminosity throughout the paper. Similarly, the term or value “error” always implies the 1σ error or standard deviation unless noted otherwise. Magnitudes are based on the Vega system.

2 Observations

2.1 Archival data

Archival data is essential for establishing the presence of an accretion burst since it constrains the source luminosity during the pre-burst state. Because G358 is located in the Galactic center region it has been covered by a wealth of surveys. For the present study, fluxes in the near-IR (NIR), mid-IR (MIR), and far-IR (FIR) as well as positions and images have been retrieved from the surveys listed in Table 1

The comparison of the HI-GAL photometry (Molinari et al. 2016) with that performed by the ATLASGAL team (Contreras et al. 2013)revealed inconsistencies at longer wavelengths and with regard to the error estimates. Therefore, we performed photometry on our own on the corresponding PACS and SPIRE images (epoch 2010 September 07), retrieved from the NASA/IPAC Infrared Science Archive (IRSA)², using the MPFit2DPeak function from the IDL Astronomy library (Landsman 1995). The point spread function was chosen to be a 2D Gaussian which yielded the best fit over other representations (Moffat, Lorentzian profiles). The background was estimated on a fixed annulus for all wavelengths to ensure consistency of sky estimates.

The target is not included in the AKARI Bright Source Catalog (Yamamura et al. 2010) but can be identified in Far-Infrared Surveyor images (FIS, Doi et al. 2015) taken with the N60 (λ_central = 65 μm) and WIDE-S (λ_central = 90 μm) filters (epoch 2007 January). Correspondingly, photometry was performed on those frames in the same way as described above. Since the WIDE-S image suffers from severe striping, the resulting flux has a substantial uncertainty and will therefore be neglected.

2.2 Near-infrared imaging

Optical and near-infrared (NIR) imaging of G358 was performed using the seven-channel Gamma-Ray Burst Optical/Near-infrared Detector (GROND, Greiner et al. 2008), using director’s discretionary time (DDT) at the MPG/ESO 2.2-m telescope at La Silla (Chile) on 2019 February 8. GROND obtains images in seven bands (optical: Sloan g’r’i’z’, NIR: JHK_s) simultaneously. The total integration time amounts to 38 min. Data processing was performed by means of the GROND pipeline (Krühler et al. 2008).

2.3 Far-infrared integral-field spectroscopy

The Field-Imaging Far-Infrared Line Spectrometer (FIFI-LS, Looney et al. 2000; Fischer et al. 2018; Colditz et al. 2018) is a far-infrared integral field spectrograph aboard SOFIA. FIFI-LS features a blue and red channel in parallel. Each channel has a field of view (FOV) consisting of five by five spatial pixels with a plate scale of 6′′ /pixel in the blue channel and 12′′/pixel in the red channel. They provide an overall wavelength coverage from 51 to 203 μm. This matches very well the range of the SED, where MYSOs emit the bulk of their energy by dust continuum radiation, and where the relative flux increase due to an accretion burst is highest (MacFarlane et al. 2019). Thus, instruments like FIFI-LS provide the best prospects to detect the luminosity increase due to accretion bursts.

A DDT proposal to perform spectro-photometry of the FIR dust emission of G358 using FIFI-LS was submitted in mid-February 2019. One hour of observing time was awarded by the Director of SOFIA Science Mission Operations. The measurements were performed on 2019 May 1, when the maser flare was still strong but already decaying (MacLeod et al., in prep.; Yonekura et al., in prep.). Several spectral bands (see Table 2) were chosen to sample the full spectral range of FIFI-LS and, at the same time, cover high rotational transitions of the CO molecule. The spectral scan length per sub-band ranged from 0.3 to 1.0 μm. A regular follow-up proposal for SOFIA Cycle 8 was submitted and accepted as well. Due to the flight suspension induced by the COVID-19 pandemic, the observations were delayed and eventually performed on 2020 August 28. The same settings and amount of observingtime as for the first epoch measurements were used. However, because of a technical issue, no data could be obtained in the blue channel. Yet, drawing conclusions on the flux evolution from the two epochs became possible.

3 Data analysis and results

As evidenced by the ALMA/SMA observations, G358 is a star forming complex harboring several massive protostars (Brogan et al. 2019). Because of its considerable distance and the large beam sizes of observing facilities at MIR and FIR wavelengths, the photometry obtained by the latter represents the total flux. The following considerations aim at identifying and characterizing the IR counterpart of MM1. This includes an approach to account for flux contributions of all other components. Eventually, conclusions on the nature of MM1 will be drawn based on radiative transfer (RT) modeling of its pre-burst and burst SEDs.

3.1 NIR imaging and source identification

As shown in Fig. 1, a NIR source is situated close to the position of G358. It is listed as 2MASS J17431001−2951460 in the 2MASS All-Sky Catalog of Point Sources (Cutri et al. 2003) but was detected only in the K_s band at 11.5 ± 0.07 mag. The deeper VISTA Variables in Via Lactea Survey (VVV, Minniti et al. 2017) which imaged the Galactic bulge and the adjacent southern plane from 2010 to 2016 yielded an H-band detection as well. The VVV K_s and H catalog magnitudes of 11.87 ± 0.01 and 15.54 ± 0.06 mag imply a color index of H − K_s = 3.67 ± 0.06 mag which indicates a very red object. In the VVV catalog, the object is flagged as variable as it has shown rapid brightness changes within a 3σ range of 0.4 mag, and a peak-to-peak variation of 0.79 mag, during five years of VVV K_s -band monitoring. Its mean brightness during the VVV monitoring amounts to K_s = 12.23 mag.

Our GROND imaging detected it in the J, H and K_s bands, but not at shorter wavelengths. The GROND Ks magnitude of 11.9 ± 0.02 indicates a brightness, increased by 0.34 mag with respect to the VVV mean, which is consistent with the source variability. Its position, within 0.′′2, is consistent with the secondary hot molecular core and dust continuum source MM3 detected by ALMA, which is located 1.′′ 09 to the southwest of the main hot molecular core MM1 (cf. Brogan et al. 2019). Therefore, the NIR source represents the IR counterpart of MM3 and is, thus, unrelated to the outburst.

The VVV and GROND NIR color composite images are shown for comparison in Fig. 1. The K_s image of the difference GROND−VVV, obtained after flux scaling and convolving to the same spatial resolution, does not provide evidence for the presence of a light echo from an accretion outburst, unlike for the case of S255IR-NIRS3 (Caratti o Garatti et al. 2017). This may be another sign for the high extinction toward the bursting source.

Fig. 1 Top: VVV JHKs color composite of the G358 region (epoch 2010 August 15). The positions of MM1 and MM3 from Brogan et al. (2019) are marked. Black plus signs denote positions from IR observations at wavelengths up to 24 μm. Bottom: GROND JHKs color composite (epoch 2019 February 8) with contours of the ALMA 0.89 mm continuum map (Brogan et al. 2019).

3.2 WISE/(NEO)WISE photometry

Due to its orbit, the WISE IR space telescope (Wright et al. 2010) visits a region in the sky twice a year. For G358 the 2019 visits occurred on March 17 and August 27. The first one almost coincided with the peak of the maser flare (MacLeod et al., in prep.; Yonekura et al., in prep.) which maximized chances for a possible MIR detection. Photometry for G358 from the WISE and subsequent (NEO)WISE (Mainzer et al. 2014) missions were retrieved from IRSA, covering observations until end of 2019. A saturation correction has been applied³ to account for a photometric bias due to warm-up of the detector which amounts to + 0.02 (W1) and + 0.33 (W2) magnitudes, respectively.

The (NEO)WISE W1 (3.4 μm) and W2 (4.6 μm) light curves are shown in Fig. 2. Because of its brightness in both filters bands, G358 led to small (W1) or mild (W2) detector saturation. In this case, the derivation of the brightness from the wings of the point spread function (cf. Sect. IV.4.c.iii of the WISE All-Sky Explanatory Supplement, Cutri et al. 2012) leads to enhanced scatter. Nevertheless, the discovery of an obvious brightening which usually accompanies enhanced accretion should have been possible. However, there is no clear-cut evidence in both light curves for a flux increase at the burst epoch or later on. From the scatter of the photometric values, a possible increase due to the burst can be constrained. Assuming a 2σ detection limit for W1 and W2 of ~ 0.5 mag, that is a joint 2.8σ limit, upperbounds for a possible flux contribution caused by the burst can be derived as 0.10 and 0.45 Jy, respectively.The failure of the burst detection suggests that, at any time, MM3 provided by far most, if not all, emission seen in the (NEO)WISE bands. Further details on the (NEO)WISE astrometry are outlined in Sect. 3.3.2.

3.3 Infrared photo-astrometry

In case of imaging an unresolved crowded region a shift of the emission centroid may be caused by differing object SEDs which would introduce a wavelength dependence and/or by variability leading to a temporal displacement. By placing upper limits on a possible centroid shift, constraints on the contribution of a single source with known position, MM1 in the present case, to the overall emission can be established. In the following this will be done using the present infrared data.

3.3.1 VVV and GROND

An upper limit to the MM1 pre-burst 2.18 μm flux can be derived from the stacked VVV K_s image. It is based on 131 exposures of 4 s each, that is a total integration time of 8.7 min. Confusion noise due to the high surface density of objects in the Galactic center neighborhood limits the detection of a possible faint NIR counterpart of MM1. It has to be brighter than K_s ~ 15.5 mag to be recognized next to 2MASS J17431001−2951460. This corresponds to a flux density of 0.42 mJy (Cohen et al. 2003). Similarly, the burst K_s magnitude can be constrained by the corresponding GROND image. Given the aperture sizes of the 2.2-m and VISTA telescopes as well as the total integration times, the GROND image should almost reach the sensitivity of the stacked VVV frame. However, inferior seeing and slightly elliptical images reduce its depth by about half a magnitude, resulting in a burst 2.15 μm upper limit of 0.66 mJy.

3.3.2 (NEO)WISE

For the following analysis, (NEO)WISE positions were retrieved within 5′′ of the MM1 location from IRSA, based on frames with the photometric quality flag A in both bands, a signal-to-noise ratio ≥ 20 and a frame quality rating of 10. The measurement quality and source reliability information flags of all frames, however, indicate contamination by the nearby 2MASS J17430939−2951517 which is situated 9.′′9 southwest of the MM3 and brighter in both (NEO)WISE bands. This results from the image full width at half maximum (FWHM) of 6.′′ 1 in the W1 band⁴.

The offsets of the individual positions of the IR source with regard to MM1 are shown in Fig. 3 where the data was divided into three groups. The pre-burst one includes all positions until the last visit before the burst (2018 August, blue), the burst group comprises those of the 2019 March visit (red) and the post-burst group represents the 2019 August visit (green). The quantitative analysis confirms the visual impression for the mean positions that both, burst and post-burst, are consistent with the pre-burst location, situated close to MM3. So most of the emission seen by (NEO)WISE arises from the MM3 NIR counterpart 2MASS J17431001−2951460, with some contamination from the nearby IR source to the southwest which causes the elongated distribution of positions. While the mean (NEO)WISE post-burst position (green) is in between MM1 and MM3, its large position error precludes drawing any conclusions from this fact on whether or not this is a late sign of the burst. These findings are consistent with the non-detection of the burst in the (NEO)WISE light curves (Sect. 3.2).

Fig. 2 (NEO)WISE W1 (blue) and W2 (red) light curves based on mean magnitudes and respective errors for each visit. The first two epochs are from the WISE mission. Vertical lines mark the dates of the flare peak (black) and the first FIFI-LS epoch (green). Horizontal lines indicate mean magnitudes (solid) and their errors (dashed).

Fig. 3 (NEO)WISE coordinate offsets relative to an origin at MM1 (from Brogan et al. 2019). Pre-burst, burst, and post-burst data are shown in blue, red and green. Thick error bars represent mean positions while individual appear thin. Positions of MM1 and MM3 (measured with ALMA) are marked by the diamond and the square, respectively. The arrow points toward the bright 2MASS J17431001−2951460, which is at a distance of 9.′′9 from the (0,0) position. The yellow circle indicates the W1 image FWHM of 6.′′ 1.

3.3.3 Spitzer MIPS

As indicated by Fig. 1 (top) the positions of the IR counterparts at wavelengths up to 24 μm are almost coincident with ALMA source MM3. Supposing about equal flux densities of MM1 and MM3 at this wavelength, the centroid of the MIPSGALimage of G358 should be located halfway between both sources. This is clearly not the case. Thus we conclude that the pre-burst MIPS 24 μm flux of MM1 was considerably smaller than the joint flux of MM1 and MM3 (4.58 ± 0.02 Jy, Hinz et al. 2009). An upper flux limit for MM1 can be derived by assuming that it leads to a detectable centroid shift of three times the positional error of 0.′′02 (Hinz et al. 2009). Since such a shift has not been detected, the contribution of MM1 to the total flux must be less than 0.42 Jy. The upper limit is valuable to constrain the pre-burst luminosity.

3.4 SOFIA FIFI-LS spectro-photometry

The FIFI-LS data was processed by the “FIFI_LS_REDUX” pipeline (version 1.7.0, Clarke et al. 2015) and downloaded from the SOFIA Data Cycle System. The observations yielded a clear detection of continuum emission from the target in all bands. G358 is the only object in the FOV. In several bands CO line emission has been detected as well which will be discussed elsewhere.

The continuum fluxes were derived from an image consisting of the pixel-wise median of a spectral data cube along the wavelength dimension, thus free from emission-line flux. For both epochs, central wavelengths and derived flux densities together with an error measure are given in Table 2. We note that our observations revealed the need for an empirical correction of the 118.9 μm flux. It was derived from SED fits of other MYSOs of our FIFI-LS data set which show a similar flux deficit as well in this band. This was also confirmed by looking at data from the flux calibrators Mars and Uranus.

The measurement errors derived from the error images provided by the FIFI_LS_REDUX pipeline do not include the calibration uncertainty. Therefore, we adopted a conservative approach by using an uncertainty of 10%, cf. Fadda et al. (2019). Since a narrow dust emission feature at 69 μm (Sturm et al. 2013) – not covered by our FIFI-LS bands – is the only one in this wavelength region, a low-order polynomial fit seems to be representative for the actual SED at the time of the observing epoch. The residuals from these fits listed in Table 2 correspond to a mean relative error which indeed equals the above uncertainty.

With the SOFIA flux densities at hand, the change of the SED due to the burst and its temporal evolution can be evaluated from Fig. 4 which shows, among others, the pre-burst SED (blue symbols), based on MIPS, AKARI, and Herschel data as well as the emission-line corrected ATLASGAL 870 μm flux from Brogan et al. (2019). The SED based on our first-epoch FIFI-LS observations (epoch 2019 May 1) and the ALMA 870 μm integrated G358 flux from Brogan et al. (2019) (epoch 2019 April 12) is shown in red. The burst SED features flux densities larger than a factor of ≳2 when compared to the pre-burst ones and a possible change of the SED shape. Thus, a luminosity increase – the prime signature of an accretion burst – has been witnessed with SOFIA. It represents the second confirmation of such an event after the S255IR-NIRS3 burst (Caratti o Garatti et al. 2017) using this unique facility. Given the non-detection in the NIR and MIR as well as the non-significant flux increase in the (sub)mm (Brogan et al. 2019) this promotes the G358 event to be the first NIR-, MIR- and (sub)mm-dark but FIR-loud MYSO accretion burst.

While the second epoch SOFIA observations suffer from the lack of the blue FIFI-LS bands, which sampled the flux density at wavelengths short-ward of the SED peak, the red channel data (green symbols in Fig. 4) nevertheless unambiguously indicate elevated flux levels during the post-burst stage compared to pre-burst and depressed compared to the bursting stages.

Fig. 4 Observed FIR/(sub)mm SEDs showing pre-burst Herschel/ AKARI data (blue) and FIFI-LS burst values (red) together with the corresponding (sub)mm fluxes from Brogan et al. (2019). The upper 24 μm MIPS limit is indicated as well. The solid lines represent the reddened graybody fits (cf. Sect. 4.1). The second epoch FIFI-LS data and the corresponding fit appear bf in green.

3.5 FIFI-LS image analysis

The angular distance of 1.′′09 between MM3 and MM1 at a position angle (PA) of 247° warrants to investigate whether both sources could be marginally resolved at the shortest FIFI-LS wavelengths, provided they have similar flux contributions. The diffraction limit for the 2.5-m mirror of SOFIA at 52.2 μm amounts to 5.′′3. However, due to various reasons (Graf et al. 2017), the FWHM of the actual point spread function (PSF) exceeds that value.

In order to address this point, we cannot rely on the absolute pointing accuracy of SOFIA but have to analyze the image morphology. So the median continuum images of the blue channel were fit by a bivariate (2D) Gaussian. If MM3 contributes a substantial flux contribution to the image, the major axis of the fit should clearly exceed the smaller one and be aligned to that direction. The corresponding quantities, axis ratio, position angle and respective errors, are given in Table 3. When judging these results, it has to be kept in mind that the FIR observations were performed in chop-nod mode. Thus, deviations from a perfect image superposition may lead to an elongated image as well. The position angle (PA) of the major axis for the 52.2 μm band is closer to that of MM3 compared to the other ones. However, the PAs of all bands are more aligned with the chopping angle of 293°. Given this fact, it remains open whether the result for the shortest wavelength of the blue channel might be interpreted in favor of a noticeable contribution from MM3.

3.6 PACS image analysis

We checked the Herschel-PACS images at 70 μm for potentialsigns of source multiplicity or source elongation. Unfortunately, the only existing data toward this region (from Hi-GAL, see Molinari et al. 2010) was obtained in the so-called fast-scan parallel-mode which comes with a quite peculiar PSF that is elongated along the two almost orthogonal scan directions. Its equivalent FWHM is typically 8.′′ 8 × 9.′′6 and is thus larger than the SOFIA PSF at the same wavelength. This impedes a study of small elongations due to source multiplicity. We nevertheless attempted a PSF photometry using the Starfinder IDL tool (Diolaiti et al. 2000), Version 1.8.2a, by using the adequate PACS PSFs (version 2.2) for the fast-scan parallel mode, provided by the instrument team⁵. The algorithm can in principle find multiple sources blended within one FWHM. While we varied several parameters for background subtraction and treatment of sub-pixel offsets, the program consistently identified just one strong source, showing a very good correlation coefficient of 0.977 with regard to the input PSF. This source is formally 1.′′ 1 away from the ALMA position of component MM1, and 0.′′7 away from ALMA component MM3. The absolute astrometry of the Herschel Spacecraft at the end of the mission came with a 3σ uncertainty of 2.′′ 7 (cf. Sánchez-Portal et al. 2014). Therefore, based on the Herschel astrometry alone, it is not possible to decide whether the revealed Herschel compact source is coincident with the location of MM1 or MM3, or somewhere in between. Based on the SED modeling, we think that the Herschel 70 μm emission is dominated by MM1. Given our PSF fitting experiments we conclude that a potential minor contribution at 70 μm from MM3 must be more than 70 times weaker than the one from MM1.

4 Graybody fits and burst parameters

4.1 Graybody fits to SEDs based on FIR and (sub)mm data

Since the observed SEDs of G358 based on FIR and (sub)mm data do not strongly differ from a Planck function, the most simple approach to approximate them is using a modified blackbody, in other words a graybody (for instance Elia & Pezzuto 2016). It comprises three parameters: dust temperature, emissivity index, and solid angle of the emission. Before performing the fits the fluxes were dereddened to account for interstellar extinction. A value of A_V = 60 mag was assumed which seems to be justified (see Sect. 5.2.2). For flux dereddening, the R_V = 5.5 dust model of Draine (2003) has been used. The fits, reddened to match the observations, are shown as solid lines in Fig. 4.

The corresponding fit to the pre-burst SED yielded a dust temperature of K which is slightly lower than the pre-burst value of 28.5 ± 1.5 K by Brogan et al. (2019). The dust emissivity index amounts to β = 1.83 ± 0.06 which is consistent with the emissivity index close to the center of spiral galaxies in the Local Group (Mattsson et al. 2014; Tabatabaei et al. 2014) and appropriate for temperature reasons (Boudet et al. 2005).

Since it seems plausible that the dust properties remain unchanged during the burst, except for the innermost region where modifications may have occurred (for example Ábrahám et al. 2009), the pre-burst emissivity index was used for the burst- and post-burst fits.

As expected, the dust temperature for the burst SED is marginally higher K as well as the solid angle of the emission which increased by a factor of 1.20 ± 0.14. However, the χ² of the burst fit is 2.5 times as large as that of the pre-burst one. This is caused by the lack of a pronounced peak in the burst SED. Since our assessment of the FIFI-LS data did not reveal any issue other than that for 118.9 μm flux, which has been corrected by our recalibration, the flat FIR burst spectrum is probably real. Future time-dependent RT modeling may show whether such a feature can be produced by a burst of short duration.

An attempt was made to come up with a graybody fit for the second FIFI-LS epoch data as well. Because of the sparse data, the pre-burst emissivity and the mean solid angle from the pre- and burst fits were adopted and only the temperature was varied. The resulting temperature amounts to K.

4.2 Empirical burst parameters

The FIR/(sub)mm luminosity increase due to the MM1 accretion burst can be determined by integrating the above graybody SED fits and taking source distance given in Sect. 1 as well as interstellar extinction into account. Here we assume that all other sources which contribute to the total G358 flux stayed constant during the pre-, burst- and post-burst epochs. The weak variability of MM3 in the NIR and MIR is of no concern here since its FIR emission is far below that of MM1 (see Sect. 5). Since the bulk of the energy is emitted in the FIR, the FIR/(sub)mm luminosity estimates only weakly depend on A_V. The impact of the distance ambiguity is considered in the following analysis.

Integration of the graybody fits yields the following FIR/(sub)mm luminosities: pre-burst , first FIFI-LS epoch , and second epoch , respectively.As emphasized, without interstellar extinction, these values shrink but only by 9%. The luminosity increase due to the accretion burst at the dates of the FIFI-LS observations amounts to and . The presence of a substantial luminosity increase during the second FIFI-LS epoch, meaning about 18 months after the peak of the maser flare, is remarkable.

For deriving major parameters of the burst, we follow the approach of Caratti o Garatti et al. (2017). In addition, for what concerns the estimate of the burst energy, the two epochs of FIFI-LS observations provide the opportunity to account for the temporal evolution of the luminosity. The simplest, yet plausible, approach is to assume a linear decrease which holds from the flare peak date to the date when the pre-burst level will be reached again. The linear flux decay approximation yields a duration Δt of 869 ± 303 days, with apre-burst-level return date of 2021 July 31.

The FIR/(sub)mm burst energy is , where is the average luminosity increase which equals to half of the peak increase for a linear drop to zero. It amounts to J. Its upper and lower bounds from the uncertainties in both duration and luminosity are 1.0 × 10³⁸ J and 2.6 × 10³⁸ J, respectively.

We emphasize that these values represent lower limits to the luminosity increase and the corresponding energy release since they are based on the FIR/(sub)mm emission only. The ultimate quantities will be derived from the results of the RT analysis below. Accretion related quantities require the knowledge of protostellar mass and radius and will be considered in Sect. 6.

5 Analysis of SEDs

5.1 SED decomposition

In order to derive representative parameters of the bursting MYSO MM1 from RT modeling, the underlying SED should be as free as possible from contributions from other objects. Thanks to the availability of NIR, MIR, and (sub)mm photometry for MM3, the contamination from this source to the overall flux can be removed by using predicted flux densities from its best RT model (see Sect. 5.2.1), assuming that its SED has not changed. This approach has been proven successful in a similar, yet less sophisticated fashion, for a study of FIR emission from W3(OH) and the neighboring W3(H₂O) (Stecklum et al. 2002).

Here we extend it by also taking into account the presumed contributions from the remaining sources detected at (sub)mm wavelengths. Since none of these has an IR counterpart we assume that they are in an early evolutionary stage similar to MM1, with SEDs that resemble the overall pre-burst SED. So for MM1 and each of those (MM2, MM4-8), the (sub)mm fluxes from Brogan et al. (2019) were fit by a graybody, using the temperature and emissivity index derived from the pre-burst SED, to obtain the individual solid angles of the emission.

For the pre-burst epoch, the following steps were performed to obtain the MM1 fluxes. First, the predicted MM3 fluxes were subtracted from the G358 total fluxes. Second, the ratio between the MM1 solid angle and the sum of all solid angles was calculated which amounts to 0.50 ± 0.05. It corresponds to the relative contribution of MM1 to the MM3-subtracted fluxes. So, by multiplying the latter with this ratio the pre-burst fluxes of MM1 were obtained.

Similarly, for the wavelengths of the burst as well as post-burst observations, the pre-burst fluxes were predicted from the MM3 modelas well as the pre-burst graybody fit. Then, the MM3 contribution was removed. Finally, the contribution of MM2+MM4-8 to the MM3-subtracted pre-burst fluxes had to be taken out from the observed burst as well as post-burst fluxes to yield those for MM1. These are listed in Tables A.2–A.4 (for MM1 pre-burst, burst, and post-burst respectively). They represent the best possible approximation of the intrinsic ones of MM1. We emphasize that the procedure to derive the MM1 fluxes is tailored in order to reproduce the total flux. Therefore, the MM1 (sub)mm fluxes exceed those given in Table 2 of Brogan et al. (2019) by a factor of about two.

Fig. 5 Pre-burst SEDs for the following sources: black symbols – total FIR/(sub)mm emission of G358, green – MM3, blue – MM1, derived from the total FIR/(sub)mm emission by removing contributions from all other sources (including MM3). For MM1 and MM3, the observed values are shown together with the ten best RT fits. At wavelengths beyond 40 μm MM1 dominates the total flux density. The seemingly gap at 10 μm in the MM1 SED is due to strong silicate absorption.

5.2 Radiative transfer analysis

Characteristic properties of YSOs, namely their luminosities as well as mass, geometry and extent of the surrounding dust, can bederived by modeling their dust continuum radiation to match the observed SEDs (for instance Wolf et al. 2012; Whitney et al. 2013). For G358 this has been done by Brogan et al. (2019) to infer the pre-burst luminosity L^pre, using the YSO model grid of Robitaille (2017).

Utilizing the same model pool, we performed an RT analysis of the SEDs of MM1 and MM3 using the Python implementation “sedfitter”⁶ of the SED fitter (Robitaille et al. 2007) which is described below. We did not fit the G358 total flux (black symbols in Fig. 5), since models with multiple sources, which would be appropriate for massive star forming regions and clearly for G358, are not included in the Robitaille (2017) model grid. Instead, we extracted the fluxes of MM1 from the total ones as described above. From the ten best fits (per epoch/source), mean values and uncertainties for all “free” parameters of the models are derived (see Appendix A). By using a weighting with the respective χ² values, we ensure that the best fits determine the corresponding mean values. For the parameters which are log-spaced in the Robitaille model pool (r_min, r_max, via R_⋆, T_⋆), we use the geometric mean, while for all other parameters, we use the arithmetic mean. Since each model is comprised in the pool with nine different inclinations from 0° to 90° to account for the inclination dependence of the SED, the ten best fits are not necessarily composed by ten different models. Instead, some models might be included more than once, but with different inclinations.

We note that these models incorporate passive disks. While this may seem inappropriate in the context of accretion bursts where active disks are often invoked, it is justified by the fact that we are primarily interested in reproducing the FIR and (sub)mm emission. Due to the strong radial dependence of their viscosity, for instance Pringle (1981), active disks differ from passive ones primarily in the very innermost region where the bulk of the dissipated energy is being released. The details of this process are not relevant for the highly reprocessed emission in the FIR and (sub)mm range, which is dominated by dust radiation from the outer regions.

Before describing the actual modeling, we emphasize that additional data and results are provided in the Appendix A. These comprise the flux tables which were used to establish the SEDs, a summary of the RT models which have been used in the analysis, tables listing the parameters of the ten best RT models for each of considered cases, and a table providing optical depths for the ten best pre-burst RT models of MM1.

5.2.1 Modeling of the dust continuum radiation of MM3

For the fitting of the MM3 SED, the combined “spubhmi” + “spubsmi” data sets from Robitaille (2017) have been used. They comprise 120 000 YSO models with nine inclinations for each model (leading to a set of 1 080 000 SEDs in total). All models employ Milky Way dust with R_V = 5.5 and a grain size distribution from Weingartner & Draine (2001).

Their designations are based on the respective model components. Each model (in both data sets) comprises the following components: star, passive circumstellar disk, bipolar cavity, Ulrich-type envelope (Ulrich 1976), and ambient medium. The “spubhmi”-setting features an inner hole in between star and disk, whereas for “spubsmi” the inner radius is governed by the dust sublimation radius (assuming T_sub = 1600 K, in accordance with Robitaille et al. 2006). We note that the Robitaille model pool includes other data sets, that are composed of less (or different) components and are thus less suited to represent the structure of a YSO at a very early evolutionary stage. For further details, we refer the reader to Robitaille (2017). A synthetic aperture size of 3′′ has been used for fitting the SED. The fluxes of the MM3 SED are listed in Table A.1.

We adopt the same distance range as for MM1 but allow for a smaller interstellar extinction. Figure 5 shows the SED (green) together with the total pre-burst SED (black symbols) and the pre-burst flux attributed to MM1 (blue), where the contribution of the other sources has been removed as described in Sect. 5.1. The ten best fits for each SED are shown with solid lines. The best fit is dark, whereas the other are slightly transparent. Our models suggest that the contribution of MM3 to the total flux at wavelengths beyond λ ≥ 40 μm is marginal. Nevertheless we take the best MM3 fit into consideration when refining the MM1 fluxes in the following analysis.

As indicated by its moderate obscuration in the NIR already: MM3 is likely the most evolved object of the G358 complex. It features a disk with a dust mass of . This may seem heavy keeping the canonical gas-to-dust mass ratio (γ) of 100 in mind. However, the latter is not applicable since γ depends on the galactocentric distance R_GC (Giannetti et al. 2017, Eq. (2)) because of the Galactic metallicity gradient. For G358, R_GC=1.6 kpc implies a value of γ = 38 ± 4 which yields a total disk mass of . The fit delivers A_V = 20 ± 5 mag and an inclination of i = 51 ± 9°. The inner radius of the disk seems to be close to the star, r_inner < 3.5 R_sub holds for each of the ten best models. The mean value of the disk outer radius amounts to au, while the maximum is as high as 1600 au. The mean luminosity of corresponds to a ZAMS star of 11 M_⊙ and 4.2 R_⊙ (Tout et al. 1996). While most of the models have luminosities below 10 000 L_⊙, one has a luminosity of 43 000 L_⊙. Such a high effective temperature would imply the presence of a compact HII region. However, MM3 escaped the detection in the radio continuum at a sensitivity level of ~ 50 μJy beam⁻¹ in the survey of Hu et al. (2016) which probably rules out this model. Recently, Bayandina et al. (in prep.) succeeded to detect faint radio continuum emission from MM3. The whole parameter set for the ten best models (including their respective χ² values) can be found (together with the weighted means and standard deviation σ) in Table A.5.

As described in Sect. 5.1 the MM3 SED fit is used for the SED decomposition. In order to check the robustness of the FIR part of the best models, two more fits were performed with reduced data sets. At first, the ALMA fluxes were omitted which yielded slightly lower FIR fluxes compared to the nominal fit described above. Then, also the WISE W4 and MIPS 24 μm fluxes were dropped. The best model for this SED shows fluxes exceeding those of the nominal fit (cf. Fig. A.2). Thus, we conclude that the FIR part of the nominal MM3 SED fit is well constrained by the incorporated fluxes longward of 20 μm.

5.2.2 Modeling of the dust continuum radiation of MM1

With the refined MM1 SEDs at hand, a comparison of the results from RT modeling for the pre-burst, post- and burst-states becomes possible. Before describing this analysis, a remark must be made concerning the interstellar extinction A_V which is a free parameter of the SED fitter. Because of the wavelength dependence of the dust optical properties, extinction is most pronounced at short wavelengths. Thus, for SEDs like that of MM3 it can be well constrained by the best models. However, for the SED of MM1 which is almost exclusively defined by FIR and (sub)mm measurements, interstellar extinction is less influential and, therefore, harder to derive. Since higher A_V implies larger source luminosities to reproduce the observed fluxes, an upper limit needs to be established to constrain L. The recalibration of a Galactic dust-reddening model based on IRAS and COBE/DIRBE results by Schlafly & Finkbeiner (2011) suggests a value of A_V = 115 ± 4.2 mag along the sight line toward G358. Since this holds for the whole path across the Galaxy, while G358 is in front of the Galactic center region, the actual A_V will in fact be lower. Therefore, we assumed an interstellar extinction range of A_V = 30−70 mag.

To begin with, we fit the pre-burst MM1 SED which represents the stationary state using the distance and interstellar extinction ranges given above. This established the ten best pre-burst fits. For this purpose, the “spubsmi” data set from Robitaille (2017) has been used which includes 40 000 models at nine inclinations (360 000 SEDs in total). For these models, the inner radius of the dust disk corresponds to the sublimation radius R_sub. This constraint seems to be justified given the high accretion rates at which massive stars form (for example Zinnecker & Yorke 2007; Kuiper et al. 2016), and in particular before and during an accretion burst. It requires a few remarks concerning the understanding of “heating” and “cooling” in this case, since the dust disk cannot get any hotter than the sublimation temperature. What happens due to a luminosity increase is that enhanced dust sublimation at the inner rim R_sub pushes the latter outward. This is accompanied by an adjustment of the temperature profile T(r) via absorption and re-emission such that, for a given radius (beyond the actual R_sub), the temperature exceeds the previous value. Conversely, for a given temperature, the growth in radius implies a larger radiating area and, therefore, leads to a flux increase. The reverse process happens once the burst ceased. R_sub will shift back inward, allowing the dust replenishment by accretion and/or dust reformation. This will lead to a flux decrease and might be considered as “cooling” but, yet, the inner rim is still at the dust sublimation temperature.

The next step was to establish a new model pool from the ten best pre-burst-SEDs which serves to fit the post-/burst SEDs. This pool contains 100 SEDs in total, where we reuse the best 10 models, with a source luminosity increasing in nine linearly spaced steps from 2 to 6 L^pre, respectively. Inclination, interstellar extinction and the distance were set to the values of the underlying model from the pre-burst-fit since none of these parameters is expected to change during the burst. Additionally, the original pre-burst-SEDs are included. Since the inner disk radius is governed by the dust sublimation temperature, assumed to be 1600 K, it shifts outward for those models with increasing luminosity. Otherwise, the system geometry is kept the same. While this ignores possible changes of the disk due to the burst, for instance in the density structure, it is the simplest approach to model the dust continuum emission due to the accretion burst, and feasible to be treated using common RT codes which are generally static. While a proper burst modeling requires time dependent RT, possibly coupled to hydrodynamics for utmost consistency, our simplified treatment nevertheless allows us to draw major conclusions.

For computing the burst models, the Hyperion code (Robitaille 2011) has been used. The so obtained database is the foundation to fit the SEDs of burst and post-burst epochs. We note that by constraining the model pool, we ensure consistency of the results, which means the best fits for pre-burst, burst, and post-burst SEDs are based on the same models. A sketch of all SEDs that are included in burst model database, can be found in the Appendix (Fig. A.1.) The range of the observed fluxes from the burst- and post-burst-epoch is well covered by the models. Only the burst observations at λ = 163 and 186 μm are outside of the flux range covered by the models. Only one model out of the ten best burst fits has the maximal luminosity increase included in the pool. Therefore, it is not necessary to include models with higher source luminosities.

Before presenting the results, a few remarks have to be made concerning the fitting. We exclude the sub-mm observations at λ > 890 μm from the SED-fit of the burst since their deviation from the stationary models is biggest at those wavelengths (see Sect. 9.1). The post-burst is observed only at wavelengths greater than 118 μm, meaning that there are no constraints in the MIR. Since the maser flux did not fall below its pre-burst level until now (MacLeod et al., in prep.; Yonekura et al., in prep.), we assume that the post-burst flux in the MIR is also not below the pre-burst level. Therefore, models fromthe post-burst fits that have MIR fluxes below those of the mean pre-burst model were excluded. An aperture size of 3′′ has been applied in for the fits. With this choice the χ²-value of the pre-burst fit becomes smallest. All obtained parameters are stable against a variation of the aperture size (they agree within their respective errors).

The results of the RT modeling of the MM1 SEDs are visualized in Fig. 6, which shows the ten best fits to the pre-burst (blue), burst (red) and post-burst SED (green). The best model is shown with the darkest line, while the nine other ones are slightly transparent. In addition to the ten best fits, the so called “mean”-model is shown in black. This model is computed using the weighted mean parameters from the combination of pre-burst, post-burst and burst fits (see below). The stellar parameters, which are obviously changing during the burst, are defined by the pre-burst models alone. The mean model is not only shown for the pre-burst, but also for source luminosities corresponding to the mean values during burst and post-burst. As expected, the mean model lies within the range, spanned by the ten best models at each epoch. The results are summarized in Table A.6.

The RT modeling indicates that during the burst the luminosity increases from pre-burst level of to . At the post-burst epoch, the luminosity isstill elevated at a level of .

This corresponds to a relative luminosity gain by a factor of during the burst and by for the post-burst epoch. The increase in luminosity during the burst amounts to which implies a total burst energyof J, where the calculation has been done similar to Sect. 4.2. The derived linear decay time of d predicts an expected return to the pre-burst luminosity in 2021 September, which agrees with the previous estimate based on the gray-body fits within the errors (see Sect. 4.2).

A comparison of the RT results to those obtained from graybody fits in Sect. 4.1 shows that the luminosity increase and energy release from the RT models are indeed higher, although the luminosities at burst and post-burst epochs agree within their respective errors. The main reasons are the inclusion of more energetic emission from hotter regions in the models which lacks in the graybody fits and the correction of the flux contribution from the other G358 sources as outlined in Sect. 5.1.

The estimated parameters from the RT modeling are given below. The interstellar extinction indicated by the fits is A_V = 60 ± 10 mag. The system has a low inclination of 22 ± 11°, which means it is seen close to pole on. This viewing geometry agrees with that of the spiral-arm accretion flow of Chen et al. (2020b) with i = 25 ± 10°. The stellar radius amounts to . It has been obtained from the pre-burst-models alone. The burst might cause an increase in the stellar radius (bloating) although it is unclear, whether the protostar will respond on time scales of a few months. The models show a considerable scatter in the derived disk properties. The derived outer radii are between 140 and 3800 au, with a mean value of au. The favored dust mass of the circumstellar disk is as low as 8.4 × 10⁻⁵ M_⊙, where the 1σ-confidence interval extends from 1.1 × 10⁻⁶ to 6.1 × 10⁻³ M_⊙. Using the appropriate gas-to-dust mass ratio (see Sect. 5.2.1) this corresponds to a most probable total disk mass of 3.2 × 10⁻³ M_⊙ within an uncertainty range of 4.4 × 10⁻⁵ to 0.24 M_⊙. Remarkably, the corresponding total disk masses are lower, and mostly much smaller, than those derived for MYSOs from SED fitting (for example Kraus et al. 2010; Johnston et al. 2015). Among the best models there is only one (90Yt0exl_03) with a total disk mass of 1.4 M_⊙ which is not that much below present estimates. In contrast to MM3, the much higher scatter of the disk mass of MM1 compared to the other log-spaced parameters indicates that it cannot be reliably estimated by the fitter. In other words, a substantial flux contribution from a disk is not necessarily required to reproduce the SEDs. Thus, while we cannot reliably estimate its mass, we may conclude that the disk is quite likely lightweight.

We note that during the burst the disk dust mass might decrease due to dust sublimation. The total mass will be unaffected, assuming that all sublimated dust just increases the gas mass. The whole parameter set for the ten best models for each of the three epochs (including their respective χ² values) can be found (together with the weighted means and standard deviation σ) in Table A.6.

The above given properties have been derived from the results of all MM1 SED-fits, pre-burst, post-burst and burst. The respective χ² -values were used to weight the obtained parameters, similar to what has been done for MM3. We normalized the sum of the respective weighting factors to unity, split up to 0.5 for the pre-burst and 0.25 for the other two epochs. With this we ensure an equal contribution of the fit to the “stationary”- (pre-burst) and “non-stationary”-system (burst and post-burst), where the “non stationary”-contribution is obtained taking into account burst and post-burst epoch equally. The χ² -values of the burst models are about 5 times higher than for the pre-burst. This might be related to the scatter in the FIFI-LS data, to a coarser sampling of the model-parameters (because of the fact that we only consider models that fit the pre-burst SED pretty well), and to differences in comparison to the static case (whereas the pre-burst-system most probably is stationary, this does not hold forthe post-/bursting stage). We note, that in the burst case the relative weights of the models are on the sameorder, whereas in the pre-burst case they differ by a factor of 1.5 at most (for the post-burst it is 1.3).

Fig. 6 Modeled pre-burst MM1 SEDs (blue), together with its burst (red) and post-burst SED (green). Triangles mark upper limits. The weighted mean-model (see text) is shown in black. It is shown three times but with different source luminosities for the respective observing epochs. Since the 870 (pre-) and 889 μm (post- and burst) fluxes are very similar, their symbols cannot easily be distinguished.

6 Stellar-dependent burst parameters

The RT modeling of the MM1 SEDs aimed at determining the luminosities of the MYSO for the various epochs to infer the total energy released by the accretion burst. For deriving the accreted mass and the mass accretion rate, protostellar mass and radius have to be known. For YSOs approaching the ZAMS, the corresponding values which match the luminosity are a good approximation. At first, we use this approach to obtain massand radius estimates which reproduce the MM1 pre-burst luminosity, assuming solar metallicity. For , these correspond to and (Tout et al. 1996), respectively. The accreted mass is inferred from E_acc = GM_*M_acc∕R_*, where G is the gravitational constant, M_* is the stellar mass, M_acc is the accreted mass, and E_acc the released energy derived in the previous section. Here it is implicitly assumed that the potential energy is released as well when the matter eventually reaches the protostar.

Finally,the mass accretion rate will be obtained from Ṁ_acc=M_acc∕Δt_acc. We have to emphasize here that, generally, the duration of the enhanced accretion Δt_acc will be shorter than that of the elevated emission Δt which has to be used to come up with an overall burst energy estimate. Since the maser excitation is due to MIR dust emission, the total maser flux can be taken as a proxy for the accretion strength (see Sect. 9.4). From the rise and fall of the maser light curve (MacLeod et al., in prep.; Yonekura et al., in prep.), an effective duration Δt_acc of about two months can be derived with a presumed uncertainty range of ± 5 days. With the above quantities we obtain , and , where the errors are dominated by the error of Δt_acc and E_acc. Using the ZAMS mass-radius relation the range of the accretion rate for a given average luminosity increase is shown in Fig. 7.

However, since MM1 is likely in an earlier evolutionary stage, preceding the ZAMS, the above assumption may not hold. A different and presumably more realistic approach is possible using the stellar radius from the RT modeling of the pre-burst SED together with the stellar mass of 12 ±3 M_⊙ derived from the kinematic model of the spiral-arm accretion flows (Chen et al. 2020b). This leads to , and . The large positive error range is mainly due to the corresponding large uncertainty of the stellar radius. To put this into perspective, during its short burst G358 MM1 consumed about 180 Earth masses. Notably, because of the small disk mass, the accreted fraction represents 16% of the total. This raises the question whether the lightweight disk is a stable or transient feature.

Fig. 7 Dependence of the derived accretion rate on the average luminosity increase ⟨ ΔLacc ⟩ and assumed ZAMS stellar mass for the parameter range bracketing MM1. The horizontal dashed line corresponds to ⟨ ΔLacc ⟩ and the lightgray region marks its 1σ uncertainty while the vertical dashed line and the darker gray region indicate the stellar mass and its uncertainty. The values of the accretion rate are indicated.

Fig. 8 Temperature distribution of the mean model in the x-z plane (innermost part of first quadrant) for the pre-burst (left), burst- (center) and post-burst epochs (right). The orange and red lines enclose the temperature range from 94 K (orange) to 120 K (red) during the pre- and burst epochs, respectively. The white contours mark gas particle volume densities of which decrease with increasing z. The vertical solid black line indicates the outer radius of the disk. The dashed black lines mark the radius of the maser ring from the first and 4th epoch of the VLBI observations (Burns et al. 2020, and in prep.). The length of the black bar corresponds to 50 mas.

7 Methanol maser relocation

It has been emphasized that methanol masers play a crucial role in identifying accretion bursts from MYSOs. The maser activity of G358 during its burst was extraordinary and unique in several aspects (Breen et al. 2019; Brogan et al. 2019; MacLeod et al. 2019). The excitation of new maser spots at larger distance has been observed for the first time during the accretion burst of S255IR-NIRS3 (Moscadelli et al. 2017). Notably, for G358 a ring-like propagation of maser spots, likely excited by the heat wave due to the burst, has been witnessed (Burns et al. 2020). Further evidence for spatial changes of the G358 maser distribution during and after the burst has been gained (Bayandina et al., in prep.).

Both methanol in the gas phase and the proper IR radiation are two of the major requirements for maser excitation. While cosmic ray sputtering of dust grain mantles can also release methanol to the gas phase (Dartois et al. 2020), thermal desorption due to grainheating will be the dominant process near the MYSO. It requires temperatures of at least 94 K (Luna et al. 2018) while the maximum desorption rate is reached at about 120–125 K (Collings et al. 2004). For such temperatures, the peak flux of the thermal IR radiation from the warm dust is in the proper wavelength range needed to excite these masers (Ostrovskii & Sobolev 2002; Cragg et al. 2005). The RT models of MM1 during the pre-burst, burst, and post-burst epochs described in Sect. 5.2.2 yielded spatial dust temperature distributions which can be used to address which regions of the circumstellar environment are potential sites of Class II 6.7 GHz methanol masers, and how they change due to the burst.

Figure 8 shows the central part of temperature distribution of the mean model in the first quadrant of the x-z plane for the pre-burst (left), burst (middle) and post-burst (right) epochs. The white lines mark the following gas densities (density decreases with z). To transform the dust densities from the RT-simulation to the above given gas densities, we use the relation , where is the molar mass of H₂, N_A is the Avogadro constant and γ is the dust to gas ratio introduced in Sect. 5.2.1. Because of its low mass, the disk, which is embedded in the rotationally flattened envelope, cannot be easily recognized. The vertical solid black line indicates the outer radius of the disk.

For gas densities exceeding 10⁸ cm⁻³, the maser brightness drops rapidly due to collisional de-excitation (Cragg et al. 2005, Fig. 2). Therefore, within the densest regions (innermost contour - disk mid-plane, envelope at the centrifugal radius) the excitation of Class II methanol masers is basically ruled out.

The orange and red lines of Fig. 8 indicate temperatures of 94 K (orange, minimum temperature for thermal methanol desorption to occur) and 120 K (red, optimum temperature for methanol desorption) during the pre-burst, burst, and post-burst, respectively. At each epoch the temperatures right of the solid orange line are too low for desorption of methanol, which means that it remains bound within icy dust grain mantles. Thus, these lines represent the methanol snow line beyond which masers cannot occur. Due to the heating by the burst, it will be shifted outward and move inward again once the disk started to cool after the burst ceased.

During the pre-burst stage (Fig. 8 left), possible maser sites are likely confined to a region which originates in the disk around ~ 300 au and stretches along the surface of the cavity wall. The situation is different for the burst epoch (Fig. 8 center). As expected, the methanol snow line moved outward and is located at about 1000 au. At the same time the region below the disk/envelope surface has become warmer and presumably got enriched with gaseous methanol. These conclusions from the RT modeling can be compared to the multi-epoch VLBI observations of the expanding maser ring (Burns et al. 2020, and in prep.). The vertical dashed black lines mark the circumference of the maser ring from the first and 4th VLBI epoch. The first epoch observations were obtained already two weeks after the flare started. It is reasonable to assume that during the very early flare rise the excitation conditions were not too far off the pre-burst case. In fact the location of the maser circumference during the first epoch is within the region where phase transition from the solid to the gaseous methanol seems to occur in the pre-burst state. The data of the fourth VLBI epoch was taken about three weeks after the first FIFI-LS observations (“burst” epoch). The maser circumference at that time is confined to the desorption temperature interval for regions not too far off the disk plane. The fact that the extent of the maser ring matches the expected maser positions for both epochs suggests that our RT modeling, although static, nevertheless describes major changes of the circumstellar environment due to the burst properly.

Moreover, Fig. 8 seems to indicate that at the first VLBI epoch, the maser sites were likely close to the interface between the outflow lobe and the disk/envelope, that is in a region of relatively low optical depth. The subliminal maser propagation speed reported by Burns et al. (2020) implies substantial optical depths. This suggests that the masers likely propagated into the disk rather than on its surface. If so it could be expected that the ring expansion slowed down over time. This conclusion is fairly speculative and its verification requires time-dependent RT simulations.

The right panel of Fig. 8 indicates the readjustment of temperature distribution after the burst. Thus, the desorption range moved inward and is basically in between those for both pre-burst and burst.

Another important parameter for the maser excitation is the specific column density, N_M ΔV⁻¹, where N_M is the methanol column density along the line of sight and ΔV the velocity range of the molecules. The low inclination i obtained from our RT modeling and by Chen et al. (2020b) implies that ΔV is only slightlylarger than thermal line-width since the radial velocity component of the orbital motion will be small. This increases the prospects for maser to occur.

When assessing these results it should be kept in mind that the circumstellar environment of G358 is almost certainly more structured than our RT model. Thus, sight lines of lower optical depth as well as density fluctuations may leave an imprint on the actual distribution of the maser spots.

8 Wavelength and time dependence of the flux variation

The RT results for G358 can be used to address the question which wavelength range is suited best to detect the flux increase induced by a MYSO burst. For this purpose, the mean model SEDs for burst and post-burst were normalized by the pre-burst model SED. It has to be kept in mind that dividing the flux values cancels the extinction correction. The result is shown in the upper panel of Fig. 9. It can be seen that the peak values of the relative flux m exceed the luminosity gains derived in the previous section. For our MYSO models, the largest relative flux increase occurs around 10 μm. This agrees with the results of MacFarlane et al. (2019) for similar luminosity gains. For their low-mass YSO models, the maximum flux ratio shifts toward the FIR for stronger accretion bursts. The least relative flux increase occurs in the (sub)mm. This is because the emission in this part of the spectrum arises from the Rayleigh-Jeans tail of the Planck function, where the spectral radiance only linearly depends on temperature. Since the temperature increase of the outer circumstellar regions is lower than for the inner ones, their flux increase is accordingly smaller.

The wavelength dependence of the relative flux increase agrees qualitatively with that of the outburst of OO Serpentis (OO Ser), as found by the first YSO outburst monitoring which covered the full wavelength range from NIR to FIR (Kóspál et al. 2007). These observations also revealed a wavelength-dependent time shift of the peak brightness which had not been witnessed for other eruptive YSOs before.

A recent study by Contreras Peña et al. (2020), based on MIR (NEOWISE W2, 4.6 μm) and (sub)mm (SCUBA-2, 850 μm) variability of deeply embedded protostars, shed more light on the wavelength dependence of the rise/drop speed of the relative flux m. From our SED fits, we can assess this issue in a “semi-empirical” manner for the whole wavelength range of the SED. The relative rise/drop speeds are obtained for the pre-burst to the burst (Pr →B) and the burst to post-burst (B →Po) epochs, respectively. The results are presented in the lower panel of Fig. 9. The two curves show m(λ) = | log(F_x(λ)∕F_y(λ)) × Δt⁻¹ | for each wavelength of the modeled SEDs. The indices x, y denote the respective epoch pairs and Δt is the corresponding epoch difference (in years). The upper curve holds for the flux rise while the lower one for the flux decrease. Since we do not know exactly when the rise started, the upper curve may be shifted while its slope is unaffected from the actual date.

With only two wavelengths at hand, Contreras Peña et al. (2020) were not able to infer the functional form of m(λ) and had to assume a proportionality such that m(4.6 μm) = η m(850 μm). From objects with the most significant variability at both wavelengths, they derived η = 5.53 ± 0.29. This can be compared tothe η values basedon m at 4.6 and 850 μm from our results (black dots in Fig. 9). These amount to η_{Pr →B} = 4.48 and η_{B →Po} = 5.88 for rise anddrop, respectively. The values from our models are similar, although η_{Pr →B} is somewhat lower.

More importantly, our approach allowed us to derive m(λ) for the whole range of the SED. The corresponding curves show that beyond λ ≳ 10 μm, log (m) depends approximately linear on log(λ). Thus, the wavelength dependence m(λ) actually resembles a power law. The black lines show the respective fits to both epoch pairs. They yielded the following relations for m(λ): m_Pr→B ∝ λ^{−0.42±0.01} and , for the rise and drop, respectively. The spectral indices of both relations agree within the errors, suggesting that the sense of the flux variability does not affect the speed of the relative flux change. These results indicate that (a) burst induced flux variations are quicker at shorter wavelengths and (b) the rise and fall of the fluxes take equal times.

Yet, a cautionary note is required here. First of all, our assessment of the correlated flux variability is based on YSO models with passive disks. The viscous energy release in an active disk leads to a temperature distribution that differs from the passive case in particular in the innermost region. This will lead to a higher value of η as shown by Contreras Peña et al. (2020). Moreover, we have to mention that the models used by these authors and us are based on static RT. Time-dependent RT shall be used to verify the above considerations. In addition, it has to be emphasized that the wavelength dependence of the fading rate as obtained by Contreras Peña et al. (2020) and in the present paper has not been found for the OO Ser event (Kóspál et al. 2007). In that case, the fading rate is only time-dependent and slows down over time. The reason for the different behavior still needs to be understood.

Fig. 9 Top: wavelength dependence of the relative flux increase for the burst (red) and post-burst (green) epochs for mean model. Markers denote grid points. Bottom: wavelength dependence of the speed of the relative flux change m for the pre-burst to burst (blue) and burst to post-burst (green) periods. The black solid lines mark linear fits to λ > 10 μm, the dashed vertical lines indicates the wavelengths used by Contreras Peña et al. (2020). The black dots mark the values for m used to compare our models to Contreras Peña et al. (2020).

9 Discussion

In the following, particular results are evaluated with regard to their credibility and impact concerning our understanding of MYSO accretion bursts. The G358 event is compared to the previous cases. Before we address individual topics, we note for completeness, that our investigation only covers the surplus dust continuum emission caused by the accretion burst. While the fraction of the energy consumed by phase transitions (sublimation of grain ice mantles, dust evaporation and dissociation of molecules) is negligible in the present context, it nevertheless seriously affects the chemistry of the YSO (Rab et al. 2017; Wiebe et al. 2019).

9.1 Misfit of the (sub)mm fluxes

Figure 6 shows that the best burst and post-burst fits predict (sub)mm fluxes that exceed the measured values. The observed values, indicated by the red and green symbols at λ = 889 μm amount to 60–70% of the respective mean burst and post-burst model. Since attempts to solve this discrepancy failed or required non-plausible assumptions, for instance on grain properties or geometry, we conclude that this mismatch indicates a deviation of the temperature distribution of the circumstellar environment from the static case. The accretion burst causes an additional inside-out heating of the circumstellar environment where inner regions try to reach a new thermal equilibrium while outer ones are still in the previous steady-state. Time-dependent RT simulations suggest that the heating timescales strongly increase with wavelength (Harries 2011; Johnstone et al. 2013). Since the MM1 SEDs are dominated by the increased FIR fluxes, their fits will, therefore, yield over-predicted (sub)mm fluxes.

For deeply embedded YSOs, the contribution of the envelope to the thermal FIR/(sub)mm emission may be significant (or even dominant), for example Contreras Peña et al. (2020). For about half of the best MM1 pre-burst models, the overall optical depth of theenvelope exceeds by far that of the disk (see Table A.7). In these cases the envelope produces the overwhelming fraction of the (sub)mm emission. Furthermore, the heating of parts of the envelope, which are in the shadow of the disk, will besuppressed or (strongly) delayed. This effect will be especially important for systems with strongly flared disks.

Not only the heating time scales have an influence on the burst-SED, but also the characteristics of the burst itself. Short/weak bursts might be not strong or long enough to heat the whole circumstellar environment which is probably the case for G358. This will lead to an overestimation of the sub-mm fluxes when using stationary models as well. Because of that, fluxes beyond 889 μm were not included to fit the MM1 burst SED in order to avoid their systematic overestimation induced by our method.

9.2 Uncertainty of the accretion parameters

It was mentioned above that the ZAMS values for R_* might not apply since the high accretion rates during the growth of a massive protostar will temporarily bloat its radius which could then be many times larger than that of a main sequence star (Hosokawa & Omukai 2009; Hosokawa et al. 2010). Recent simulations suggest that bloating is even intensified by accretion bursts (Meyer et al. 2019a). As a consequence of bloating, more mass needs to be accreted to produce the observed luminosity increase. In the cold disk accretion model of Hosokawa et al. (2010), the protostellar swelling occurs in the mass range between 5–9 M_⊙. Since the mass of MM1 is likely of that order, the supposition of being bloated is perhaps valid. If so, the accretion rate may be elevated as it is proportional to the stellar radius. Unfortunately, due to the high extinction of deeply embedded protostars, it is hard to derive their surface gravity from photospheric absorption lines because of strong veiling. Only recently this has been achieved for the first time. The MYSO G015.1288–00.6717 shows an A-type spectrum with a relatively low surface gravity, pointing to a bloated protostar (Pomohaci et al. 2017).

For that reason, caution is required for what concerns the comparison of stellar-related burst properties. For a MYSO in a relatively evolved stage, like S255IR-NIRS3 which shows a flat pre-burst SED, the derivation of an accretion rate using the ZAMS approach seems to be justified. However, its application to MYSOs in earlier stages, like G358 or NGC6334I-MM1, is questionable. In addition, the ZAMS approach has its own challenges because, unlike low-mass stars, massive stars arrive on the ZAMS while accreting. Therefore, their pre-burst luminosity includes a likely small but unknown fraction of accretion luminosity. For low- and intermediate mass-stars, accretion rates can be derived from emission-line tracers (for example Mendigutía et al. 2011) which, however, might be only exceptionally seen in scattered light from MYSOs (Neckel & Staude 1995).

In any case the burst energy E_acc provides a good proxy on the energetics of the event. In this regard, the values for G358-MM1, S255IR-NIRS3, and NGC6334I-MM1 of J (see Sect. 5.2.2), 1.2 ± 0.4 × 10³⁹ J (Caratti o Garatti et al. 2017), and 8 × 10³⁸ J (Hunter et al. 2017) already indicate a range of about an order of magnitude. Since the luminosities of S255IR-NIRS3 and NGC6334I-MM1 were still elevated when the aforementioned papers appeared, their final numbers will be even higher.

9.3 Burst strength bias

It is well established that young massive stars have a high stellar multiplicity (Zinnecker & Yorke 2007; Pomohaci et al. 2019). Given their distances, high spatial resolution obtained by interferometric observations is required to resolve the individual components of MYSOs like G358 or to prove their single nature. For less embedded objects, this can be achieved with the VLTI in the thermal IR, while for deeply embedded sources (sub)mm interferometers like NOEMA or ALMA have to be used. FIR interferometry, which requires groups of satellites in space, does not yet exist, although it has been envisaged (Linz et al. 2020). Keeping in mind that MYSOs emit the bulk of their energy in the FIR, the coarse spatial resolution in the FIR implies that in most if not all cases, their luminosities consist of individual contributions. Thus, relating the luminosity increase caused by an accretion burst of single source to the total pre-burst luminosity only yields a lower limit to the burst strength. In case of G358 the luminosity gain derived from the graybody fits of pre- and burst total fluxes amounts to 2.5 while the RT results based on the SED decomposition (see Sect. 5.1) yielded a value of 4.7 for MM1.

9.4 Burst duration and temporal evolution of the IR emission

Class II methanol masers are excited by MIR thermal emission from warm dust grains with temperatures exceeding 100 K (Ostrovskii & Sobolev 2002; Cragg et al. 2005). Regardless of the details of the enhanced accretion, the increase of released energy will in all likelihood cause the evaporation of dust grains at the inner boundary of the disk and thus shift its rim outward. Since the radiative transfer will adjust the temperature profile accordingly, the MIR-emitting region moves toward a larger radius. Due to its larger surface area the emission increases and that of the masers as well. Thus, the rise and fall of the maser emission is coupled to the accretion variability through the dust MIR continuum emission. First direct evidence for this relation has been found by us for the S255IR-NIRS3 burst (Stecklum et al. 2018a), see Durjasz et al. (2019) and Olech et al. (2020) for other cases.

Since the maser flare of G358 lasted only for a few months (MacLeod et al., in prep.; Yonekura et al., in prep.), it can be concluded that the increased accretion rate dropped shortly before the end of the flare. So the question arises how does it come that elevated FIR emission is still present after 18 months past the flare peak? The current notion is that dust cooling should be quick because of the low optical depth at those wavelengths. The bulk of the FIR/(sub)mm emission comes from outer regions of the circumstellar environment. Thus, their optical depth toward the observer is low indeed. However, this is not necessarily the case for the optical depth toward the protostar where the energy is being released. An effective optical depth of a few will slow down the heating of the outer YSO environment and “trap” burst energy in the dusty medium. The latter will be radiated away and cause elevated FIR emission for quite some time after the burst terminated. Evidence of post-burst elevated FIR emission has been found for the OO Ser event (Kóspál et al. 2007) and by us for S255IR-NIRS3 (Stecklum et al., in prep.). G358 is another example which shows this behavior as well.

The projected linear separations between nearest neighbors of the G358 complex of up to 10⁴ au (Brogan et al. 2019) imply individual object sizes with corresponding light travel times of up to two months. These can easily stretch to one year or more for moderate efficient optical depths. Furthermore, it has to be kept in mind that there is a variety of sight lines with different optical depths along which outer regions will be heated. The different propagation speeds along those lines incur a broadening of the heating duration. A subluminal speed (0.04…0.08 c) of the heat wave during the burst has been witnessed in G358 by tracing, for the first time, the outward motion of maser spots (Burns et al. 2020). Subluminal speed for the relocation of maser spots was found for S255IR-NIRS3 (Stecklum et al. 2018b) as well. This illustrates the imprint of the optical depth on the radiative transfer of the burst energy throughout the YSO.

9.5 Accretion burst drivers

YSO accretion bursts are thought to be caused by enhanced mass transfer through the circumstellar disk, triggered by various reasons, which may drive the inner disk to become active and self-luminous (Hartmann & Kenyon 1996). Among the possible causes, erratic infall from a protostellar envelope (cf. Vorobyov & Basu 2015; Meyer et al. 2017) seems to be favorable for a MYSO in an evolutionary stage as early as G358. Interestingly, recent high-resolution radio imaging of MM1 suggests the presence of two spiral-arm accretion flows, winding toward the protostar as traced by methanol masers (Brogan et al. 2019; Chen et al. 2020b). This is in line with very high-resolution ALMA imaging of high-mass protostars that provided evidence for filamentary streamers pointing onto the central sources (Goddi et al. 2020; Johnston et al. 2020). These might represent multi-directional accretion channels which possibly inhibit the formation of a large, steady disk at the very early stages of massive star formation. Thus, the small disk mass suggested by our best RT SED fits may be seen as a hint for the extreme youth of MM1. Possibly, its environment is not relaxed enough to allow the formation of a larger sustainable circumstellar disk. Such a disk might be prone to gravitational instability, considered to be another cause for episodic accretion (Vorobyov & Basu 2015). Faint radio continuum emission from MM1 has been detected during the burst as well (Bayandina et al., in prep.). Whether it is variable and perhaps related to an active disk needs to be explored.

9.6 G358 in the accretion burst context

The G358 burst is the second one from a MYSO for which the accretion luminosity could be measured from the SED changes. This allows us to compare the derived quantities with the results obtained for S255IR-NIRS3 (Caratti o Garatti et al. 2017), which is thefirst step toward a statistics of MYSO burst properties. Assuming a ZAMS star, an accreted mass of ~ 3.4 × 10⁻³ M_⊙ has been estimated for S255IR-NIRS3. The corresponding accretion rate amounts to 5 ± 2 × 10⁻³ M_⊙ yr⁻¹. While these values will likely be revised once the comprehensive monitoring data set has been analyzed, we can assume that their order of magnitude is correct. The accretion rate of the G358 burst almost corresponds to that of the S255IR-NIRS3 event. However, the accreted mass is about one order of magnitude smaller due to the shorter duration of the burst. Thus, compared to S255IR-NIRS3, the G358 burst was a minor one. Simulations of MYSO accretion (Meyer et al. 2019b) suggest that minor bursts similar to that of G358 are much more frequent compared to major ones, and occur predominantly at very early stages of protostellar evolution. We emphasize that the above-mentioned disk accretion rates during the burst are at least an order of magnitude larger than those derived from SED modeling of non-bursting high-mass protostellar objects (Fazal et al. 2008; Grave & Kumar 2009).

While the MYSO accretion burst sample is still small, it will grow in the mid-term by following up methanol maser alerts on a regular basis. Recently, the suggestion by Proven-Adzri et al. (2019) that the periodicity of methanol masers in G323.46–0.08 could be due to an accretion burst was confirmed (Wolf et al., in prep.) using (NEO)WISE and VVV data. While this is another event which was identified a posteriori, it raised the total number of known MYSO bursts to five, including V723 Car. So a thorough comparison of MYSO burst properties with corresponding models (for example Meyer et al. 2019b), should become possible in the foreseeable future.

10 Conclusions

The SOFIA FIFI-LS observations of G358 during two epochs, supplemented by NIR imaging and supplementary archival data, and their RT analysis yield the following major results:

• Increased FIR fluxes measured with FIFI-LS confirmed the accretion burst of G358 which was alerted by the maser flare. Since no excess emission associated with the burst was detected in NIR, MIR or (sub)mm observations, it represents the first NIR-, FIR- and (sub)mm-dark but FIR-loud MYSO accretion burst.

• By means of the SED decomposition approach actual luminosity estimates of the driving source of the accretion burst could be obtained which yielded more representative burst parameters.

• From the change of the luminosities at the two FIFI-LS epochs the decay-time time of the FIR excess emission could be determined. This yielded more reliable estimates of the burst energy and the derived parameters. Moreover, wavelength-dependent rates for the rise and fall of the relative fluxes could be obtained.

• The circumstellar disk of MM1 is less massive than those found for other MYSOs and possibly transient. This may be an indication of the extreme youth of the source.

• A comparison with previous MYSO bursts indicates a considerable range of burst characteristics. The burst of G358 is the least energetic of the small sample.

• The RT modeling allowed us to draw conclusions on the possible sites of maser excitation and their relocation due to the burst which seem to be supported by VLBI radio observations.

• As already demonstrated for the S255IR-NIRS3 event, the G358 results underline once more that the observational capabilities of SOFIA in the FIR are of utmost importance to study the erratic growth of stars.

As soon as a maser parallax for G358 becomes available, a reanalysis should be performed to narrow down the uncertainties of the derived parameters which will increase their credibility.

Acknowledgements

We thank the anonymous referee for accurate review and highly appreciated comments and suggestions which helped to improve the quality of this paper. V.W. is supported by the by the German Aerospace Center (DLR) under grant number 50OR1718. A.C.G. has received funding from the European Research Council (ERC) under the European Union’s Horizon2020 research and innovation programme (grant agreement No. 743029). A.S. acknowledges support by the Ministry of Science and Higher Education of the Russian Federation under the grant 075-15-2020-780 (contract 780-10). A.V. acknowledges the support from the Government of Russian Federation and Ministry of Science and Higher Education of the Russian Federation (grant N13.1902.21.0039) in the part of the data analysis. Part of the funding for GROND (both hardware and personnel) was generously granted by the Leibniz-Prize to G. Hasinger (DFG grant HA 1850/28-1) and by the Thüringer Landessternwarte Tautenburg. Based on observations made with the NASA/DLR Stratospheric Observatoryfor Infrared Astronomy (SOFIA). SOFIA is jointly operated by the Universities Space Research Association,Inc. (USRA), under NASA contract NNA17BF53C, and the Deutsches SOFIA Institut (DSI) under DLR contract 50 OK 0901 to the University of Stuttgart. This research has made use of the NASA/IPAC Infrared Science Archive, which is funded by the National Aeronautics and Space Administration and operated by the California Institute of Technology. This publication makes use of data products from the Near-Earth Object Wide-field Infrared Survey Explorer ((NEO)WISE), which is a joint project of the Jet Propulsion Laboratory/California Institute of Technology and the University of Arizona. (NEO)WISE is funded by the National Aeronautics and Space Administration. The ATLASGAL project is a collaboration between the Max-Planck-Gesellschaft, the European Southern Observatory (ESO) and the Universidad de Chile. It includes projects E-181.C-0885, E-78.F-9040(A), M-079.C-9501(A), M-081.C-9501(A) plus Chilean data. This research has made use of NASA’s Astrophysics Data System.

Appendix A Additional information

Fig. A.1 Model pool (in gray) showing all SEDs, used to fit burst and post-burst epochs. The observed fluxes from both epochs are plotted in red (burst) and green (post-burst). They are well covered by the model SEDs (except for the burst data points at 163 and 186 μm, as discussed when introducing the model pool in Sect. 5.2.2). The dashed black line indicates the mean model, as obtained from the pre-burst fit alone. The NIR/MIR fluxes of that model have been used to set a lower limit to the post-burst, which was only observed in the red channel of FIFI (at λ ≥ 118 μm).

Fig. A.2 MM3 RT models illustrating constraints on the FIR emission. Filled circles indicate measured fluxes. The model shown in green is the best one (cf. Fig. 5). The best fit with the (sub)mm fluxes (green) omitted is shown in red. When also omitting the WISE W4 and MIPS 24 μm measurements (yellow), the best fit (blue) lies above the nominal one. This confirms that the FIR emission of MM3 and, thus, its contribution to the total FIR fluxes, is well constrained by the WISE, MIPS and ALMA measurements.

References

All Tables

All Figures

	Fig. 1 Top: VVV JHKs color composite of the G358 region (epoch 2010 August 15). The positions of MM1 and MM3 from Brogan et al. (2019) are marked. Black plus signs denote positions from IR observations at wavelengths up to 24 μm. Bottom: GROND JHKs color composite (epoch 2019 February 8) with contours of the ALMA 0.89 mm continuum map (Brogan et al. 2019).
In the text

	Fig. 2 (NEO)WISE W1 (blue) and W2 (red) light curves based on mean magnitudes and respective errors for each visit. The first two epochs are from the WISE mission. Vertical lines mark the dates of the flare peak (black) and the first FIFI-LS epoch (green). Horizontal lines indicate mean magnitudes (solid) and their errors (dashed).
In the text

Fig. 3 (NEO)WISE coordinate offsets relative to an origin at MM1 (from Brogan et al. 2019). Pre-burst, burst, and post-burst data are shown in blue, red and green. Thick error bars represent mean positions while individual appear thin. Positions of MM1 and MM3 (measured with ALMA) are marked by the diamond and the square, respectively. The arrow points toward the bright 2MASS J17431001−2951460, which is at a distance of 9.′′9 from the (0,0) position. The yellow circle indicates the W1 image FWHM of 6.′′ 1.

In the text

	Fig. 4 Observed FIR/(sub)mm SEDs showing pre-burst Herschel/ AKARI data (blue) and FIFI-LS burst values (red) together with the corresponding (sub)mm fluxes from Brogan et al. (2019). The upper 24 μm MIPS limit is indicated as well. The solid lines represent the reddened graybody fits (cf. Sect. 4.1). The second epoch FIFI-LS data and the corresponding fit appear bf in green.
In the text

Fig. 5 Pre-burst SEDs for the following sources: black symbols – total FIR/(sub)mm emission of G358, green – MM3, blue – MM1, derived from the total FIR/(sub)mm emission by removing contributions from all other sources (including MM3). For MM1 and MM3, the observed values are shown together with the ten best RT fits. At wavelengths beyond 40 μm MM1 dominates the total flux density. The seemingly gap at 10 μm in the MM1 SED is due to strong silicate absorption.

In the text

Fig. 6 Modeled pre-burst MM1 SEDs (blue), together with its burst (red) and post-burst SED (green). Triangles mark upper limits. The weighted mean-model (see text) is shown in black. It is shown three times but with different source luminosities for the respective observing epochs. Since the 870 (pre-) and 889 μm (post- and burst) fluxes are very similar, their symbols cannot easily be distinguished.

In the text

Fig. 7 Dependence of the derived accretion rate on the average luminosity increase ⟨ ΔLacc ⟩ and assumed ZAMS stellar mass for the parameter range bracketing MM1. The horizontal dashed line corresponds to ⟨ ΔLacc ⟩ and the lightgray region marks its 1σ uncertainty while the vertical dashed line and the darker gray region indicate the stellar mass and its uncertainty. The values of the accretion rate are indicated.

In the text

Fig. 8 Temperature distribution of the mean model in the x-z plane (innermost part of first quadrant) for the pre-burst (left), burst- (center) and post-burst epochs (right). The orange and red lines enclose the temperature range from 94 K (orange) to 120 K (red) during the pre- and burst epochs, respectively. The white contours mark gas particle volume densities of which decrease with increasing z. The vertical solid black line indicates the outer radius of the disk. The dashed black lines mark the radius of the maser ring from the first and 4th epoch of the VLBI observations (Burns et al. 2020, and in prep.). The length of the black bar corresponds to 50 mas.

In the text

Fig. 9 Top: wavelength dependence of the relative flux increase for the burst (red) and post-burst (green) epochs for mean model. Markers denote grid points. Bottom: wavelength dependence of the speed of the relative flux change m for the pre-burst to burst (blue) and burst to post-burst (green) periods. The black solid lines mark linear fits to λ > 10 μm, the dashed vertical lines indicates the wavelengths used by Contreras Peña et al. (2020). The black dots mark the values for m used to compare our models to Contreras Peña et al. (2020).

In the text

Fig. A.1 Model pool (in gray) showing all SEDs, used to fit burst and post-burst epochs. The observed fluxes from both epochs are plotted in red (burst) and green (post-burst). They are well covered by the model SEDs (except for the burst data points at 163 and 186 μm, as discussed when introducing the model pool in Sect. 5.2.2). The dashed black line indicates the mean model, as obtained from the pre-burst fit alone. The NIR/MIR fluxes of that model have been used to set a lower limit to the post-burst, which was only observed in the red channel of FIFI (at λ ≥ 118 μm).

In the text

Fig. A.2 MM3 RT models illustrating constraints on the FIR emission. Filled circles indicate measured fluxes. The model shown in green is the best one (cf. Fig. 5). The best fit with the (sub)mm fluxes (green) omitted is shown in red. When also omitting the WISE W4 and MIPS 24 μm measurements (yellow), the best fit (blue) lies above the nominal one. This confirms that the FIR emission of MM3 and, thus, its contribution to the total FIR fluxes, is well constrained by the WISE, MIPS and ALMA measurements.

In the text