M17 MIR: A Massive Protostar with Multiple Accretion Outbursts ¹¹1Based on observations collected at the European Southern Observatory, Chile, Programs 073.C-0170(A+B), and 077.C-0174(A).

The wide-field $JHK_{s}$- and $L^{\prime}$-band images were obtained with ISAAC at the Very Large Telescope (VLT) in 2002 September and 2004 May to September, respectively. The pixel scale is $0\farcs 148$ and $0\farcs 07$, and the FWHM is $0\farcs 6$ and $0\farcs 5$ for the $JHK_{s}$- and $L^{\prime}$-band images, respectively (see more details of observations in Hoffmeister et al., 2008). These $JHK_{s}L^{\prime}$-band images and photometric data were published by Hoffmeister et al. (2008).

The MIR images were obtained toward the M17 SW photodissociation region (PDR) (e.g., Crété et al., 1999; Pérez-Beaupuits et al., 2012) with the VLT imager and spectrometer for mid-infrared (VISIR; Lagage et al., 2004) at the VLT in 2006 April through the SiC filter ($\lambda=11.8\,\mu$m, $\Delta\lambda=2.1\,\mu$m). The data were taken with the VISIR_img_obs_AutoChopNod template. The photometric standard star was HD 178345 during this observation. The data were reduced and calibrated using the EsoReflex workflows (Freudling et al., 2013) and the VISIR pipeline recipes 4.3.7²²2http://www.eso.org/sci/software/pipelines/. The image is of good quality with FWHM $\approx 0\farcs 4$, reaching the diffraction limit.

The archival imaging data of ISO (Kessler et al., 1996) were obtained for the M17 SW PDR from the ISO Data Archive⁵⁵5https://www.cosmos.esa.int/web/iso/access-the-archive. The ISO observations for the M17 SW PDR were carried out in 1997 March and published in Crété et al. (1999); they aim at the interstellar medium inside the PDR. The calibrated ISO/ISOCAM images in the three broad filters (LW1, LW4, and LW6) with central wavelengths of 4.5, 6.0, and 7.7 $\mu$m are used in this paper.

The Galactic Legacy Infrared Mid-Plane Survey Extraordinaire (GLIMPSE; Benjamin et al., 2003; Churchwell et al., 2009) is a legacy science program of the Spitzer Space Telescope that surveyed the inner Galactic plane in the 3.6, 4.5, 5.8 and $8.0\,\mu$m bands with the IR camera IRAC. The angular resolution in the four bands (I1,I2,I3,and I4) is about $2\farcs 5$. The GLIMPSE mosaicked images of M17 were obtained from the NASA/IPAC Infrared Science Archive⁶⁶6http://irsa.ipac.caltech.edu/frontpage/; the median observation date was 2005 September. The GLIMPSE point sources were retrieved from the GLIMPSE I Spring 07 Archive⁷⁷7See more descriptions in http://irsa.ipac.caltech.edu/Missions/spitzer.html.

M17 was revisited by Spizter in the warm phase in 2017 July. During these observations, only the two bluest bands ($3.6$ and $4.5\,\mu$m) were available, and the flux-calibrated images (program ID: 13117, observer: Robert Benjamin) were obtained from the Spitzer Heritage Archive⁸⁸8Hosted at https://irsa.ipac.caltech.edu/Missions/spitzer.html.

In order to cover the far-IR regime, we incorporated the publicly available imaging observations from the Herschel Space Observatory (Herschel) and its PACS instrument at 70, 100, and 160 $\mu$m (Poglitsch et al., 2010). The Herschel data used in this work are Level 3.0 products (eight observations with IDs from 1342192767 to 1342192774) provided by the Herschel Science Archive⁹⁹9http://archives.esac.esa.int/hsa/whsa/. The pixel scales are $1\farcs 6$, $1\farcs 6$, and $3\farcs 2$ at 70, 100, and 160 $\mu$m, and the measured resolutions are $8\farcs 0$, $8\farcs 0$, and $12\farcs 0$, respectively.

The Wide-field Infrared Survey Explorer (WISE, Wright et al., 2010) is a 40 cm telescope in a low earth orbit that surveyed the entire sky in 2010 using four IR bands at 3.4, 4.6, 12, and 22 $\mu$m. The angular resolution in the four bands (W1, W2, W3, and W4) are $6\farcs 1$, $6\farcs 4$, $6\farcs 5$, and $12\farcs 0$, respectively. Only the WISE W1 and W2 data of M17 are considered in this work, because the M17 H II region is saturated in the W3 and W4 bands. With the primary aim of studying near-Earth Objects, NEOWISE observations resumed in 2013 December and continue to date, after the telescope’s cryogen tanks were depleted. Due to the orbit of the telescope, WISE/NEOWISE observations visit M17 twice a year, in March and September. From the NEOWISE 2020 data release containing W1 and W2 observations from 2013 December until 2019 December, the single-exposure data of M17 collected over 6 yr (12 epochs) are analyzed in this work, in addition to the WISE all-sky survey image of M17 (two epochs).

The single-exposure images of the WISE/NEOWISE mission toward M17 were retrieved from the WISE data hosted on IRSA ¹⁰¹⁰10https://irsa.ipac.caltech.edu/Missions/wise.html, with the service of WISE/NEOWISE Coadder. The single exposures were coadded to obtain higher-quality W1/W2 images with enhanced resolution of $0\farcs 6875/\mathrm{pix}$. Finally WISE/NEOWISE images from 14 epochs of M17 were used in this work. Aperture photometry with a fixed radius of 8.5 pixel ($=5\farcs 8$) is applied to all WISE/NEOWISE images with a field of view of $3\arcmin\times 3\arcmin$. Because the WISE/NEOWISE W2 band is very similar to the Spitzer I2 band, the photometric results from multi epochs are checked against the isolated star G015.0198-00.6768 located $33\arcsec$ south of the image center. G015.0198-00.6768 is not cataloged in the ALLWISE catalog (Cutri et al., 2021) or in the unWISE catalog (Schlafly et al., 2019), but has a cataloged Spitzer I2 magnitude of 8.55 (Spitzer Science, 2009).

