© The Authors 2024

1. Introduction

The mass-gain history of young stellar objects (YSOs) is characterized by variable accretion rates (Audard et al. 2014; Fischer et al. 2023). Compared to steady-state accretion scenarios, variable accretion introduces complexities into our understanding of protostellar jets and envelope chemistry. For example, episodic accretion may be the cause of the clumpiness in high-velocity protostellar jets (Ellerbroek et al. 2014; Vorobyov et al. 2018; Garufi et al. 2019). Additionally, the periodic release of accretion energy by low-mass YSOs (LMYSOs, with M_* < 2 M_⊙) is known to affect the chemical abundances, ice, and dust properties of the surrounding environment (Ábrahám et al. 2009; Taquet et al. 2016; Vorobyov et al. 2022). It has been argued that protostars of all masses grow by disk accretion (Beltrán & de Wit 2016; Caratti o Garatti et al. 2017). However, the mechanisms and environmental effects of episodic accretion by LMYSOs are much better understood than for their intermediate-mass (M_* = 2 − 8 M_⊙) and high-mass (M_* > 8 M_⊙) counterparts (Elbakyan et al. 2021, 2023; Fischer et al. 2023). The larger mass gain of high-mass YSO (HMYSO) bursts may also lead to stronger effects on the chemistry of the surrounding environment than is expected for LMYSOs (Dunham & Vorobyov 2012; Guadarrama et al. 2024). Hydrodynamic simulations have suggested that HMYSOs gain up to 50% of their zero-age main-sequence mass in the bursting stage of accretion, which constitutes < 2% of the stellar formation times (Meyer et al. 2019a, 2021). On the other hand, LMYSOs only gain 5%−35% of their final mass in bursts (Dunham & Vorobyov 2012). To understand how episodic accretion in HMYSOs impacts the broader theory of star formation, we must investigate its rates, underlying mechanisms, and effects on the surrounding environment. There are theoretical predictions to the number and magnitude of accretion bursting events for different HMYSO masses, where the modelling is done without magnetic fields, and for elementary initial geometries (Meyer et al. 2021). Only about four accretion-bursting HMYSOs have been studied in depth during the onset of their accretion bursts: S255IR-NIRS3, NGC6334I-MM1B, G358.93−0.03 and G24.33+0.14 (Hunter et al. 2017; Caratti o Garatti et al. 2017; Burns et al. 2020; Hirota et al. 2022). There are also a few HMYSOs with suspected accretion bursts, based on archival maser or infrared observations: V723 Car, G323.46−0.08 and M17 MIR (Tapia et al. 2015; Proven-Adzri et al. 2019; Chen et al. 2021). Disks around massive stars are difficult to observe, and it has been theoretically challenging to reproduce the rise times, accretion rates and accretion luminosities of the observed HMYSO bursts with mechanisms proposed for LMYSO bursts (Elbakyan et al. 2021). Finally, the consequences of observed intermittent matter and energy ejections into the surrounding environment need to be investigated (Cesaroni et al. 2018; Burns et al. 2020).

All of the HMYSO bursts were detected by follow-up observations prompted by the flaring of the 6.7 GHz methanol maser transition (Fujisawa et al. 2015; MacLeod et al. 2018; Sugiyama et al. 2019). A coordinated effort to detect and investigate bursting HMYSOs is being led by the Maser Monitoring Organization (M2O, Burns et al. 2022). Single-dish monitoring of multiple maser lines has shown to be effective in identifying accretion bursts. The radiatively pumped 6.7 GHz methanol maser transition occurs in hot (T_dust > 150 K) and moderately dense (10⁵ < n_H2 < 10⁹ cm⁻³) gas (Cragg et al. 2002, 2005). The correlation of 6.7 GHz methanol masers with the dust temperature implies that these masers could serve as an indicator and a ‘timer’ tracking the evolution of the radiation for accretion bursts in HMYSOs (e.g. MacLeod et al. 2018). It should be noted that a 6.7 GHz maser flare on its own is not enough to confirm an accretion burst and contemporaneous infrared and/or millimetre increase should also be detected (Caratti o Garatti et al. 2017; Hunter et al. 2021; Stecklum et al. 2021). Another bright and common maser species that has been observed to flare during accretion bursts is the 22 GHz water maser transition (e.g. MacLeod et al. 2018; Brogan et al. 2018; Bayandina et al. 2022a). In high-mass star-forming regions (HMSFRs), the 22 GHz water maser transition is typically understood to be collisionally pumped in dense (10⁸ < n_H2 < 10¹⁰ cm⁻³) shocked gas. The masers have different properties based on whether they originate in warm (T_kin ∼ 400 K) post-shock gas from a high-velocity (v_shock > 40 km s⁻¹) dissociative J-shock or in hot (T_kin ∼ 1000 K) gas within a low-velocity (v_shock < 40 km s⁻¹) non-dissociative C-shock (Elitzur et al. 1992; Kaufman & Neufeld 1996; Hollenbach et al. 2013).

The most detailed models of water masers in J-shocks and C-shocks, as well as simpler broad parameter searches, all neglect the effects of external ultraviolet radiation (Kaufman & Neufeld 1996; Hollenbach et al. 2013; Gray et al. 2016). Preliminary theoretical calculations of water masers in irradiated gas have suggested that increased stellar radiation dampens the effectiveness of the 22 GHz water maser pump (Nesterenok 2013). Observationally, an inverse correlation between 6.7 GHz methanol masers and 22 GHz water masers have been found in G107.298+5.639, G025.65+1.05 and the masers < 1000 au from the bursting source NGC6334I-MM1B (Szymczak et al. 2016; Brogan et al. 2018; Sobolev et al. 2019). On the other hand, it was also observed that in NGC6334I some velocity components of 22 GHz water and 6.7 GHz methanol masers flared within 0.2 yr of each other (Figure 9 of MacLeod et al. 2018). A lag of two months for the onset of 6.7 GHz methanol and 22 GHz water maser flares can be significant, as in G358.93−0.03 MM1 (Bayandina et al. 2022b). However, in the case of NGC6334I, large random variations in the 22 GHz water maser time series do not allow a determination of the flare onset with two-month precision. The contemporaneously flaring masers were at a similar velocity, but interferometric observations have revealed that they are likely not co-spatial. The 22 GHz water masers and the 6.7 GHz methanol masers are around 2750 au and 2000 au north of the bursting source MM1B, respectively (Hunter et al. 2018; Brogan et al. 2018). Nonetheless, the precise mechanisms leading to the contemporaneous flaring of two maser species in NGC6334I that have contrary pumping mechanisms have not yet been adequately addressed. Previous studies have proposed that the flaring of the 22 GHz water masers in sources such as NGC6334I, S255IR-NIRS3, and G358.93−0.03 was caused by amplified radiation fields (MacLeod et al. 2018; Hirota et al. 2021; Bayandina et al. 2022a). This explanation has invoked both collisional pumping as well as radiative pumping of 22 GHz water masers in conditions with T_dust > 1400 K, T_kin < 500 K and a narrower density range than collisionally pumped water masers (Brogan et al. 2018; Gray et al. 2022).

Apart from the 22 GHz transition, water also has many other maser lines such as 183 GHz, 321 GHz and 325 GHz. The most up-to-date list of water maser detections can be found at www.maserdb.net (Ladeyschikov et al. 2022a) For predictions of maser lines, we refer to Gray et al. (2016). Simultaneous observations of multiple co-spatial maser lines, combined with model predictions, can allow us to constrain hydrogen densities, kinetic temperatures, dust temperatures, and water abundances. The 10_2, 9 → 9_3, 6 321.225677 GHz (hereafter 321 GHz) water maser line can be combined with the 22 GHz water maser to distinguish between different types of shocks. It is rare compared to the 22 GHz transition, but has been detected with single-dish (e.g. W3(OH), W49N and W51 IRS 2) and interferometric observations (e.g. Cepheus A and Orion I; Menten et al. 1990; Patel et al. 2007; Hirota et al. 2014a). The 321 GHz water line has an upper energy level at E_u/k = 1861 K, which implies that it forms at high kinetic temperatures (T ∼ 1000 K). Suppose that at 321 GHz and 22 GHz, the water masers in a particular source are co-spatial, then the emissivity ratio R = L₂₁/L₃₂₁ can be used to discriminate between excitation from different shock types under the assumption of equal beaming for the two lines. For R ∼ 1, it implies that the water maser is probably excited in C-shocks, where T_kin ∼ 1000 K (Kaufman & Neufeld 1996). Otherwise, for R ∼ 10 it implies J-shocks, where the masers form in T_kin = 400 K gas (Hollenbach et al. 2013). The line ratios can also be combined with maser models to constrain other maser physical conditions (e.g. Yates et al. 1997; Gray et al. 2016). The first detection of 321 GHz masers in an accretion-bursting source was in S255IR-NIRS 3 by Hirota et al. (2021). They found regions with both high R ∼ 10 and low R ∼ 1. This shows a variety of shock types that excite water masers in star-forming regions, which might respond differently in the time-variable radiation environments of accretion-bursting sources.

Our target, NGC6334I (d = 1.3 kpc Chibueze et al. 2014; Wu et al. 2014), is a complex HMSFR. It contains nine millimetre sources MM1-9 (Hunter et al. 2006; Brogan et al. 2016). Figure 1 shows an overview of the millimetre and centimetre continuum of the source. The largest star-forming clump, MM1, is a hot multi-core, containing seven dust cores with T_dust > 100 K and high column densities > 1 × 10²⁵ cm⁻² (Brogan et al. 2016). The most compact core is MM1B; it has a radius smaller than 155 au, and started undergoing an accretion burst in January 2015 (Hunter et al. 2017; MacLeod et al. 2018). Infrared observations with SOFIA combined with radiative transfer modelling have found that MM1B has a mass of 6.7 M_⊙ and is undergoing a burst with an accretion rate of Ṁ ≥ 2.3 × 10⁻³ M_⊙ yr⁻¹ that will last for 40 − 130 yr (Hunter et al. 2021).

