11institutetext:  Thüringer Landessternwarte Tautenburg, Sternwarte 5, 07778
Tautenburg, Germany, 11email<EMAIL_ADDRESS>22institutetext: Max
Planck Institute for Astronomy, Königstuhl 17, 69117 Heidelberg, Germany
33institutetext: Dublin Institute for Advanced Studies, 31 Fitzwilliam Place,
D02 XF86, Dublin, Ireland 44institutetext: Deutsches SOFIA Institut,
University of Stuttgart, 70569 Stuttgart, Germany 55institutetext: National
Radio Astronomy Observatory, 520 Edgemont Road, Charlottesville, VA 22903, USA
66institutetext: Mizusawa VLBI Observatory, National Astronomical Observatory
of Japan, Osawa 2-21-1, Mitaka, Tokyo 181-8588, Japan 77institutetext: Joint
Institute for VLBI ERIC, Oude Hoogeveensedijk 4, 7991 PD Dwingeloo, The
Netherlands 88institutetext: Astro Space Center, P.N. Lebedev Physical
Institute of RAS, 84/32 Profsoyuznaya st., Moscow, 117997, Russia
99institutetext: SUPA, School of Physics and Astronomy, University of St.
Andrews, North Haugh, St. Andrews KY16 9SS, UK 1010institutetext: Center for
Astrophysics — Harvard Smithsonian, Cambridge, MA 02138, USA
1111institutetext: Academia Sinica Inst of Astronomy & Astrophysics ASIAA PO
Box 23-141. Taipei 106. China 1212institutetext: Korea Astronomy and Space
Science Institute, 776 Daedeokdae-ro, Yuseong-gu, Daejeon, 34055, Republic of
Korea 1313institutetext: Hartebeesthoek Radio Astronomy Observatory, P.O. Box
443, Krugersdorp 1740, South Africa ; The University of Western Ontario, 1151
Richmond Street, London, ON N6A 3K7, Canada 1414institutetext: Max-Plank-
Institut für Radioastronomie, Auf dem Hügel 69, 53121 Bonn, Germany
1515institutetext: Institute of Astronomy, Faculty of Physics, Astronomy and
Informatics, Nicolaus Copernicus University, Grudziadzka 5, 87-100 Torun,
Poland 1616institutetext: Xinjiang Astronomical Observatory, Chinese Academy
of Sciences, Urumqi, Xinjiang, China 1717institutetext: Astronomical
Observatory, Institute for Natural Sciences and Mathematics, Ural Federal
University, 19 Mira street, Ekaterinburg 620002, Russia) 1818institutetext:
INAF-Osservatorio Astronomico di Cagliari, Via della Scienza 5, 09047,
Selargius, CA, Italy 1919institutetext: Space Research Unit, Physics
Department, North West University, Potchefstroom 2520, South Africa
2020institutetext: Radio Astronomy Laboratory of Crimean Astrophysical
Observatory RAS, Katsively, RT-22 Crimea 2121institutetext: Center for
Astronomy, Ibaraki University, 2-1-1 Bunkyo, Mito, Ibaraki 310-8512, Japan

# Infrared observations of the flaring maser source G358.93$-$0.03††thanks:
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.

SOFIA confirms an accretion burst from a massive young stellar object

B. Stecklum 11 V. Wolf 11 H. Linz 22 A. Caratti o Garatti 33 S. Schmidl 11 S.
Klose 11 J. Eislöffel 11 Ch. Fischer 44 C. Brogan 55 R. Burns 66 O. Bayandina
77 C. Cyganowski 99 M. Gurwell 1010 T. Hunter 44 N. Hirano 1111 K.-T. Kim 1212
G. MacLeod 1313 K. M. Menten 1414 M. Olech 1515 G. Orosz 1616 A. Sobolev 1717
T. K. Sridharan 1010 G. Surcis 1818 K. Sugiyama 66 J. van der Walt 1919 A.
Volvach 2020 Y. Yonekura 2121

(Received October 11, 2020; accepted December 21, 2020)

###### Abstract

Context.  Class II methanol masers are signposts of massive young stellar
objects (MYSOs). Recent evidence shows that flares of these masers are driven
by MYSO accretion bursts. Thus, maser monitoring can be used to identify such
bursts which are hard to discover otherwise. Infrared observations reveal
burst-induced changes in the spectral energy distribution - first and foremost
a luminosity increase - which provide valuable information on a very intense
phase of high-mass star formation.

Aims.  In mid-January 2019, flaring of the 6.7 GHz CH3OH maser (hereafter
maser) of the MYSO G358.93-0.03 (hereafter G358) was reported. The
international maser community (M2O) initiated an extensive observational
campaign which revealed extraordinary maser activity and yielded the detection
of numerous new masering transitions. Interferometric imaging with the Atacama
Large Millimeter/submillimeter Array (ALMA) and the Submillimeter Array (SMA)
resolved the maser emitting core of the star forming region and proved the
association of the masers with the brightest continuum source (MM1), which
hosts a hot molecular core. These observations, however, failed to detect a
significant rise in the (sub)millimeter dust continuum emission. Therefore, we
performed near-infrared (NIR) and far-infrared (FIR) observations to prove or
disprove whether the CH3OH flare was driven by an accretion burst.

Methods.  NIR imaging with the Gamma-Ray Burst Optical/Near-infrared Detector
(GROND) has been acquired and integral-field spectroscopy with the Field-
Imaging Far-Infrared Line Spectrometer (FIFI-LS) aboard the Stratospheric
Observatory for Infrared Astronomy (SOFIA) was carried out on two occasions to
detect possible counterparts to the (sub)millimeter sources and compare their
photometry to archival measurements. The comparison of pre-burst and burst
spectral energy distributions is of crucial importance to judge whether a
substantial luminosity increase, caused by an accretion burst, is present and
if it triggered the maser flare. Radiative transfer modeling of the spectral
energy distribution (SED) of the dust continuum emission at multiple epochs
provides valuable information on the bursting MYSO.

Results.  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 radiative transfer modeling yielded major burst
parameters and suggests that the MYSO features a circumstellar disk which
might be transient. From the pre-burst, burst, and post-burst SEDs,
conclusions on heating and cooling time-scales could be drawn. Circumstances
of the burst-induced maser relocation have been explored.

Conclusions.  The verification of the accretion burst from G358 is another
confirmation that Class II methanol maser flares represent an alert for such
events. Thus, monitoring of these masers greatly enhances the chances of
identifying MYSOs during periods of intense growth. 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 accretion burst known to date. According to
models, bursts of this kind occur most often.

###### Key Words.:

Accretion, accretion disks – Stars: formation – Stars: protostars – Stars:
individual objects (G358.93-0.03) – Radiative transfer

##  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_{\star}{\gtrsim}8\,{\rm M}_{\odot}$) 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\,{\times}\,10^{3}\,{\rm L}_{\odot}$ (Povich et al., 2011) to
${\sim}\,4\,{\times}\,10^{3}\,{\rm L}_{\odot}$ (Tapia et al., 2015) which
correspond to a mass of 8 – 9 ${\rm M}_{\odot}$ (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., 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)111See M2O website at http://MaserMonitoring.org to coordinate single-
dish monitoring of masers and interferometric follow-up measurements.