The SOFIA¹¹¹¹11SOFIA is jointly operated by the Universities Space Research Association, Inc. (USRA), under NASA contract NAS2-97001, and the Deutsches SOFIA Institut (DSI) under DLR contract 50 OK 0901 to the University of Stuttgart. (Young et al., 2012) observations toward the M17 SW PDR were carried out with the FORCAST instrument (Herter et al., 2013) at $19\,\mu$m ($\lambda_{\mathrm{eff}}=19.7\,\mu$m; $\Delta\lambda=5.5\,\mu$m) and $37\,\mu$m ($\lambda_{\mathrm{eff}}=37.1\,\mu$m; $\Delta\lambda=5.5\,\mu$m) on 2017 August 2 (PI: James M. De Buizer). The calibrated Level 3 images at 19 and $37\,\mu$m were downloaded from the SOFIA Data Cycle System ¹²¹²12https://dcs.arc.nasa.gov/. The resolutions at 19 and $37\,\mu$m are $3\farcs 3$ and $3\farcs 5$, respectively (De Buizer et al., 2017). The SOFIA/FORCAST images fill the gap between the VLT/VISIR and Herschel images in wavelength coverage.

The James Clerk Maxwell Telescope (JCMT) SCUBA2 submillimeter observations were taken at $450\,\mu$m and $850\,\mu$m on 2016 April 10 for M17, as a part of the JCMT large program “SCOPE: SCUBA-2 Continuum Observations of Pre-protostellar Evolution” (JCMT program code: M16AL003). The pipeline-reduced and calibrated maps of M17 at $450\,\mu$m and $850\,\mu$m were retrieved from Canadian Astronomy Data Centre ¹³¹³13http://www.cadc-ccda.hia-iha.nrc-cnrc.gc.ca/en/jcmt/.

The earliest detection of M17 MIR was at $3.9\,\mu$m $L$-band by the Canada–France–Hawaii Telescope (CFHT) observations in 1992 December and 1993 May with a spatial resolution of $1\arcsec$ (Giard et al., 1994). This object, named as source 2 by the authors, and six other extremely red objects are suggested to be young stars in very early evolutionary stages (see Table 1 in Giard et al., 1994). Figure 1 presents the multiwavelength images centered on M17 MIR. A summary of the IR flux densities obtained by the various facilities except for the earliest CFHT observation is presented in Table 1.

The source sizes in the various IR images are in proportional to FWHMs of the corresponding images. In the $11.8\micron$ image of $0\farcs 4$ resolution ($=800\,\mathrm{au}$), M17 MIR is still a single source. M17 MIR appears as a compact point source at $3-37\,\micron$. At Herschel $70\micron$ and $160\micron$, the far-IR emission is dominated by the contamination from M17 IRS5 and M17 UC1 to the northeast. M17 MIR is not resolved at its IR location of M17 MIR, likely due to the strong ambient far-IR emission. For M17 MIR, we estimated an upper limit $<77\,\mathrm{Jy}$ at $70\,\micron$.

The faint IR brightness of M17 MIR suggests that it is a low-luminosity object or the extinction to M17 MIR is quite high. From a quick examination of the IR images of M17 MIR, we think high extinction is the more likely reason. M17 MIR’s brightness shows a steep rise in the $3-8\,\micron$ range (see also Figure 8). However, M17 MIR is not visible in the ISO broad filter image at $9.6\,\micron$ (see Figure 2a–f; Crété et al., 1999). We know that the silicate absorption feature at $9.7\,\micron$ can cause high and broad absorption to a highly obscured object at around this wavelength. The high extinction to M17 MIR can explain the nondetection at $9.6\,\micron$.

The IR flux densities of M17 MIR are at least two orders of magnitude lower than those of the IR-bright high-mass protostar IRS5A (Chen et al., 2015) within the same M17 SW PDR. On the other hand, the outwards motion of the CJO16 H₂O masers very likely traces the shocked gas accompanying the protostellar outflow driven by M17 MIR. In the earliest stage of high-mass protostars, the mass loss might evolve from a highly uncollimated expanding shell to a collimated outflow (e.g., Carrasco-González et al., 2015). The dynamical age of the expanding shell driven by M17 MIR is 12.5 yr (CJO16), younger than the dynamical age ($\sim 25$ yr) of W75N(B)-VLA 2 (Carrasco-González et al., 2015). The faint IR brightness of M17 MIR and the young dynamical age of the CJO16 H₂O masers suggest that M17 MIR is a promising candidate for a high-mass protostar in the earliest stage. The analysis of the spectral energy distribution (SED) constructed from Table 1 will help to better constrain the fundamental properties of M17 MIR (Sect. 5.3).