Fig. 1. Overview of the source. Left: ALMA 2015.6 1.3 mm dust continuum from Brogan et al. (2016). The main regions MM1-4 are labelled in white. The peak positions of the bursting source, MM1B, as well as MM2 and MM4 are marked with crosses. The location of the millimetre cavity is marked in red. Right: JVLA 2015.9 5 cm continuum with 2017.8 22 GHz water masers overlaid in blue (Brogan et al. 2018). The positions and names of the water maser associations introduced in Section 1 are in green. The positions and orientations of the N-S and NW jets (solid and dashed-dotted respectively, from Brogan et al. 2018) and NW-SE jet (dashed line, from Chibueze et al. 2021) as well as the 0.5 pc NE-SW outflow (dotted line, from Qiu et al. 2011) are marked in red.

NGC6334I-MM1 harbours a large-scale wide-angle outflow, and at least three jets (Qiu et al. 2011; Brogan et al. 2018; Chibueze et al. 2021). The right panel of Figure 1 shows the positions, orientations, and sizes of the NE-SW outflow and the N-S, NW-SE, and MM1B-NW jets. The positions of 22 GHz water masers observed with the JVLA in 2017.8 and the names of the maser associations relative to the outflows and jets from Brogan et al. (2018) are also shown. The NW-SE jet is associated with the CM2-W2 water maser association at the northern edge and UCHII-W3 at the southern edge. The NW-SE jet has a small inclination angle of −6° in the plane of the sky (Chibueze et al. 2021). On the other hand, the southern edge of the N-S jet is associated with UCHII-W1 and UCHII-W2 (Figures 8 and 9 of Brogan et al. 2018). The maser association MM1-W1 is very close to the bursting source MM1B and might be affected by multiple outflows (Chibueze et al. 2021).

We wanted to characterize the effect of the accretion burst in NGC6334I on 22 GHz water maser variability and sought to determine its underlying mechanisms. We conducted VLBI monitoring observations of 22 GHz water masers in NGC6334I before, during and after the onset of the accretion burst, using the VLBI Exploration of Radio Astrometry array (VERA, VERA Collaboration 2020). Additionally, we report the first detection of 321 GHz water masers in NGC6334I with Atacama Large Millimeter/submillimeter Array (ALMA) Cycle 6 observations. Our observations are co-temporal with the burst onset, allowing us to track changes in water masers as the burst begins, in contrast to previous studies by Brogan et al. (2018). This work also extends the proper motion study of Chibueze et al. (2021) by covering a longer time frame and considering a broader set of observations and theoretical backgrounds.

2. Observations and data reduction

2.1. VERA observations

We did seven epochs of observations of the 22.23508 GHz water maser transition with VERA between 2014 and 2016. These epochs cover the time before, during and after the onset of the accretion burst (see Table 1). NRAO530 was used as the phase and bandpass calibrator while J1713−3418 was used as the flux calibrator. The data were correlated with a software correlator in the NAOJ Mizusawa campus (Oyama et al. 2016). The spectral resolution was 0.42 km s⁻¹ for all epochs and the beam size was 2.2 × 0.8 mas with a position angle of −13°. The data were reduced with standard data-reduction procedures with the Astronomical Image Processing System (AIPS) (van Moorsel et al. 1996). For each epoch, self-calibration was applied to the bright −7.6 km s⁻¹ feature to improve the signal-to-noise ratio of the images. After self-calibration, the line-free RMS noise, σ_RMS, was around 0.8 Jy beam⁻¹ for each epoch, but could be worse in channels with bright maser emission. Phase referencing was not successful in any observing session, due to sensitivity constraints from the low elevation of NGC6334I (δ = −34°) observed with VERA at a latitude 30° −40°. The loss of precision in absolute position does not affect any of our results. The absolute position on the sky was estimated to sufficient precision by finding the position of our −7.6 km s⁻¹ self-calibration reference feature in the more sensitive contemporaneous observations by Chibueze et al. (2021). We accounted for the differences in the velocity reference frame between the KaVA and VERA observations, as well as the proper motion of the feature (see Section 2.1.2). In Appendix A we compare our observations to those of Chibueze et al. (2021).

2.1.1. Maser spots and features

The maser maps were produced by fitting 2D Gaussian functions on the emission peaks in each channel with the SAD AIPS task. SAD erroneously identified imaging artefacts, such as sidelobe emission, as detections. Sidelobe emission was identified as faint maser spots spatially symmetric around a bright feature. Sidelobe spots also displayed non-Gaussian spectral profiles. For our observations, we found that grouping spots that extended over at least four contiguous channels and showed a Gaussian spectral profile showed a good balance of minimizing false detections and maximizing real detections. Maser spots are individual Gaussian fits to single channels, while maser features are physically grounded, have a full-width at half maximum of 0.5 − 2 km s⁻¹ and are spread across multiple channels for observations which spectrally resolve the line. Maser features have a characteristic size of 1 au (e.g. Gwinn 1994a; Torrelles et al. 2001; Marvel et al. 2008). The grouping of spots into features was done with the flux-weighted mean of maser spots, given by:

$$\begin{array}{c}\hfill {\mathit{x}}_{\mathrm{f}}=\frac{\sum I_i{\mathit{x}}_i}{\sum I_i}\end{array}$$(1)

with x_i and I_i the spot spatial coordinates and intensity, respectively. The uncertainties on the position of the feature were calculated with 10⁵ Monte Carlo (MC) simulations for each feature. In each MC iteration, a position was sampled from a normal distribution with a width equal to the astrometric uncertainty of the maser spot, and Equation (1) was applied. The uncertainties reported in Table 2 are 3σ_RMS of the 10⁵ MC runs for each feature. The uncertainties for maser features are 30 − 200 μas in Right Ascension and 70 − 500 μas in Declination. The S/N of the maser feature determined its positional uncertainty. The peak intensity, I_peak, radial velocity V_LSR and full-width half-maximum (FWHM) of the maser feature was calculated with a Gaussian fit to the spectral profile of the maser feature.

2.1.2. Proper motions

We calculated relative proper motions for persistent maser features over consecutive epochs and with V_LSR within 0.5 km s⁻¹. Due to burst-induced variability, we also split the epochs into pre-burst (epochs 1−4) and burst (epochs 5−7). The accretion burst introduced variability in all epochs after 2015.0, but certain features persisted until epoch 4 (2015.3). Relative proper motions were calculated with a linear fit to the position over time of the feature positions x_f. The uncertainties were uncertainties on the fit.

Relative proper motions are not physical, as they are in the rest frame of the reference feature. Further, without absolute phase referencing, the proper motions cannot be put in absolute coordinates in the sky. We shifted the proper motions to an absolute position and velocity frame.

The absolute position was obtained by shifting the reference frame to the reference feature in UCHII-W1 at (α,  δ) = (${17}^{\mathrm{h}}{20}^{\mathrm{m}}52\stackrel{\text{s}}{.}600$, $-35°46\prime 50\stackrel{″}{.}508$) of 2015.95−2016.01 KaVA observations Chibueze et al. (2021), that corresponds to the time between epoch five (2015.88) and six (2016.11) of our observations. This was done by setting the pre-burst and burst proper motions of our VERA reference feature equal to the proper motions of the same maser feature in the KaVA observations (proper motion IDs: 70 − 73 of Table 1 in Chibueze et al. 2021). The change in velocity frame was done by changing the spot maps, not the proper motions. It was done with the equation

$$\begin{array}{c}\hfill {\mathit{x}}_{\mathrm{shifted}}={\mathit{x}}_{\mathrm{unshifted}}+{\mathit{\mu }}_{\mathrm{T}}{t}_{\mathrm{obs}}+{\mathit{P}}_{\mathrm{T}}.\end{array}$$(2)

Here x_(un)shifted is the (un)shifted position of the maser spots. The proper motion transformation vector μ_T in mas yr⁻¹ is the difference in the proper motion measurements of the same feature between our VERA observations and the contemporaneous KaVA observations. The time after the first epoch is t_obs in years. Lastly, P_T is the change in offset to make the KaVA reference feature the offset zero-point. For the pre-burst proper motions μ_T = ( − 1.9 ± 0.3,  17.3 ± 0.4) mas yr⁻¹ and for the burst proper motions μ_T = ( − 2.0 ± 0.3,  16.9 ± 0.4) mas yr⁻¹. The position shift was the same for the pre-burst and burst epochs with P_T = ( − 0.68921, 5.28800) arcseconds. The uncertainty in the proper motions due to the velocity frameshift changes in the following way:

$$\begin{array}{c}\hfill \mathrm{\Delta }\mathit{\mu }=\sqrt{{(\mathrm{\Delta }{\mathit{\mu }}_{\mathrm{V}})}^2+{(\mathrm{\Delta }{\mathit{\mu }}_{\mathrm{K}})}^2}\end{array}$$(3)

with Δμ the uncertainty after the frame shift, Δμ_V mas yr⁻¹ the uncertainty of each VERA proper motion, and Δμ_K the uncertainty of the same feature measured in Chibueze et al. (2021). After this shift, the proper motions are absolute in position, but not yet in velocity.

The proper motions were shifted to an absolute velocity frame by assuming that all the maser features in CM2-W2 have a mean speed equal to the upper limit of 3.8 mas yr⁻¹ set by the JVLA observations of Brogan et al. (2018). In our case, the absolute velocity frame could not be calculated by assuming both poles of a bipolar outflow to have the same velocity (as in e.g. Burns et al. 2016). If we assume the masers are moving at a position angle of −79.4° (Chibueze et al. 2021), it implies the upper limits (μ_α < −0.7 mas yr⁻¹, μ_δ < 3.735 mas yr⁻¹). These correspond to a proper motion shift of (0.7, −15.292) mas yr⁻¹ for the pre-burst proper motions and (0.7, −15.301) mas yr⁻¹ for the burst proper motions. Our proper motions are upper limits for northward proper motions and lower limits for southward proper motions. None of our conclusions are affected by uncertainties included by the shift to an absolute position and velocity frame.

2.2. ALMA band 7 observations

