A multi-transition methanol maser study of the accretion burst source G358.93-0.03-MM1

We present the most complete to date interferometric study of the centimeter wavelength methanol masers detected in G358.93-0.03 at the burst and post-burst epochs. A unique, NIR/(sub)mm-dark and FIR-loud MYSO accretion burst was recently discovered in G358.93-0.03. The event was accompanied by flares of an unprecedented number of rare methanol maser transitions. The first images of three of the newly-discovered methanol masers at 6.18, 12.23, and 20.97 GHz are presented in this work. The spatial structure evolution of the methanol masers at 6.67, 12.18, and 23.12 GHz is studied at two epochs. The maser emission in all detected transitions resides in a region of $\sim$0.2$^{\prime\prime}$ around the bursting source and shows a clear velocity gradient in the north-south direction, with red-shifted features to the north and blue-shifted features to the south. A drastic change in the spatial morphology of the masing region is found: a dense and compact"spiral"cluster detected at epoch I evolved into a disperse,"round"structure at epoch II. During the transition from the first epoch to the second, the region traced by masers expanded. The comparison of our results with the complementary VLA, VLBI, SMA, and ALMA maser data is conducted. The obtained methanol maser data support the hypothesis of the presence of spiral-arm structures within the accretion disk, which was suggested in previous studies of the source.