The warm area close to the central object emits radiation mostly at MIR wavelengths, while the surrounding envelope of much larger size is expected to radiate at submillimeter to millimeter wavelengths due to a much lower effective temperature. Figure 2 presents maps of the $850\mu$m dust continuum emission and the HCO+ $J=9\rightarrow 8$ line integrated intensity of M17 SW PDR. M17 MIR, represented by the H₂O maser source 2 detected by VLA observations (Johnson et al., 1998, hereafter JDG98), is located exactly in the secondary peak of the 850 $\mu$m dust continuum emission. This dust continuum emission peak, along with the other two peaks, make up the northern condensation in M17 SW PDR that was first noticed by Stutzki & Guesten (1990). The three main components in the northern condensation in M17 SW PDR were studied in detail through submillimeter and millimeter observations. Their physical properties were estimated by Hobson et al. (1993) (their Table 4). M17 MIR is located within the FIR3 source in Hobson et al. (1993). On the other hand, the map of integrated line intensity (range: $0-40\,\mathrm{km\,s^{-1}}$, peak velocity $\sim+20\,\mathrm{km\,s^{-1}}$) of the HCO+ $J=9\rightarrow 8$ line shows clumpy structures similar to the dust continuum emission. A dense core whose size is comparable to the beam size of the observations ($\sim 8\arcsec$, or 0.08 pc) is coinciding with M17 MIR. Given the critical density $n_{\mathrm{crit}}=9\times 10^{7}\,\mathrm{cm}^{-3}$ and an excitation temperature $E_{\mathrm{up}}=192.58\,\mathrm{K}$ of the $802.458\,\mathrm{GHz}$ HCO+ $J=9\rightarrow 8$ line (Pérez-Beaupuits et al., 2015), the dense core associated with M17 MIR is likely a candidate for a hot molecular core (HMC; Cesaroni, 2005). HMCs commonly corresponds to a hydrogen column density $N_{\mathrm{H}}\gtrsim 10^{23}$ or even higher (e.g. Bonfand et al., 2017). Although a high-mass protostar embedded in an HMC is undetectable at short wavelengths due to high extinction, the detection of the central protostar becomes possible at MIR wavelengths. For instance, the IR counterpart F4 of the HMC G9.62+0.19-F appears at $3.8\,\mu$m in the $L^{\prime}$ band and at longer wavelengths (Linz et al., 2005).

The multiwavelength data of M17 MIR and the submillimeter to millimeter data of its surrounding medium indicate that M17 MIR is a very early high-mass protostar embedded in an HMC and is also driving a protostellar outflow.

For M17 MIR, the photometric magnitude in the 1993 CFHT observation was $10.2\pm 0.2$ (source 2 in Table 2; Giard et al., 1994), while the magnitude in the 2004 VLT observation was $12.78\pm 0.14$ (Hoffmeister et al., 2008). The filter profiles in the two observations are almost identical. Table 2 lists the $L$-band magnitudes of the embedded YSOs in the two observations. Interestingly, M17 MIR is about 10 times fainter in the 2004.09 epoch than in the 1993.03 epoch, whereas the other sources are fairly stable.

The VLT $L^{\prime}$-band observations in 2004.09 were one year earlier than the Spitzer observations in 2005.09. By interpolating the Spitzer flux densities at I2 and I3 to that at $L^{\prime}$ under the assumption of a blackbody radiation, the flux density of M17 MIR in the $L^{\prime}$ band is about 2.2 mJy in the 2005.09 epoch. This predicted flux density in 2005.09 is close to the $1.91\pm 0.17$ mJy measured in 2004.09. We therefore speculate that M17 MIR’s luminosity is stable during 2004.09–2005.09.

M17 MIR shows variabilities in the Spitzer I1 and I2 images from 2005.09 and 2017.07 (see Figure 3). Comparing the flux densities at I1 in the two epochs for an aperture radius of $1\farcs 6$, we find an increase by a factor of 2.9. Similarly, the flux densities at I2 show an increase by a factor of 3.3 in an aperture radius of $3\farcs 2$. In contrast, the other sources are stable in the two bands (see Table 3).

The MIR variabilities of M17 MIR found for different epochs motivate us to compare its flux densities from the ISO observations in 1997.03 and the Spitzer observations in 2005.09 and 2017.07. Table 4 shows the flux densities of M17 MIR and other nonvariable sources in the ISO/LW1 and Spitzer/I2 bands, which have almost identical filter profiles centered at $4.5\,\mu$m. Figure 4 compares the ISO LW1 and Spitzer I2 images obtained in 1997.03 and 2005.09, respectively. Direct comparison from Figure 4 clearly indicates that M17 MIR is significantly brighter than sources 3 and 5 in the ISO LW1 image in the 1997.03 epoch, while it becomes fainter than sources 3 and 5 in the Spitzer I2 image in the 2005.09 epoch. Given that sources 3 and 5 are unchanged (see context above), M17 MIR is varying between 1997.03 and 2005.09. The aperture photometry of M17 MIR in the $4.5\,\mu$m band shows variability from 1997 to 2017.

For M17, WISE/NEOWISE images in 14 epochs, including 2010.03 and 2010.09, and March and September in 2014 to 2019, were analyzed in this work. M17 MIR is only detected in the W2 band. The WISE W1 flux densities, extrapolated from its flux densities in the Spitzer I1 and I2 bands, are about 0.6 and 1.8 mJy during the 2005.09 and 2017.07 epochs, respectively. The extreme faintness at W1 is consistent with the nondetection of M17 MIR in the corresponding W1 images.

Figure 5 presents the WISE/NEOWISE W2 images centered on M17 MIR in 14 epochs. The pixel resolution of $0\farcs 6875$ clearly separates M17 MIR from the nearby IR source Anon 3 to the southwest. The photometric results of M17 MIR, source 7, and G015.0198-00.6768 are tabulated in Table 5 for 14 epochs. In 2010.03, the W2 flux density of M17 MIR was $17.7\pm 4.4$ mJy, which is very close to the value $18.4\pm 2.0$ in the Spitzer I2 band in 2005.09. Given the fact that the WISE/NEOWISE W2 band is very similar to the Spitzer I2 band, M17 MIR was roughly stable in brightness from 2004 to 2010. However, six months later M17 MIR was twice as bright in 2010.09, with a flux density of $33.7\pm 4.4$ mJy. It is very likely that the rebrightening of M17 MIR started between 2010.03 and 2010.09. In the observations of NEOWISE, M17 MIR generally shows a steady rise in brightness, but with a jump and fall motion epoch by epoch in its light curve (see Figure 6). Between 2016 and 2017, M17 MIR underwent a steep rise in brightness with an amplitude of 1.7. In 2018 M17 MIR’s brightness dropped back to the level of 2016. and was followed by a rise with a mean amplitude of 1.5 between 2018 and 2019.