We observed NGC6334I with ALMA in Cycle 6 (PI: T. Hirota, 2018.1.00024.S) in band 7. The observations were carried out on 1 September 2019 toward (α,  δ) = (17^h20^m53$\stackrel{\text{s}}{.}$44, −35°47′02$\stackrel{″}{.}$20) in configuration C43−7 with 48 antennas. The amplitude and bandpass calibrator was J1924−2914, while the phase calibrator was J1717−3342. The on-source time was 81.2 minutes. The pre-imaging calibration was done by the East-Asian ALMA Regional Center (EA-ARC) calibration pipeline. We then carried out phase-only self-calibration and imaging with the Common Astronomy Software Applications (CASA, McMullin et al. 2007). We only report the 321.2 GHz narrow-band spectral window relevant to the 321.225677 GHz water maser for this work. Line data were obtained by subtracting the continuum in the UV domain. The spectral resolution was 0.245 MHz corresponding to a velocity resolution of 0.23 km s⁻¹. We regridded the spectra to a velocity resolution of 0.25 km s⁻¹. The beam size was 0.14 × 0.10 arcsec with a position angle of −78°. The CASA task tclean produced the spectral cubes with a robust clean parameter of 0.5 was used with a cell size of 0.0125′. After self-calibration, the line-free RMS noise was 27 mJy beam⁻¹.

3. Results

The overall distribution of water masers in the region was consistent with the bipolar structure reported by Brogan et al. (2018). We detected water masers in six distinct regions, previously identified and labelled as CM2-W1, CM2-W2, MM1-W1, UCHII-W1, UCHII-W2, and UCHII-W3¹. Figure 2 shows the positions and radial velocities of these maser associations across the entire field. Detailed information about individual maser features, including their associations, positions, radial velocities, peak flux densities, and FWHM for epochs 1−7, is provided in Table 2.

Fig. 2. Positions and radial velocities of 22 GHz water masers at 2014.7−2016.2 as observed with VERA. The names of the water maser associations are annotated in the leftmost panel. The linear scale for d = 1.3 kpc is shown at the top left of each panel. The accretion burst started at 2015.0 (after epoch 2), and peaked in 2015.6 (between epochs 4 and 5).

According to MacLeod et al. (2018), the onset of the accretion burst was 2015.01 and the methanol maser flare peak was 2015.62. Our observations allow us to see the evolution of the water maser emission before (2014.7−2015.3) and after the onset of the burst (2015.9−2016.2). Time dependence in spatial distribution and flux density was found in all the water masers, over scales from 6000 au (5″) down to 1 au (submilliarcsecond) scale. Figure B.1 shows zoom-ins on the positions and radial velocities of the water masers in the different associations over the seven epochs. Figure 3 shows the velocities over time for the maser features in each association. Figure 4 shows the total line area (velocity integrated flux density) of each association over time. The bottom panel of Figure 4 also shows their velocity range over time. The maser association CM2-W2 had many more maser features than the other associations.

Fig. 3. Centre velocities and peak intensities of maser features in each association over time. The marker size is proportional to the square of the feature peak intensity ($\propto I_{\mathrm{peak}}^{0.5}$). A scale for the marker size is shown at the bottom left of each panel. The markers are colour-coded according to velocity. The grey horizontal lines indicate the dates of our observations. The date of the onset and peak of the accretion burst traced by 6.7 GHz methanol masers is shown in red and black, respectively (MacLeod et al. 2018).

Fig. 4. Time dependence of total velocity integrated intensity (line area) and velocity spread for each association from our VERA observations. The first and second panel shows respectively the line area for associations containing increasing and decreasing trends at the onset of the burst. The third panel shows the velocity spread for each association over time. The colours for the third panel are the same as in the legends in the first two panels. The red and black lines indicate the date of the onset and peak of the burst, respectively.

3.1. Maser variability in CM2-W2

On a large scale, before the burst in 2014.7, CM2-W2 was only a few maser features around the systemic velocity of −7.6 km s⁻¹ spread over 100 au (Figures 2 and B.1). In the following epochs, the water masers progressively traced a bow shape with a size of ∼400 au. CM2-W2 also brightened over time, peaked at 2015.3, and stayed flaring for all epochs (Figure 4).

There was also a velocity-dependent response in the maser features as the epochs progressed. As a whole, the velocity extent of the region increased linearly with time (Figures 3 and 4). Figure 5 shows the time dependence of maser features in different velocity bins over time. The size of the markers is proportional to the peak flux density of the maser features. The majority of the bright masers were detected at between −15 km s⁻¹ and 0 km s⁻¹. After the burst, these masers traced a linear structure spread over 0.27″ (356 au). The population of masers with V_LSR ∈ [ − 30, −15] km s⁻¹ were only detected after the onset of the burst (2015.08) reaching a peak in intensity at 2015.88, forming a 0.16″ (210 au) linear structure with a clear NE-SW velocity gradient. In the epochs after 2015.88, these features reduced in brightness and linear extent, while retaining the velocity gradient. The V_LSR ∈ [ − 30, −15] km s⁻¹ reached their peak one epoch (0.6 yr) later than the 6.7 GHz methanol masers and −7.6 km s⁻¹ water masers in CM2-W2. Lastly, two high-velocity features at −40 km s⁻¹ and −50 km s⁻¹ were only detected after 2015.9. Velocity-dependent time evolution was also seen in 2017.0, and 2017.8 (Brogan et al. 2018).

Fig. 5. Velocity-dependent response of the water maser features in CM2-W2. The radius of the spots is proportional to the peak flux density. The rows indicate different velocity bins as indicated in the panels in the first column. Each column is the observation shown in the first row. Maser features that are outside of the respective velocity bin are shown in grey. The spatial scale for all panels is the same. The offset is in terms of the reference maser feature.

There were also variations in the water masers in CM2-W2 on scales of 1 au. Figure 6 shows the position and spectra of the reference maser feature. In 2014.7, the reference feature is a single-peaked feature with an extent of 460 μas (0.6 au), which had been detected with JVLA in 2011 (Brogan et al. 2018). The maser feature is unresolved, but the size of the maser feature is defined by the distance between the most red- and blue-shifted spots (e.g. Gwinn 1994a; Marvel et al. 2008). In 2014.9 and 2015.1, the feature was double peaked, while having the same extent on the sky. Between 2015.1 and 2015.3, the maser spots were over a larger extent of 1060 μas (1.4 au) and flared by a factor of 4. The maser retained the large size and flaring state as the epochs continued. The spatial precision which allows us to detect the increase in the size of the masers is due to the high signal-to-noise (S/N = 700) of the maser observations in the flaring state.

Fig. 6. Time variation of the brightest maser feature in CM2-W2. Top panels: Offsets and VLSR of maser spots in the brightest feature of CM2-W2 for each epoch. The position, uncertainty and VLSR of each spot is shown. We note the difference in spatial scale between right ascension and declination. Bottom panels: Spectral profile of the maser spots in this feature, with Gaussian fits shown in black. In 2014.9 and 2015.1 the feature is well approximated by a double Gaussian.

3.2. Maser variability in other associations

CM2-W1. This association was detected in the first four epochs at −9 km s⁻¹ (2014-7−2015.3, Figures 2 and 3). It peaked in flux density in 2015.3 (Figure 4), but was not detected after that.

MM1-W1. This association is close to the bursting source. In epochs 2014.7−2015.1 it formed a linear structure of ∼15 au (Figure B.1). After 2015.3, the masers were displaced and then had a smaller linear scale of 7 au. The association monotonically decreased in line area and velocity range as the epochs progressed (middle panel of Figure 4).

UCHII-W1. This association consisted of multiple maser associations spread over 1000 au at the southern point of the N-S molecular outflow detected by Brogan et al. (2018). As a whole, the region flared in flux density after the onset of the burst, peaking in 2015.3, after which it dramatically dropped in line area (Figure 4). The region also dramatically decreased in velocity spread from 30 km s⁻¹ in 2014.7 to 10 km s⁻¹ in 2016.2 (Figures 3 and 4). Figure 7 shows the four sub-associations, W1-i−iv, numbered from south to north. The southernmost association, W1-i, with two features at −21.2 km s⁻¹ and −6.5 km s⁻¹, disappeared in 2015.3. The second sub-association W1-ii had a single maser feature at −34.4 km s⁻¹ which disappeared after the onset of the burst. A new maser feature at −16.2 km s⁻¹ was detected in 2016.1 and 2016.2. The third sub-association, W1-iii, was detected in all epochs around −10 km s⁻¹, but monotonically decreased in flux density as the epochs continued. The last sub-association, W1-iv, was also detected in all epochs, with monotonically decreasing flux densities.

Fig. 7. Velocity-dependent response of the water maser features in UCHII-W1. The area of the spots is proportional to the square of the peak flux density. The subregions of UCHII-W1, i−iv are marked in the first panel.

UCHII-W2. This association consisted of a single maser feature at −30 km s⁻¹ before the burst. This feature monotonically decreased in line area from 2014.7 until disappearing in 2015.3 (middle panel of Figure 4). It was not detected in 2015.9, with a new double-peaked feature detected in 2016.1−2016.2 at −31 km s⁻¹ and −26 km s⁻¹ displaced by ∼30 au. The masers in UCHII-W2 were also slightly red-shifted in centre velocity as the epochs progressed (Figure 3).

UCHII-W3. This association was first detected after the onset of the burst at 2015.1−2015.3 at −31 km s⁻¹, but it disappeared with two new features detected in 2016 at −36 km s⁻¹ and −49 km s⁻¹ separated by about 50 au. The features in this association showed a blue-shifting trend over time (Figure 3).

3.3. Pre-burst and burst proper motions

