key: cord-0490267-hrhlomuy
authors: Stecklum, B.; Wolf, V.; Linz, H.; Garatti, A. Caratti o; Schmidl, S.; Klose, S.; Eisloffel, J.; Fischer, Ch.; Brogan, C.; Burns, R.; Bayandina, O.; Cyganowski, C.; Gurwell, M.; Hunter, T.; Hirano, N.; Kim, K.-T.; MacLeod, G.; Menten, K. M.; Olech, M.; Orosz, G.; Sobolev, A.; Sridharan, T. K.; Surcis, G.; Sugiyama, K.; Walt, J. van der; Volvach, A.; Yonekura, Y.
title: Infrared observations of the flaring maser source G358.93-0.03 -- SOFIA confirms an accretion burst from a massive young stellar object
date: 2021-01-05
journal: nan
DOI: nan
sha: 3d4230bbecf328aebb93e096b07abf9b69a1f4f5
doc_id: 490267
cord_uid: hrhlomuy

Class II methanol masers are signs of massive young stellar objects (MYSOs). Recent findings show that MYSO accretion bursts cause flares of these masers. Thus, maser monitoring can be used to identify such bursts. Burst-induced SED changes provide valuable information on a very intense phase of high-mass star formation. In mid-January 2019, a maser flare of the MYSO G358.93-0.03 was reported. ALMA and SMA imaging resolved the core of the star forming region and proved the association of the masers with the brightest continuum source MM1. However, no significant flux rise of the (sub)mm dust continuum was found. Thus, we performed NIR imaging with GROND and IFU spectroscopy with FIFI-LS aboard SOFIA to detect possible counterparts to the (sub)mm sources, and compare their photometry to archival measurements. The comparison of pre-burst and burst SEDs is of crucial importance to judge whether a luminosity increase due to the burst is present and if it triggered the maser flare. The FIR fluxes of MM1 measured with FIFI-LS exceed those from Herschel significantly, which clearly confirms the presence of an accretion burst. The second epoch data, taken about 16 months later, still show increased fluxes. Our RT modeling yielded major burst parameters and suggests that the MYSO features a circumstellar disk which might be transient. From the multi-epoch SEDs, conclusions on heating and cooling time-scales could be drawn. Circumstances of the burst-induced maser relocation have been explored. The verification of the accretion burst from G358 is another confirmation that Class II methanol maser flares represent an alert for such events. The few events known to date already indicate that there is a broad range in burst strength and duration as well as environmental characteristics. The G358 event is the shortest and least luminous MYSO accretion burst so far.

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 Based on observations collected at the European Organization for Astronomical Research in the Southern Hemisphere under ESO program 0103.C-9033(A). Based on observations made with the NASA/DLR Stratospheric Observatory for Infrared Astronomy (SOFIA) under Proposal IDs 75_0037 and 08_0163. 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 3 L (Povich et al. 2011) to ∼ 4 × 10 3 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 Article number, page 1 of 22 arXiv:2101.01812v1 [astro-ph.SR] 5 Jan 2021 A&A proofs: manuscript no. G358-SOFIA 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 (Hunter et al. , 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. 2018b; 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. 2018a; 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) 1 to coordinate single-dish monitoring of masers and interferometric follow-up measurements.

G358.93-0.03 (hereafter G358, RA: 17 h 43 m 10. s 02, 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 + − 0.37 0.68 kpc 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 3 OH maser line (Menten 1991b) in G358 was announced . 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 . 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 , 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. 1 See M2O website at http://MaserMonitoring.org 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 was made 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.

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 Tab. 1 The comparison of the HI-GAL photometry (Molinari et al. 2016) with that performed by the ATLASGAL team ) 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) 2 , using the MP-Fit2DPeak 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. 

Optical and near-infrared (NIR) imaging of G358 was performed using the seven-channel Gamma-Ray Burst Optical/Nearinfrared Detector GROND , 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 minutes. Data processing was performed by means of the GROND pipeline (Krühler et al. 2008 ).

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 µm 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 observing time 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.

As evidenced by the ALMA/SMA observations, G358 is a star forming complex harboring several massive protostars . 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.

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 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. 

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 3 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 (c.f. section 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, upper bounds for a possible flux contribution caused by the burst can be derived as 0.10 Jy and 0.45 Jy, respectively. The failure of the burst detection suggests that, at any time, MM3 provided 3 See http://wise2.ipac.caltech.edu/docs/release/ neowise/expsup/sec2_1civa.html 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.

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.

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.

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 4 .

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 

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 MIPSGAL image 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.

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 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 ) 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 preburst and depressed compared to the bursting stages.

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.