The MIR flux measurements of M17 MIR span several wavelengths and observation epochs. The $3.8\,\mu$m $L^{\prime}$ filter is closest to the $L$ filter of the 1993.03 epoch observation (Giard et al., 1994), while the observations of other epochs provided flux densities in at least two filters. For a direct comparison from 1993 to 2017, the flux density in the $3.8\,\mu$m $L^{\prime}$ filter is estimated for 1993.03, 1997.03, 2005.09, and 2017.07. The multiwavelength flux densities of M17 MIR in 1997.03, 2005.09, and 2017.07 epochs yield a reasonable estimate for the continuum in the range $3-8\,\mu$m, assuming blackbody emission. The flux density of M17 MIR in the $3.8\,\mu$m $L^{\prime}$ filter is the continuum emission inside the transmission curve of the $L^{\prime}$ filter. Note that the blackbody temperatures for the three different epochs are different. For the $L$-band observations in 1993, the three blackbody temperatures are used for estimating the $L^{\prime}$ flux density of M17 MIR. We find that the largest difference among the three values is less than 15%, thus the mean value is used as the $L^{\prime}$ filter flux density. For 1993.03, 1997.03, 2005.09, and 2017.07 we obtain $L^{\prime}$ flux densities; for 2004.09 there is direct $L^{\prime}$-band photometry for M17 MIR. Figure 6 presents the combined light curve of M17 MIR in the $L^{\prime}$ and W2 bands. The time span of the multiepoch data is 26.5 yr.

M17 MIR shows decreasing, quiescent and rebrightening phases in its MIR light curve, as shown in Figure 6. Due to the lack of $3-4\,\mu$m observations before 1993 and during 1997–2004, the duration of the decreasing phase is poorly constrained. From the available data, the quiescent phase between the decreasing and rebrightening phase might last about 6 yr, from mid-2004 to mid-2010. Because M17 MIR seems to continue the trend of rebrightening, the timescale of the rebrightening phase is $\gtrsim 9$ yr. The continuing NEOWISE observations will provide vital information on the rebrightening timescale.

Table 6 summarizes the 22 GHz H₂O maser observations from 1984 to 2008. The earliest H₂O maser detection associated with M17 MIR was achieved at the VLA in C configuration in 1984 June (Forster & Caswell, 1989, hereafter FC89). The complete catalog of H₂O masers found by Forster & Caswell (1989) was compiled by Forster & Caswell (1999). Three FC89 H₂O masers are associated with M17 MIR, and show two velocity components of $16-18\,\mathrm{km\,s^{-1}}$ and $23.3\,\mathrm{km\,s^{-1}}$. One H₂O maser coinciding with M17 MIR was detected with the VLA in D configuration and shows a flux density of 18.57 Jy (Massi et al., 1988, hereafter MCF88). One JDG98 H₂O maser (source 2 in JDG98) with two velocity components is detected around the same source. The velocities of the FC89 H₂O masers are consistent with those of the two components of the JDG98 H₂O maser source 2. Both the FC89 and JDG98 H₂O masers are likely pumped by the interaction between the outflow from M17 MIR and the surrounding gas. Interestingly, the flux density of the JDG98 H₂O maser source 2 in the 1996.09 epoch is ten times higher than that of the FC89 H₂O maser in 1984.06.

Several years after the last VLA observations, three water maser observations were taken by ATCA during three epochs. A fast fading of the water maser intensity from 197 Jy in 2003.10 (Breen et al., 2010, hereafter B10), through 67 Jy in 2004.07 (B10), to 4.2 Jy in 2008.08 (Breen & Ellingsen, 2011, hereafter BE11) was observed. In the earlier two epochs, the spectra of water maser emission show two components– a stronger blueshifted one at $\sim 19-20\mathrm{km\,s^{-1}}$ and a weaker redshifted one at $24\mathrm{km\,s^{-1}}$, which are consistent with the VLA observations in the earlier epochs. In the 2008.08 epoch, two components were observed in the spectrum – a weaker blueshifted one at $\sim 18\,\mathrm{km\,s^{-1}}$ and a stronger redshifted one at $21\,\mathrm{km\,s^{-1}}$.

The latest water maser observations were conducted by very long baseline interferometry (VLBI) with VLBI Exploration of Radio Astrometry (VERA) in 2009.10–2010.12 and consist of six epochs (CJO16). The VERA observations detected 203 water maser spots in these observations. The dominant blueshifted maser spots ($V_{lsr}\lesssim\,22\,\mathrm{km\,s^{-1}}$) are mostly distributed to the south of M17 MIR, while a few redshifted maser spots ($V_{lsr}\gtrsim\,22\,\mathrm{km\,s^{-1}}$) are distributed north to the north of it (see Figure 7). This feature is consistent with the stronger blue-shifted component and the weaker red-shifted component observed by the VLA and ATCA.

The IR variability of M17 MIR overlaps with the surrounding H₂O maser variability in time. The connection between them helps to clarify the physical origin of these variabilities. The flux density of the JDG98 H₂O maser source 2 in 1996.09 is ten times higher than that of the FC89 H₂O maser source in 1984.06. This H₂O maser variation from 1984 to 1996 and the decreasing phase in the IR light curve during 1993 and 2004 suggest that the brightness peak likely occurred between 1984 and 1993. Moreover, the fast fading of the H₂O maser emission from 2003 to 2008 and the quiescent phase in the IR light curve during 2004 and 2010 together constrain the occurrence of the quiescent phase to between mid-2004 and mid-2010, characterized by both low IR brightness and H₂O maser emission.