We detected a total of 20 and 22 relative proper motions in the pre-burst (Table 3, 2014.72 − 2015.28) and burst (Table 4, 2015.88 − 2016.19) epochs, respectively. Proper motions were detected in all associations where maser features were detected. Figures 8 and 9 show the proper motions for the northern regions (CM2-W2, MM1-W1) and southern regions (UCHII-W1, UCHII-W2 and UCHII-W3) respectively. There was only one proper motion in CM2-W1 that disappeared after the burst. The proper motions can be converted to transverse velocity with the relation v_t = 4.74μ ⋅ d with v_t in km s⁻¹, μ in mas yr⁻¹, d in kpc and the factor of 4.74 arising from unit conversions. The mean speeds did not change with ${\overline{v}}_{\mathrm{pre}}=50\pm 40$ km s⁻¹ and ${\overline{v}}_{\mathrm{burst}}=54\pm 42$ km s⁻¹. After the burst, blue-shifted proper motions with V_LSR ≈ −50 km s⁻¹ were detected. The mean speed of CM2-W2 was 33 km s⁻¹ before and after the burst, with only the velocity dispersion increasing after the burst. This is expected as the reference feature was set to 30 km s⁻¹ in all epochs to set the absolute velocity frame (see Section 2.1.2). The proper motions in CM2-W2 point in the NW direction, the same as the molecular outflow seen by Brogan et al. (2018) and through velocity-variance covariance analysis of proper motions (Chibueze et al. 2021). In MM1-W1 the mean speeds increased from v_MM1, pre = 11 ± 4 km s⁻¹ to v_MM1, burst = 28 ± 17 km s⁻¹, and the centroid position of the proper motions shifted by 15 au. This is consistent with either maser dropping below detection limits and new masers being excited (as seen in Brogan et al. 2018, between 2017.0 and 2017.8) or the same maser cloud undergoing a bulk displacement of 120 km s⁻¹ between 2015.3 and 2015.9. In the southern regions, almost all proper motions detected before the burst were not detected afterwards. In UCHII-W1, four associations named W1B-i−iv were detected. All of these associations had high speeds v_pre ∼ 90 km s⁻¹, except W1B-iii with v = 17 km s⁻¹. After the burst, only two high-velocity associations W1A-i (corresponding to W1B-iii) and W1A-ii (corresponding to W1B-iv) were detected. The association W1A-i was displaced by 12 au to the northeast, similar to MM1-W1. In UCHII-W2, there were southward pointing proper motions with V_LSR = −30.5 km s⁻¹ and v_{UCHII − W2, pre} = 70 km s⁻¹ and a new feature at V_LSR = −28.3 km s⁻¹ with v_{UCHII − W2, burst} = 76 km s⁻¹ displaced 6 au to the north-east. Finally, UCHII-W3 consisted of a single feature before the burst with V_LSR = −31.8 and v_{UCHII − W3, pre} = 147 km s⁻¹. After the burst, two new high-speed (v = 160 km s⁻¹) features at V_LSR = −48 km s⁻¹ and −36.8 km s⁻¹ were detected, 61 au and 110 au displaced from the pre-burst feature.

Fig. 8. Proper motions of 22 GHz water masers in the northern regions of NGC6334I, before and during the accretion burst. The positions, orientations and colours of the arrows show the positions, proper motions and VLSR of the maser features, respectively. The colour scale is on the right. The left panels show the pre-burst epochs (2014.72−2015.28) and the right panels show the burst epochs (2015.88−2016.19). The linear scales in position and space are shown in the bottom left corner of each panel. In the bottom panel, the black dashed lines show the position and orientation of the linear features shown in the maser features, and the black arrow shows the displacement of the maser features after the accretion burst.

Fig. 9. As in Figure 8, but for the southern regions.

In summary, the proper motions showed variability in all associations. High-velocity > 70 km s⁻¹ proper motions were detected in the southern regions, likely connected to either the N-S or NW-SE jets originating in MM1 (Brogan et al. 2018; Chibueze et al. 2021). The speeds for the southern regions are lower limits, as the mean velocity of CM2-W2, v = 30 km s⁻¹, is close to the upper limit of 23 km s⁻¹ from JVLA bulk motion between 2017.0 and 2017.8 (Brogan et al. 2018).

3.4. 321 GHz water masers

We detected a single 321 GHz water maser feature at the southwestern end of CM2-W2. This is the first detection of 321 GHz water masers in NGC6334I. Figure 10 shows the integrated intensity of the 2019.7 321 GHz water maser, with the positions and velocities of 2017.0 and 2017.8 22 GHz water masers detected with JVLA (Brogan et al. 2018). The 22 GHz water masers in CM2-W2 in 2011 were more localized than those observed in 2017.0, and not co-spatial with the 2019.7 321 GHz maser (Brogan et al. 2018). The right panel of Figure 10 shows the spectrum of the maser, which peaks at −7.75 km s⁻¹ with an integrated flux of 27.3 Jy. The small peaks in the spectrum are thermal HCOOCH₃ lines close to 321 GHz (Hirota et al. 2014b). Table 5 shows the position, V_LSR, integrated flux, peak intensity and deconvolved size from channel-by-channel 2D Gaussian fits to the maser feature. The maser feature is unresolved, showing the same elongation and size as the ALMA synthesized beam. Using the peak intensity and the upper limit for size reported in Table 5 the lower limit of the brightness temperature for the 321 GHz emission is 3.2 × 10⁴ K, implying maser emission. The comparison with the closest published 22 GHz water maser observations with JVLA shows some maser spots from 2017.8 which are co-spatial and at the same velocity as the 2019.7 321 GHz water maser. It is interesting that in 2017.0 there were no water masers in the same position and frequency as the 321 GHz masers. As 22 GHz water masers are bright in a larger parameter space than 321 GHz water masers (Gray et al. 2016), it could be that the 321 GHz water masers were only detectable between 2017.0 and 2019.7. We consider the significance of the 321 GHz maser detection and the possible effect of time variability in CM2-W2 in Section 4.3.

Fig. 10. Detection of the 321 GHz water maser with ALMA. Left: Integrated intensity map of 2019.7 321 GHz water masers as detected with ALMA in greyscale. The coloured crosses and triangles show the peak positions of 22 GHz water masers detected in 2017.0 and 2017.8, respectively, while the purple contour is the 5 cm continuum point source (Brogan et al. 2018). The ALMA beam and linear scale are shown at the bottom left and bottom right, respectively. Right: Spatially integrated spectrum of the 321 GHz water masers. The small peaks are HCOOCH3 thermal line emission, which is close to 321 GHz.

4. Discussion

Our VLBI monitoring observations contemporaneous with the onset of the accretion burst give some clues to the mechanisms behind the water maser flare in NGC6334I. In Section 4.1, we discuss the time variation of water masers in our VLBI observations taking the geometry of the various outflows and the neighbouring UCHII region into account. In Section 4.2 we use a toy model to distinguish three different sources of the variability and propose a statistical test on single-dish monitoring observations to distinguish mechanisms of water maser variability. Lastly, in Section 4.3, we aim to find the best explanation for the unexpected 22 GHz water maser flare in the bow-shock CM2-W2. In Appendix A we compare our proper motion calculations to previous estimates in NGC6334I and consider the conditions of reliability of proper motion measurements generally.

4.1. Burst-induced water maser evolution

4.1.1. Source radiation and outflow geometry of MM1

Before we consider the effects of the accretion burst on the water masers, we first consider the radiation spectrum from the accretion bursting source and the geometry of the complex outflows in MM1. The radiation spectrum for the protostar, disk, and cocoon environment of a bursting HMYSO is a complex open question. We can qualitatively consider the main radiation components of a protostar disk system to understand how episodic accretion drives diverse types of maser variability. An accretion-bursting source is expected to have infrared, optical, ultraviolet and possibly X-ray radiation components. Infrared radiation originates from the accretion disk and dusty natal cocoon (Hosokawa et al. 2010; Beltrán & de Wit 2016). Photospheric ultraviolet and optical radiation are expected to be bright, where low-optical depth cavities have been carved by jets and outflows, but faint when obscured by the disk which is often perpendicular to the jet axis (Rosen & Krumholz 2020). During the burst, the photospheric SED reddens but becomes more luminous due to stellar bloating (Meyer et al. 2019b). It is also possible that bursting HMYSOs have an X-ray emission component. Some FU Orionis sources have much higher X-ray emission than quiescent T Tauri stars, but the effect is not universal (Skinner et al. 2009; Kuhn & Hillenbrand 2019). Further, intermediate- to high-mass YSOs might have less X-ray emission, due to weak, complex magnetic fields (Augustson et al. 2016; Villebrun et al. 2019). In summary, for most angles relative to the jet axis, infrared emission dominates the SED, but for some viewing angles, the high-energy component should not be neglected.

We consider the geometry of the natal envelope MM1 and its corresponding outflows to estimate which radiation components are incident on each maser association. To do this, we must first establish the relationship between MM3-UCHII, MM1, and the various outflows. This understanding is crucial because the large UCHII region can deflect outflowing material, but it is not immediately obvious whether MM3-UCHII and MM1 are at the same distance. Additionally, we consider the matter along the path from the bursting source to the maser association. High dust content in the intervening path will reprocess any high-energy radiation from the bursting source to infrared wavelengths. We can make qualitative inferences by considering the morphology of the outflows, jets and maser associations relative to MM3-UCHII. These are introduced in Figure 1. Figures 8 and 9 of Brogan et al. (2018) are also instructive.

Firstly, MM3-UCHII must be either co-distant or behind MM1 for the free-free emission of the UCHII to provide seed photons for the maser associations UCHII-W1−3. We do see some influence of the UCHII region in deflecting the large-scale red-shifted lobe of NE-SW outflow (Figure 9 of Brogan et al. 2018). The undetected SE redshifted lobe of the NW-only jet might also be due to deflection from the UCHII region (Brogan et al. 2018).

The N-S and NW-SE jets provide crucial information about the system’s geometry. The NW-SE jet, with a small inclination angle of −6° (Chibueze et al. 2021), is associated with CM2-W2 and UCHII-W3 masers. Its SE lobe terminates at the bright edge of MM3-UCHII (Figure 1). The N-S jet’s maser associations and thermal emission extend into MM3-UCHII’s peak, suggesting MM3-UCHII might be behind UCHII-W1 and W2. This orientation could explain the lack of apparent deflection in the N-S jet. However, the effect of MM3-UCHII on the N-S jet remains tentative, with proper motion interpretations requiring careful consideration of velocity reference frames (see Appendix A).