INTRODUCTION
Is there a common mechanism of star formation across the entire stellar mass spectrum? While star formation theory is well established for low-mass stars (e.g., McKee & Ostriker 2007), much less is known about high-mass star formation. Recent studies such as Caratti o Garatti et al. (2017) and Hunter et al. (2017) provide direct evidence that massive young stellar objects (MYSOs), analogously to low-mass stellar objects, form via disk-mediated accretion accompanied by episodic accretion bursts. These bursts may be a result of disk fragmentation (Ahmadi et al. 2019;Meyer et al. 2019b). Such a mechanism provides a means to overcome radiation pressure (Hosokawa et al. 2016), and it is thought that up to half of the final stellar mass may be acquired in these accretion events. But the accretion process itself is poorly understood -largely due to scarce observational evidence. Presently, only three cases of accretion bursts in MYSOs have been reported: S255IR (Caratti o Garatti et al. 2017), NGC6334I (Hunter et al. 2017), and -the topic of this paper -G358.93−0.03 (Stecklum et al. 2021). A potential fourth accretion-burst source, G323.46−0.08, may have recently undergone such an event (Proven-Adzri et al. 2019).
A larger sample is clearly needed, but identifying MYSOs undergoing accretion bursts is difficult. It is in this respect that masers have proven to be a powerful probe of the mechanisms of massive star formation. Masing occurs only within certain ranges of physical conditions of the gas and radiation field (e.g., see the reviews by Ellingsen et al. 2007 andBreen et al. 2019). Thus the spatial distribution of masers can reveal the distribution of temperature, density and radiation enhancements in the region, while the kinematics of maser spots can trace gas motions. All of these properties are essential for understanding the episodic accretion phenomenon.
If the accretion burst augments the local radiation field, the resulting increase in incident photons will cause all masers in the foreground of the continuum emission to increase in flux. Masers, either compact or extended, and covering a wide range of velocities may be involved -as was seen in S255IR (Szymczak et al. 2018), NGC6334I (MacLeod et al. 2018), and recently in G358.93−0.03 Brogan et al. 2019;Breen et al. 2019). Moreover, in all these sources, multiple maser transitions were seen to flare 1 and rare maser transitions appeared. The masers then weakened or disappeared in a span of weeks to months for the rarer masers but longer for the more common transitions (MacLeod et al. 2018).
1.1. G358.93-0.03 On January 14, 2019, a rapid rise of the 6.67 GHz class II methanol maser flux density, accompanied by the appearance of several new velocity features, was detected in G358.93−0.03 ) by the Ibaraki 6.7 GHz Methanol Maser Monitor (iMet) program (Yonekura et al. 2016). The flare of 6.67 GHz maser emission reached a first peak on February 15, 2019, when the flux density of the maser feature at −15.6 km s −1 reached ∼660 Jy (Figure 1). A second peak occurred on March 12, when the velocity feature at −17.2 km s −1 reached a flux density of ∼900 Jy. A daily spectrum of the 6.67 GHz methanol maser in G358.93−0.03 is published regularly on the iMet project website 2 . Intensive follow-up observations were made by the Maser Monitoring Organization 3 (M2O), a self-organized collaboration of maser monitoring observatories. Monitoring with the 26-m telescope of the Hartebeesthoek Radio Astronomy Observatory (HartRAO, Republic of South Africa; Figure 1) confirmed that the flaring 6.67 GHz maser flux density had been rapidly increasing since February of 2019 and in ∼1.5 months reached ∼1000 Jy -about 200× greater than the stable flux density of ∼5 Jy that had been reported by observations spanning the previous decade (Caswell et al. 2010;Chambers et al. 2014;Hu et al. 2016;Rickert et al. 2019).
In addition to the 6.67 GHz methanol maser flare, the appearance of multiple new maser transitions was observed in G358.93−0.03. A rare 23.12 GHz class II methanol maser (only the fourth known occurrence of this transition) was detected. A 12.18 GHz class II methanol maser, previously undetected in this source despite several epochs of observations (Breen et al. 2012), was found to have a higher flux density than the flaring 6.67 GHz maser. The ATCA observations of G358.93−0.03 conducted in March, 2019 detected emission from 12 additional methanol transitions, six of which were previously undetected class II methanol masers; among these were the first known torsionally excited methanol masers (Breen et al. 2019). In January-March, 2019, Volvach et al. (2020) reported a short-lived (44 days) 19.97 GHz methanol maser emission. More transitions may have been missed owing to insufficient frequency coverage. The latest observations with the SMA (March, 2019) and ALMA (April, 2019) yielded an additional 14 methanol maser discoveries in the (sub)millimeter range, primarily from transitions in excited torsional states (Brogan et al. 2019).
G358.93−0.03 was a maser-quiet source at the pre-flare epoch, with only two maser transitions detected: the 6.67 GHz methanol maser (Caswell et al. 2010;Chambers et al. 2014;Hu et al. 2016;Rickert et al. 2019) and 22 GHz water maser (Titmarsh et al. 2016). Moreover, G358.93−0.03 was not as well-studied as S255IR or NGC6334I. For example, no maser proper motion data is available to provide a parallax distance. According to the BeSSeL Revised Kinematic Distance Calculator 4 (Reid et al. 2014), the near kinematic distance to the source is 6.75 +0.37 −0.68 kpc. The SMA and ALMA (sub)millimeter imaging resolved G358.93−0.03 into a cluster of eight continuum sources (Brogan et al. 2019). Two of these sources, MM1 and MM3, were found to host hot molecular cores. All of the newlydiscovered methanol masers were found to be associated with the brightest continuum source, MM1 (e.g. Brogan et al. (2019); Burns et al. (2020)).
The flaring class II methanol maser transitions detected in G358.93−0.03 require similar physical conditions according to theoretical models, albeit over a wider range of values (Sobolev et al. 1997a,b;Cragg et al. 2004Cragg et al. , 2005. Thus, the simultaneous flaring of these masers indicates a sudden change in the local physical conditions -probably related to the radiation field -which could provoke the necessary conditions for maser amplification. The multi-epoch Long Baseline Array (LBA) observations, performed during the 6.67 GHz methanol maser flare, allowed the imaging of the subluminal propagation of a thermal radiation "heatwave" emanating from the accreting high-mass protostar (Burns et al. 2020).
The SMA and ALMA images revealed a partial elliptical ring of the methanol masers, which was interpreted as a coherent physical structure illuminated by a radiative event from the central object (Brogan et al. 2019).
Remarkably, the VLA observations of Chen et al. (2020a,b), conducted during the maser flare epoch in March, 2019, revealed that the discovered rare maser transitions of HDO, HNCO, and 13 CH 3 OH appear to trace a spiral-arm structure around the bursting source (Chen et al. 2020b). This finding was the first observational evidence of a link between accretion bursts and disk substructures (Chen et al. 2020b).
The accretion burst in G358.93−0.03 was decisively confirmed by multi-epoch SOFIA observations (Stecklum et al. 2021). The event is found to be the first NIR/(sub)mm-dark and FIR-loud MYSO accretion burst, showing an increase in the source flux only in the FIR, and not in the NIR or (sub)mm bands (Stecklum et al. 2021).
The excellent coordination of the M2O observatories ensured that monitoring of G358.93−0.03 was done in a timely manner. As a result, we had the unique opportunity to observe a stage of the flare that was missed in the cases of S255 and NGC6334I. In this paper, we present the most complete to date interferometric study of the methanol masers detected in G358.93−0.03 at the burst and post-burst epochs.

OBSERVATIONS AND DATA REDUCTION
Two observing sessions of G358.93−0.03 were carried out with the Karl G. Jansky Very Large Array (VLA) on February 25, 2019 (epoch I) and on June 4, 2019 (epoch II) as the Target of Opportunity programs 19A-448 and 19A-476. The first session was a follow-up observation in response to the rapid rise of the 6.67 GHz methanol maser flux density detected in the source in January 2019 . The second session was made in response to the 22 GHz water maser flare of April 2019 (the M2O monitoring data, Fig. 1).
The epoch I observation was made during a C → B reconfiguration which led to an asymmetrical beam, highly elongated in the N-S direction. At the first epoch, the priority was to catch the flaring methanol masers before they faded, so promptness was the driving factor more than resolution. The epoch II observation was made in Aconfiguration.
In this paper, we present the C-, Ku-, and K-bands spectral line and K-band continuum observations. The observing time for each session was two hours. The pointing coordinates for G358.93−0.03 were RA (J2000) = 17 h 43 m 10.02 s and Dec (J2000) = −29 • 51 45.8 , with an LSR velocity of −15.55 km s −1 . The same calibration sources were used at both epochs: 3C286 was used as a flux density, bandpass, and delay calibrator; J1744−3116 was the complex gain calibrator (with an angular separation from the target source of 4.6 • ).
The maser lines were observed in narrow spectral windows (1, 2 and 4 MHz at C-, Ku-, and K-bands, respectively) with 512 or 1024 channels. Continuum emission was observed in 31 spectral windows with 128 1-MHz channels. Observation parameters, including the synthesized beam size and the rms noise level, for spectral line and continuum data, are presented in Tables 1 and 2. Table 1 contains the list of detected maser transitions (detection of maser emission in each band and epoch is marked with 'Y', non-detection is marked with 'N').
All steps of the post-correlation data reduction were done with the Common Astronomy Software Applications (CASA, McMullin et al. (2007)). For basic flagging and data calibration, we used the VLA CASA Calibration Pipeline. Nevertheless, there were two issues that required special, additional treatment. First, during the observations we used the VLA catalog coordinates for the phase calibrator J1744−3116. Post-observation, we learned that the VLBI Source Position Catalog 5 lists a more reliable position. The discrepancy from the VLA catalog position is 0.313 arcsec. Prior to calibration and imaging we corrected the calibrator position to the VLBI catalog values with CASA task fixvis. Second, the observations were carried out during a multi-frequency maser flare in the source, which led to the detection of a number of methanol maser lines in the wide-band VLA continuum window (for the list of the detected maser lines, see the note to Table 2). The maser lines were not flagged by the pipeline, and their data appeared in the continuum images as a false detection of continuum peaks. Additional manual flagging was done to avoid this contamination.
Calibrated data were imaged with the CASA task clean. Briggs weighting was used for maser data and natural weighting was used for continuum data. Gaussian fitting of the images was performed with the CASA task imfit. A two-dimensional Gaussian brightness distribution was fit to all maser and continuum emission peaks with a flux density above the 3σ level (see Tables 1 and 2 for 1σ-levels) to determine their positions and flux densities. Further in the text, we list and discuss so-called "maser spots" -maser emission detected in a single velocity channel of a data cube.
To estimate the uncertainties in the absolute positions of the masers detected with the VLA, we compared our positions with those obtained by other observations made recently by the M2O collaboration. The 6.7 GHz methanol maser spots detected with the VLA at the first epoch (February 25, 2019) were superposed with the 6.7 GHz data obtained with the LBA on February 28, 2019 from Burns et al. (2020). A N-S position difference was noted and estimated to be ∼0.07 arcsec. This value is within the accuracy of the VLA position measurement and is within one pixel of our C-band VLA images. Absolute coordinates measured with VLA can be affected by uncompensated ionospheric propagation delays, and the effect is most notorious for declinations ≤30 • (Argon et al. 2000). To estimate the displacements for other types of masers, we (1) measured the median positions of the centres of the maser clusters, (2) estimated shifts between centres of the maser clusters and MM1 (the central source of our maps (Brogan et al. 2019)). The estimated offsets are indicated in the notes of Tables 5-9. The offsets were introduced to our data and used in the data analysis and preparation of Figures 3-7, 10-11, and 13-15. Note that the absolute positional uncertainty of MM1, with which we compare positions of the detected masers, is ∼0.03 arcsec (Brogan et al. 2019). The absolute position shift between the first and second VLA epochs appeared to be of ∼0.04 arcsec. a The adopted rest frequency (with errors in the last digit); references: (1)   b The source is not considered further in the present study.

Methanol Maser Emission
Maser emission is detected in all frequency bands: there are not only well-known methanol masers at 6.67 and 12.18 GHz, but also rare and recently discovered ones. The four newly discovered methanol maser transitions at 6. 18,12.23,20.35,and 20.97 GHz,detected Tables 4-9. These tables, in machine-readable format, are presented in the Appendix. Images of the brightest emission spots detected at a particular frequency, as well as spectra of these maser transitions and their spatial distribution, are presented in Note that the structures seen in the maser spot maps are much smaller than the angular size of the VLA synthesized beam in each band (Table 1). The high signal-to-noise achieved for the bright maser emission allows us to fit the maser spot positions with sub-beamsize accuracy. Nevertheless, if there is more than one maser spot in a velocity channel, the spatial and velocity structure will be dominated by the brighter components (i.e., we obtain 'centroid mapping'; see Section 4.2 for the further discussion on the issue). This is especially true for the crowded maser regions detected at the flare epoch.
6.18 GHz. The 6.18 GHz class II methanol maser is one of the torsionally excited lines discovered with the ATCA during the burst in G358.93−0.03 (Breen et al. 2019). By the time of the VLA epoch II observations, the flux density of the 6.18 GHz maser had dropped significantly and remained at the ∼0.1 Jy level (Figure 2b), compared to the flaring flux density of ∼300 Jy. At the VLA epoch II, only the −18.6 km s −1 feature remained detectable. The ATCA spectrum showed emission from −18.8 to −15.5 km s −1 with a peak emission at −16.2 km s −1 . However, both the ATCA and VLA positions of the 6.18 GHz maser coincide within their positional uncertainties.   Table 4. 6.67 GHz. The first epoch of the VLA observations was conducted soon after the first 6.67 GHz maser flare in G358.93−0.03. On February 25, 2019 (VLA epoch I), the single-dish flux density of the velocity feature at ∼ −16 km s −1 had decreased from ∼600 Jy to ∼450 Jy, according to iMet data 6 . We note that February 25 was also the date when the dominant feature at V LSR ∼ −16 km s −1 was overtaken by the feature at V LSR ∼ −17 km s −1 .
The brightest flare of the 6.67 GHz maser occurred on March 12, when the ∼ −17 km s −1 velocity feature reached 900 Jy. By the second VLA epoch (June 4, 2019), the single-dish flux density of the 6.67 GHz methanol maser had decreased to ∼90 Jy, with both lines -at V LSR ∼ −16 km/s (dominant in the first flare) and ∼ −17 km/s (dominant in the second flare) -having about the same amplitude (see Fig. 1 and 3).
The methanol maser emission at 6.67 GHz is detected at both VLA epochs, coming from an area of ∼ 0.2 (∼1000 AU, for an assumed distance of 6.75 kpc) around the MM1 position (Fig. 3). Several fainter (<1 Jy) maser spots were detected outside of the main cluster. These "distant" maser spots are located to the west of the main cluster ( Fig.  3(a,b)). At the first VLA epoch, the cluster of the 6.67 GHz maser spots was elongated in the NE-SW direction, with most of the blue-shifted masers to the NE and the red-shifted features to the SW. Here, 'blue' and 'red' refer to the centroid of the velocity range, of about −17 km s −1 . In contrast, the 6.67 GHz maser spots detected at the second epoch are arranged in a bow-shaped structure, again with the blue-shifted velocity features to the north and red-shifted features to the south. The overall velocity pattern of the spectrum did not change greatly between epochs, although the weaker velocity features at V LSR ∼ −18 -−19.5 km s −1 detected in the first epoch were not detected at the second epoch.  Right panels: (c) superimposed 6.67 GHz methanol maser spot maps and spectra (the markers on the spectra correspond to the maser spots on the map) from the epochs I and II. Plots are color-coded by radial velocity (see colorbar for color scale). The error bars indicate the position fitting errors from Table 5. 12.18 GHz. The methanol maser emission at 12.18 GHz is found in a region of ∼0.2 around MM1 ( Figure 4). This maser transition shows a behaviour quite similar to that of the 6.67 GHz maser. At the VLA epoch I, the 12.18 GHz maser emission is also elongated in the NE-SW direction with blue-shifted velocity features to the north-east and red-shifted features to the south-west. At the VLA epoch II, there is a wide bow-shaped structure with red-shifted masers to the north and a less-ordered structure to the south, comprised of blue-shifted masers.
At the post-flare epoch (VLA epoch II), both maser transitions faded to about the same flux density with the 6.67 GHz maser slightly stronger: the peak flux densities of the 6.67 GHz and 12.18 GHz emission were ∼84 Jy and ∼79 Jy, respectively.     23.12 GHz. During the burst, the 23.12 GHz methanol masers associated with G358.93−0.03 showed the brightest and the most complex spectrum (MacLeod et al. 2019) of this rare maser transition detected to date (see Galván-Madrid et al. 2010 and references therein). The 23.12 GHz emission covers the same velocity range as of the 6.67 and 12.18 GHz masers, but with about 10× lower flux density. By the VLA epoch II, the 23.12 GHz maser had faded to ∼2 Jy, but continued to have a complex, multi-component spectrum.
Spatially, the 23.12 GHz maser emission is found around MM1 (Figure 7). At the flare epoch, the 23.12 GHz maser cluster appeared to be more compact than the 6.67/12.18 GHz masers, with a size of ∼0.15 . The flaring 23.12 GHz emission consists of two elongated sub-clusters, with the northern cluster hosting blue-shifted velocity features and the southern cluster with red-shifted features. At VLA epoch II, the 23.12 GHz emission is found in an expanded region of ∼0.2 (note that the lower flux densities at epoch II result in higher positional uncertainties, which could affect the apparent size of the region). The structure of the masing region at the post-flare epoch is less ordered, but shows the same north-south velocity gradient.

Continuum Emission
The VLA observations toward G358.93−0.03 detected two continuum sources ( Table 2). The hot core MM3 showed continuum emission in all three frequency bands at both epochs. Nevertheless, MM3 has no evident association with the bursting source and the detected methanol masers, and is not considered further in the present study. The bursting source MM1 was detected in K-band at the first epoch of the VLA observations. No K-band emission above the 3σ noise level was found at the second epoch.
In At the post-burst epoch, on April 12, 2019 in ALMA observations, the flux density was ∼282 mJy at 337 GHz (Brogan et al. 2019). Note, however, that the observations were made at different stages of the source's activity, and, therefore, the results of the observations cannot be directly compared.
The K-band continuum images of MM1 at both epochs are shown in Figure 8 and a summary of the parameters of the detected continuum peak is presented in Table 3. The offset between the ALMA and VLA positions of MM1 is ∼0.2 ( Fig. 9) which is similar to the absolute positional uncertainty reported in Brogan et al. (2019).  . Three of these, the maser transitions at 6.18, 12.23, and 20.35 GHz, were not even predicted by maser pumping calculations performed up to now. Others had been predicted in theoretical works (e.g., Sobolev et al. 1997b;Cragg et al. 2005), but had not been detected previously. The reasons for their rarity are not entirely clear.
Several of the newly detected methanol masers were imaged for the first time in the VLA observations presented here. The high sensitivity and moderate resolution of the VLA allow us to study the spatial structure of very weak masers, that are not accessible to VLBI imaging, but are readily detectable by the VLA. We were able to image and precisely locate positions of faint masers at 6.18 and 12.23 GHz with flux densities of ∼0.1 Jy (Figures 2 and 5). Notably, not only do the imaged methanol masers, including the newly discovered ones, arise in the same region, but they also trace the same structures ( Figure 10 Pre-flare VLA data from early 2012 exist and were published by Hu et al. (2016). In the quiescent state, the peak 6.67 GHz maser flux density was ∼5 Jy (Hu et al. 2016). The data of Hu et al. (2016) were obtained with the ∼ 4 beam of the VLA C-configuration, i.e. with much lower resolution than our VLA data, which precludes a comparison of the spot maps.
Our VLA data indicate that all of the flaring methanol masers originate from a single region of size ∼0.2 ( Figures  10 and 11) which coincides with the brightest mm source, MM1, detected with ALMA (Brogan et al. 2019). The close spatial association of the various flaring maser transitions is expected if caused by an accretion burst.
In order to analyse the maser density distribution at two VLA epochs and its evolution between the observations, we calculated the number of detected maser spots around MM1 -see histograms in Figure 11. The blue markers on the map and blue histogram bars correspond to the VLA epoch I, and the red markers/bars correspond to the epoch II. The bin width of the histograms is 0.02 arcsec which is one tenth of the area occupied by the maser emission and corresponds to the median position fitting error of our data (see Tables 4-9).   At the first VLA epoch, the flaring methanol maser emission consists of two linked sub-clusters, where the northern cluster hosts blue-shifted velocity spots while the southern cluster has red-shifted ones (Figure 10a). The analyses of the density of the maser spot distribution around MM1 (Figure 11), shows two density peaks located almost symmetrically around the source at the separations of ∼0.05 arcsec. One 'knot' is located to NE and another to SW of the MM1 position. Following the idea of Chen et al. (2020b), we interpret the VLA maser spot map obtained during the flare epoch for the 6.67, 12.18 and 23.12 class II methanol masers as a two-arm spiral pattern. The most crowded parts of the maser clusters could be explained as turning points of the gas motion in the arms (Meyer et al. 2018).
At the post-flare epoch, all detected methanol masers (including the newly-discovered and rare transitions) trace a bow-shaped structure extending eastward from the MM1 position (Figure 10b). While the velocity pattern persisted (the blue-shifted velocity spots are found to the north and the red-shifted to the south), the region affected by the maser emission expanded. The histograms in Figure 11 show that the maximums of the maser spot density shifted outwards to the radius of ∼0.09 arcsec from the position of the central source.
To estimate the evolution of the maser distribution in the region, we compare the positions of the maser emission density "knots" evaluated from the histograms in Figure 11. The density peaks were found at the diameter of ∼0.1 arcsec around MM1 at epoch I, and of ∼0.2 arcsec at epoch II. Thus, the region of physical conditions suitable to sustain maser emission doubled in size.
The drastic change between flare and post-flare epochs affected the entire masing region of ∼1000 AU (adopting the BeSSeL estimate of a near kinematic distance of 6.75 kpc) and happened in a period of about three months. Assuming these parameters, the triggering event may have propagated at a speed of ∼0.06c from the bursting source.  Such a high speed is almost certainly not caused by any physical movement of material. Hence, the flare and spatial rearrangement of the maser emission must be caused by a radiative event from the protostar that propagated through regions of more favorable physical conditions, more distant from the protostar.

Comparison of the VLA data with other observations
A thermal radiation "heatwave" emanating from an accreting high-mass protostar in G358.93−0.03, propagating at subluminal speed, has been inferred from multi-epoch LBA observations of 6.67 GHz methanol masers (Burns et al. 2020). The lower spatial resolution of the VLA prohibits such a precise analysis of the spatial structure of the maser emission in the region. On the other hand, it provides a unique insight into an extended component of maser emission.
The comparison of the VLA cross-correlation spectra (with a synthesized beam size of ∼3.5 arcsec, see Table 1) and single-dish spectra obtained with the Hitachi 32-m radio telescope as a part of the Ibaraki 6.7-GHz Methanol Maser Monitor (iMet) program shows that 100% of the single-dish flux density is recovered on the VLA baselines ( Figure  12). In contrast, only ∼10% of the maser flux density was recovered with the milli-arcsecond synthesized beam of the LBA (Burns et al. 2020). Based on these results and in accordance with the core/halo hypothesis of the structure of class II methanol maser masers (e.g. Minier et al. 2002), we conclude that with the VLA baselines we detect both compact (core) and extended (halo) emission while with the LBA baselines we detect only the compact core emission. The comparison of the 6.7 GHz methanol maser spot maps obtained in the VLBI observations of Burns et al. (2020) and in our VLA observations (Figure 13), shows that compact and extended components of the maser emission trace the same region around MM1, but highlight rather different structures. Both compact (VLBI data) and extended (VLA data) emission can be divided in two clusters: an upper (northern) cluster consisting of blue-shifted velocity spots and a lower (southern) cluster consisting of red-shifted velocity spots. However, the northern cluster of extended VLA maser emission is smaller than the southern cluster, while for the compact VLBI emission, the two clusters have about the same size. Also noteworthy is that both extended and compact maser emission trace roughly the same turning points in the both the northern and southern clusters.
Nevertheless, the LBA (Burns et al. 2020) and VLA observations were carried out with a three-day interval from each other during the period of high activity of the source. Therefore, it is possible that not all differences in the structure of the spot maps can be attributed to the extended/compact emission duality. For example, the most southern, linearly elongated maser cluster detected with the LBA has no counterpart in the VLA map. The iMet monitoring data indicated that the source shows daily/intraday variability. If we assume that a ∼0.05 arcsec cluster was ignited by the propagation of the heat wave over three days between the observations, the wave speed must be of ∼0.7c. On the one hand, a wave propagating through the low density material is supposed to have a speed of light, thus the speed of ∼0.7c is achievable, especially in the lower density outer regions. However, the absence of the elongated southern substructure is most likely due to the effect of centroid mapping. The most distant spot from MM1 with a particular velocity will be averaged in the synthesised beam with any other spots with the same velocity in the region, on average that will bring the centroid map position to be closer to center of the red cluster. It is also found that the cm-wavelength methanol masers detected with the VLA (this paper) and newly discovered high-frequency methanol masers detected with the SMA and ALMA (Brogan et al. 2019) trace similar spatial patterns (see Figure 14). Here we present a comparison of our VLA data with the 199.575 GHz methanol maser -the strongest (sub)millimeter maser detected in Brogan et al. (2019). Our analysis showed that the VLA epoch I data (February, 2019) resembles the SMA 199.575 GHz maser data obtained on March 14, 2019, while the VLA epoch II data (June, 2019) is similar to the spatial distribution of the 199.575 GHz methanol maser detected on April, 16 with the ALMA. In February-March, 2019, during the peak activity of the methanol maser flare (see Figure 1), both the common (VLA data) and (sub)millimeter, torsionally excited (SMA and ALMA data) methanol maser transitions are found in two sub-clusters separated by location and velocity, with the northern cluster containing the blue-shifted maser spots and the southern cluster predominated by the red-shifted velocity spots (Figure 14a). At the post-flare epoch, in April-June, 2019, all detected methanol masers trace a wider and less-dense/more-elongated formation, which can be interpreted as a spiral-arm structure (Figure 14b).
The cross-matching of our data and the results from Chen et al. (2020b) showed that the spatial distribution of the methanol masers detected at the second VLA epoch closely reassembles that which Chen et al. (2020b) reported for the 13 CH 3 OH, HDO, and HNCO masers in G358.93−0.03 (Figure 15b). Note that despite the fact that the observations   of Chen et al. (2020b) were conducted with the VLA on April 4, 2019, i.e. closer to our VLA epoch I (38 day apart -see Figure 15a), the data of the VLA epoch II (61 day apart) appeared to be more fitting (Figure 15b). Chen et al. (2020b) argue that the rare masers detected in their VLA observation trace a two-armed spiral structure; in this scenario the arms reveal accretion flows falling onto the central star in spiral trajectories. According to theory (e.g. Meyer et al. 2019a;Jankovic et al. 2019), spiral arms, as well as accretion burst, are products of gravitational instability of a disk around MYSOs. Our VLA data seems to suggest that the spiral-arm features of the accretion disk are traced not only by rare maser species detected by Chen et al. (2020b), but also by more common methanol masers at 6.67 and 12.18 GHz. Note that the observations of Chen et al. (2020b) were performed still at the flare epoch (April 4), while our VLA observations took place at the post-flare epoch (June 4) -see Figure 1. In the period between observations, no drastic changes in the structure of the spectrum were observed for the methanol maser at 6.7 GHz, for example; the only striking difference was the gradual decrease in the maser flux density from ∼500 Jy to ∼100 Jy 7 . The two-arm spiral infall gas flow structure seems to still exist at the post-flare period, and the physical conditions (e.g., temperature) in these flows change with the flare propagation, making different type of masers (with different pumping conditions) gradually arise or disappear in the flows. The resemblance of the maser clusters detected in Chen et al. (2020b) and in our observations at epoch II ( Figure 15b) indicates that the spatial distribution of the methanol masers remained roughly constant after April 2019. Note that the most recent VLBI/ALMA results obtained by the M2O indicate the presence of a more complex spatial structure of the regions traced by the methanol masers than that found with the VLA (Burns et al. in prep., Brogan et al. in prep). For example, VLBI observations suggest a 4-arm system (Burns et al. in prep.), not the two-arm system that was discovered with the VLA (Chen et al. (2020b), this paper). Nevertheless, the VLBI and VLA results do not contradict each other, but rather indicate the limitations of spatially unresolved maser spot fitting when there is complex underlying spatial morphology in each channel. It is possible that the two arms in the blue side of the disk are averaged within the VLA synthesised beam into one arm, and the same for the two arms on the red side of the disk. Since the disk is almost face-on, the velocity differences across the disk are small which worsens this effect as the centroid mapping technique derives a single spot per velocity channel and many maser spots exist in different parts of the disk (in different arms) at the same velocity. A more detailed comparison of the VLA and VLBI/ALMA results obtained for the source will be presented in the future M2O publications.

Cm-continuum emission of MM1
The strong millimeter continuum emission of MM1 Brogan et al. (2019), coupled with the lack of significant centimeter emission (Figure 8) indicates that the source is probably in a precursor state of ultracompact HII region (Churchwell 2002), as it was indicated in Brogan et al. (2019).
Our continuum observations suggest a moderate decrease in the flux densities of MM1 in the period between the flare and post-flare epochs. As we indicated in Section 3, the ALMA mm-observations and VLA cm-observations were performed during different epochs of the source activity. The analysis is particularly difficult because the presumptive bursting source MM1 was detected only in K-band and at one epoch of observation. Moreover, in contrast to S255 and NGC6334I, cm-continuum emission in G358.93−0.03 shows much lower flux densities.
Nevertheless, a decrease of the cm-continuum flux density of a source can happen during accretion event. In the case of NGC6334I, for example, Hunter et al. (2017) reported a 4-fold increase in the dust continuum emission, while the free-free emission (at 1.5 cm) fell by about the same factor (Brogan & Hunter 2018;Hunter 2019). Such a decrease would naturally occur if the ionizing photon flux from the young stellar object is reduced due to the bloating of the star. At the high densities of these MSFR, the recombination timescale for the ionized gas can be of the order of days to weeks. Hence, the free-free continuum emission can track the same accretion events, but show a decrease while the IR luminosity increases (Hunter 2019).

SUMMARY
A multi-epoch and multi-frequency VLA imaging of maser and continuum emission in C, Ku, and K-bands was performed for the massive star-forming region G358.93−0.03. Two observing sessions were conducted, at the maser flare epoch (VLA epoch I) and at a post-flare epoch (VLA epoch II). The main outcomes and scientific insights obtained from the observations are summarized as follows: (1) Maser emission is detected and imaged in several methanol transitions. Spatial structure evolution is studied for methanol masers at 6.67, 12.18, and 23.12 GHz at two observational epochs.
(2) The first interferometric images are obtained for the new methanol maser transitions at 6.18, 12.23, and 20.97 GHz (the masers were discovered in G358.93−0.03 in single-dish and ATCA observations Brogan et al. 2019)).
(3) Methanol maser emission in all detected transitions and at both epochs is found in a region with a diameter of ∼0.2 around the MM1 position and shows a velocity gradient in the NS direction.
(4) A drastic change in the spatial distribution of the detected methanol masers is found. At the flare epoch, the 6.67, 12.18, and 23.12 GHz methanol masers are found in elongated regions, aligned in a NE-SW direction. At the post-flare epoch, the methanol masers trace bow-shaped structures extending eastward from the MM1 position. During the transition from the first epoch to the second, the region traced by masers expanded while the velocity gradient decreased.
(5) The obtained data suggests that the methanol masers detected with the VLA trace the spiral-arm structures within the accretion disk which were first discovered in rare maser lines in Chen et al. (2020b).  Note-(1) Table 5 is published in its entirety in the machine-readable format. A portion is shown here for guidance regarding its form and content.
(2) The positional shifts of ∆RA=0.010, ∆DEC=0.079 (epoch I) and ∆RA=0.017, ∆DEC=0.017 (epoch II) were introduced to the data to prepare Figures 3 and 10-11, and 13-15 (see section 2). Note-(1) Table 6 is published in its entirety in the machine-readable format. A portion is shown here for guidance regarding its form and content.  Note-(1) Table 7 is published in its entirety in the machine-readable format. A portion is shown here for guidance regarding its form and content.
(2) The positional shift of ∆RA=-0.1, ∆DEC=-0.016 (epoch II) was introduced to the data to prepare Figures 5 and 10-11, and 13-15 (see section 2). Note-(1) Table 8 is published in its entirety in the machine-readable format. A portion is shown here for guidance regarding its form and content.
(2) The positional shift of ∆RA=-0.055, ∆DEC=0.022 (epoch II) was introduced to the data to prepare Figures 6 and 10-11, and 13-15 (see section 2). Note-(1) Table 9 is published in its entirety in the machine-readable format. A portion is shown here for guidance regarding its form and content.