We checked the Herschel-PACS images at 70 µm for potential signs 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 ( . 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.

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 Sec. 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 T pre d = 26.4 ± 0.5 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 T burst d = 30.1 ± 0.5 K as well as the solid angle of the emission which increased by a factor of 1.20 ± 0.14. However, the χ 2 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 preand burst fits were adopted and only the temperature was varied. The resulting temperature amounts to T post d = 28.9 ± 0.6 K.

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 Sec. 1 as well as in-terstellar 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 Sec. 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 L 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 a pre-burst-level return date of 2021 July 31.

The FIR/(sub)mm burst energy is E acc FIR = <∆L acc FIR > × ∆t, where <∆L acc FIR > is the average luminosity increase which equals to half of the peak increase for a linear drop to zero. It amounts to E acc FIR = 1.9 × 10 38 J. Its upper and lower bounds from the uncertainties in both duration and luminosity are 1.0 × 10 38 J and 2.6 × 10 38 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 Sec. 6.

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 Sec. 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 2 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 temper- 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. ature 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 postburst observations, the pre-burst fluxes were predicted from the MM3 model as 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 Table A .2, A.3 and 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 Tab. 2 of Brogan et al. (2019) by a factor of about two.

Characteristic properties of YSOs, namely their luminosities as well as mass, geometry and extent of the surrounding dust, can be derived 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" 6 of the SED fitter (Robitaille et al. 2007 ) which is described below. We did not fit the G358 total flux (black sym-6 https://zenodo.org/record/235786 bols 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 χ 2 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 , M dust disk , L 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.

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 Tab. 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 Sec. 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 M dust disk = 0.068 + − 0.11 0.043 M . 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 2.6 + − 6.8

1.0 M . 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 510 + − 450 240 au, while the maximum is as high as 1600 au. The mean luminosity of 7540 + − 5340 3130 L 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 −1 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 χ 2 values) can be found (together with the weighted means and standard deviation σ) in Table A.5.

As described in Sec. 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.

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 preburst 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 Sec. 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 from the 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 χ 2 -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 preburst, 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 L pre = 5000 + − A comparison of the RT results to those obtained from graybody fits in Sec. 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 Sec. 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 geom-etry agrees with that of the spiral-arm accretion flow of Chen et al. (2020b) with i = 25 ± 10 • . The stellar radius amounts to 8.4 + − 15.7 5.5 R . 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 950 + − 1500 580 au. The favored dust mass of the circumstellar disk is as low as 8.4 × 10 −5 M , where the 1σ-confidence interval extends from 1.1 × 10 −6 to 6.1 × 10 −3 M . Using the appropriate gas-to-dust mass ratio (see Sec. 5.2.1) this corresponds to a most probable total disk mass of 3.2 × 10 −3 M within an uncertainty range of 4.4 × 10 −5 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 χ 2 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 χ 2 -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 χ 2 -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 preburst-system most probably is stationary, this does not hold for the post-/bursting stage). We note, that in the burst case the relative weights of the models are on the same order, whereas in the pre-burst case they differ by a factor of 1.5 at most (for the post-burst it is 1.3).

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 mass and radius estimates which reproduce the MM1 pre-burst luminosity, (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 froṁ M 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 Sec. 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 M acc = 3.1 + −

where the errors are dominated by the error of ∆t acc and E acc . Using the ZAMS massradius 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 M acc = 5. 3.0 × 10 −3 M yr −1 . 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.

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 grain heating 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 Sec. 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 n H 2 = [1, 0.5, 0.3, 0.2] × 10 8 cm −3 (density decreases with z). To transform the dust densities from the RT-simulation to the above given gas densities, we use the relation ρ = m H 2 N −1 A × n H 2 γ −1 , where m H 2 is the molar mass of H 2 , N A is the Avogadro constant and γ is the dust to gas ratio introduced in Sec. 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 8 cm −3 , 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; Burns et al., 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 preburst and burst.

Another important parameter for the maser excitation is the specific column density, N M ∆V −1 , 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 slightly larger 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.

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 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 "semiempirical" 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 −1 | 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 to the η values based on 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 and drop, 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 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. Moreover, 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.

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 ). 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 the envelope 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 be suppressed 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.

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 2.9 + − (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.

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 Sec. 5.1) yielded a value of 4.7 for MM1.

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 4 au 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.

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 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. 2018; Johnston et al. 2020) . These might represent multidirectional 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.

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 the first step toward a statistics of MYSO burst properties. Assuming a ZAMS star, an accreted mass of ∼ 3.4 × 10 −3 M has been estimated for S255IR-NIRS3. The corresponding accretion rate amounts to 5 ± 2 × 10 −3 M yr −1 . 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.

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 G358is 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.