Intervening dust and gas between MM1B and the maser associations is crucial for understanding the effects of accretion bursts on maser variability. A potential cavity SE of MM1B, visible in the 1.3 mm continuum (Figure 1), may allow a direct radiative connection between the bursting source and UCHII-W3. The N-S jet’s cavities are not visible (Figure 1), but this might be the case if the N-S jet has a large inclination angle or if it has a different driving source than MM1B. This geometry could allow UCHII-W1 and W2 to receive only thermal and infrared components from the accretion burst due to dust reprocessing of high-energy radiation. However, these conclusions depend on the accuracy of our source geometry interpretations.

4.1.2. Responses of water masers to the burst in NGC6334I

This section examines the response of water masers to an accretion burst in NGC6334I, focusing on several key maser associations. We analyse the behaviour of CM2-W2, MM1-W1, UCHII-W1, UCHII-W2, and UCHII-W3, comparing their responses and exploring the underlying physical mechanisms. Table 6 shows a summary of the observed trends in each maser association, predictions on the types of shocks exciting the water masers and estimates of each association’s incident energy components.

The bow-shock CM2-W2 dominates the water maser emission in NGC6334I. Single-dish observations of the 22 GHz water masers with HartRAO by MacLeod et al. (2018) indicate that > 80% of the post-burst flux density is between −5 km s⁻¹ and −12 km s⁻¹, which is dominated by CM2-W2 in our VLBI maps and the JVLA maps (Brogan et al. 2018). The water maser features in CM2-W2 gradually increased in number, spatial extent and flux density from the onset of the burst (Figures 2 and 5). This process of gradually tracing more features has been observed up to two years after the burst (2017.8, Brogan et al. 2018). This indicates that whatever process is causing the maser flare in CM2-W2 at large (> 50 au) scales is continuous, starting at the onset of the burst, rather than a short (< 1 yr) single pulse ‘heatwave’ as in G358.93−0.03 (Burns et al. 2020; Bayandina et al. 2022a). Further, our VLBI observations confirm what was shown by Brogan et al. (2018) that the large maser flares in the −5 km s⁻¹ to −12 km s⁻¹ velocity features seen by single-dish monitoring (MacLeod et al. 2018) are the combination of multiple maser features spread spatially over > 300 au. New persistent −7 km s⁻¹ features were excited south of CM2-W2 between 2015.88 and 2016.11 (top panel of Figure 5). Different velocity features in CM2-W2 respond to the radiative changes due to the accretion burst at different times. The −15 km s⁻¹ to 0 km s⁻¹ features (top panel of Figure 5) respond immediately to the accretion burst, while the −30 km s⁻¹ to −15 km s⁻¹ features peak 0.6 yr later (middle panel of Figure 5). The <  − 30 km s⁻¹ features flare even later, almost a year after the start of the burst (bottom panel of Figure 5). Individual features show jumps in flux density and size, but overall there’s a velocity-dependent continual ‘lighting up’ of features (Figures 5 and 6). This process continues at least 2.8 years after the accretion burst started (Brogan et al. 2018). This behaviour could be explained by an irradiated surface layer which radiatively heats up immediately at the onset of the burst and then mixes with the rest of the gas at the shock front over multiple months.

The shock type in CM2-W2 can be inferred from observations of both 22 GHz and 321 GHz water masers. In the 2017.8 water masers, a maser spot was detected co-spatial and at the same velocity as the 321 GHz water maser reported in this work (ID: 911 of Table 4 of Brogan et al. 2018). The 22 GHz spot has a reported integrated flux of 56 Jy, implying R = F₂₂/F₃₂₁ = 2.1, which is more consistent with high-temperature T ∼ 1000 K maser emission from a C-shock (Yates et al. 1997; Patel et al. 2007). This finding supports the interpretation of CM2-W2 as a C-shock region, which is crucial for understanding its unique response to the accretion burst. In Sections 4.2 and 4.3, we lay out a scenario in which variable radiation fields might explain the unique maser behaviour in CM2-W2.

The behaviour seen in CM2-W2 differs from that seen in other associations. In MM1-W1, Brogan et al. (2018) identified a dampening of the water masers that was caused by the high dust temperature in the protocluster, which does not allow cold dust to be a sink and reduces the maser efficiency. Our observations of MM1-W1 during the accretion burst are consistent with the picture put forth by Brogan et al. (2018). We find that the fading of the water masers happened < 0.1 yr after the onset of the burst (Figure 2). In UCHII-W2, a similar behaviour was observed.

The behaviour of UCHII-W3 differed significantly from CM2-W2, despite both being radiatively connected to MM1B (Section 4.1.1). One explanation could be some unknown effects due to interactions with MM3-UCHII. Another could be the difference in shock types between UCHII-W3 and CM2-W2. UCHII-W3 likely features J-shocks (high-velocity proper motions of 150 km s⁻¹), while CM2-W2 shows characteristics of C-shocks (lower proper motions and R = 2.1). These two shock types have different chemical pathways through which water is produced and population inversion is sustained (Kaufman & Neufeld 1996; Hollenbach et al. 2013). The effects of UV irradiation on water maser emission for these shocks have not yet been explicitly modelled, but it has been argued that 22 GHz water masers from C-shocks are generally brighter and might have stronger responses to thermal and high-energy energy from the burst (Hollenbach et al. 2013; Gray et al. 2024a).

Why did CM2-W2 brighten and trace a larger velocity range after the burst, while UCHII-W1 brightened initially, but had reduced flux density and velocity range after that? The differences may lie in the shock type and source geometry. If there is no low optical depth cavity, then UCHII-W1 only sees the thermal and infrared heat waves. The infrared heatwave is known to dampen masers rather than amplify them. On the other hand, UCHII-W1 also has high-velocity (> 70 km s⁻¹) proper motions that indicate maser emission due to J-shocks.

We can compare the water maser variability in NGC6334I to that of other sources. In S255IR-NIRS3, it was not clear that variability associated with the accretion burst was identified, as no 22 GHz maser monitoring was done before and after the accretion burst (Hirota et al. 2021). The case of G358.93−0.03 is a noteworthy comparison with NGC6334I. In both cases, there was dampening of the masers close (< 1500 au) to the bursting source, as expected from maser theory as explained in Section 1 (Brogan et al. 2018; Bayandina et al. 2022a; Miao et al. 2024). Between the two sources, it is clear that the large scale (> 10 000 au) density structure of the natal envelope modulates water maser variability. The density structure modulates the motion of the thermal heatwave, which is a pulse that can propagate at up to ∼4% of the speed of light (Burns et al. 2020). For G358.93−0.03, water maser associations showed directionally dependent variability, which was attributed to inhomogeneous density distributions around the bursting source (Miao et al. 2024).

Our analysis of NGC6334I reveals three fundamental considerations about water maser response to accretion bursts. The effect of thermal, infrared and high-energy components from accretion bursts affect water maser response separately. Then the intervening density distribution between the bursting source and the maser association modulates the burst energy components. Finally, water masers in C-shocks and J-shocks respond differently to these energetic components. In the next section, we consider different mechanisms of water maser variability.

4.2. Variability mechanisms in 22 GHz water masers

The inverse problem of calculating physical variables (H₂ density and temperature, H₂O abundance, shock velocity, maser path length, incident radiation) from 22 GHz water maser observations (position, radial velocity, flux density and time variability) is non-trivial. How do you know that your maser observations are informative about the underlying physical processes? For water maser observations to be informative, considering the observational degeneracy mentioned above, guidelines for judging the difference between fast approximately random and slow systematic variability effects need to be considered.

4.2.1. Three sources of variability

To illustrate the three sources of variability in the collisionally pumped maser time series, we adopt a simple toy model as an illustration. Consider a cylinder of velocity-coherent H₂O gas with a length s, a radius a, and a large aspect ratio of s ≫ a. The gas has an optical depth τ = κ_νs at 22 GHz. The opacity, κ_ν, is negative for masers, which amplify background radiation. The value of κ_ν for a molecule and volume element is dependent on the response of the level populations to the macroscopic variables of the gas. The maser saturates at a length s_sat and amplifies a background intensity I₀. For s < s_sat, the maser amplification is exponential and for s > s_sat the amplification is linear. In the toy model, we take into account that κ_ν may behave differently in the saturated and unsaturated regimes. We also add random variations in the path length due to chaotic fluid flow in which collisionally pumped masers form, δs = ϵ(t)s, where |ϵ(t)| < ϵ_max is a uniformly random value that changes with time. The distribution of ϵ_max < 1 is unknown. The observed intensity I_ν(s) for s > s_sat + δs(t) is then given by:

$$\begin{array}{c}\hfill I(s,t)=I_0(t)e^{{\kappa }_{\nu },\mathrm{sat}(t)·s_{\mathrm{sat}}}[{\kappa }_{\nu }(t)·(s-s_{\mathrm{sat}}+\delta s(t))].\end{array}$$(4)

The toy model in Equation (4) does not aim to be physically precise, but a few helpful conclusions can be drawn from it. For a maser column with a mean length s, there are three sources of variability, with their corresponding timescales. The first is background variations I₀(t) (hereafter background variability), the second is variations in the maser pump efficiency, measured by κ_ν (hereafter, pump variability), the third is the hydrodynamic component δs(t) which changes the path length of the maser (hereafter hydrodynamic variability, Gwinn 1994a; Ladeyschikov et al. 2022b). Theoretical work has studied maser time series in the case of rotation, line-of-sight overlap and shocks, which all fall under the umbrella of hydrodynamic variability (Gray et al. 2019, 2024a,b). For 22 GHz water masers, positive or negative changes in κ_ν can be induced by changes in the external radiation field, kinetic temperature or density changes of the collisional partner H₂ (Gray et al. 2016, 2022). If we take into account the fact that variations in δs(t) can be faster than the separations between two observations at times t₁ and t₂, we cannot directly use the ratio of the two observations I(t₁)/I(t₂) to draw any conclusions on the source of variability, due to the random component. But, if we have constraints on the timescales of I₀(t), t_background or of κ_ν(t), t_pump and we know that either t_background ≫ t_hydro or t_pump ≫ t_hydro, combined with multiple observations before and after the change in I₀(t) or κ_ν(t), then the changes in the dispersion σ_burst/σ_pre would probe changes in the physical conditions if ⟨δs(t)⟩ = 0. Pump and hydrodynamic variability are not independent in terms of how they change I_ν. According to Equation (4), pump variability can amplify or dampen the effects of hydrodynamic variability. We have reason to suspect that background variability was not the cause of the maser flare in NGC6334I (see Section 4.3.1), so we looked for changes in the dispersion before and after the accretion burst.