VERA observations in the earliest two epochs (2009.10 and 2010.01) were made in the quiescent phase, while the later observations in the epochs 2010.04, 2010.09, 2010.11, and 2010.12 took place in the rebrightening phase. The proper motions of the H₂O masers spots are traced in more than three epochs (J. O. Chibueze 2021, private communication). The proper motions of the H₂O maser spots are measured in the transient phase between quiescent and rebrightening phases from 2009.10 to 2010.12. The 3D motions of the H₂O maser spots detected by the VERA observations show an expanding bubble originating from M17 MIR. Together with the rebrightening, we suggest that both phenomena are due to the accretion/outflow activity of M17 MIR.

Similar flaring of the H₂O masers associated with the outflow from the high-mass protostar NGC 6334I-MM1B shows a mean factor of 6.5 before and after the burst of MM1B in the submillimeter range (Brogan et al., 2018). Those authors attribute the enhanced accretion rate to the coincident flaring in the submillimeter continuum and the H₂O maser emission of NGC 6334I–MM1B. For M17 MIR, we suggest an accretion outburst to explain both the MIR variations and the correlated H₂O maser variations. Accompanied by the expanding motions of H₂O maser spots measured during the transient phase from quiescent to rebrightening and the close positive relation between mass loss and accretion rate, the most promising explanation for the outburst in MIR brightness since 2010.03 is an enhanced accretion rate.

The proper motions and radial velocities of the CJO16 H₂O masers are overlaid on the VLT $L^{\prime}$ image, as shown in Figure 7. The 3D motions of the H₂O maser spots can be used to estimate the momentum rate of the gas flow, given the assumption that all the momentum of the ejected gas is transferred to the surrounding molecular environment. The momentum rate $\dot{P}$ is then estimated from the relation (Goddi et al., 2011)

where $V_{10}$ is the mean maser velocity in units of $10\,\mathrm{km\,s^{-1}}$, $R_{100}$ is the average distance of H₂O masers from the driving source in units of 100 au, $\Omega$ is the solid angle of the gas flow, and $n_{8}$ is the gas volume density in units of $10^{8}\,\mathrm{cm^{-3}}$. $R_{100}$ is calculated as 0.85 by averaging the distances of the H₂O masers from the dynamical center, which is determined at the position of $\mathrm{R.A.(J2000)}=18^{\mathrm{h}}20^{\mathrm{m}}23\fs 017,\mathrm{decl.(J2000)}=-16\arcdeg 11\arcmin 47\farcs 98$ (CJO16). The gas density is about $10^{8}\,\mathrm{cm^{-3}}$ according to the discussion in Sect. 3.1. $\Omega$ is about $4\pi$ for an expanding bubble. The above physical properties from observations yield a momentum rate $\dot{P}\sim 4\times 10^{-3}\,M_{\odot}\,\mathrm{yr^{-1}\,km\,s^{-1}}$ for the protostellar outflow of M17 MIR. This is comparable to outflows of $10^{-3}-10^{0}\,M_{\odot}\,\mathrm{yr^{-1}\,km\,s^{-1}}$ driven by MYSOs (Moscadelli et al., 2016), and one magnitude higher than the outflow from a low-mass YSO in IRAS 20231+3440 (Ogbodo et al., 2017).

The protostellar outflow provides an indirect measure of the disk accretion rate. The outflow mass rate is derived by dividing the outflow momentum rate $\dot{P}$ by the outflow velocity. Given $\dot{P}\sim 4\times 10^{-3}\,M_{\odot}\,\mathrm{yr^{-1}\,km\,s^{-1}}$ estimated in this work and $V_{10}\approx 1.9\,\mathrm{km\,s^{-1}}$ determined by CJO16, the outflow mass rate $\dot{M}_{out}$ is $\sim 2\times 10^{-4}\,M_{\odot}\,\mathrm{yr^{-1}}$. Assuming a ratio of outflow mass rate to disk accretion rate of 30% (Beltrán & de Wit, 2016), the disk accretion rate onto M17 MIR is $\dot{M}_{\mathrm{acc}}\sim 7\times 10^{-4}\,M_{\odot}\,\mathrm{yr^{-1}}$.

This disk accretion rate is estimated from the H₂O maser observations during 2009 to 2010, when M17 MIR just started to rebrighten. The MIR emission between 2010.09 and 2010.03 varies by a factor of about 2 in the WISE W2 band. Johnstone et al. (2013) simulate the MIR flux variation induced by an enhanced disk accretion rate and find that the MIR flux variation is proportional to the disk accretion rate. The W2 flux between 2010 and 2019 varies by an average factor of 3. We speculate that the disk accretion rate in 2019 is likely 3 times the rate estimated above, reaching $\sim 2\times 10^{-3}\,M_{\odot}\,\mathrm{yr^{-1}}$.

During the quiescent phase from mid-2004 to mid-2010 (hereafter the 2005 epoch), we collect IR flux densities at six wavelengths to construct the SED of M17 MIR. During the rebrightening phase, the flux densities obtained in 2017.07 (Spitzer I1 and I2 bands) and 2017.08 (SOFIA/FORCAST $19\,\mu$m and $37\,\mu$m bands) are separated by 14 days. The NEOWISE W2 flux density was obtained at the end of 2017 September, nearly two months later than the Spitzer and SOFIA observations. We collect only the Spitzer and SOFIA observations in 2017.07/08 to construct the SED of M17 MIR during the outburst phase (hereafter the 2017 epoch). The flux densities of M17 MIR in the 2005 and 2017 epochs are tabulated in Table 7. The $3\sigma$ upper limit of the Herschel $70\mu$m image is 77 Jy for M17 MIR, which is considered as the upper limit in both the 2005 and 2017 epochs. The additional input parameters besides the flux densities are distance and foreground extinction $A_{V}$. The distance to M17 MIR is fixed to 2.0 kpc. We assume a wide range of 10–200 mag for the foreground extinction $A_{V}$.