A&A proofs: manuscript no. G358-SOFIA Table A .1. MM3-SED. The fluxes at λ ≥ 889 µm, are MM3-fluxes, whereas the fluxes at λ ≤ 24 µm are total fluxes. The contribution of all the other sources, which are less evolved than MM3 at NIR wavelengths, can be neglected (as indicated in Fig. 5 ). For the fit of MM3 (see Sec. 5.2.1), we use the total fluxes at short wavelengths. Only in the FIR and (sub)mm regimes, where the contribution of the other sources is not negligible, we use the MM3-fluxes. Facilities and instruments are given behind corresponding references.

Flux Ref. [µm] [erg cm −2 s −1 ]

1.63 1.16 ± 0.07 × 10 −12 1 2. 13 1.67 ± 0.02 × 10 −11 1 3.55

3.64 ± 0.02 × 10 −10 2 4.49

6.10 ± 0.03 × 10 −10 2 5.73

9.73 ± 0.05 × 10 −10 2 7.0 6.90 ± 0.49 × 10 −10 3 7.87 5.30 ± 0.04 × 10 −10 2 11. 6 3.54 ± 0.07 × 10 −10 4 15.0

6.46 ± 0.20 × 10 −10 3 22. 1 5.13 ± 0.13 × 10 −10 4 23.7

4.49 ± 0.13 × 10 −10 5 889

1.34 ± 0.04 × 10 −13 6 1282

2.43 ± 0.10 × 10 −14 6 1420

1.50 ± 0.19 × 10 −14 6 1532

1.06 ± 0.08 × 10 −14 6 2.93 ± 0.30 × 10 −10 1 52.0 5.51 ± 0.56 × 10 −9 1 4.85 ± 0.49 × 10 −9 54. 8 6.94 ± 0.70 × 10 −9 1 6.12 ± 0.62 × 10 −9 60.7

6.20 ± 0.62 × 10 −9 1 5.04 ± 0.51 × 10 −9 87.2 7.48 ± 0.75 × 10 −9 1 5.35 ± 0.54 × 10 −9 118. 6 6.45 ± 0.65 × 10 −9 1 4.59 ± 0.46 × 10 −9 124.2 8.97 ± 0.90 × 10 −9 1 7.22 ± 0.72 × 10 −9 142.2 5.26 ± 0.53 × 10 −9 1 5.25 ± 0.53 × 10 −9 153. 3 5.29 ± 0.53 × 10 −9 1 4.08 ± 0.41 × 10 −9 162. 8 5.47 ± 0.55 × 10 −9 1 4.41 ± 0.45 × 10 −9 186. 4 4.58 ± 0.46 × 10 −9 1 3.82 ± 0.39 × 10 −9 889

3.81 ± 0.11 × 10 −12 2 1.72 ± 0.07 × 10 −12 1282

4.21 ± 0.03 × 10 −13 2 1420

2.74 ± 0.28 × 10 −13 2 1532

1.88 ± 0.01 × 10 −13 2

References. (1) 118.6 5.17 ± 0.52 × 10 −9 1 3.55 ± 0.36 × 10 −9 124.2 7.70 ± 0.77 × 10 −9 1 5.55 ± 0.56 × 10 −9 142.2 4.74 ± 0.48 × 10 −9 1 3.28 ± 0.33 × 10 −9 153. 3 4.35 ± 0.44 × 10 −9 1 3.09 ± 0.31 × 10 −9 162. 8 3.01 ± 0.31 × 10 −9 1 1.90 ± 0.19 × 10 −9 186.4

2.53 ± 0.26 × 10 −9 1 1.73 ± 0.18 × 10 −9 889

3.81 ± 0.11 × 10 −12 2 1.72 ± 0.07 × 10 −12 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.