4.2.2. Testing variance increase in water maser single-dish time series

We used HartRAO monitoring data of 22 GHz water masers F_ti, Vi = F(t_i, V_i), where F_ν is the flux density in Jy, t_i is the time in days and V_i is the velocity channel. We segmented the data to pre-burst MJD 56500 < t_i < MJD 57080 and burst MJD 57700 < t_i < MJD 58700 data. We ignored the rise time at the onset of the accretion burst between MJD 57080 and MJD 57300. These data were binned into two arrays F_pre and F_burst. The first panel of Figure 11 shows the −7.32 km s⁻¹ HartRAO water maser time series, as well as the time windows used to construct F_pre and F_burst (for details on the observations and data reduction, see MacLeod et al. 2018). We applied a non-parametric permutation test outlined in Tibshirani & Efron (1993) to see whether there was a statistically significant change in the dispersion before and after the accretion burst. The test assumes the data to be uncorrelated over time, and there to be no long-term trend influencing the data variance. In the case of only observing hydrodynamic variability in 22 GHz water masers, we can treat the data as uncorrelated as the path length changes in the masering gas are essentially random. To eliminate variance due to long-term trends, we take data in which the 6.7 GHz methanol masers do not change over time. We also fit a linear function of F_pre and F_burst versus time separately for each velocity channel to remove possible long-term (0.5 yr) linear trends. We stress that the residuals after the subtraction of the linear function are not instrumental noise, as we have S/N = 20−700 in most cases, and the residuals are the intrinsic short-term variability in the water masers. We only are interested in short-term hydrodynamic variability, so removing the linear trends removes variance increase which might be due to pump variability with long timescales. The test statistic t_stat is given by

$$\begin{array}{c}\hfill t_{\mathrm{stat}}={\log }_{10}\left(\frac{{\sigma }_{\mathrm{post}}^2}{{\sigma }_{\mathrm{pre}}^2}\right)·\end{array}$$(5)

Fig. 11. Permutation test of 22 GHz water maser long-term time series taken with HartRAO. Top panel: Example time series of the −7.32 km s−1 velocity time series. The pre-burst and burst time windows are annotated. Second panel: Distribution of tstat for 104 iterations where the array is shuffled for the −7.32 km s−1 time series. The value of the unshuffled time series is shown as the black line. Third panel: Value of the ratio of the unshuffled tstat over the standard deviation of the shuffled distribution. The velocity channels that have a statistically significant increased and decreased dispersion are shown in blue and red, respectively. Bottom panel: Dispersion ratio for channels where there were statistically significant changes. The blue points are σburst/σpre where the dispersion increased, while the red points are σpre/σburst where the dispersion decreased.

The value of t_stat will be positive if the variance increases between the pre-burst and post-burst cases, and negative if the variance decreases. If σ_pre² = σ_post² then t_stat = 0. The middle panel of Figure 11 shows an example of this test for the −7.32 km s⁻¹ time series where there is a clear detection (13σ_t) of increased dispersion. The bottom panel of Figure 11 shows the results of this test for each velocity channel using the same time bins. This test does not suggest significant changes in variance for most of the velocity channels, which gives some confidence to its robustness. The bottom panel of Figure 11 shows the ratio of amplification or dampening between the pre-burst and post-burst periods. We report individual velocity channels as there can be multiple spatially independent maser features in a single velocity feature of a single-dish observation. The velocity channels at −44 km s⁻¹, −7 km s⁻¹ and −4 km s⁻¹ showed an order of magnitude increase in σ_burst/σ_pre. On the other hand, velocity channels at −33 km s⁻¹, −22 km s⁻¹ and −1 km s⁻¹ showed respective dampening of a factor of 26, 15, and 7 in their dispersions. Using current and previous interferometric maps of 22 GHz water masers, we identified the dominant features contributing to σ_burst/σ_pre changes in each velocity range. Figure 3 and Table 2 are helpful for this comparison. The increase in σ_burst/σ_pre for the −7 km s⁻¹ and −4 km s⁻¹ can be attributed to the flaring of CM2-W2, which is the dominant maser feature in the field after 2015. The increase in the variance for the −44 km s⁻¹ might be associated with the high-velocity features in CM2-W2 after 2015.5. On the other hand, the decrease in the variance of the features around −33 km s⁻¹, −22 km s⁻¹ and −1 km s⁻¹ can be attributed to the dampening of features in UCHII-W3, UCHII-W1, and MM1-W1, respectively (Figure 3). These comparisons show that variance changes before and during an accretion burst in 22 GHz water maser monitoring data can be a useful probe of the excitation and destruction of water maser features. Testing on theoretically modelled time series, and long-term monitoring data for other sources help improve our understanding of the strengths and weaknesses of this analysis. These results support our earlier classification of water maser variability types. Pump efficiency variations can amplify or dampen random flux density fluctuations caused by chaotic fluid motion in shocked regions outlined by Gwinn (1994a,b). This occurs through changes in maser depth per path length. This analysis and the VLBI observations reported in this work imply that the variable radiation fields due to the accretion burst caused pump or background variability (dampening and amplification) in the water masers in NGC6334I. In the case of a water maser flare, chance line of sight overlap can be distinguished from pump variability by using the method in this section to check whether or not there are changes in the dispersion of the maser time series. Although the dispersion test possibly indicates pump variability, it does not explain the cause of the variability. We consider the possible sources of non-hydrodynamic variability in Section 4.3.3.

These results also show that single water maser observations in a time series cannot be predicted, but that their statistics hold information on the underlying physics. Time-dependent background radiation and changes in the maser depth are important components, but the stochastic hydrodynamic variability cannot be ignored. The distribution of timescales for hydrodynamic variability would need to be determined from simulations of the complex fluid flow in the shock fronts in which these masers appear. Further, it should be considered that single-dish time series rarely only observe the radiation from a single maser feature and that single maser features have been seen to drift in velocity. Characterizing these effects could enhance the analysis of long-term (> 10 year) single-dish water maser monitoring.

4.3. How did the accretion burst cause the water masers to flare in CM2-W2?

How does the accretion burst lead to a positive correlation between the radiation field strength and flux density of water masers in CM2-W2? Here we consider possible explanations: Amplification of background radiation, radiatively pumped 22 GHz water masers, or radiative heating of the collisional partner.

4.3.1. Change in background continuum

The change in the background continuum can be tested with pre-burst and burst centimetre continuum measurements of CM2-W2. The pre-burst 22 GHz (1.5 cm) continuum was not detected at the 3σ = 0.30 mJy beam⁻¹ level with JVLA (Brogan et al. 2016), and the post-burst 1.5 cm continuum was not reported. On the other hand, for synchrotron continuum emission, the flux density at 5 cm and 1.5 cm are proportional. The 5 cm continuum was detected, and it decreased from 0.36 mJy beam⁻¹ in 2011 to 0.19 mJy beam⁻¹ in 2016.9 (Brogan et al. 2016; Hunter et al. 2018). It does not seem that an increase in the background continuum is responsible for the maser flare in CM2-W2.

4.3.2. Radiatively pumped water masers

It has been proposed that strong infrared radiation could radiatively pump bright 22 GHz water masers (Gray et al. 2022). In this scheme, although water masers are dampened by increased dust temperature in most densities and kinetic temperatures, there is a part of the parameter space, 10⁴ < n_H2O < 10⁶ cm⁻³, T_kin < 500 K and T_dust > 900 K, where maser optical depth scales very quickly with T_dust. Radiatively pumped 22 GHz water masers then correspond to ‘cold-gas hot-dust’ regions. We consider this proposal of radiatively pumped water masers to be insufficient in our case for two reasons. First, CM2-W2 is at the edge of a protostellar jet, which has high T_kin compared to T_dust, as seen in maser morphology, proper motion direction and thermal line emission (Brogan et al. 2018; Chibueze et al. 2021; Hunter et al. 2021). Second, this proposal is valid if there was a significant change in the dust temperature from the pre-burst to burst epochs, but the changes in millimetre dust emission were not observed in CM2-W2 (Hunter et al. 2017).

4.3.3. Radiative heating of H₂

We also consider the possibility that the water masers stayed collisionally pumped, but that the accretion burst increased the maser pump efficiency (κ_ν in Equation 4) enough to cause the maser flare. Analysis of the effect of the accretion burst on masers can be broken up into three parts. The first is the spectrum of the bursting YSO and its radiative transfer to the masering gas. Followed by the changes in the physical conditions of the masering gas, such as heating and chemical effects. Lastly, one can consider the maser properties in the irradiated gas. We suggest what might be the most reasonable explanation given the observations of NGC6334I-MM1 and water maser models. Our analysis here is relatively simple, and we leave more detailed calculations on each of the steps above for future work.

The accretion burst led to a bolometric luminosity increase by a factor of 16 ± 4.4 in MM1 (Hunter et al. 2021). All continuum observations are in the far-infrared to submillimetre wavelengths which only observed reprocessed dust emission (Hunter et al. 2017, 2021; Brogan et al. 2018). We assume that the radiation that terminates at CM2-W2 is weakly attenuated by dust and preserves its original spectral shape (Section 4.1.1). The main spectral components are the photosphere (mainly UV for T_eff = 16 000 K, Hunter et al. 2021) and the accretion disk (infrared and non-thermal X-rays, e.g. de Sá et al. 2019). We rule out infrared radiation as an explanation as it would rather dim the water masers in CM2-W2 as in NGC6334I-MM1 and G107.298+5.639 (Szymczak et al. 2016; Brogan et al. 2018).