We employed the SED model grid for YSOs that spans a wide range of evolutionary stages, from the youngest deeply embedded protostars to pre-main-sequence stars with debris or no disk, and which covers a wide and uniform region of parameter space (Robitaille, 2017a, hereafter R17). We fit the SED of M17 MIR using several R17 model sets of different model components, using the Python-based fitting tool (1.3 version; Robitaille, 2017b). Based on the prior knowledge of M17 MIR, it is a highly embedded massive protostar with an accretion disk. The model sets we used are sp-s-i,sp-smi,spu-smi,spubsmi. All four model sets have a passive disk with an inner radius equal to the dust sublimation radius. sp-s-i is the simplest model set including only the central star and passive disk. The other three model sets include ambient gas (sp-smi), envelope+ambient gas (spu-smi), and envelope+ambient+cavity (spubsmi).

Fitting the SED of M17 MIR to the individual models of the four model sets returns a $\chi^{2}$ value for each individual model. The criterion $\chi^{2}-\chi_{\mathrm{best}}^{2}<3n_{\mathrm{data}}$ can generally split good and bad fits. We used this criterion to select good fits from the four model sets. Among these good fits, we use the following criteria to pick reasonable models.

Disk mass $\geqslant\,0.01M_{\sun}$. The disk accretion of M17 MIR estimated from H₂O maser motions is $7\times 10^{-4}\,M_{\sun}\,\mathrm{yr}^{-1}$. To sustain such a high accretion rate, we require this criterion.

Disk inclination angle $\geqslant 30\arcdeg$. Because the blue- and redshifted H₂O maser spots are distinctly separated in the vicinity of M17 MIR, the disk inclination angle should not be small.

Exclude unphysical models. Some models are unphysical with too small radii but too high temperatures. We compute the Stefan–Boltzmann luminosity ($4\pi R_{\ast}^{2}\sigma T_{\ast}^{4}$) from the radius and temperature of each model. These unphysical models are below the ZAMS track in the H-R diagram. We keep only the models that are lying on/above the ZAMS track in the H-R diagram.

The above criteria for the four model sets return reasonably good fits for the SEDs of M17 MIR in the 2005 and 2017 epochs. Because of the degeneracy of the four model sets, one has to score these sets not only based on $\chi^{2}$ values. In order to assess the model sets in comparison with each other in a statistically robust way, Robitaille (2017a) suggested calculating $P(D|M)\varpropto N_{\mathrm{good}}/N$, where $N_{\mathrm{good}}$ is the total number of models from a given model set, and $N$ is the total number of models in that set. Towner et al. (2019) used the R17 model sets to find the best YSO models for 12 MYSOs based on the overall assessment of $\chi^{2}$ and $P(D|M)$ scores. For the good fits of M17 MIR, we show the $\chi_{\mathrm{best}}^{2}$ and $P(D|M)$ scores of the four model sets in Table 8. In the 2005 and 2017 epochs, sp-s-i is the ”best representation” model set according to both the $\chi_{best}^{2}$ and $P(D|M)$ scores.

Among the sp-s-i model set, 44 and 13 models remain for the SEDs of M17 MIR in the 2005 and 2017 epochs, respectively. In the 2017 epoch, the SOFIA observations at $19.7\micron$ and $37.1\micron$ are very useful in distinguishing the various models of a given model set (Figure 8). Among the 13 models in the 2017 epoch (Table 9), two models correctly reproduce the SED slope between $19.7\micron$ and $37.1\micron$. Indeed, the $\chi^{2}_{\mathrm{best}}$ model is different from the second fit only in the disk inclination angle. Aided by the observations at $19.7\micron$ and $37.1\micron$ in the 2017 epoch, the $\chi^{2}_{\mathrm{best}}$ model is distinctly better than the other sp-s-i models.

In the 2005 epoch, the 44 sp-s-i models show a degeneracy in the long-wavelength SED beyond $12\micron$, due to the lack of observations at long wavelengths. Povich et al. (2011) found similar SED degeneracy at long wavelengths for a deeply embedded YSO with hard X-ray emission (PCYC source 699 from Povich et al. (2011)) in the Carina nebula region, when they fitted the SED of this object to YSO model SEDs from Robitaille et al. (2006). In Figure 9 of Povich et al. (2011), the distributions of the well-fit YSO models show two distinct groups: intermediate-mass YSOs with relatively high circumstellar extinction and low-mass YSOs with lower circumstellar extinction. The group with high circumstellar extinction is consistent with the high absorption of the X-ray source (Povich et al., 2011). For the 44 sp-s-i models of M17 MIR in the 2005 epoch, two distinct families are found in the plot of infrared luminosity integrated between 12 and $500\micron$ versus disk maximum radius $R_{\mathrm{max}}^{\mathrm{disk}}$: high infrared luminosity with $R_{\mathrm{max}}^{\mathrm{disk}}>1000\,\mathrm{au}$ and lower infrared luminosity with $R_{\mathrm{max}}^{\mathrm{disk}}<1000\,\mathrm{au}$. The $R_{\mathrm{max}}^{\mathrm{disk}}>1000\,\mathrm{au}$ family is more consistent with the general understanding of MYSOs, which are luminous at long infrared wavelengths and have large circumstellar disks. Moreover, the $\chi^{2}_{\mathrm{best}}$ sp-s-i model in the 2017 epoch also suggests a large circumstellar disk for M17 MIR. The $R_{\mathrm{max}}^{\mathrm{disk}}>1000\,\mathrm{au}$ family of the sp-s-i models is preferred for the SED of M17 MIR in the 2005 epoch. Among the nine sp-s-i models with $R_{\mathrm{max}}^{\mathrm{disk}}>1000\,\mathrm{au}$ in Table 9, two models with $T_{eff}=24,540\,\mathrm{K}$ likely produce significant 22 GHz continuum emission that could be detected by the ATCA observations during 2003–2008 (B10,BE11). However, no 22 GHz continuum source was reported around M17 MIR according to the ATCA observations. We excluded the two very hot sp-s-i models in the following analysis. On the other hand, stellar luminosity in the 2005 epoch should be much lower than that in the 2017 epoch. This constraint leads to a final collection of only three reasonable sp-s-i models. For the three sp-s-i models, we calculated the mean $\langle T_{\mathrm{eff}}\rangle$, $\langle R_{\ast}\rangle$,$\langle L_{\ast}\rangle$, and foreground $\langle A_{V}\rangle$, weighted by $\frac{1}{\chi^{2}}$. These values are listed in Table 10 for the 2005 epoch.