G358.93-0.03 (hereafter G358, RA: $\rm 17^{h}43^{m}10.\\!^{s}02{}$,
$\rm{DEC\\!:\leavevmode\nobreak\ }{-}29\degr
51^{\prime\prime}45\aas@@fstack{\prime\prime}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$\,{\mathbin{\smash{\raisebox{2.79857pt}{$\underset{\raisebox{0.21529pt}{$\smash{-}$}}{\smash{+}}$}}}}\,^{\,0.37}_{\,0.68}$
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
$\lesssim$ 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 CH3OH 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 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$\sigma$ 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 Tab.1

Table 1: Archival resources Survey |  Spectral range |  Facility |  Reference
---|---|---|---
2MASS |  NIR |  1.3-m telescopes at Fred Lawrence Whipple Observatory and Cerro Tololo Inter-American Observatory |  Cutri et al. (2003)
VVV |  NIR |  4-m VISTA telescope at ESO Cerro Paranal |  Minniti et al. (2010)
ISOGAL |  MIR |  0.6-m ISO space telescope, ISOCAM camera |  Omont et al. (2003)
GLIMPSE |  MIR |  0.85-m Spitzer space telescope, IRAC camera |  Spitzer Science (2009)
(NEO)WISE |  MIR |  0.4-m WISE space telescope |  Mainzer et al. (2014)
MIPSGAL |  MIR and FIR |  0.85-m Spitzer space telescope, MIPS camera |  Gutermuth & Heyer (2015)
FIS |  FIR |  0.67-m AKARI space telescope, FIS scanning photometer |  Doi et al. (2015)
HI-GAL |  FIR |  3.5-m Herschel space telescope, PACS and SPIRE cameras |  Molinari et al. (2016)
ATLASGAL |  (sub)mm |  12-m APEX telescope at Mount Chajnantor |  Schuller et al. (2009)

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)222https://irsa.ipac.caltech.edu, 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 ($\lambda_{\rm central}\,{=65}\,\mu$m) and
WIDE-S ($\lambda_{\rm central}\,{=90}\,\mu$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_{\rm 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).

###  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 $\mu$m to 203 $\mu$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 $\mu$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.

##  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_{\rm s}$ band at
11.5$\,{\pm}\,$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_{\rm s}$ and $H$ catalog magnitudes of $11.87\,{\pm}\,0.01$ and
$15.54\,{\pm}\,0.06$ mag imply a color index of $H-K_{\rm
s}\,{=}\,3.67\,{\pm}\,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$\sigma$ range of 0.4 mag, and a peak-to-peak variation of
0.79 mag, during five years of VVV $K_{\rm s}$-band monitoring. Its mean
brightness during the VVV monitoring amounts to $K_{\rm s}\,{=}\,$12.23 mag.

Our GROND imaging detected it in the $J,H$ and $K_{\rm s}$ bands, but not at
shorter wavelengths. The GROND Ks magnitude of $11.9\,{\pm}\,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$\aas@@fstack{\prime\prime}$2, is consistent with the secondary hot molecular
core and dust continuum source MM3 detected by ALMA, which is located
1$\aas@@fstack{\prime\prime}$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_{\rm 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.

Figure 1: top: VVV $JHK_{s}$ 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
$\mu$m. bottom: GROND $JHK_{s}$ 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

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

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 applied333See
http://wise2.ipac.caltech.edu/docs/release/neowise/expsup/sec2_1civa.html 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 $\mu$m) and W2 (4.6 $\mu$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$\sigma$ detection limit for W1 and W2 of $\sim$ 0.5 mag, that is
a joint 2.8$\sigma$ 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 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 $\mu$m flux can be derived from the
stacked VVV $K_{\rm 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_{\rm s}\,{\sim}\,$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_{\rm 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 $\mu$m upper limit of 0.66 mJy.

####  3.3.2 (NEO)WISE

Figure 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$\aas@@fstack{\prime\prime}$9 from the (0,0) position. The yellow circle
indicates the W1 image FWHM of 6$\aas@@fstack{\prime\prime}$1\.

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 ${\geq}$ 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$\aas@@fstack{\prime\prime}$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$\aas@@fstack{\prime\prime}$1 in the W1
band444See
https://wise2.ipac.caltech.edu/docs/release/allwise/expsup/sec1_2.html.

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 (Sec. 3.2).

####  3.3.3 Spitzer MIPS

As indicated by Fig. 1 (top) the positions of the IR counterparts at
wavelengths up to 24 $\mu$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 $\mu$m flux of MM1 was considerably smaller than the joint flux of MM1
and MM3 (4.58$\,{\pm}\,$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$\aas@@fstack{\prime\prime}$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
$\mu$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 $\mu$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 $\mu$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 $\mu$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 $\gtrsim$ 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.

Table 2: FIFI-LS photometry: Left first epoch, right second epoch

Central | Flux | —Residual— | Flux | —Residual—
---|---|---|---|---
wavelength | density |  | density |
$[\mu m]$ | [Jy] | [Jy] | [Jy] | [Jy]
52.0 |  95.6 |  7.7 |  |
54.8 | 126.8 | 15.0 |  |
60.7 | 125.6 |  5.7 |  |
87.2 | 217.6 |  4.4 |  |
118.6 | 255.1* | 37.0 | 216.5* | 42.6
124.2 | 371.5 | 72.8 | 304.5 | 52.4
142.2 | 315.4 |  8.1 | 224.3 |  2.2
153.3 | 270.4 | 34.5 | 221.4 | 12.8
162.8 | 297.1 |  1.9 | 162.3 | 29.7
186.0 | 284.0 |  9.4 | 155.9 |  9.3

  * *

recalibrated value

###  3.5 FIFI-LS image analysis

The angular distance of 1$\aas@@fstack{\prime\prime}$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 $\mu$m amounts to 5$\aas@@fstack{\prime\prime}$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.

Table 3: FIFI-LS bivariate Gaussfit results Effective wavelength | Axis ratio | Major axis PA
---|---|---
$[\mu m]$ |  | [°]
52.0 | $1.12\pm 0.03$ | $310.1\pm 6.1$
54.8 | $1.22\pm 0.02$ | $347.4\pm 1.9$
60.7 | $1.11\pm 0.02$ | $336.4\pm 3.6$
87.2 | $1.10\pm 0.01$ | $330.6\pm 1.1$

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 $\mu$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 $\mu$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$\aas@@fstack{\prime\prime}$8 ${\times}$
9$\aas@@fstack{\prime\prime}$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
team555https://www.cosmos.esa.int/documents/12133/996891/PACS+photometer+point+spread+function.
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$\aas@@fstack{\prime\prime}$1 away from the ALMA
position of component MM1, and 0$\aas@@fstack{\prime\prime}$7 away from ALMA
component MM3. The absolute astrometry of the Herschel Spacecraft at the end
of the mission came with a 3$\sigma$ uncertainty of
2$\aas@@fstack{\prime\prime}$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 $\mu$m emission is dominated by MM1. Given our PSF fitting
experiments we conclude that a potential minor contribution at 70 $\mu$m from
MM3 must be more than 70 times weaker than the one from MM1.

Figure 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 $\mu$m MIPS limit is indicated
as well. The solid lines represent the reddened graybody fits (cf. 4.1). The
second epoch FIFI-LS data and the corresponding fit appear bf in green.

##  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 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_{\rm d}^{\rm pre}{=}\,26.4\,{\pm}\,0.5$ K which is slightly lower than the
pre-burst value of $28.5\,{\pm}\,1.5$ K by Brogan et al. (2019). The dust
emissivity index amounts to $\beta\,{=}\,1.83\,{\pm}\,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_{\rm d}^{\rm burst}{=}\,30.1\,{\pm}\,0.5$ K as well as the solid angle of
the emission which increased by a factor of 1.20$\,{\pm}\,0.14$. However, the
$\raisebox{2.0pt}{$\chi$}^{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 $\mu$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 $T_{\rm d}^{\rm
post}{=}\,28.9\,{\pm}\,0.6$ K.

Figure 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\,\mu$m MM1 dominates the total
flux density. The seemingly gap at $10\,\mu$m in the MM1 SED is due to strong
silicate absorption.

###  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 Sec. 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 Sec. 5). Since the bulk of the energy
is emitted in the FIR, the FIR/(sub)mm luminosity estimates only weakly depend
on $A_{\rm 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^{\rm pre}_{\rm
FIR}\,{=}\,7600\,{\mathbin{\smash{\raisebox{2.79857pt}{$\underset{\raisebox{0.21529pt}{$\smash{-}$}}{\smash{+}}$}}}}\,^{\leavevmode\nobreak\
800}_{1400}\,{\rm L}_{\odot}$, first FIFI-LS epoch $L^{\rm burst}_{\rm
FIR}\,{=}\,19\,300\,{\mathbin{\smash{\raisebox{2.79857pt}{$\underset{\raisebox{0.21529pt}{$\smash{-}$}}{\smash{+}}$}}}}\,^{2200}_{3700}\,{\rm
L}_{\odot}$, and second epoch $L^{\rm post}_{\rm
FIR}\,{=}\,12\,700\,{\mathbin{\smash{\raisebox{2.79857pt}{$\underset{\raisebox{0.21529pt}{$\smash{-}$}}{\smash{+}}$}}}}\,^{1500}_{2600}\,{\rm
L}_{\odot}$, 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 $\Delta
L^{{\rm burst}}_{\rm
FIR}\,{=}\,11\,700\,\mathbin{\smash{\raisebox{2.79857pt}{$\underset{\raisebox{0.21529pt}{$\smash{-}$}}{\smash{+}}$}}}^{\,2300}_{\,3900}\,{\rm
L}_{\odot}$ and $\Delta L^{{\rm post}}_{\rm
FIR}\,{=}\,5100\,{\mathbin{\smash{\raisebox{2.79857pt}{$\underset{\raisebox{0.21529pt}{$\smash{-}$}}{\smash{+}}$}}}}\,^{\,1700}_{\,2900}\,{\rm
L}_{\odot}$. 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 $\Delta t$ of
$869\,{\pm}\,303$ days, with a pre-burst-level return date of 2021 July 31.

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

##  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 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(H2O) (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\,{\pm}\,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 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 5, 6 and 7 (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.

###  5.2 Radiative transfer analysis

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^{\rm 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”666https://zenodo.org/record/235786 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
${\raisebox{2.0pt}{$\chi$}}^{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_{\rm min},\,r_{\rm max}$, $M_{\rm
disk}^{\rm dust},\,L$ via $R_{\star},\,T_{\star}$), 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_{\rm 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_{\rm sub}\,{=}\,1600\,{\rm 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\arcsec{}$ has been used for fitting
the SED. The fluxes of the MM3 SED are listed in Tab. 4.

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 $\lambda\,{\geq}\,40\,\mu$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_{\rm disk}^{\rm
dust}\,{=}\,0.068\,{\mathbin{\smash{\raisebox{2.79857pt}{$\underset{\raisebox{0.21529pt}{$\smash{-}$}}{\smash{+}}$}}}}\,_{0.043}^{0.11}\,{\rm
M}_{\odot}$. This may seem heavy keeping the canonical gas-to-dust mass ratio
($\gamma$) of 100 in mind. However, the latter is not applicable since
$\gamma$ depends on the galactocentric distance $R_{\rm GC}$ (Giannetti et
al., 2017, Eq. 2,) because of the Galactic metallicity gradient. For G358,
$R_{\rm GC}{=}$1.6 kpc implies a value of $\gamma\,{=}\,38\,{\pm}\,4$ which
yields a total disk mass of
$2.6\,{\mathbin{\smash{\raisebox{2.79857pt}{$\underset{\raisebox{0.21529pt}{$\smash{-}$}}{\smash{+}}$}}}}\,_{1.0}^{6.8}\,{\rm
M}_{\odot}$. The fit delivers $A_{V}\,{=}\,20\pm 5\,$mag and an inclination of
$i\,{=}\,51\,{\pm}\,9$°. The inner radius of the disk seems to be close to the
star, $r_{\rm inner}\,{<}\,3.5\,R_{\rm sub}$ holds for each of the ten best
models . The mean value of the disk outer radius amounts to
$510\,{\mathbin{\smash{\raisebox{2.79857pt}{$\underset{\raisebox{0.21529pt}{$\smash{-}$}}{\smash{+}}$}}}}\,_{240}^{450}\,$au,
while the maximum is as high as $1600\,$au. The mean luminosity of
$7540\,{\mathbin{\smash{\raisebox{2.79857pt}{$\underset{\raisebox{0.21529pt}{$\smash{-}$}}{\smash{+}}$}}}}\,_{3130}^{5340}\,{\rm
L}_{\odot}$ corresponds to a ZAMS star of $11\,{\rm M}_{\odot}$ and $4.2\,{\rm
R}_{\odot}$ (Tout et al., 1996). While most of the models have luminosities
below $10\,000\,{\rm L}_{\odot}$, one has a luminosity of $43\,000\,{\rm
L}_{\odot}$. 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 ${\sim}\,50\,\mu$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
$\raisebox{2.0pt}{$\chi$}^{2}$ values) can be found (together with the
weighted means and standard deviation $\sigma$) in Table 8.

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 $\mu$m fluxes were
dropped. The best model for this SED shows fluxes exceeding those of the
nominal fit (cf. Fig. 11). Thus, we conclude that the FIR part of the nominal
MM3 SED fit is well constrained by the incorporated fluxes longward of 20
$\mu$m.

####  5.2.2 Modeling of the dust continuum radiation of MM1

Figure 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\,\mu$m (post- and burst) fluxes are very similar, their
symbols cannot easily be distinguished.

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\,{\pm}\,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_{\rm 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_{\rm 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_{\rm 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_{\rm
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^{\rm 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. 10.)
The range of the observed fluxes from the burst- and post-burst-epoch is well
covered by the models. Only the burst observations at $\lambda\,{=}\,163$ and
$186\,\mu$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 $\lambda\,{>}\,890\,\mu$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\,\mu$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\arcsec{}$ has been
applied in for the fits. With this choice the
$\raisebox{2.0pt}{$\chi$}^{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 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 9.

The RT modeling indicates that during the burst the luminosity increases from
pre-burst level of $L^{\rm
pre}\,{=}\,5000\,{\mathbin{\smash{\raisebox{2.79857pt}{$\underset{\raisebox{0.21529pt}{$\smash{-}$}}{\smash{+}}$}}}}\,_{\,900}^{1100}\,{\rm
L}_{\odot}$ to $L^{\rm
burst}\,{=}\,23\,400\,{\mathbin{\smash{\raisebox{2.79857pt}{$\underset{\raisebox{0.21529pt}{$\smash{-}$}}{\smash{+}}$}}}}\,_{3700}^{4400}\,{\rm
L}_{\odot}$. At the post-burst epoch, the luminosity is still elevated at a
level of $L^{\rm
post}\,{=}\,12\,400\,{\mathbin{\smash{\raisebox{2.79857pt}{$\underset{\raisebox{0.21529pt}{$\smash{-}$}}{\smash{+}}$}}}}\,_{1700}^{2000}\,{\rm
L}_{\odot}$.

This corresponds to a relative luminosity gain by a factor of
$4.7\,{\mathbin{\smash{\raisebox{2.79857pt}{$\underset{\raisebox{0.21529pt}{$\smash{-}$}}{\smash{+}}$}}}}\,_{1.5}^{2.1}$
during the burst and by
$2.5\,{\mathbin{\smash{\raisebox{2.79857pt}{$\underset{\raisebox{0.21529pt}{$\smash{-}$}}{\smash{+}}$}}}}\,_{0.8}^{1.1}$
for the post-burst epoch. The increase in luminosity during the burst amounts
to $\Delta L^{\rm
burst}\,{=}\,18\,000\,{\mathbin{\smash{\raisebox{2.79857pt}{$\underset{\raisebox{0.21529pt}{$\smash{-}$}}{\smash{+}}$}}}}\,_{4900}^{6100}\,{\rm
L}_{\odot}$ which implies a total burst energy of $E_{\rm
acc}\,{=}\,2.9\mathbin{\smash{\raisebox{2.79857pt}{$\underset{\raisebox{0.21529pt}{$\smash{-}$}}{\smash{+}}$}}}_{0.8}^{1.1}\,{\times}\,10^{38}$
J, where the calculation has been done similar to Sec. 4.2. The derived linear
decay time of $\Delta
t\,{=}\,907\,{\mathbin{\smash{\raisebox{2.79857pt}{$\underset{\raisebox{0.21529pt}{$\smash{-}$}}{\smash{+}}$}}}}\,_{236}^{318}$
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 Sec. 4.2).

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\,{\pm}\,10$
mag. The system has a low inclination of $22\,{\pm}\,11\degr$, 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\,{\pm}\,10$°. The stellar radius amounts to
$8.4\,{\mathbin{\smash{\raisebox{2.79857pt}{$\underset{\raisebox{0.21529pt}{$\smash{-}$}}{\smash{+}}$}}}}\,_{5.5}^{15.7}{\,\rm
R}_{\odot}$. 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\,{\mathbin{\smash{\raisebox{2.79857pt}{$\underset{\raisebox{0.21529pt}{$\smash{-}$}}{\smash{+}}$}}}}\,_{580}^{1500}$
au. The favored dust mass of the circumstellar disk is as low as
$8.4\,{\times}\,10^{-5}\,{\rm M}_{\odot}$, where the $1\sigma$-confidence
interval extends from $1.1\,{\times}\,10^{-6}$ to
$6.1\,{\times}\,10^{-3}\,{\rm M}_{\odot}$. 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\,{\times}\,10^{-3}\,{\rm M}_{\odot}$ within an uncertainty range
of $4.4\,{\times}\,10^{-5}$ to $0.24\,{\rm M}_{\odot}$. 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\,{\rm M}_{\odot}$ 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
$\raisebox{2.0pt}{$\chi$}^{2}$ values) can be found (together with the
weighted means and standard deviation $\sigma$) in Table 9.

The above given properties have been derived from the results of all MM1 SED-
fits, pre-burst, post-burst and burst. The respective
$\raisebox{2.0pt}{$\chi$}^{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
$\raisebox{2.0pt}{$\chi$}^{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 pre-burst-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).

##  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 mass and radius estimates which
reproduce the MM1 pre-burst luminosity, assuming solar metallicity. For
$L^{\rm
pre}\,{=}\,5000\,{\mathbin{\smash{\raisebox{2.79857pt}{$\underset{\raisebox{0.21529pt}{$\smash{-}$}}{\smash{+}}$}}}}\,_{\,900}^{1100}\,{\rm
L}_{\odot}$, these correspond to
$9.7\,{\mathbin{\smash{\raisebox{2.79857pt}{$\underset{\raisebox{0.21529pt}{$\smash{-}$}}{\smash{+}}$}}}}\,_{0.6}^{0.3}\,{\rm
M}_{\odot}$ and
$3.9\,{\mathbin{\smash{\raisebox{2.79857pt}{$\underset{\raisebox{0.21529pt}{$\smash{-}$}}{\smash{+}}$}}}}\,_{0.2}^{0.1}\,{\rm
R}_{\odot}$ (Tout et al., 1996), respectively. The accreted mass is inferred
from $E_{\rm acc}\,{=}\,GM_{*}M_{\rm acc}/R_{*}$, where G is the gravitational
constant, $M_{*}$ is the stellar mass, $M_{\rm acc}$ is the accreted mass, and
$E_{\rm 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.

Figure 7: Dependence of the derived accretion rate on the average luminosity
increase ${<}\Delta L^{\rm acc}{>}$ and assumed ZAMS stellar mass for the
parameter range bracketing MM1. The horizontal dashed line corresponds to
${<}\Delta L^{\rm acc}{>}$ and the light gray region marks its 1$\sigma$
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.

Finally, the mass accretion rate will be obtained from $\dot{M}_{\rm
acc}{=}M_{\rm acc}/\Delta t_{\rm acc}$. We have to emphasize here that,
generally, the duration of the enhanced accretion $\Delta t_{\rm acc}$ will be
shorter than that of the elevated emission $\Delta 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 $\Delta t_{\rm acc}$ of about two months can be derived with a
presumed uncertainty range of ${\pm}$5 days. With the above quantities we
obtain $M_{\rm
acc}\,{=}\,3.1\mathbin{\smash{\raisebox{2.79857pt}{$\underset{\raisebox{0.21529pt}{$\smash{-}$}}{\smash{+}}$}}}^{1.2}_{0.9}\,{\times}\,10^{-4}\,{\rm
M}_{\odot}$, and $\dot{M}_{\rm
acc}{=}1.8\mathbin{\smash{\raisebox{2.79857pt}{$\underset{\raisebox{0.21529pt}{$\smash{-}$}}{\smash{+}}$}}}^{1.2}_{1.1}\,{\times}\,10^{-3}\,{\rm
M}_{\odot}\rm yr^{-1}$, where the errors are dominated by the error of $\Delta
t_{\rm acc}$ and $E_{\rm 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\,{\pm}3\,{\rm
M}_{\odot}$ derived from the kinematic model of the spiral-arm accretion flows
(Chen et al., 2020b). This leads to $M_{\rm
acc}\,{=}\,5.3\mathbin{\smash{\raisebox{2.79857pt}{$\underset{\raisebox{0.21529pt}{$\smash{-}$}}{\smash{+}}$}}}_{4.4}^{11.1}\,{\times}\,10^{-4}\,{\rm
M}_{\odot}$, and $\dot{M}_{\rm
acc}\,{=}\,3.2\mathbin{\smash{\raisebox{2.79857pt}{$\underset{\raisebox{0.21529pt}{$\smash{-}$}}{\smash{+}}$}}}^{5.4}_{3.0}\,{\times}\,10^{-3}\,{\rm
M}_{\odot}\rm 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.

Figure 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
$n_{\rm H_{2}}\,{=}\,[0.2,0.3,0.5,1]\,{\times}\,10^{8}\,\rm cm^{-3}$ 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,
Burns et al., 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 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_{\rm H_{2}}\,{=}\,[1,0.5,0.3,0.2]\,{\times}\,10^{8}\,\rm 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
$\rho\,{=}\,{m_{\rm H_{2}}}{N_{A}^{-1}}\,{\times}\,n_{\rm H_{2}}\gamma^{-1}$,
where $m_{\rm H_{2}}$ is the molar mass of H2, $N_{A}$ is the Avogadro
constant and $\gamma$ 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 $\rm 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 $\rm 94\,K$
(orange, minimum temperature for thermal methanol desorption to occur) and
$\rm 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 ${\sim}\,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 pre-burst and burst.

Another important parameter for the maser excitation is the specific column
density, $N_{\rm M}{\Delta V}^{-1}$, where $N_{\rm M}$ is the methanol column
density along the line of sight and ${\Delta V}$ the velocity range of the
molecules. The low inclination $i$ obtained from our RT modeling and by Chen
et al. (2020b) implies that $\Delta 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.

##  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 $\mu$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.

Figure 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 $\lambda\,{>}\,10\mu 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).

A recent study by Contreras Peña et al. (2020), based on MIR (NEOWISE W2, 4.6
$\mu$m) and (sub)mm (SCUBA-2, 850 $\mu$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$\to$B) and the burst to post-burst (B$\to$Po) epochs, respectively. The
results are presented in the lower panel of Fig. 9. The two curves show
$m(\lambda)=|\,log(F_{\rm x}(\lambda)/F_{\rm y}(\lambda))\,{\times}\,{\Delta
t}^{-1}\,|$ for each wavelength of the modeled SEDs. The indices $x,y$ denote
the respective epoch pairs and $\Delta 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(\lambda)$ and had to assume a
proportionality such that $m({\rm 4.6\mu m})\,{=}\,\eta\,m({\rm 850\mu m})$.
From objects with the most significant variability at both wavelengths, they
derived $\eta\,{=}\,5.53\,{\pm}\,0.29$. This can be compared to the $\eta$
values based on $m$ at 4.6 and 850 $\mu$m from our results (black dots in Fig.
9). These amount to $\eta_{\rm Pr\to B}\,{=}\,4.48$ and $\eta_{\rm B\to
Po}\,{=}\,5.88$ for rise and drop, respectively. The values from our models
are similar, although $\eta_{\rm Pr\to B}$ is somewhat lower.

More importantly, our approach allowed us to derive $m(\lambda)$ for the whole
range of the SED. The corresponding curves show that beyond
$\lambda\,{\gtrsim}\,10\,\mu$m, $log(m)$ depends approximately linear on
$log(\lambda)$. Thus, the wavelength dependence $m(\lambda)$ actually
resembles a power law. The black lines show the respective fits to both epoch
pairs. They yielded the following relations for $m(\lambda)$: $m_{\rm Pr\to
B}\,{\propto}\,\lambda^{-0.42\pm 0.01}$ and $m_{\rm B\to
Po}\,{=}\,(1.16\mathbin{\smash{\raisebox{2.79857pt}{$\underset{\raisebox{0.21529pt}{$\smash{-}$}}{\smash{+}}$}}}^{0.04}_{0.03})\times\,\lambda^{-0.43\pm
0.01}$, 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 $\eta$ 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.

##  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 $\lambda\,{=}\,889\,\mu$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 10). 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\,\mu$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_{\rm 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\mathbin{\smash{\raisebox{2.79857pt}{$\underset{\raisebox{0.21529pt}{$\smash{-}$}}{\smash{+}}$}}}_{0.8}^{1.1}\,{\times}\,10^{38}$
J (see Sec. 5.2.2), $1.2\pm{0.4}\,{\times}\,10^{39}$ J (Caratti o Garatti et
al., 2017), and $8\,{\times}\,10^{38}$ 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 Sec. 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 104 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$\dots 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., 2018; 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 the first step toward a statistics of MYSO
burst properties. Assuming a ZAMS star, an accreted mass of
${\sim}\,3.4\,{\times}\,$10-3 M⊙ has been estimated for S255IR-NIRS3. The
corresponding accretion rate amounts to 5$\,{\pm}\,$2 ${\times}\,$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.

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

###### 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
(grant agreement No. 743029). AS acknowledges support by the Ministry of
Science and Higher Education of the Russian Federation under the grant
075-15-2020-780 (contract 780-10). AV 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 Observatory for 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.

## References

  * Ábrahám et al. (2009) Ábrahám, P., Juhász, A., Dullemond, C. P., et al. 2009, Nature, 459, 224
  * Audard et al. (2014) Audard, M., Ábrahám, P., Dunham, M. M., et al. 2014, in Protostars and Planets VI, ed. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning, 387
  * Boudet et al. (2005) Boudet, N., Mutschke, H., Nayral, C., et al. 2005, ApJ, 633, 272
  * Breen et al. (2013) Breen, S. L., Ellingsen, S. P., Contreras, Y., et al. 2013, MNRAS, 435, 524
  * Breen et al. (2019) Breen, S. L., Sobolev, A. M., Kaczmarek, J. F., et al. 2019, ApJ, 876, L25
  * Brogan et al. (2019) Brogan, C. L., Hunter, T. R., Towner, A. P. M., et al. 2019, ApJ, 881, L39
  * Burns et al. (2020) Burns, R. A., Sugiyama, K., Hirota, T., et al. 2020, Nature Astronomy, 4, 506
  * Caratti o Garatti et al. (2017) Caratti o Garatti, A., Stecklum, B., Garcia Lopez, R., et al. 2017, Nature Physics, 13, 276
  * Chen et al. (2020a) Chen, X., Sobolev, A. M., Breen, S. L., et al. 2020a, ApJ, 890, L22
  * Chen et al. (2020b) Chen, X., Sobolev, A. M., Ren, Z.-Y., et al. 2020b, Nature Astronomy
  * Clarke et al. (2015) Clarke, M., Vacca, W. D., & Shuping, R. Y. 2015, in Astronomical Society of the Pacific Conference Series, Vol. 495, Astronomical Data Analysis Software an Systems XXIV (ADASS XXIV), ed. A. R. Taylor & E. Rosolowsky, 355
  * Cohen et al. (2003) Cohen, M., Wheaton, W. A., & Megeath, S. T. 2003, AJ, 126, 1090
  * Colditz et al. (2018) Colditz, S., Beckmann, S., Bryant, A., et al. 2018, Journal of Astronomical Instrumentation, 7, 1840004
  * Collings et al. (2004) Collings, M. P., Anderson, M. A., Chen, R., et al. 2004, MNRAS, 354, 1133
  * Contreras et al. (2013) Contreras, Y., Schuller, F., Urquhart, J. S., et al. 2013, A&A, 549, A45
  * Contreras Peña et al. (2020) Contreras Peña, C., Johnstone, D., Baek, G., et al. 2020, MNRAS, 495, 3614
  * Cragg et al. (2005) Cragg, D. M., Sobolev, A. M., & Godfrey, P. D. 2005, MNRAS, 360, 533
  * Cutri & et al. (2014) Cutri, R. M. & et al. 2014, VizieR Online Data Catalog, II/328
  * Cutri et al. (2003) Cutri, R. M., Skrutskie, M. F., van Dyk, S., et al. 2003, VizieR Online Data Catalog, II/246
  * Cutri et al. (2012) Cutri, R. M., Wright, E. L., Conrow, T., et al. 2012, Explanatory Supplement to the WISE All-Sky Data Release Products, 1
  * Dartois et al. (2020) Dartois, E., Chabot, M., Bacmann, A., et al. 2020, A&A, 634, A103
  * Diolaiti et al. (2000) Diolaiti, E., Bendinelli, O., Bonaccini, D., et al. 2000, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 4007, Adaptive Optical Systems Technology, ed. P. L. Wizinowich, 879–888
  * Doi et al. (2015) Doi, Y., Takita, S., Ootsubo, T., et al. 2015, PASJ, 67, 50
  * Draine (2003) Draine, B. T. 2003, ARA&A, 41, 241
  * Durjasz et al. (2019) Durjasz, M., Szymczak, M., & Olech, M. 2019, MNRAS, 485, 777
  * Elia & Pezzuto (2016) Elia, D. & Pezzuto, S. 2016, MNRAS, 461, 1328
  * Erickson & Davidson (1993) Erickson, E. F. & Davidson, J. A. 1993, Advances in Space Research, 13, 549
  * Fadda et al. (2019) Fadda, D., Vacca, W. D., Fischer, C., et al. 2019, in American Astronomical Society Meeting Abstracts, Vol. 51, American Astronomical Society Meeting Abstracts #234, 208.05
  * Fazal et al. (2008) Fazal, F. M., Sridharan, T. K., Qiu, K., et al. 2008, ApJ, 688, L41
  * Fischer et al. (2018) Fischer, C., Beckmann, S., Bryant, A., et al. 2018, Journal of Astronomical Instrumentation, 7, 1840003
  * Fujisawa et al. (2015) Fujisawa, K., Yonekura, Y., Sugiyama, K., et al. 2015, The Astronomer’s Telegram, 8286, 1
  * Giannetti et al. (2017) Giannetti, A., Leurini, S., König, C., et al. 2017, A&A, 606, L12
  * Goddi et al. (2018) Goddi, C., Ginsburg, A., Maud, L., Zhang, Q., & Zapata, L. 2018, arXiv e-prints, arXiv:1805.05364
  * Graf et al. (2017) Graf, F., Reinacher, A., Spohr, D., Jakob, H., & Fasoulas, S. 2017, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 10401, Proc. SPIE, 1040112
  * Grave & Kumar (2009) Grave, J. M. C. & Kumar, M. S. N. 2009, A&A, 498, 147
  * Greiner et al. (2008) Greiner, J., Bornemann, W., Clemens, C., et al. 2008, PASP, 120, 405
  * Gutermuth & Heyer (2015) Gutermuth, R. A. & Heyer, M. 2015, AJ, 149, 64
  * Harries (2011) Harries, T. J. 2011, MNRAS, 416, 1500
  * Hartmann & Kenyon (1985) Hartmann, L. & Kenyon, S. J. 1985, ApJ, 299, 462
  * Hartmann & Kenyon (1996) Hartmann, L. & Kenyon, S. J. 1996, ARA&A, 34, 207
  * Herbig (1966) Herbig, G. H. 1966, Vistas in Astronomy, 8, 109
  * Hinz et al. (2009) Hinz, J. L., Rieke, G. H., Yusef-Zadeh, F., et al. 2009, ApJS, 181, 227
  * Hosokawa & Omukai (2009) Hosokawa, T. & Omukai, K. 2009, ApJ, 691, 823
  * Hosokawa et al. (2010) Hosokawa, T., Yorke, H. W., & Omukai, K. 2010, ApJ, 721, 478
  * Hu et al. (2016) Hu, B., Menten, K. M., Wu, Y., et al. 2016, ApJ, 833, 18
  * Hunter et al. (2017) Hunter, T. R., Brogan, C. L., MacLeod, G., et al. 2017, ApJ, 837, L29
  * Hunter et al. (2018) Hunter, T. R., Brogan, C. L., MacLeod, G. C., et al. 2018, ApJ, 854, 170
  * Johnston et al. (2020) Johnston, K. G., Hoare, M. G., Beuther, H., et al. 2020, A&A, 634, L11
  * Johnston et al. (2015) Johnston, K. G., Robitaille, T. P., Beuther, H., et al. 2015, ApJ, 813, L19
  * Johnstone et al. (2013) Johnstone, D., Hendricks, B., Herczeg, G. J., & Bruderer, S. 2013, ApJ, 765, 133
  * Kóspál et al. (2007) Kóspál, Á., Ábrahám, P., Prusti, T., et al. 2007, A&A, 470, 211
  * Kraus et al. (2010) Kraus, S., Hofmann, K.-H., Menten, K. M., et al. 2010, Nature, 466, 339
  * Krühler et al. (2008) Krühler, T., Küpcü Yoldaş, A., Greiner, J., et al. 2008, ApJ, 685, 376
  * Kuiper et al. (2016) Kuiper, R., Turner, N. J., & Yorke, H. W. 2016, ApJ, 832, 40
  * Landsman (1995) Landsman, W. B. 1995, in Astronomical Society of the Pacific Conference Series, Vol. 77, Astronomical Data Analysis Software and Systems IV, ed. R. A. Shaw, H. E. Payne, & J. J. E. Hayes, 437
  * Linz et al. (2020) Linz, H., Bhatia, D., Buinhas, L., et al. 2020, Advances in Space Research, 65, 831
  * Liu et al. (2018) Liu, S.-Y., Su, Y.-N., Zinchenko, I., Wang, K.-S., & Wang, Y. 2018, ApJ, 863, L12
  * Looney et al. (2000) Looney, L. W., Geis, N., Genzel, R., et al. 2000, in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 4014, Airborne Telescope Systems, ed. R. K. Melugin & H.-P. Röser, 14–22
  * Luna et al. (2018) Luna, R., Satorre, M. Á., Domingo, M., et al. 2018, MNRAS, 473, 1967
  * MacFarlane et al. (2019) MacFarlane, B., Stamatellos, D., Johnstone, D., et al. 2019, MNRAS, 487, 4465
  * MacLeod et al. (2018) MacLeod, G. C., Smits, D. P., Goedhart, S., et al. 2018, MNRAS, 478, 1077
  * MacLeod et al. (2019) MacLeod, G. C., Sugiyama, K., Hunter, T. R., et al. 2019, MNRAS, 489, 3981
  * Mainzer et al. (2014) Mainzer, A., Bauer, J., Cutri, R. M., et al. 2014, ApJ, 792, 30
  * Mattsson et al. (2014) Mattsson, L., Gomez, H. L., Andersen, A. C., et al. 2014, MNRAS, 444, 797
  * Mendigutía et al. (2011) Mendigutía, I., Calvet, N., Montesinos, B., et al. 2011, A&A, 535, A99
  * Menten (1991a) Menten, K. M. 1991a, in Astronomical Society of the Pacific Conference Series, Vol. 16, Atoms, Ions and Molecules: New Results in Spectral Line Astrophysics, ed. A. D. Haschick & P. T. P. Ho, 119–136
  * Menten (1991b) Menten, K. M. 1991b, ApJ, 380, L75
  * Meyer et al. (2019a) Meyer, D. M. A., Haemmerlé, L., & Vorobyov, E. I. 2019a, MNRAS, 484, 2482
  * Meyer et al. (2019b) Meyer, D. M. A., Vorobyov, E. I., Elbakyan, V. G., et al. 2019b, MNRAS, 482, 5459
  * Meyer et al. (2017) Meyer, D. M. A., Vorobyov, E. I., Kuiper, R., & Kley, W. 2017, MNRAS, 464, L90
  * Minniti et al. (2017) Minniti, D., Lucas, P., & VVV Team. 2017, VizieR Online Data Catalog, 2348
  * Minniti et al. (2010) Minniti, D., Lucas, P. W., Emerson, J. P., et al. 2010, New A, 15, 433
  * Molinari et al. (2016) Molinari, S., Schisano, E., Elia, D., et al. 2016, VizieR Online Data Catalog, J/A+A/591/A149
  * Molinari et al. (2010) Molinari, S., Swinyard, B., Bally, J., et al. 2010, PASP, 122, 314
  * Moscadelli et al. (2017) Moscadelli, L., Sanna, A., Goddi, C., et al. 2017, A&A, 600, L8
  * Neckel & Staude (1995) Neckel, T. & Staude, H. J. 1995, ApJ, 448, 832
  * Olech et al. (2020) Olech, M., Szymczak, M., Wolak, P., Gérard, E., & Bartkiewicz, A. 2020, A&A, 634, A41
  * Omont et al. (2003) Omont, A., Gilmore, G. F., Alard, C., et al. 2003, A&A, 403, 975
  * Ostrovskii & Sobolev (2002) Ostrovskii, A. B. & Sobolev, A. M. 2002, in IAU Symposium, Vol. 206, Cosmic Masers: From Proto-Stars to Black Holes, ed. V. Migenes & M. J. Reid, 183
  * Parsons et al. (2018) Parsons, H., Dempsey, J. T., Thomas, H. S., et al. 2018, ApJS, 234, 22
  * Pomohaci et al. (2019) Pomohaci, R., Oudmaijer, R. D., & Goodwin, S. P. 2019, MNRAS, 484, 226
  * Pomohaci et al. (2017) Pomohaci, R., Oudmaijer, R. D., Lumsden, S. L., Hoare, M. G., & Mendigutía, I. 2017, MNRAS, 472, 3624
  * Povich et al. (2011) Povich, M. S., Smith, N., Majewski, S. R., et al. 2011, ApJS, 194, 14
  * Pringle (1981) Pringle, J. E. 1981, ARA&A, 19, 137
  * Proven-Adzri et al. (2019) Proven-Adzri, E., MacLeod, G. C., Heever, S. P. v. d., et al. 2019, MNRAS, 487, 2407
  * Rab et al. (2017) Rab, C., Elbakyan, V., Vorobyov, E., et al. 2017, A&A, 604, A15
  * Ramírez et al. (2008) Ramírez, S. V., Arendt, R. G., Sellgren, K., et al. 2008, ApJS, 175, 147
  * Robitaille (2011) Robitaille, T. P. 2011, A&A, 536, A79
  * Robitaille (2017) Robitaille, T. P. 2017, A&A, 600, A11
  * Robitaille et al. (2007) Robitaille, T. P., Whitney, B. A., Indebetouw, R., & Wood, K. 2007, ApJS, 169, 328
  * Robitaille et al. (2006) Robitaille, T. P., Whitney, B. A., Indebetouw, R., Wood, K., & Denzmore, P. 2006, ApJS, 167, 256
  * Sánchez-Portal et al. (2014) Sánchez-Portal, M., Marston, A., Altieri, B., et al. 2014, Experimental Astronomy, 37, 453
  * Schlafly & Finkbeiner (2011) Schlafly, E. F. & Finkbeiner, D. P. 2011, ApJ, 737, 103
  * Schuller et al. (2009) Schuller, F., Menten, K. M., Contreras, Y., et al. 2009, A&A, 504, 415
  * Sobolev et al. (1997) Sobolev, A. M., Cragg, D. M., & Godfrey, P. D. 1997, A&A, 324, 211
  * Spitzer Science (2009) Spitzer Science, C. 2009, VizieR Online Data Catalog, II/293
  * Stecklum et al. (2002) Stecklum, B., Brandl, B., Henning, T., et al. 2002, A&A, 392, 1025
  * Stecklum et al. (2018a) Stecklum, B., Caratti o Garatti, A., Cardenas, M. C., et al. 2018a, in The wonders of star formation, online proceedings, https://www-archive.ph.ed.ac.uk/star-formation/sites/default/files/attachments/stecklum_update.pdf
  * Stecklum et al. (2016) Stecklum, B., Caratti o Garatti, A., Cardenas, M. C., et al. 2016, The Astronomer’s Telegram, 8732, 1
  * Stecklum et al. (2018b) Stecklum, B., Caratti o Garatti, A., Hodapp, K., et al. 2018b, in Astrophysical Masers: Unlocking the Mysteries of the Universe, ed. A. Tarchi, M. J. Reid, & P. Castangia, Vol. 336, 37–40
  * Sturm et al. (2013) Sturm, B., Bouwman, J., Henning, T., et al. 2013, A&A, 553, A5
  * Sugiyama et al. (2019) Sugiyama, K., Saito, Y., Yonekura, Y., & Momose, M. 2019, The Astronomer’s Telegram, 12446, 1
  * Szymczak et al. (2018a) Szymczak, M., Olech, M., Sarniak, R., Wolak, P., & Bartkiewicz, A. 2018a, MNRAS, 474, 219
  * Szymczak et al. (2018b) Szymczak, M., Olech, M., Wolak, P., Gérard, E., & Bartkiewicz, A. 2018b, A&A, 617, A80
  * Tabatabaei et al. (2014) Tabatabaei, F. S., Braine, J., Xilouris, E. M., et al. 2014, A&A, 561, A95
  * Tapia et al. (2015) Tapia, M., Roth, M., & Persi, P. 2015, MNRAS, 446, 4088
  * Tout et al. (1996) Tout, C. A., Pols, O. R., Eggleton, P. P., & Han, Z. 1996, MNRAS, 281, 257
  * Ulrich (1976) Ulrich, R. K. 1976, ApJ, 210, 377
  * Vorobyov & Basu (2015) Vorobyov, E. I. & Basu, S. 2015, ApJ, 805, 115
  * Weingartner & Draine (2001) Weingartner, J. C. & Draine, B. T. 2001, ApJ, 548, 296
  * Wenger et al. (2000) Wenger, M., Ochsenbein, F., Egret, D., et al. 2000, A&AS, 143, 9
  * Whitney et al. (2013) Whitney, B. A., Robitaille, T. P., Bjorkman, J. E., et al. 2013, ApJS, 207, 30
  * Wiebe et al. (2019) Wiebe, D. S., Molyarova, T. S., Akimkin, V. V., Vorobyov, E. I., & Semenov, D. A. 2019, MNRAS, 485, 1843
  * Wolf et al. (2012) Wolf, S., Henning, T., & Stecklum, B. 2012, MC3D: Monte-Carlo 3D Radiative Transfer Code
  * Wright et al. (2010) Wright, E. L., Eisenhardt, P. R. M., Mainzer, A. K., et al. 2010, AJ, 140, 1868
  * Yamamura et al. (2010) Yamamura, I., Makiuti, S., Ikeda, N., et al. 2010, VizieR Online Data Catalog, II/298
  * Young et al. (2012) Young, E. T., Becklin, E. E., Marcum, P. M., et al. 2012, ApJ, 749, L17
  * Zinnecker & Yorke (2007) Zinnecker, H. & Yorke, H. W. 2007, ARA&A, 45, 481

##  Appendix A Additional information

Table 4: MM3-SED. The fluxes at $\lambda\,{\geq}\,889\,\mu$m, are MM3-fluxes,
whereas the fluxes at $\lambda\,{\leq}\,24\,\mu$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.

Wavelength | Flux | Ref.
---|---|---
$[\mu$m] | ${[\rm erg\,cm^{-2}\,s^{-1}}]$ |
$\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ 1.63$ | $1.16\pm 0.07\,{\times}\,10^{-12}$ | 1
$\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ 2.13$ | $1.67\pm 0.02\,{\times}\,10^{-11}$ | 1
$\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ 3.55$ | $3.64\pm 0.02\,{\times}\,10^{-10}$ | 2
$\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ 4.49$ | $6.10\pm 0.03\,{\times}\,10^{-10}$ | 2
$\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ 5.73$ | $9.73\pm 0.05\,{\times}\,10^{-10}$ | 2
$\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ 7.0\leavevmode\nobreak\ $ | $6.90\pm 0.49\,{\times}\,10^{-10}$ | 3
$\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ 7.87$ | $5.30\pm 0.04\,{\times}\,10^{-10}$ | 2
$11.6$ | $3.54\pm 0.07\,{\times}\,10^{-10}$ | 4
$15.0$ | $6.46\pm 0.20\,{\times}\,10^{-10}$ | 3
$22.1$ | $5.13\pm 0.13\,{\times}\,10^{-10}$ | 4
$23.7$ | $4.49\pm 0.13\,{\times}\,10^{-10}$ | 5
$889\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ $ | $1.34\pm 0.04\,{\times}\,10^{-13}$ | 6
$1282\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ $ | $2.43\pm 0.10\,{\times}\,10^{-14}$ | 6
$1420\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ $ | $1.50\pm 0.19\,{\times}\,10^{-14}$ | 6
$1532\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ $ | $1.06\pm 0.08\,{\times}\,10^{-14}$ | 6

(1) Minniti et al. (2017) (VISTA/VIRCAM); (2) Ramírez et al. (2008)
(Spitzer/IRAC); (3) Omont et al. (2003) (ISO/ISOCAM); (4) Cutri & et al.
(2014) (ALLWISE); (5) Gutermuth & Heyer (2015) (Spitzer/MIPS); (6) Brogan et
al. (2019) (ALMA, SMA)

Table 5: Total- and MM1-pre-burst fluxes. We obtained the MM1 flux densities
used for the SED-fit in Sec. 5.2.2 by removing the contribution from all other
sources according to Sec. 5.1. Facilities and instruments are given behind
corresponding references.

Wavelength | Total flux | Ref. | MM1 flux
---|---|---|---
$[\mu$m] | $[{\rm erg\,cm^{-2}\,s^{-1}}]$ |  | $[{\rm erg\,cm^{-2}\,s^{-1}}]$
$2.15$**$*$Upper limit | $5.86\pm 0.59\,{\times}\,10^{-13}$ | 1 |
$24$**$*$Upper limit | $5.25\,{\pm}0.53\,{\times}\,10^{-11}$ | 1 |
$65$ | $3.10\pm 0.06\,{\times}\,10^{-9}$ | 1 | $1.55\pm 0.03\,{\times}\,10^{-9}$
$70$ | $3.13\pm 0.38\,{\times}\,10^{-9}$ | 1 | $1.57\pm 0.06\,{\times}\,10^{-9}$
$160$ | $1.90\pm 0.21\,{\times}\,10^{-9}$ | 1 | $9.50\pm 0.52\,{\times}\,10^{-10}$
$250$ | $7.82\pm 1.04\,{\times}\,10^{-10}$ | 1 | $3.91\pm 0.27\,{\times}\,10^{-10}$
$350$ | $2.34\pm 0.52\,{\times}\,10^{-10}$ | 1 | $1.17\pm 0.26\,{\times}\,10^{-10}$
$500$ | $4.98\pm 2.82\,{\times}\,10^{-11}$ | 1 | $2.49\pm 0.33\,{\times}\,10^{-11}$
$850$ | $4.73\pm 0.15\,{\times}\,10^{-12}$ | 2 | $2.36\pm 0.08\,{\times}\,10^{-12}$
$870$ | $4.03\pm 0.07\,{\times}\,10^{-12}$ | 3 | $2.00\pm 0.35\,{\times}\,10^{-12}$

(1) present paper (VISTA/VIRCAM, Spitzer/MIPS, SOFIA/FIFI-LS); (2) Parsons et
al. (2018) (JCMT/SCUBA-2); (3) Brogan et al. (2019) (APEX/LABOCA) 777

Table 6: Total- and MM1-burst fluxes: Similar to Table 5 but for the burst
epoch. We assume that during the burst only the luminosity of MM1 increased
while all other sources remained constant. We note that we do not provide MM1
fluxes for $\lambda>890\mu m$, since the data points at theses wavelengths are
discarded for the burst SED fit as discussed in Sec. 9.1. Facilities and
instruments are given behind corresponding references.

wavelength | total flux | Ref. | MM1 flux
---|---|---|---
$[\mu$m] | $[{\rm erg\,cm^{-2}\,s^{-1}}]$ |  | $[{\rm erg\,cm^{-2}\,s^{-1}}]$
$\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ 2.15$**$*$Upper limit | $9.20\pm 0.92\,{\times}\,10^{-13}$ | 1 |
$\leavevmode\nobreak\ \leavevmode\nobreak\ 3.4\leavevmode\nobreak\ $**$*$Upper limit | $8.82\pm 0.89\,{\times}\,10^{-11}$ | 1 |
$\leavevmode\nobreak\ \leavevmode\nobreak\ 4.6\leavevmode\nobreak\ $**$*$Upper limit | $2.93\pm 0.30\,{\times}\,10^{-10}$ | 1 |
$52.0$ | $5.51\pm 0.56\,{\times}\,10^{-9}$ | 1 | $4.85\pm 0.49\,{\times}\,10^{-9}$
$54.8$ | $6.94\pm 0.70\,{\times}\,10^{-9}$ | 1 | $6.12\pm 0.62\,{\times}\,10^{-9}$
$60.7$ | $6.20\pm 0.62\,{\times}\,10^{-9}$ | 1 | $5.04\pm 0.51\,{\times}\,10^{-9}$
$87.2$ | $7.48\pm 0.75\,{\times}\,10^{-9}$ | 1 | $5.35\pm 0.54\,{\times}\,10^{-9}$
$118.6\leavevmode\nobreak\ $ | $6.45\pm 0.65\,{\times}\,10^{-9}$ | 1 | $4.59\pm 0.46\,{\times}\,10^{-9}$
$124.2\leavevmode\nobreak\ $ | $8.97\pm 0.90\,{\times}\,10^{-9}$ | 1 | $7.22\pm 0.72\,{\times}\,10^{-9}$
$142.2\leavevmode\nobreak\ $ | $5.26\pm 0.53\,{\times}\,10^{-9}$ | 1 | $5.25\pm 0.53\,{\times}\,10^{-9}$
$153.3\leavevmode\nobreak\ $ | $5.29\pm 0.53\,{\times}\,10^{-9}$ | 1 | $4.08\pm 0.41\,{\times}\,10^{-9}$
$162.8\leavevmode\nobreak\ $ | $5.47\pm 0.55\,{\times}\,10^{-9}$ | 1 | $4.41\pm 0.45\,{\times}\,10^{-9}$
$186.4\leavevmode\nobreak\ $ | $4.58\pm 0.46\,{\times}\,10^{-9}$ | 1 | $3.82\pm 0.39\,{\times}\,10^{-9}$
$889\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ $ | $3.81\pm 0.11\,{\times}\,10^{-12}$ | 2 | $1.72\pm 0.07\,{\times}\,10^{-12}$
$1282\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ $ | $4.21\pm 0.03\,{\times}\,10^{-13}$ | 2 |
$1420\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ $ | $2.74\pm 0.28\,{\times}\,10^{-13}$ | 2 |
$1532\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ $ | $1.88\pm 0.01\,{\times}\,10^{-13}$ | 2 |

(1) present paper (2.2-m MPG/ESO telescope/GROND, NEOWISE, SOFIA/FIFI-LS); (2)
Brogan et al. (2019) (ALMA, SMA) 888

Table 7: Total- and MM1-post-burst fluxes: Similar to Table 5 but for the
post-burst. We assume that only the luminosity of MM1 has changed while all
other sources remained constant. Facilities and instruments are given behind
corresponding references.

wavelength | total flux | Ref. | MM1 flux
---|---|---|---
$[\mu$m] | $[{\rm erg\,cm^{-2}\,s^{-1}}]$ |  | $[{\rm erg\,cm^{-2}\,s^{-1}}]$
$118.6$ | $5.17\pm 0.52\,{\times}\,10^{-9}$ | 1 | $3.55\pm 0.36\,{\times}\,10^{-9}$
$124.2$ | $7.70\pm 0.77\,{\times}\,10^{-9}$ | 1 | $5.55\pm 0.56\,{\times}\,10^{-9}$
$142.2$ | $4.74\pm 0.48\,{\times}\,10^{-9}$ | 1 | $3.28\pm 0.33\,{\times}\,10^{-9}$
$153.3$ | $4.35\pm 0.44\,{\times}\,10^{-9}$ | 1 | $3.09\pm 0.31\,{\times}\,10^{-9}$
$162.8$ | $3.01\pm 0.31\,{\times}\,10^{-9}$ | 1 | $1.90\pm 0.19\,{\times}\,10^{-9}$
$186.4$ | $2.53\pm 0.26\,{\times}\,10^{-9}$ | 1 | $1.73\pm 0.18\,{\times}\,10^{-9}$
$889\leavevmode\nobreak\ \leavevmode\nobreak\ $ | $3.81\pm 0.11\,{\times}\,10^{-12}$ | 2 | $1.72\pm 0.07\,{\times}\,10^{-12}$

(1) present paper (SOFIA/FIFI-LS); (2) Brogan et al. (2019) (ALMA)

Figure 10: 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 $\mu$m, as discussed when introducing
the model pool in Sec. 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 $\lambda\,{\geq}\,118\,\mu$m). Figure 11: 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 $\mu$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.

Table 8: MM3 SED fit results: The table contains the full parameter set as well as the derived luminosities (see Sec. 5.2) for the ten best models, together wit their respective $\raisebox{2.0pt}{$\chi$}^{2}$-values. The weighted mean and standard deviation $\sigma$ is given in the bottom. The confidence interval for log-sampled values extends from $\frac{x}{\sigma}$ to $x\,{\times}\,\sigma$, where for lin-spaced values it is $x\pm\sigma$ as usual. All masses (and densities) are given in unit dust mass, where we use the Weingartner & Draine (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_{disk}^{min}$, its centrifugal radius $r_{C}$ is set to $r_{disk}^{max}$. All models extends to an outer radius, where density and temperature have dropped to the ambient level.  Grid spacing |  | lin | log | log | log | log | log | lin | lin | log | log | lin | lin | log | log | lin | log
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
Model | $\raisebox{2.0pt}{$\chi$}^{2}$ | $\rm A_{V}$ | d | $\rm R_{*}$ | $\rm T_{*}$ | $\rm m_{\rm disk}$ | $\rm r_{\rm disk}^{\rm max}$ | $\beta_{\rm disk}$ | $p_{\rm disk}$ | $h_{100}^{\rm disk}$ | $\rho_{0}^{\rm env}$ | $p_{\rm cav}$ | $\theta_{0}^{\rm cav}$ | $\rho_{0}^{\rm cav}$ | $r_{\rm disk}^{\rm min}$ | inc | $\rm L_{*}$
|  | mag | kpc | ${\rm R}_{\odot}$ | K | ${\rm M}_{\odot}$ | au |  |  | au | g/cc |  | $\degr$ | g/cc | $R_{\rm sub}$ | $\degr$ | ${\rm L}_{\odot}$
aer04B9J | $5842$ | $16.5$ | $6.75$ | $11.5$ | $15720$ | $0.0981$ | $440$ | $1.06$ | $-0.205$ | $9.29$ | $7.28\,{\times}\,10^{-21}$ | $1.16$ | $40.4$ | $1.09\,{\times}\,10^{-23}$ | $3.35$ | $60.5$ | $7161$
aer04B9J | $6799$ | $17.1$ | $6.75$ | $11.5$ | $15720$ | $0.0981$ | $440$ | $1.06$ | $-0.205$ | $9.29$ | $7.28\,{\times}\,10^{-21}$ | $1.16$ | $40.4$ | $1.09\,{\times}\,10^{-23}$ | $3.35$ | $52$ | $7161$
rlHDFVxN | $7251$ | $16.6$ | $6.07$ | $9.32$ | $15680$ | $0.0321$ | $1216$ | $1.13$ | $-1.32$ | $11.8$ | $9.75\,{\times}\,10^{-23}$ | $1.64$ | $32.8$ | $8.59\,{\times}\,10^{-22}$ | $3.05$ | $59.6$ | $4639$