Article number, page 19 of 22

A&A proofs: manuscript no. G358-SOFIA (2001) Milky Way grain size distribution A for R V = 5.5. We note that, some models appear more than ones, but for different inclinations. In our analysis we thread them as different models. The envelope inner radius is set to r min disk , its centrifugal radius r C is set to r max disk . All models extends to an outer radius,

where density and temperature have dropped to the ambient level. The table features the same structure as Table A .5, but additionally the burst-fit is included. The burst names are composed of the original model-name and the adapted luminosity increase. Consequently the weighted mean and σ are obtained taking into account both, the pre-and burst-results (see text). Only the parameters, which refers to the source are taken from the pre-burst alone (they change due to the burst). Again, the confidence interval for log-sampled values extends from x σ to x × σ.

We note that the inner radius of the disk is not a free parameter (unlike for MM3). Instead, it is governed by the sublimation radius (at T sub = 1600K). An outward shift of the sublimation radius may lower the dust masses for the post-/burst. 

Protostars and Planets VI

Astronomical Data Analysis Software an Systems XXIV (ADASS XXIV)

VizieR Online Data Catalog

VizieR Online Data Catalog

Explanatory Supplement to the WISE All-Sky Data Release Products

Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series

American Astronomical Society Meeting Abstracts

The Astronomer's Telegram

Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series

Vistas in Astronomy

Infrared observations of the flaring maser source G358

Astronomical Society of the Pacific Conference Series

Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series

Atoms, Ions and Molecules: New Results in Spectral Line Astrophysics

VizieR Online Data Catalog

VizieR Online Data Catalog

Cosmic Masers: From Proto-Stars to Black Holes

VizieR Online Data Catalog

The wonders of star formation, online proceedings

The Astronomer's Telegram

Astrophysical Masers: Unlocking the Mysteries of the Universe

The Astronomer's Telegram

VizieR Online Data Catalog

Sternwarte 5, 07778 Tautenburg

Dublin Institute for Advanced Studies, 31 Fitzwilliam Place, D02 XF86

Germany 5 National Radio Astronomy Observatory, 520 Edgemont Road, Charlottesville

Lebedev Physical Institute of RAS, 84/32 Profsoyuznaya st

South Africa ; The University of Western Ontario

Germany 15 Institute of Astronomy, Faculty of Physics

Institute for Natural Sciences and Mathematics

Via della Scienza 5, 09047

South Africa 20 Radio Astronomy Laboratory of Crimean Astrophysical Observatory RAS, Katsively, RT-22 Crimea 21 Center for Astronomy

ISO/ISOCAM)

APEX/LABOCA) Notes. (*) Upper limit References. (1) present paper (SOFIA/FIFI-LS)

Acknowledgements. We thank the anonymous referee for accurate review and highly appreciated comments and suggestions which helped to improve the quality of this paper. VW is supported by the by the German Aerospace Center (DLR) under grant number 50OR1718. ACG has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme ( 

3.10 ± 0.06 × 10 −9 1 1.55 ± 0.03 × 10 −9 703.13 ± 0.38 × 10 −9 1 1.57 ± 0.06 × 10 −9 1601.90 ± 0.21 × 10 −9 1 9.50 ± 0.52 × 10 −10 2507.82 ± 1.04 × 10 −10 1 3.91 ± 0.27 × 10 −10 3502.34 ± 0.52 × 10 −10 1 1.17 ± 0.26 × 10 −10 5004.98 ± 2.82 × 10 −11 1 2.49 ± 0.33 × 10 −11 8504.73 ± 0.15 × 10 −12 2 2.36 ± 0.08 × 10 −12 8704.03 ± 0.07 × 10 −12 3 2.00 ± 0.35 × 10 −12 A&A proofs: manuscript no. G358-SOFIA Table A .7. Optical depth in the V-band for the ten best MM1 pre-burst models along the mid-plane and along the line of sight toward the center. The intrinsic extinction dominates the total extinction (including the interstellar extinction) by far, as indicated by the huge values in the Vband. For τ mid−plane the contribution of the disk is given separately. For some models, τ disk almost covers the total extinction (along the line of sight), whereas for some models the contribution of the envelope (which is the same as τ mid−plane total − τ mid−plane disk ) is the dominant one. We note that the optical depths are given for the pre-burst, during the burst the optical depths may be slightly lower because of dust sublimation.