In the above analysis for the R17 models satisfying $\chi^{2}-\chi^{2}_{\mathrm{best}}<3n_{\mathrm{data}}$, no constraint on the foreground extinction $A_{V}$ was applied. The good SED models in the 2005 and 2017 epochs converge on a narrow range of foreground extinction $A_{V}=120-150$ mag, despite a broad $A_{V}$ range of $10-200$ mag used in the SED fitting procedure. This narrow $A_{V}$ range still holds for the other “bad” SED models. The observations at $7.8\micron$ and $11.8\micron$ in the 2005 epoch that cover the broad absorption feature at $9.7\micron$ yield a strong constrain on the foreground extinction. This is consistent with the deep absorption feature at $9.7\,\mu$m caused by the large absorption of silicate grains, which is commonly observed for deeply embedded high-mass protostars (e.g., Morales et al., 2009; Pitts et al., 2018).

The difference in luminosity between the 2005 and 2017 epochs indicates the luminosity burst of $\Delta L=7600\pm 900L_{\sun}$ during an accretion outburst of M17 MIR. If $\Delta L$ is entirely attributed to increased accretion, we can predict the flux density of M17 MIR in the Spitzer I1 and I2 bands during the outburst. We assume that the blackbody temperature of $\Delta L$ during outburst is $1000\,\mathrm{K}$. Reddenning the blackbody radiation from $\Delta L$ with a certain amount of extinction essentially matches the flux density in the Spitzer I1 and I2 bands. We assume a total extinction to M17 MIR of 150 mag, which is the sum of a foreground extinction of 130 mag according to the good SED models for M17 MIR and an additional extinction of 20 mag from the circumstellar disk. The blackbody radiation from $\Delta L=7600\pm 900\,L_{\sun}$ suffers from a total extinction of 150 mag, equivalent to $0.066\times 150$ mag and $0.052\times 150$ mag in the Spitzer I1 and I2 bands. This toy model yields flux density of 5.8 mJy and 52.5 mJy in the Spitzer I1 and I2 bands, respectively, which are comparable to the observed values $3.8\pm 0.7$ mJy and $59.7\pm 5.9$ mJy in the 2017 epoch, differing by a factor of less than 2.

The luminosities of M17 MIR in Table 10 are indeed the total luminosites, i.e., the combination of photospheric luminosity ($L_{\ast}$) and accretion luminosity ($L_{\mathrm{acc}}$),

In the 2005 epoch, we assume that half of $L_{\mathrm{2005}}=1400\,L_{\sun}$ arises from a ZAMS photosphere with a solar metalicity, and the progenitor’s properties are: mass $M_{\ast}=5.4\,M_{\sun}$, radius $R_{\ast}=2.7\,R_{\sun}$, and effective temperature $T_{\mathrm{eff}}=18000\,\mathrm{K}$ (Tout et al., 1996), which lies between spectral type B2.5V and B3V (Pecaut & Mamajek, 2013). The accretion rate $\dot{M}_{\mathrm{acc}}$ in the 2005 epoch needed to produce the half of $L_{\mathrm{2005}}$ is $\sim 1.1\times 10^{-5}\,M_{\sun}\,\mathrm{yr^{-1}}$, comparable to the steady accretion rates for the bursting MYSO of Hunter et al. (2021) and for the $6\,M_{\sun}$ formation model of Haemmerlé et al. (2013). In the 2017 epoch, we assume $L_{\mathrm{acc}}\approx 8300\,L_{\sun}$, which yields $\dot{M}_{\mathrm{acc}}\approx 1.7\times 10^{-3}\,M_{\sun}\,\mathrm{yr^{-1}}$ in the 2017 epoch, where $M_{\ast}=5.4\,M_{\sun}$ and $R_{\ast}=34.1\,R_{\sun}$ are assumed. The value of $\dot{M}_{\mathrm{acc}}$ in the 2017 epoch estimated from the stellar parameters agrees with the value of $\sim 2\times 10^{-3}\,M_{\sun}\,\mathrm{yr^{-1}}$ obtained from the proper motions of H₂O maser spots (Sect. 5.2). In the 2017 epoch, a larger stellar radius is required to produce an accretion rate comparable to the value estimated from observations. This fact is consistent with the simulation prediction that a large accretion rate may cause the protostellar photosphere to bloat, resulting in larger radius (Yorke & Bodenheimer, 2008; Hosokawa & Omukai, 2009; Meyer et al., 2019a).