High-energy radiation would affect the chemistry and ionization state of the masering gas. Maser models for C-shocks and J-shocks neglect external ultraviolet radiation, so it is unclear what higher-order effects this radiation would have on the masers (Kaufman & Neufeld 1996; Hollenbach et al. 2013). We can see how the radiation might affect T_kin of the masering gas. We use simple energy arguments to make an order of magnitude estimate for ΔT_kin due to the accretion burst. The energy balance equation is

$$\begin{array}{c}\hfill \frac{\mathrm{d}U}{\mathrm{d}t}=\mathrm{\Gamma }-\mathrm{\Lambda }\end{array}$$(6)

where U is the internal energy of the gas in ergs while Γ and Λ are the heating and cooling rates. For simplicity, we assume Γ ≫ Λ. In C-shocks with water masers, the main source of cooling is water maser emission, the dust and H₂ thermal line cooling being weaker by an order of magnitude (Kaufman & Neufeld 1996). We can integrate Equation (6) to get

$$\begin{array}{c}\hfill U=U_0+{\displaystyle {\int }_{t_{\mathrm{pre}}}^{t_{\mathrm{post}}}}\mathrm{\Gamma }(t)\mathrm{d}t\end{array}$$(7)

where U₀ is the energy of the masering gas before the burst, t_pre the time before the burst, t_post the time of the peak 6.7 GHz methanol and t_rise = t_post − t_pre. If we assume that the masering gas has an optical depth so that some constant fraction η of the incident radiation is absorbed as heat, then Γ = ηFdS. The incident F is from the YSO at MM1B, and dS is the area of the maser exposed to F. The brightest masers are perpendicular to the shock propagation direction (Hollenbach et al. 2013). For a cylindrical maser perpendicular to the YSO, dS = ar_maser² where r_maser is the maser radius (∼1 au, Gwinn 1994a) and a is the maser aspect ratio (10 − 100, Hollenbach et al. 2013).

The YSO flux is given by F = L_*/4πR², with L_* the accretion and photospheric bolometric luminosity of the YSO, and R is the distance from the YSO and the maser (2800 au). If we assume an ideal equation of state U = mC_VT_kin, where m = n_H2V = μ_H2n_H2ar_maser³ is the mass of the maser. With μ_H2 the mass of a hydrogen molecule and n_H2 ∼ 10^8 − 10 cm⁻³ the hydrogen number density. The heat capacity per unit mass at constant volume, C_V, can be obtained from the heat capacity per mole at constant pressure (C_P) from Chase (1985) with γ = 7/5 for an ideal diatomic gas and T the kinetic temperature, Equation (7) becomes

$$\begin{array}{c}\hfill \mathrm{\Delta }{T}_{\mathrm{kin}}=\frac{\eta \gamma N_{\mathrm{A}}}{4\pi R^2·n_{{\mathrm{H}}_2}{r}_{\mathrm{maser}}{C}_{\mathrm{P}}}{\displaystyle {\int }_{t_{\mathrm{pre}}}^{t_{\mathrm{post}}}}{L}_{\ast }(t)\mathrm{d}t.\end{array}$$(8)

We can assume that L_*(t) has a time dependence similar to the −7.26 km s⁻¹ 6.7 GHz methanol maser time series (MacLeod et al. 2018). The pre-burst and burst luminosity has been estimated with L_*(t_pre) = 4300 ± 900 L_⊙ and L_*(t_post) = 49 000 ± 8000 L_⊙. For a step function in luminosity, the integral in Equation 8 becomes t_riseΔL_*. Keeping n_H2 and η as free parameters Equation (8) becomes

$$\begin{array}{c}\hfill \mathrm{\Delta }{T}_{\mathrm{kin}}=2.8\times {10}^5\frac{\eta }{n_8}\mathrm{K},\end{array}$$(9)

where n₈ is n_H2 in units of 10⁸ cm⁻³. This shows that even if the (non-maser) cooling Λ is half of the heating Λ = 0.5Γ, and the heating efficiency η is less than one percent η < 0.01, the gas in CM2-W2 would heat by 10³ K.

To estimate η, the spectral shape of the incident radiation is essential as H₂ has wavelength-dependent opacity. Modelling of interstellar gas heated with X-rays suggests that 1 keV X-rays can have heating efficiencies of up to 10% (Glassgold et al. 2012). On the other hand H₂ has a non-negligible opacity to ultraviolet at the densities (n_H2 ∼ 10^8 − 10 cm⁻³) in which water masers are found (Sternberg & Dalgarno 1989). The spectrum of an accretion-bursting HMYSO is still an open question. However, we consider it possible that a bursting HMYSO could radiate enough high-energy radiation to produce the heating efficiency required for a significant increase in T_kin.

There is empirical reason to think that CM2-W2 is being heated up more than the other regions, and that is the detection of 321 GHz water masers at the south of the bow-shock (Figure 10). These masers require high gas temperatures T_kin ∼ 1000 K.

After we consider the change in T_kin, we need to consider the response of the water maser. There are differences in the model predictions to the maser response to gas heating. The J-shock model of Hollenbach et al. (2013) has a weak dependence of maser brightness temperature to gas kinetic temperature. In this model, to have a maser flare in a single feature with a factor of 4, the gas temperature needs to increase by ∼300 K. On the other hand, in the model of Gray et al. (2016) the maser depth responds differently to temperature changes at different densities. For hydrogen number densities 3 × 10⁸ cm⁻³, the maser depth is very weakly dependent on T_kin, while at very high densities between 3 × 10⁸ cm⁻³ and 3 × 10¹⁰ cm⁻³, maser depth changes rapidly between T_kin = 200 K and T_kin = 500 K. This can be most clearly seen in Figure 1 of Gray et al. (2022), where for T_dust = 50 K, maser depth changes by a factor of 2 between T_kin = 250 K and 350 K. As maser depth goes at least linearly (for saturated masers) or exponentially (for unsaturated masers), small temperature changes between 200 K and 400 K at hydrogen number densities between 10⁸ cm⁻³ and 10¹⁰ cm⁻³ can adequately explain the maser flare in the unsaturated case. The calculations of Yates et al. (1997) (independent of shock type) indicate that water masers have a sharp rise in maser depth between 500 and 1000 K, but that additional increase in T_kin after 1000 K would not help the maser action.

These model calculations show that under negligible non-maser cooling and relatively low heating efficiencies, the accretion burst can cause temperature jumps that would amplify the maser emission in CM2-W2 to observed levels. Our explanation for water maser flares in accretion-bursting sources is similar to that of Gray et al. (2024a), although their argument follows different details and they focus on the burst in G358.93-0.03-MM1.

5. Conclusion

We did multi-epoch VLBI observations of 22 GHz water masers in 2014−2016 with VERA. These observations cover the epochs before, during, and after the onset of the accretion burst in NGC6334I. We also report the first detection and imaging of 321 GHz water masers in this source with ALMA. These multi-epoch observations allowed us to track maser evolution throughout the accretion burst, providing unique insights into how the burst affects different maser environments within the complex source structure. Our findings can be summarized as follows:

1. In CM2-W2, associated with a bow shock at the northern edge of the NW-SE jet, we found large-scale (> 50 au) and small-scale (∼1 au) variations in the water masers. The large-scale size increase of the maser region and the velocity-dependent temporal variations suggest a propagating effect from the accretion burst, likely due to increased energy input into the surrounding medium. This expansion and the detection of new maser features at different velocities over time provide direct evidence of the burst’s impact on the maser environment.

2. The masers in MM1-W1 close to the bursting source and UCHII-W1-3 at the southern edge of two protostellar jets show complex variability patterns that reflect the intricate source structure and different shock types. MM1-W1 exhibited negative pump variability consistent with high dust temperatures. UCHII-W3, despite being radiatively connected to MM1B, did not show similar flaring to CM2-W2, possibly due to differences in shock types (J-shocks vs. C-shocks). UCHII-W2 showed gradual dimming, while UCHII-W1 displayed complex behaviour. These varied water maser responses highlight the importance of source geometry, shock types, and the different effects of thermal, infrared, and high-energy components from the accretion burst.

3. We estimated absolute water maser proper motions in the pre-burst cases (2014.72 − 2015.28) and burst cases (2015.88 − 2016.19). Low-velocity v < 40 km s⁻¹ proper motions were detected in CM2-W2 and MM1-W1, likely from non-dissociative C-shocks. Alternatively, the proper motions in UCHII-W1, UCHII-W2, and UCHII-W3 were consistent with high-velocity v > 70 km s⁻¹ dissociative J-shocks. No significant changes were found in the mean velocity of the proper motions. The consistency in mean proper motion velocities before and during the burst provides further evidence that observed changes in maser distribution are primarily due to excitation effects rather than physical motions. The presence of excitation effects underscores the importance of careful interpretation of maser proper motions, especially in dynamic environments experiencing accretion bursts.

4. We used a toy maser model to introduce the concepts of background, hydrodynamic, and pump variability. Our VLBI observations of individual maser features complement this analysis by providing high-resolution spatial information, allowing us to connect changes in maser variance to specific regions within NGC6334I. With long-term single-dish monitoring observations from HartRAO, we identified statistically significant changes in the variance of the time-dependent flux density in certain velocity channels. Variance changes after the burst are associated with changes in the maser pump efficiency, which amplifies hydrodynamic variability. We propose the changes in the variance of a collisionally pumped maser flux density time series as a way to identify pump variability in water masers.

5. We considered possible explanations for the contemporaneous flaring of collisionally pumped 22 GHz water masers and 6.7 GHz methanol masers in NGC6334I. The spatial distribution and intensity changes of masers observed in our VLBI monitoring support the hypothesis of radiative heating from high-energy radiation propagating through cavities with low optical depth. These observations highlight that the maser response of irradiated shocks deserve further theoretical investigation.