The amplitude of MIR flux variation is directly proportional to the enhanced $L_{\mathrm{acc}}$. Nevertheless, the enhanced $L_{\mathrm{acc}}$ is a reduced version of enhanced $\dot{M}_{\mathrm{acc}}$, due to the significant inflation of the central object during the burst accretion phase. In the burst accretion phase, the $\dot{M}_{\mathrm{acc}}$ is amplified by a factor of hundreds compared with that in the quiescent phase. M17 MIR even shows a WISE/NEOWISE W2 flux variation in a cadence of six months, despite the overall trend of becoming continuously brighter. The overall trend is suggested to originate from the increasing $\dot{M}_{\mathrm{acc}}$. On the other hand, the smaller variation in the cadence of six months might imply unsteady $\dot{M}_{\mathrm{acc}}$ with a smaller amplitude. For instance, the variation in the W2 flux density is $-42\%$ between 2017.09 and 2018.03. This variation corresponds to $\Delta\dot{M}_{\mathrm{acc}}=-42\%$, if we assume the stellar radius to remain unchanged. This is the largest variation in $\dot{M}_{\mathrm{acc}}$ during the burst accretion phase.

Considering the MIR variation during the decreasing phase 1993.03–2004.09, M17 MIR is about three times brighter in 1993.03 than in 2017.08. We speculate that $\dot{M}_{\mathrm{acc}}$ in 1993.03 is $\sim 5\times 10^{-3}\,M_{\sun}\,\mathrm{yr}^{-1}$, three times that in the rebrightening phase. However, $\dot{M}_{\mathrm{acc}}$ in 1993.03 is probably even higher, since M17 MIR may expand to even larger size than that in the rebrightening phase.

The MIR light curve of M17 MIR exhibits three phases characterized by M17 MIR’s MIR variation in flux density over 26 yr. From the above results and analyses, we conclude that the decreasing and rebrightening phase of M17 MIR indeed reflect the two accretion bursts, which are separated by the quiescent phase lasting about 6 yr. The first accretion burst likely began between 1984 and 1993, and lasted longer than 11 yr but less than 20 yr. Thus we use a middle value of 15 yr. We fortunately caught the initial stage of the second accretion burst, aided by the multiepoch observations of the Spitzer and WISE/NEOWISE space telescopes. The duration of the second accretion burst is definitely longer than 9 yr. It is very likely that M17 MIR has not yet reached to its peak brightness, and the duration of the second accretion burst is similar to that of the first epoch. M17 MIR is the first MYSO showing multiple accretion bursts with durations of tens of years and a burst magnitude of about 2 mag. The observed properties of the accretion bursts of M17 MIR agree with the prediction from theoretical simulation that 2 mag bursts have a mean duration of about 20 yr (Meyer et al., 2019b).

We assume a bursting accretion timescale of 15 yr, which is comparable to the duration of the decreasing phase in the MIR light curve (somewhere between 1984 and 1993 to mid-2004). If M17 MIR can sustain this enhanced disk accretion rate of $\sim 2\times 10^{-3}\,M_{\odot}\,\mathrm{yr^{-1}}$, it could gain about $0.03\,M_{\odot}$ during one accretion burst. This amount of accreted mass is an order of magnitude higher than the mass accretion onto S255IR–NIRS3 (Caratti o Garatti et al., 2017) and two orders higher than that onto G358.93–MM1 (Stecklum et al., 2021). The quiescent phase of M17 MIR lasts about 6 yr. During the time span of 26 yr, the fraction of time in accretion burst is about 83%. The main accretion phase of a high-mass protostellar object in an HMC lasts a few times $10^{4}$ yr (Furuya et al., 2005; Bonfand et al., 2017). A crude estimate, assuming $t_{\mathrm{acc}}=1\times 10^{4}\,\mathrm{yr}$, is that this fraction of time in accretion burst yields $\sim 14\,M_{\odot}$ mass that M17 MIR could accrete from its circumstellar disk.

The sample of MYSO accretion bursts is still small. Six MYSOs show a luminosity increase due to an accretion burst: V723 Car (Tapia et al., 2015), S255IR–NIRS3 (Caratti o Garatti et al., 2017; Liu et al., 2018), NGC 6334I–MM1 (Hunter et al., 2017, 2021), G323.46–0.08 (Proven-Adzri et al., 2019), G358.93–0.03 MM1 (Stecklum et al., 2021), and M17 MIR in the present work. However, these MYSOs are in different pre-main-sequence (PMS) evolutionary stages. The durations of accretion bursts are poorly constrained except for M17 MIR. To qualitatively describe the relation between the PMS evolutionary stage and an accretion burst, their key parameters are summarized in Table 11. Relatively evolved MYSOs are commonly brighter in the near-IR $K$ band, while the earliest MYSOs are very faint or even invisible in the same band. Generally the $K$-band brightness of these MYSOs is representative of the PMS evolutionary stage. Among the six MYSOs tabulated in Table 11, four sources show variabilities in the $K$ band, indicating relatively late PMS stages. In particular, G323.46–0.08 was the first one that showed periodic 6.7 GHz methanol maser emission with a period of about 90 days (Proven-Adzri et al., 2019). Later – on the basis of multiepoch IR data – it was confirmed as an accretion burst MYSO (see Stecklum et al., 2021, and references therein).

For V723 Car, S255IR–NIRS3, and G358.93–MM1 the duration of the accretion burst $\Delta t$ was constrained from their IR luminosity changes. The durations of their luminosity increases are of order of years and thus comparable to each other. The two earliest MYSOs, NGC 6334I–MM1 and M17 MIR, are both embedded in HMCs. The progenitor of NGC 6334I-MM1 in its pre-burst phase might be able to produce a hypercompact region MM1B (Brogan et al., 2016), implying that NGC 6334I-MM1 is more massive and at a later evolutionary stage than M17 MIR. Acoording to the schematic evolutionary sequence of a massive protostar (Motte et al., 2018), M17 MIR might be the earliest MYSO with accretion outbursts. Considering the extreme youth of M17 MIR, its multiple accretion bursts suggest that minor accretion bursts are frequent in the very early stages of massive star formation.