In conclusion, our multi-epoch VLBI observations, combined with theoretical analysis, reveal the complex interplay between accretion bursts, source structure, and water maser behaviour in NGC 6334I. Future work in the theoretical investigation of masers in irradiated shocks would be useful for the interpretation of water masers in future accretion bursts. Additionally, the contemporaneous flaring of collisionally pumped and radiatively pumped transitions in bursting sources prompts observational investigation in infrared and high-energy regimes.

Data availability

Full Table 2 is available at the CDS via anonymous ftp to cdsarc.cds.unistra.fr (130.79.128.5) or via https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/691/A157

Acknowledgments

We thank the anonymous referee for thoughtful comments which improved the quality and clarity of this manuscript. JMV acknowledges funding from the National Research Foundation (UID: 134192). JMV thanks Alexander Rawlings for valuable input on statistical methods. JMV and MJ acknowledge support from the Research Council of Finland grant No. 348342. T.H. is financially supported by the MEXT/JSPS KAKENHI grant Nos. 17K05398, 18H05222, and 20H05845. E.I.V. and A.M.S. acknowledge support from the Russian Science Foundation, project No. 23-12-00258. JOC acknowledges financial support from the South African Department of Science and Innovation’s National Research Foundation under the ISARP RADIOMAP Joint Research Scheme (DSI-NRF Grant Number 150551). This paper makes use of the following ALMA data: ADS/JAO.ALMA#2018.1.00024. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), NSTC and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. HartRAO forms part of SARAO, which is a National facility operating under the auspices of the National Research Foundation (NRF), South Africa.

References

Appendix A: Proper motion considerations unrelated to the accretion burst

A.1. Comparison with Chibueze et al. (2021)

The proper motions that we estimated show a characteristic bipolar outflow pattern. Relative proper motions for this source were also estimated by Chibueze et al. (2021) with KaVA. Their observations had a higher angular and spectral resolution and better sensitivity but only covered three epochs a year after the 2015.0 onset of the accretion burst (spanning 2015.89 − 2016.01). On the other hand, our seven observations span times from before the burst (2014.7) until a year after the burst (2016.2). The largest systematic difference between our results and those of Chibueze et al. (2021) is that they did not shift their proper motions to an absolute velocity reference frame. They reported a large (∼120 km s⁻¹) northward transverse speed for CM2-W2 and northward pointing proper motions in UCHII-W1 contrary to our observations and Brogan et al. (2018). This led them to suspect that the MM3-UCHII might be affecting the proper motions of masers in UCHII-W2 and UCHII-W3. This is due to their proper motions being in the reference frame of a bright maser in UCHII-W1, and therefore their relative proper motions should be shifted to an absolute velocity frame to be physically interpreted. They also calculated proper motions for spots (individual velocity channels), rather than maser features (see Section 2.1.1). Their Velocity-Variance Covariance Matrix (VVCM Bloemhof 1993, 2000) estimate of the outflow and inclination angle of the NW-SE outflow is independent of velocity frame, and is therefore robust. Our high-velocity (∼100 km s⁻¹) proper motions for UCHII-W1 are consistent with the bulk motion observed by Brogan et al. (2018) between 2011 and 2017.8.

A.2. Reliability of proper motions

We do not expect the underlying gas motion to be affected by the accretion burst (see the timing argument at the end of Section 4.1.2). Changes in the proper motions should rather be due to changes in pumping conditions, rather than real gas motion. This is analogous but with a different physical mechanism than the propagation of 6.7 GHz methanol masers at high velocities seen by Burns et al. (2020). The water masers’ proper motions in CM2-W2, MM1-W1 and UCHII-W1 show more variation in their vector direction as those in some other sources (e.g. Burns et al. 2016). It has been argued by Goddi et al. (2006) that water maser proper motions reliably trace gas motion. Their argument is empirical, and based on the ordered motion of water masers in G24.78+0.08 over three epochs. Their argument only shows that some (one source) water masers trace gas motion, but not all. The empirical (see Section 3) and theoretical (see Section 4.2) considerations in this work indicate that great care should be taken when interpreting water maser proper motions as physical gas motion. For water maser proper motions to be reliable observable, robust inter-epoch identification methods for maser features should be considered. It should be considered whether the three types of variability (background, pump and hydrodynamic, defined in Section 4.2) are on longer or shorter timescales than the observations.

Appendix B: Supporting Figures

Fig. B.1. Water maser spot positions epochs 2014.7−2016.2. Each horizontal panel shows a different association over time, with the name of the region in the leftmost panel. The spatial scale is constant over epochs for each association. The dates of observation are shown above the top panel. A linear scale is shown in each panel. The colour scale is the same as Figure 2.

All Tables

All Figures

Fig. 1. Overview of the source. Left: ALMA 2015.6 1.3 mm dust continuum from Brogan et al. (2016). The main regions MM1-4 are labelled in white. The peak positions of the bursting source, MM1B, as well as MM2 and MM4 are marked with crosses. The location of the millimetre cavity is marked in red. Right: JVLA 2015.9 5 cm continuum with 2017.8 22 GHz water masers overlaid in blue (Brogan et al. 2018). The positions and names of the water maser associations introduced in Section 1 are in green. The positions and orientations of the N-S and NW jets (solid and dashed-dotted respectively, from Brogan et al. 2018) and NW-SE jet (dashed line, from Chibueze et al. 2021) as well as the 0.5 pc NE-SW outflow (dotted line, from Qiu et al. 2011) are marked in red.

In the text

	Fig. 2. Positions and radial velocities of 22 GHz water masers at 2014.7−2016.2 as observed with VERA. The names of the water maser associations are annotated in the leftmost panel. The linear scale for d = 1.3 kpc is shown at the top left of each panel. The accretion burst started at 2015.0 (after epoch 2), and peaked in 2015.6 (between epochs 4 and 5).
In the text

Fig. 3. Centre velocities and peak intensities of maser features in each association over time. The marker size is proportional to the square of the feature peak intensity ($\propto I_{\mathrm{peak}}^{0.5}$). A scale for the marker size is shown at the bottom left of each panel. The markers are colour-coded according to velocity. The grey horizontal lines indicate the dates of our observations. The date of the onset and peak of the accretion burst traced by 6.7 GHz methanol masers is shown in red and black, respectively (MacLeod et al. 2018).

In the text

Fig. 4. Time dependence of total velocity integrated intensity (line area) and velocity spread for each association from our VERA observations. The first and second panel shows respectively the line area for associations containing increasing and decreasing trends at the onset of the burst. The third panel shows the velocity spread for each association over time. The colours for the third panel are the same as in the legends in the first two panels. The red and black lines indicate the date of the onset and peak of the burst, respectively.

In the text

Fig. 5. Velocity-dependent response of the water maser features in CM2-W2. The radius of the spots is proportional to the peak flux density. The rows indicate different velocity bins as indicated in the panels in the first column. Each column is the observation shown in the first row. Maser features that are outside of the respective velocity bin are shown in grey. The spatial scale for all panels is the same. The offset is in terms of the reference maser feature.

In the text

Fig. 6. Time variation of the brightest maser feature in CM2-W2. Top panels: Offsets and VLSR of maser spots in the brightest feature of CM2-W2 for each epoch. The position, uncertainty and VLSR of each spot is shown. We note the difference in spatial scale between right ascension and declination. Bottom panels: Spectral profile of the maser spots in this feature, with Gaussian fits shown in black. In 2014.9 and 2015.1 the feature is well approximated by a double Gaussian.

In the text

	Fig. 7. Velocity-dependent response of the water maser features in UCHII-W1. The area of the spots is proportional to the square of the peak flux density. The subregions of UCHII-W1, i−iv are marked in the first panel.
In the text

Fig. 8. Proper motions of 22 GHz water masers in the northern regions of NGC6334I, before and during the accretion burst. The positions, orientations and colours of the arrows show the positions, proper motions and VLSR of the maser features, respectively. The colour scale is on the right. The left panels show the pre-burst epochs (2014.72−2015.28) and the right panels show the burst epochs (2015.88−2016.19). The linear scales in position and space are shown in the bottom left corner of each panel. In the bottom panel, the black dashed lines show the position and orientation of the linear features shown in the maser features, and the black arrow shows the displacement of the maser features after the accretion burst.

In the text

	Fig. 9. As in Figure 8, but for the southern regions.
In the text

Fig. 10. Detection of the 321 GHz water maser with ALMA. Left: Integrated intensity map of 2019.7 321 GHz water masers as detected with ALMA in greyscale. The coloured crosses and triangles show the peak positions of 22 GHz water masers detected in 2017.0 and 2017.8, respectively, while the purple contour is the 5 cm continuum point source (Brogan et al. 2018). The ALMA beam and linear scale are shown at the bottom left and bottom right, respectively. Right: Spatially integrated spectrum of the 321 GHz water masers. The small peaks are HCOOCH3 thermal line emission, which is close to 321 GHz.

In the text

Fig. 11. Permutation test of 22 GHz water maser long-term time series taken with HartRAO. Top panel: Example time series of the −7.32 km s−1 velocity time series. The pre-burst and burst time windows are annotated. Second panel: Distribution of tstat for 104 iterations where the array is shuffled for the −7.32 km s−1 time series. The value of the unshuffled time series is shown as the black line. Third panel: Value of the ratio of the unshuffled tstat over the standard deviation of the shuffled distribution. The velocity channels that have a statistically significant increased and decreased dispersion are shown in blue and red, respectively. Bottom panel: Dispersion ratio for channels where there were statistically significant changes. The blue points are σburst/σpre where the dispersion increased, while the red points are σpre/σburst where the dispersion decreased.

In the text

	Fig. B.1. Water maser spot positions epochs 2014.7−2016.2. Each horizontal panel shows a different association over time, with the name of the region in the leftmost panel. The spatial scale is constant over epochs for each association. The dates of observation are shown above the top panel. A linear scale is shown in each panel. The colour scale is the same as Figure 2.
In the text