Free Access
Issue
A&A
Volume 647, March 2021
Article Number A23
Number of page(s) 18
Section Interstellar and circumstellar matter
DOI https://doi.org/10.1051/0004-6361/202039798
Published online 02 March 2021

© ESO 2021

1 Introduction

In the last few years, mass accretion events in high-mass young stellar objects (HMYSOs), which have masses and luminosities of M* > 8 M and L* > 5 × 103 L, respectively, have been discovered via the sudden increase of continuum emission from radio to infrared (IR) wavelengths. One of the first examples was an outburst from the 20 M HMYSO NIRS 3 in the high-mass star-forming region S255 IR (Caratti o Garatti et al. 2017) located at a distance of 1.78 kpc (Burns et al. 2016). At almost the same time, a similar accretion burst event was identified through the millimeter continuum emission from another HMYSO, NGC6334I-MM1 (Hunter et al. 2017). These discoveries suggest that mass accretion burst events are common in HMYSOs, as theoretically predicted for the formation of both low- and high-mass stars (Acosta-Pulido et al. 2007; Caratti o Garatti et al. 2011; Contreras Peña et al. 2017) through disk-mediated accretion (Krumholz et al. 2009; Kuiper et al. 2011).

These accretion bursts triggered a large increase in the flux density of the associated Class II methanol (CH3 OH) masers at 6.7 GHz. For S255 IR, a CH3OH maser flare was reported in November 2015 (Fujisawa et al. 2015; Moscadelli et al. 2017; Szymczak et al. 2018), through pumping by mid-infrared (MIR) radiation from dust grains. This maser flare could be attributed to the increase of the accretion luminosity, estimated to have begun in mid-June 2015 (Caratti o Garatti et al. 2017). Follow-up near-infrared (NIR) observations showed a brightening of the source of ΔK ≃ 2.9 mag and Δ H ≃ 3.5 mag, and an increase in luminosity up to 1.3 × 105 L, which corresponds to a mass accretion rate of 5 × 10−3 M yr−1 (Caratti o Garatti et al. 2017; Moscadelli et al. 2017; Cesaroni et al. 2018). Subsequent submillimeter (Liu et al. 2018) and NIR (Uchiyama et al. 2020) observations also reported a significant increase of the corresponding flux densities. A similar 6.7 GHz CH3 OH maser flare was also detected in NGC6334I-MM1 (Hunter et al. 2018; MacLeod et al. 2018), along with a 22 GHz water (H2O) maser burst (Brogan et al. 2018). Another HMYSO G358.93-0.03-MM1 (Sugiyama et al. 2019; Burns et al. 2020) showed a flare of the 6.7 GHz CH3 OH masers and other new methanol transitions (Breen et al. 2019; Brogan et al. 2019; MacLeod et al. 2019; Chen et al. 2020), a result of therecent extensive observations by the global collaboration network, the Maser Monitoring Organization (M2O; MaserMonitoring.org). The outbursts from these HMYSOs represent a unique opportunity to investigate the accretion, and, in turn, ejection processes in these objects; and to derive useful parameters, such as the mass accretion and ejection rates, infall velocity, and duration of the events.

According to subsequent monitoring observations of S255 NIRS 3, the outburst was also detected at centimeter wavelengths with the Karl Jansky Very Large Array (VLA), although with a delay of ~1 yr (Cesaroni et al. 2018). The slope of the continuum spectrum (~0.78) is suggestive of emission from a thermal radio jet (Reynolds 1986), of which the mass loss rate has been boosted by the accretion burst. The delay between the IR and the radio burst (Cesaroni et al. 2018; Uchiyama et al. 2020) is not surprising, as the former propagates at the speed of light and the latter at the typical speed of the ejecta (500–1000 km s−1). Indeed, Cesaroni et al. (2018) detected possible evidence of the thermal shock produced by the ejecta/wind as the jet recollimates. Assuming an escape velocity of 700 km s−1, the observed time delay corresponds to a traveling distance of ~160–180 au from the source (or ~100 mas).

Alongside the continuum emission observations with the VLA, we have been simultaneously monitoring the H2O maser emission, which is also expected to undergo a variation due to the change of physical conditions in the shock ahead of the jet (e.g., for NGC6334I-MM1; Brogan et al. 2018), where the masers are located (Goddi et al. 2007; Burns et al. 2016). One can note that the richest H2O maser clusters are separated from the continuum peak (the putative YSO position) by 100–200 mas, which is comparable to the expected distance traveled by the material ejected at the time of the IR burst. The “new” ejecta could now be reaching and interacting with the H2O maser environment.

In the present paper, we aim to discuss the relation between the accretion burst in S255 NIRS 3 and the H2O masers associated with this HMYSO. We describe the observations and data analysis in Sect. 2. In particular, we devote Sects. 2.12.3, respectively, to the observations made with the Japanese very long baseline interferometer (VLBI) VLBI Exploration of Radio Astrometry (VERA), the Atacama Large Millimeter/submillimeter Array (ALMA), and the VLA. In Sect. 3, we present the results from the observations with VERA single-dish and VLBI monitoring to investigate possible changes in flux density (Sect. 3.1) and structure (Sect. 3.2) of the H2O maser features, which are compared with the VLA data and previous VLBI observations (Goddi et al. 2007; Burns et al. 2016). The observations prior to the burst are compared with ours to verify a possible change in the H2O maser spatial and velocity structure. At the same time, we also report maser absolute proper motions, and thus their 3D expansion velocities, for which technical information is provided in Appendix A. Furthermore, we give ALMA results for the submillimeter continuum and H2O masers at Band 7 to study the physical properties very close to the central HMYSO (Sects. 3.3 and 3.4). A possible origin of the maser variability and the maser physical properties are discussed in Sect. 4. Finally, we summarize our conclusions in Sect. 5.

2 Observations and data analysis

2.1 VERA

The VLBI monitoring observations of the H2O (61,6 –52,3) line at 22.235080 GHz were carried out with VERA in 2017, as summarized in Table 1. VERA consists of four 20-m radio telescopes located in Japan, with a maximum baseline length of 2270 km, providing a typical beam size of 1.2 mas. The pointing and phase-tracking center position of S255 NIRS 3 was

which is the same as in Burns et al. (2016). All observations lasted nine hours from horizon to horizon, and the typical on-source time was ~6.4 h. Left-handed circular polarization data were recorded on hard disk devices at 1 Gbps. The total bandwidths were 16 and 240 MHz for S255 NIRS 3 and the calibrator, J0613+1708, respectively. The spectral resolution for the maser line was chosen equal to 15.625 kHz, corresponding to a velocity resolution of 0.21 km s−1.

We employed the standard dual beam observation mode (VERA collaboration 2020), as in previous observations of S255 NIRS 3, in which the target and a nearby calibrator (J0613+1708) were observed simultaneously with a separation angle of 0.87 degrees (Burns et al. 2016). The latter was used as a phase or position reference (Reid & Honma 2014) to determine the absolute positions of the H2O masers. During the observations, artificial wideband radio signals were injected into both dual-beam receivers to calibrate the phase difference between them (Honma et al. 2008a). Delay and bandpass calibration were done by observing a strong continuum source, DA193. Then, residual phase calibration was done using J0613+1708. This was detected at a sufficiently high signal-to-noise ratio, with a total flux density of 0.27–0.33 Jy for all epochs.

Correlation processing was done using the software correlator developed at the National Astronomical Observatory of Japan (NAOJ) Mizusawa campus. After correlation, we applied the correction of an a priori delay-tracking model to reduce residual phase offsets during the following calibration procedure (Honma et al. 2008b). Standard calibration and synthesis imaging was performed using the National Radio Astronomy Observatory (NRAO) Astronomical Image Processing System (AIPS) software package (Greisen et al. 2003). We conducted a phase-referencing analysis to determine the absolute position of the reference maser spot1 for each epoch (e.g., Burns et al. 2016; VERA collaboration 2020). The resultant astrometric accuracy is estimated to be 0.5 mas in the worst case according to the absolute proper motion fitting, and this is discussed in the next section. The rms noise levels were 40–60 mJy beam−1 in the line-free channels.

The H2O masers in S255 NIRS 3 have been monitored with the single-dish radio antennas of VERA since October 20, 2016. The typical monitoring interval was two months, each time using in most cases one out of the four VERA antennas. Observations were done in the conventional position switching mode, and amplitude calibrations were done by means of the chopper-wheel method. The spectral resolution of a digital spectrometer was 15.625 or 31.25 kHz.

Table 1

Summary of VERA observations.

2.2 ALMA

We used ALMA cycle 5 data (ADS/JAO.ALMA #2017.1.00178.S) to target the continuum emission in S255 NIR3 and identify the powering source of the 22 GHz H2O masers. Observations were carried out at Band 7 (in the range of 320–340 GHz) on December 02, 2017, and the on-source time was ten minutes. The array consisted of 48 12-m antennas in the C43-7 configuration with baseline lengths of 41–6855 m. Amplitude and bandpass calibrations were done with J0510+1800, and we observed J0625+1440 for the phase calibration.

We used the pipeline-calibrated visibility data delivered by the East-Asian ALMA Regional Center (EA-ARC) to carry out synthesis imaging using the Common Astronomy Software Applications (CASA) package (McMullin et al. 2007). We selected line-free channels from four spectral windows, three with a bandwidth of 1875 MHz and one 469 MHz wide, centered at 335.9, 334.0, 322.5 (broader bands), and 321.2 GHz (narrower band), respectively. The effective bandwidth integrated over the line-free channels was 3593 MHz. Using these line-free channels, a continuum image was produced by means of the clean procedure with a Briggs robust weighting parameter of 0.5. Self-calibration was done using the continuum image to solve only the phase terms. The flux density of the continuum emission increased by a factor of 1.5. The resultant rms noise level was 0.33 mJy beam−1. The synthesized beam size of the continuum emission was 0.073′′ × 0.059′′, with a position angle of − 77 degrees. Pure line data were obtained by subtracting the continuum from the original line and continuum data in the UV domain.

We assigned one of the spectral windows to the H2O 102,9 –93,6 line at 321.225656 GHz (Chen et al. 2000, hereafter referred to as 321 GHz H2O maser). The lower state energy level of the 321 GHz H2O maser transition is 1846 K. The spectral resolution of the 321 GHz H2O maser map was 0.244 MHz, which corresponds to a velocity resolution of 0.23 km s−1. We made a channel map with a velocity resolution of 0.25 km s−1 by regridding the spectral channels. The rms noise level of the line-free channel data was 0.01 Jy beam−1. The synthesized beam size for the 321 GHz H2O maser line was 0.082′′ × 0.070′′, with a position angle of − 73 degrees.

2.3 VLA

The object S255 NIRS 3 was also observed using the VLA of the NRAO during five epochs in 2016. Array configurations were in C-configuration (March 11, 2016), B-configuration (July 10 and August 01, 2016), and A-configuration (October 15 and November 20, 2016). All observations were done at frequencies of 6.0, 10.0, 22.2, and 45.5 GHz to image the continuum emission as described in Cesaroni et al. (2018). In this paper, we analyze only the 22 GHz H2O maser data. The WIDAR correlator was used to record dual circular polarization with a spectral resolution of 15.625 kHz for the C- and B-configurations, and 331.25 kHz for the A-configuration.

3 Results

3.1 The 22 GHz H2O maser spectra

Figure 1 shows the total power spectra of the 22 GHz H2O masers in S255 NIRS 3 observed with the VERA 20-m single-dish antennas. Some of the epochs have multiple independent data from different antennas/dates. In such cases, we plot the spectra that have higher frequency resolution or higher sensitivity in Fig. 1. By comparing the data acquired with different VERA antennas, we conservatively estimate an uncertainty on the flux measurements of ~30%, due to calibration and/or pointing errors. The brightest spectral components can be seen at 6.4 km s−1, as observed during VLBI monitoring (see below), which shows a slight offset from the systemic velocity of 5.2 km s−1 (Wang et al. 2011; Zinchenko et al. 2015). A significant flux increase is seen on April 24, 2017, which implies that such an increase began between February and April, 2017, namely about 21–23 months after the 6.7 GHz CH3 OH maser flare detected on July 07, 2015 (Fujisawa et al. 2015; Moscadelli et al. 2017; Szymczak et al. 2018) and seven-to-nine months after the flux increase of the radio continuum emission in July 2016 (Cesaroni et al. 2018). Other velocity components in the velocity range of 0 to 20 km s−1 are clearly detected and are highly variable.

The H2O maser spectra also exhibit significant changes in their peak velocities. In particular, the redshifted components at ~15 km s−1 show drasticchanges in their spectral profiles. Such variations were observed in a long-term monitoring of this source (Felli et al. 2007). In addition, the 6.4 km s−1 feature seems to drift in velocity in the interval between February 09, 2018 and December 05, 2018. Similar velocity drifts were found in thehigh-mass protostar IRAS 20126+4106, and in that case the motion was interpreted as deceleration (Fig. 5 in Moscadelli et al. 2005).

The auto-correlation and cross-correlation spectra of the 22 GHz H2O masers taken with VERA are shown in Fig. 2. By comparing the flux densities of two nearby epochs, February 23 and 25, 2017, we estimate a flux calibration error of ~30% for the cross-correlated spectrum. However, the spectra of two other epochs, April 19 and 21, 2017 show excellent agreement with each other, within a few %. The cross-correlated flux density is increased from 19 to 56 Jy at the spectral peak around 6.4 km s−1. The total power fluxes obtained from the auto-correlation data are 37–88 Jy at 6.4 km s−1. The ratio of the peak cross-correlation flux density to the total power flux is 60%, and it is almost unchanged during the VLBI monitoring observations for this 6.4 km s−1 feature. On the other hand, the higher and lower velocity spectral features are more dominant in the total-power spectra, while the cross-correlated spectra show marginal detection in some of the velocity components. This implies that they are spatially more extended than the synthesized beam size, and hence resolved out with the VERA baselines. We discuss this issue further in Sect. 4.3.

Figure 3 shows the integrated spectra observed with the VLA in C-configuration (March 11, 2016), B-configuration (July 10 and August 01, 2016), and A-configuration (October 15 and November 20, 2016). The spectra shown in this plot were extracted from the VLA image cubes. The velocity of the strongest feature observed with the VLA is ~5 km s−1, which is significantly different from those observed during the VERA monitoring in 2017–2018. The maser emission associated with NIRS 3 was brightest in the first epoch, on March 11, 2016, and it decreased after July 10, 2016. This decreasing trend in flux densities could be due either to variability of the maser or to the different array configurations used at the different epochs. In the latter case, there could be significant missing flux emitted from spatially extended components, in particular for the A-configuration data.

The picture is further complicated by the presence of a secondary source of maser emission that we discovered about 5′′ northwest (NW) of NIRS 3 at the local-standard-of-rest (LSR) velocity of 0 km s−1. The position of this NW feature is RA(J2000) = 06h12m53.7711s and Dec = +17°59′26.208′′. This feature increased in flux density during the monitoring and became brightest on August 01, 2016. At the maximum phase, the flux density was comparable to that of the main feature around NIRS 3. A more detailed discussion of the flux variability can be found in Sect. 4.2.

thumbnail Fig. 1

Total-power spectra from the VERA 20-m single-dish monitoring observations. A vertical dashed line indicates the LSR velocityof 6.4 km s−1, corresponding to the brightest components during the VERA monitoring (Fig. 2).

thumbnail Fig. 2

Auto-correlation total-power (solid red lines, AC) and scalar-averaged cross-power (blue solid lines, XC) spectra of the H2O masers. Thespectra were computed by averaging all the antennas/baselines and time ranges. Baselines were subtracted by fitting the line-free channels with a polynomial function. A vertical dashed line indicates the LSR velocity of 6.4 km s−1, corresponding to the brightest components.

3.2 VLBI mapping of the 22 GHz H2O masers

The positions of the maser features were determined through VERA observations for seven epochs in 2017, as listed in Table 1. Spatial distributions and radial velocities of the H2O masers detected in at least one of the observed epochs are plotted in the left panel of Fig. 4. Details of the identified maser features are given in Table 2 and Appendix A (see discussion below). The H2O maser features are aligned along the northeast-southwest (NE-SW) direction. The most blueshifted masers tend to be located at the SW part, while the systemic and redshifted components are distributed on the opposite side, toward NE. The brightest feature at 6.4 km s−1, identified as feature ID13, is located in the NE lobe. Around the midway, or slightly south of the midpoint between these two main maser clusters, there are two smaller clusters including the most redshifted components. The position angle of the H2O maser distribution is slightly different from that of the VLA K-band continuum emission (Cesaroni et al. 2018), which traces the radio jet driven by the central protostar.

The proper motions of all maser features were measured by fitting their position offsets as a function of time. In the present study, the monitoring period of about three months was too short to determine a reliable value of the annual parallax. Furthermore, a detailed analysis of absolute motions and annual parallax measurements for the H2O masers was already made by Burns et al. (2016). Thus, we only focus on the internal proper motions of the H2O masers, which trace outflow motions, as we discuss later. For this purpose, we followed the method described in Burns et al. (2016). Details of the adopted procedure are given in Appendix A.

We successfully produced phase-referenced images for five epochs: February 23, February 25, April 19, May 14, and June 09, 2017. The absolute coordinates of the reference maser spot ID19 derived from the Gaussian fitting to the spot map are

for the first epoch, that is, February 23, 2017. The proper motions of all maser features with respect to the barycenter of S255 NIRS 3 are listed in Table 2. Due to the short monitoring period of three months, some proper motions have large uncertainties. However, we can reveal systematic proper motions in the outflow, as shown in Fig. 5. The less ordered proper motions in the NE lobe could be due to interaction with the burst-regenerated radio jet (see Fig. 4a).

Figure 5 shows the maser features detected in all previous works (Goddi et al. 2007; Burns et al. 2016), and in the present study. These maps are registered with respect to the absolute coordinate frame with uncertainties of ~0.5 mas (see Sect. 2.1). The spatial and velocity structures traced by the H2O masers are in good agreement with those of previous observations in 2005 made with the Very Long Baseline Array (VLBA; Goddi et al. 2007) and conducted with VERA in 2010 (Burns et al. 2016), as mentioned above. The features are aligned along the NE-SW direction with a position angle of 43.5 degrees obtained by fitting a line to their positions (solid line in Fig. 5), as explained in Sect. 4.1.

thumbnail Fig. 3

Velocity-integrated flux densities of the 22 GHz H2O masers observed with the VLA. A vertical dashed line indicates the LSR velocity of 6.4 km s−1, corresponding to the brightest components during the VERA monitoring (Fig. 2). The red and blue lines correspond to the spectrum integrated around the central source NIRS 3 and that of the newly found feature offset by 5′′ northwest (NW) of NIRS 3, respectively.

thumbnail Fig. 4

(a) Distribution of the 22 GHz H2O maser features (colored symbols; Table 2) superposed on the VLA K-band continuum observed on December 27, 2016 (purple contours; Cesaroni et al. 2018) and the moment 0 map of the 321 GHz H2O maser emission (gray scale). Contour levels are 4, 8, 16, ... times the rms noise level of the 0.13 mJy beam−1 and the 32 mJy beam−1 km s−1 for the VLA K-band continuum and the moment 0 map of the 321 GHz H2O, respectively. The (0, 0) position is the phase tracking center, RA(J2000) = 06h12m54.006s and Dec = +17°59′22.96′′. The synthesized beam size of the ALMA Band 7 observation is shown in the bottom left corner of each panel. (b) Distribution of peak centroids of the 321 GHz H2O maser features (colored symbols) superposed on the ALMA Band 7 continuum (magenta contours) and the moment 0 map of the 321 GHz H2O maser emission (gray scale). Contour levels for the ALMA Band 7 continuum are 32, 64, and 128 times the rms noise level of 0.33 mJy beam−1, showing only the brightest part around the continuum peak.

Table 2

Positions and proper motions of all features.

thumbnail Fig. 5

Distributions of the 22 GHz H2O masers indicated by the plus symbols with the proper motion vectors measured with respect to the barycenter. The color code indicates the LSR velocities of the maser features. The 22 GHz H2O maser features observed in previous observations by Goddi et al. (2007) in 2005 and Burns et al. (2016) in 2010 are also plottedwith open circles and filled squares, respectively. The solid line indicates the outflow axis derived from a linear fit to the positions of all maser features, including previous results.

3.3 The ALMA Band 7 continuum emission

As shown in the right panel of Fig. 4, the peak of the submillimeter continuum emission is located on the outflow axis and is close to the central redshifted maser features. The submillimeter continuum peak is consistent with that observed previously with the Submillimeter Array (SMA) and ALMA, identified as SMA1 (Wang et al. 2011; Zinchenko et al. 2015; Liu et al. 2018). It is also coincident with the radio continuum emission peak observed with the VLA at various wavelengths (see the left panel of Fig. 4; Moscadelli et al. 2017; Cesaroni et al. 2018). The intensity, flux density, peak position, and size of the submillimeter emission are determined by 2D Gaussian fitting on the image, as listed in Table 3. The size is slightly larger than the synthesized beam, and the deconvolved size is 57.9 mas × 53.6 mas with a position angle of 89 degrees. The deconvolved size is not elongated in the direction of the outflow, unlike the radio continuum emission tracing the newly ejected jet (Cesaroni et al. 2018). Thus, the submillimeter continuum source could trace the dust emission from an envelope and/or disk associated with the driving source of the NE-SW outflow.

The flux density of the submillimeter source at ~330 GHz was measured on December 02, 2017, after previous observations with the VLA, SMA, and ALMA (Zinchenko et al. 2015; Moscadelli et al. 2017; Cesaroni et al. 2018; Liu et al. 2018). The peak intensity and integrated flux density of the previous observations at 335.4 GHz were 0.14 Jy beam−1 and 0.41 Jy, respectively, on July 20, 2017 (Liu et al. 2018). The peak intensity of our measurement at almost the same frequency is ~2 times lower than that of Liu et al. (2018). It is most likely that the twice higher resolution ( ~0.066′′) with respect to the previous ALMA observation (0.14′′) can justify the different flux values. In fact, the integrated flux density over the 0.5′′ aperture is measured to be 440 mJy for our data, consistent with that from the previous observations (Liu et al. 2018). In order to compare these results, we also produced a continuum image with the same deconvolved beam size of 0.14′′ as the previous observations by Liu et al. (2018). In this case, the peak intensity and integrated flux density are 127.4 ± 0.8 mJy beam−1 and 230.4 ± 2.1 mJy, respectively. The integrated flux density within the 0.5′′ aperture (420 mJy) does not depend on the resolution. The larger value of the integrated flux density than that of the Gaussian fitting (Table 3) is attributed to the diffuse extended emission around the compact continuum peak. Given the flux calibration accuracy and different UV coverage, our results suggest that the submillimeter continuum emission does not show any significant variation between July 20 (Liu et al. 2018) and December 02, 2017 (present study).

Table 3

Summary of ALMA Band 7 continuum data.

3.4 The 321 GHz H2O masers

We also observed the 321 GHz H2O masers with ALMA. Although this strong H2O maser emission was identified in single-dish radio observations 30 yr ago (Menten et al. 1990), this type of maser has been observed with interferometers in only three other star-forming regions: Cepheus A (Patel et al. 2007), W49A (Kramer et al. 2017), and Source I in Orion KL (Hirota et al. 2014a). S255 NIRS 3 is the fourth star-forming region for which the 321 GHz H2O masers have been imaged. Figure 6 shows the spectrum of the 321 GHz H2O maser integrated over a 0.5′′ × 0.5′′ region. There are multiple spike-like velocity components in the range − 5 to 25 km s−1 covering a slightly larger velocity range than that of the 22 GHz H2O maser spectra (Fig. 2). We note that the broader feature from − 10 to +20 km s−1 could be due to multiple HCOOCH3 lines (Hirota et al. 2014b).

The right panel of Fig. 4 shows the distribution of the 321 GHz H2O maser features superposed on the moment 0 map of the 321 GHz H2O maser line, and the continuum emission at 1 mm obtained from the ALMA Band 7 observation. The centroid positions of all features were determined via Gaussian fitting of the emission at the peak velocities. The moment 0 map shows three peaks: one at the continuum position, and one at each of the NE and SW lobes, which are also detected in the 22 GHz H2O maser line. Therefore, the 321 GHz H2O masers trace the outflow, as suggested for other HMYSOs (Patel et al. 2007; Hirota et al. 2014a). It is worth noting that there seems to be a velocity gradient in the 321 GHz H2O maser spots close to the central continuum peak. The direction is perpendicular to the outflow axis, and it is consistent with those detected in other thermal molecular lines (Zinchenko et al. 2015). This gradient could trace a rotating disk/envelope or the base of the outflow (Hirota et al. 2014a). The 321 GHz maser detections are given in Table 4, where we list positions, velocities, and flux densities for the spots at the peak velocity in individual features. Altogether, we identified 11 features at different positions and velocities. Among them, five features correspond to the positions of the 22 GHz H2O maser features detected on May 14, 2017. The other six features, including the most redshifted and blueshifted spots, have no counterpart in the 22 GHz H2O maser features, as in the case of Cepheus A (Patel et al. 2007). The peak brightness temperature Tb of the 321 GHz H2O emission is close to 104 K for the brightest components (see Table 4). This value is to be taken as a lower limit if the emission is unresolved and clearly proves that we are observing maser emission in this water line, although we cannot rule out the possibility of thermal emission for the features with the lowest peak brightness temperatures (Tb = 942 K; Table 4) and for the broad/weak spectral components (Fig. 6).

Table 4

Summary of the 321 GHz and 22 GHz H2O maser features.

thumbnail Fig. 6

Cross-power spectrum of the 321 GHz H2O maser integrated over the 0.5′′ × 0.5′′ region.

4 Discussion

4.1 Overall distribution of the 22 GHz H2O masers

The position angle of the NE-SW outflow traced by the 22 GHz H2O maser features was derived to be 43.5 ± 1.1 degrees (1σ error) by fitting a line to the maser positions on the map. This is slightly different from that of the radio jet observed on December 27, 2016 (Cesaroni et al. 2018). On the other hand, the 22 GHz H2O maser distribution coincides well with the 321 GHz H2O masers. The difference between the radio jet and maser orientations could be interpreted as the result of different outflow/jet directions in multiple/episodic ejection events; the masers could be excited in the shock front away from the central protostar produced by a previous ejection episode. One could even speculate the mass accretion burst could have been responsible for changing the disk orientation.

We find that the post-burst water maser positions and proper motions described in this work are basically the same as those in the pre-burst epochs (Goddi et al. 2007; Burns et al. 2016). The SW bow shock is persistent from the first VLBI observations in 2005 to the present observations from 2017 (Fig. 5). As found by Burns et al. (2016), the shock front is moving toward the SW away from the central HMYSO. Figure 7 shows the distribution of the maser features rotated by 43.5 degrees from north to east, where the abscissa corresponds to the outflow axis fit on the image. The distribution of maser features in the SW lobe extends over 78, 87, and 110 mas along the outflow axis in 2005 (Goddi et al. 2007), 2010 (Burns et al. 2016), and 2017 (present study), respectively (horizontal bars in Fig. 7). In the present study, there are redshifted features in the SW lobe that only appeared in 2017 at the trailing part of the lobe. If these two features are excluded, the distribution of the maser features in the SW lobe is smaller in 2017 than in previous epochs. On the other hand, the positions of the head of the bow shock measured as the minimum value along the outflow axis are moving significantly from the center. The position offsets between two consecutive epochs are 13 and 25 mas for 2005–2010 and 2010–2017, respectively. These offsets are considerably larger than the astrometric accuracy of VERA and the VLBA. The proper motion of the head of the bow shock is 3.26 ± 0.27 mas yr−1 (1σ) toward the SW, which corresponds to 28 km s−1. This value is in good agreement with those of the maser features given in Table 2. If the bow shock is propagating from the continuum peak position without acceleration or deceleration, the dynamical timescale of the outflow traced by the H2O masers is about 60 yr. Burns et al. (2016) estimated a dynamical timescale of <130 yr for the NE-SW outflow, based on the proper motion measurements of the H2O maser features in two epochs between 2005 and 2010. We improve on their result using three epochs of data from 2005 to 2017, and our new value is consistent with the upper limit reported by them.

thumbnail Fig. 7

(a) Distributions of the 22 GHz H2O maser features obtained after rotating the plot in Fig. 5 by 43.5 degrees. The abscissa is parallel to the outflow as indicated by a solid line in Fig. 5. Horizontal bars at the bottom right corner indicate the minimum, maximum, and mean positions of the SW lobealong the outflow axis for each epoch (2005, 2010, and 2017 from bottom to top). (b) Same as (a), where the displacement along the Y axis has been replaced by the observing epoch. The black solid line represents the movement of the head of the bow shock along the outflow axis.

4.2 Maser flux variability

Flux variations are very common in H2O maser sources, and S255 NIRS 3 is no exception to this rule (Felli et al. 2007; Goddi et al. 2007). This is also evident in previous VLBI observations in the pre-burst phase. For instance, Burns et al. (2016) reported that one of the redshifted features labeled J, located close to the continuum position, shows a steady increase from 4 to 82 Jy, while feature G in the SW blueshifted bow shock exhibits a gradual weakening in the observations with VERA from 2008 to 2010. The other maser features also show substantial variations (see Fig. 1 in Burns et al. 2016). In addition, the maximum flux densities of 15 and 20 Jy observed, respectively, in the NE and SW lobes by Goddi et al. (2007) are significantly different from those of the present study, although such a difference could be due to the different angular resolutions.

Flux variations of all maser features detected with VERA are shown in Appendix A. Although there is no overall trend in variability, about half (11) of the identified features increase their flux densities, including the brightest feature ID13.

We plot the light curves of the H2O masers in Fig. 8. For the single-dish monitoring data, we simply integrated them over the whole velocity range where emission is detected, namely from − 15 to + 25 km s−1. We note that the newly found NW component (see Fig. 3) could contribute to the total flux by less than 10%, if we assume that the 0 km s−1 spectral feature belongs to this NW component. Furthermore, the 0 km s−1 feature is not detected in the VERA monitoring (Fig. 2), and thus cannot affect our variability study of the masers in NIRS3.

The VERA single-dish monitoring suggests a gradual increase in the total flux density after the end of the VLBI monitoring period in 2017, but subsequently, during 2018, the total flux remains almost constant (see Fig. 8a). During VLBI monitoring in the first half of 2017, integrated spectra obtained from the VLBI image cubes follow a trend similar to that of the single-dish monitoring. These results suggest a presence of spatially extended emission components resolved out with VERA. We note that the total flux of both the extended and unresolved compact features increase with time. The existence and variability of an extended maser component is further discussed in Sects. 4.3 and 4.4.

In Fig. 9, we plot the integrated flux densities extracted for the NE and SW outflow lobes and from the central region close to the continuum peak, of which the spectra are shown in Fig. 10. The flux densities are obtained by integrating the VERA spectra over the entire velocity range in which emission is detected (i.e., from − 15 to + 25 km s−1), with 30% uncertainties as discussed in Sect. 3.1. The NE lobe clearly dominates the total flux, whereas the central region and the SW lobe contribute to only about 30% and 10% of the total emission, respectively. This could imply that the interaction between the radio jet and the surrounding environment is occurring mostly in the NE lobe. Significant flaring activity is detected in no regions, and in all cases a moderate flux increase, at most by a factor of 2, is seen during the VERA monitoring period. The auto-correlation spectra show the same trend.

Prior to the VERA monitoring, we observed the H2O maser spectra with the VLA in 2016. When compared with the VERA results (Fig. 2), the spectral profiles observed with the VLA (Fig. 3) look significantly different, as we can see plotted in Fig. 11. In the VLA spectra, only weak emission at 6.4 km s−1 is detected, whereas at this velocity one finds the brightest emission during the VERA observations. The 5 km s−1 feature becomes less prominent at the end of the VLA monitoring and it can hardly be seen in the VERA spectra, although there may be a corresponding feature slightly more blueshifted at a few km s−1. The most redshifted 15 km s−1 feature is gradually increasing its flux, as is clearly seen during the VERA monitoring and the total-power spectra especially.

Figure 8a seems to show a flux increase of the H2O masers between mid 2016 (July and August) and late 2016 (October and November). In the latter two epochs, the VLA was in a more extended configuration (A-configuration) than in the former (B-configuration), which rules out the possibility that this increase could be due to part of the flux being resolved out in the most extended configuration. However, the increase is not very prominent, and the flux appears to decrease between October and November. Furthermore, the VLA spectra in October and November are partly resolved out compared to the corresponding single-dish measurements (Fig. 11), and the total flux density of the VLA C-configuration data on March 11, 2016 is almost at the same level as those observed with the VERA single-dish monitoring in late 2016. We conclude that the total flux has probably undergone some variation between March and November 2016, but a true persistent increase was only observed during 2017, as proven by our single-dish monitoring (blue line in Fig. 8a).

thumbnail Fig. 8

Light curves of maser lines and continuum emissions. A red vertical line indicates the start of the 6.7 GHz methanol maser flare (July 07, 2015) reported by Fujisawa et al. (2015). Dashed lines indicate the dates of January 01, from 2015 to 2019. (a) Flux densities of the H2O masers integrated between − 15 and + 25 km s−1, from the VLA and VERA observations. The VERA results include the single-dish monitoring and the integrated spectra extracted from the VLBI image cubes. The VLA data were extracted from the image cubes, excluding the NW component. (b) Flux densities of representative 6.7 GHz methanol maser features observed with the Hitachi 32-m radio telescope (Yonekura et al. 2016), showing different behaviors among the light curves. (c) IR light curve at Ks band observed with the 1.5-m Kanata telescope in Hiroshima, Japan (Uchiyama et al. 2020). (d) Variation of the radio continuum emission at 22 GHz observed with the VLA (Cesaroni et al. 2018) and ALMA Band 7 at ~330 GHz (Liu et al. 2018, present paper). For the 330 GHz continuum emission on December 02, 2017 measured in the present study, we used the value obtained by integrating over a 0.5′′ aperture (420 mJy) rather than the result of the Gaussian fit (see discussion in Sect. 3.3).

thumbnail Fig. 9

Flux densities of the H2O maser features in the whole region around NIRS 3, the NE lobe, central (continuum) position, and the SW lobe observed with VERA (the corresponding spectra are shown in Fig. 10). The flux densities were computed by integrating the emission over the velocity range between − 15 and + 25 km s−1. The error bars represent the 30% flux uncertainty. A dashed vertical line indicates January 01, 2017.

thumbnail Fig. 10

Spectra of the H2O maser features integrated over different regions in the image cubes from the VERA monitoring data: (a) the whole emission region around NIRS 3; (b) the NE lobe; (c) the central (continuum) position; and (d) the SW lobe.A vertical dashed line indicates the LSR velocity of 6.4 km s−1, corresponding to the brightest components. Panel b: one can see the profiles of the weakest components at ~15 km s−1, while the brightest features at 6.4 km s−1 are better seen in panel a.

thumbnail Fig. 11

Comparison between VERA single-dish and VLA spectra observed at almost the same time in October 2016 and November (VLA)/December(VERA) 2016. Black, red, and blue solid lines show, respectively, the total-power spectrum from the VERA single-dish monitoring, the integrated VLA spectrum of NIRS 3, and that of the newly found maser source NW of NIRS 3.

4.3 Extended H2O maser emission

We already highlighted a difference between the total flux recovered by the VERA interferometer and the corresponding single-dish flux. This finding hints at the existence of extended H2O maser emission. Our VLA observations seem to confirm this idea, because the flux measured with the VLA decreases significantly as the array configuration becomes more and more extended (from C to B, to A), as one can see in Fig. 3. While such an effect could be the result of intrinsic variability of the maser emission, the comparison between the VLA spectrum obtained on October 15 and November 20, 2016 with the VERA single-dish spectra observed at almost the same time (on October 20 and December 19, 2016, see Fig. 11) proves that 40–50% of the maser emission is indeed resolved out by the A-array configuration of the VLA. This implies that such an extended emission must arise from a region of at least ~1′′ (1780 au). To our knowledge, this is the first time that water maser emission has been resolved out with the VLA in a star-forming region. This finding raises questions about the actual origin of the emission, as collisions in shock fronts appear inadequate to pump the 22 GHz line over such a large scale. We discuss this issue and propose a possible explanation in Sect. 4.4.

4.4 Origin of the maser flux increase in relation to the accretion burst

While our observations have revealed no strong flare in the 22 GHz H2O maser line, Fig. 8a clearly shows that a steady increase of the maser flux was present for almost a year, starting from the beginning of 2017. In the following, we attempt to explain the origin of such an increase.

In Figs. 8c and d, we plot the total flux emitted by S255 NIRS 3 at NIR (Uchiyama et al. 2020), submillimeter (Liu et al. 2018, and present results), and radio wavelengths (Cesaroni et al. 2018) on the same timescale used in Fig. 8a for the H2O masers. Thecomparison between the different light curves is very informative. As noted by Cesaroni et al. (2018), the time lag between the NIR and radio bursts indicates that the latter cannot be directly related to the former, due to the long time lag (≳7 months) between the peak of the NIR burst and the beginning of the radio brightening. Cesaroni et al. (2018) demonstrated that the increase of the centimeter radio flux can be due to the expansion of a thermal jet from NIRS 3. Vice versa, the submillimeter flux increase occurs prior to that at radio wavelengths, and the delay of <4 months with respect to the NIR peak is consistent with the heating of the dusty core enshrouding NIRS3 with IR photons.

It is suggested that the 22 GHz H2O masers show anticorrelation with the CH3OH masers. For instance, an intermediate-mass protostar in G107.298+5.639 exhibits alternating flares of the H2O and CH3OH masers with a period of 34.4 days, even though these two masers spatially coexist around the same protostar within 360 au (Szymczak et al. 2016). Although the detailed mechanism of the anticorrelation between the H2O and CH3OH masers in G107.298+5.639 is still unclear, it could be caused by the different physical conditions required for these masers. One possible explanation is that the increase of the NIR emission has opposing effects on the intensities of the H2O and CH3OH masers. The 22 GHz H2O maser could be quenched by the NIR radiation, while the 6.7 GHz CH3OH maser could be radiatively pumped by NIR emission, as suggested by Szymczak et al. (2016). However, this is probably not the case for S255 NIRS 3 because the flare of the CH3OH masers does not occur at the same time as the fading of the H2O masers (Moscadelli et al. 2017), and these two different masers could reflect different physical conditions in different locations. In fact, the CH3OH maser flux is not always anticorrelated with that of the H2O maser, as shown in Figs. 8a and b. Furthermore, it would be difficult to discuss any anticorrelation of these two masers in detail, such as the dimming of the H2O maser at the CH3OH flares, because of a lack of simultaneous observational data around the maximum of the CH3OH maser in late 2015. Future simultaneous monitoring of the CH3OH and H2O masers will be key to revealing a possible anticorrelation of these two masers in the maser flare events.

Before discussing the origin of the increase of the H2O maser flux, it is important to remember that this is mostly due to features lying in the NE lobe (feature ID13 at 6.40 km s−1 and ID22 at 15.96 km s−1s) at a distance of a few 100 au from the central star (see Table 2, Figs. A.9 and A.15). This implies that the flux of these features could be affected by the expansion of the radio jet once the latter has reached the location of the NE masers at approximately 300 au from the star. According to Cesaroni et al. (2018), this should occur at the end of 2016, which is consistent with the beginning of the increase of the maser flux. Despite the temporal coincidence between the two events, their relationship is not fully established, as one should explain why only a moderate and slow increase in the maser intensity is seen, instead of a strong flare. Another problem is that the H2O masers do not lie in front of the expanding radio jet imaged by Cesaroni et al. (2018, as plotted by purple contours in Fig. 5a), as one would expect if such masers were excited through collisions in a bow shock.

The previous considerations apply to the compact maser features, but it is also interesting to find out what happens to the extended maser emission resolved out in our VERA and (to some extent) VLA observations. To study the variation of this component, we computed the difference between the integrated flux of the autocorrelation spectrum and that of the spectrum extracted from the VERA image. The result is plotted as a function of time in Fig. 12 and clearly shows that the flux increase is also seen in the extended maser component. On one hand, this result suggests that both the extended and compact maser emissions are excited through the same mechanism. On the other hand, it makes it difficult to believe that the increase of maser intensity can be due to the radio jet impinging on the H2O gas, becausethe jet is too collimated to trigger maser emission over an extended region of a few arcseconds or ~ 104 au.

In conclusion, the lack of a true maser flare and the presence of extended maser emission undergoing the same increase as the compact maser features cast some doubt on the idea that the maser flux increase can be directly triggered by collisions in the radio jet. While it is possible that the water maser variation observed in our monitoring is not related at all to the outburst observed in other tracers, we propose a solution of the conundrum in the wake of a similar event observed in another HMYSO. Brogan et al. (2018) observed a strong maser flare associated with an accretion outburst in NGC6334I-MM1, with a delay of two years after the submillimeter burst. They proposed that the maser flare is caused by the radiation of the outburst propagating along the outflow cavity and terminating at the bow shock described by the masers. A similar scenario might apply to our case, where the radio jet could be piercing the molecular core, thus opening a free path along which the photons of the (already declining) IR outburst could escape and reach the shock where the masers are excited. In this scenario, the interplay between the expanding jet and the declining IR emission (see Fig. 5) could justify the lack of a sudden and strong increase of the maser intensity. At the same time, maser excitation through radiation could also explain the presence of maser emission over an extended region. A possibility of radiative excitation of the H2O masers is supported theoretically by the model of Gray et al. (2016), in which the 22 GHz H2O maser is both collisionally and radiatively pumped.

thumbnail Fig. 12

Difference between auto-correlation spectra and those from VLA/VERA images. Left and right panels: flux densities integrated over all velocity components and differential spectra. Left panel: error bars correspond to 30% uncertainty, and the dotted vertical line indicates January 01, 2017. Right panel: dotted vertical line marks the systemic velocity.

4.5 Physical properties of the shocked region

The physicalproperties of the H2O maser-emitting region can be estimated from the ratio between the luminosity of the 22 GHz maser and that of the 321 GHz masers, RL22L321, for maser features with the same position and velocity (Neufeld et al. 1990; Patel et al. 2007). Assuming that the beaming angle is the same for both maser transitions, R is also equal to the flux density ratio, R = F22F321. Table 4 compares the positions, peak velocities, and peak flux densities of the five features with both 321 and 22 GHz H2O maser emission. The variability of the 22 GHz H2O masers was established on the basis of Figs. A.1A.16. From these, we estimate a maximum uncertainty on the flux density, due to variability, of about a factor 2. We note that the faint features ID23 and ID25 were detected only in one epoch, and, consequently, for these two features we can only derive an upper limit for the ratio R. The ratio R ranges from0.5 to 14.2 for the five features, where the lowest value is found for the redshifted ( ~17 km s−1) feature (J) located close to the continuum peak. The features in the blueshifted SW lobe (D) and the redshifted feature located to the SW of the continuum peak (H) have a ratio of ~2. On the other hand, the other two features in the NE lobe (E) and closer to the continuum peak (F) have the highest value of R ~14.

According to Fig. 3 of Neufeld et al. (1990), features with low R values have higher temperatures of >1000 K. This implies that the continuum peak and SW bow shock should have higher temperatures ( >400 K) than that of the NE lobe. More recent theoretical models suggest that the 321 GHz H2O masers tend to be excited most efficiently at a hydrogen density of ~ 108 cm−3, which is slightly higher than that needed to excite the 22 GHz masers (Fig. 8 in Richards et al. 2014). Similar physical conditions were predicted by Gray et al. (2016), in which the 22 GHz H2O maser is also pumped in a lower H2O density region of ≤104 cm−3.

The long-term VLBI monitoring of more than a decade suggests that the SW maser features are more stable than the others. This behavior is the opposite to that of the radio jet, where the structure of the NE free-free knot is more stable than that of the SW knot (see Fig. 1 of Cesaroni et al. 2018). We speculate that the inhomogeneous density structure of the envelope and/or outflow cavity might affect the shape of the smaller scale outflow/jet. The higher density in the SW direction could produce a stronger maser emission at the bow shock. This scenario seems consistent with the lower values of R. In contrast, the outflow cavity in the NE lobe could favor a more prominent radio jet structure due to the lower density envelope or outflow cavity, possibly created by a previous outburst. Although the NE masers are also stable, the main difference with respect to the SW lobe consists of more scattered positions and line-of-sight velocities. This could be due to the ambient medium in the NE lobe being less homogeneous, although not necessarily less dense.

There are several 321 GHz H2O maser features not associated with the 22 GHz masers, with R values close to zero ( <1). It has been reported that the majority of the 321 GHz maser features in Cepheus A do not coexist with 22 GHz features (Patel et al. 2007). In the case of the red supergiant VY CMa, long-baseline observations with ALMA revealed that the 321 and 22 GHz masers do not always overlap (Richards et al. 2014). In these cases, the discussion based on the emissivity ratio R (Neufeld et al. 1990) cannot be applied. We cannot rule out the possibility that in our case the non-detection of the 22 GHz H2O masers and the low value of R are due to the variability of the 22 GHz masers and/or the higher spatial resolution of the VLBI data compared to the ALMA images, which could resolve the 22 GHz H2O maser features. In fact, possible signatures of velocity components other than 3, 6, 13, and 18 km s−1 appear in the single-dish spectra from October 2017 to February 2018, showing significant variation (Fig. 1). Further coordinated observations with ALMA and VLBI (and the VLA as well) are needed to address this issue.

5 Summary

We carried out monitoring observations of the 22 GHz H2O masers withboth single-dish telescopes (VERA 20-m antennas) and interferometers (VLA and VERA) after the mass accretion burst event in the high-mass protostar S255 NIRS 3. The H2O maser flux density shows a gradually increasing trend in 2017, about two years after the outburst in mid 2015, and it becomes almost stable in 2018. The flux increase is mainly seen in the maser features associated with the NE outflow lobe or the continuum source.

We found that a part of the H2O maser emission is resolved out even with the A-configuration of the VLA. To our knowledge, this is the first time that extended (> 1700 au) maser emission has been detected in a star-forming region. Both the compact features and the extended component show a similar flux variation over time, suggesting a common origin and excitation mechanism, which we propose to be radiative excitation due to the combined effect of the IR outburst and the expanding jet.

The H2O masers associated with the SW bow shocks seem to be unchanged compared with the previous VLBI monitoring (Goddi et al. 2007; Burns et al. 2016). The measured proper motions of the H2O masers in this SW bow shock suggests a dynamical timescale of 60 yr, which is much longer than the current outburst.

In addition, we conducted followup observations of the submillimeter continuum and the 321 GHz H2O masers with ALMA at Band 7. The outflow structure is also traced by the 321 GHz H2O masers and is slightly tilted compared with the radio jet powered by the mass accretion burst (Cesaroni et al. 2018). The ALMA Band 7 continuum does not show significant variation compared with the previous observations (Liu et al. 2018).

Considering all the above characteristics, we suggest that the gradual increase of the 22 GHz H2O maser flux in S255 NIRS 3 could be caused by the expanding radio jet. The change in the flux density seems to be more moderate compared with another case of mass accretion burst events in NGC6334I-MM1 (Brogan et al. 2018). As already suggested (e.g., Burns et al. 2020), mass accretion bursts in high-mass star-forming regions appear to be quite different from each other. Further surveying and monitoring of these kinds of sources/events, such as the M2O program, are needed to shed light on the mass accretion processes in high-mass star-formation.

Acknowledgement

This paper makes use of the following ALMA data: ADS/JAO.ALMA#2017.1.00178.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), NSC 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. T.H. is financially supported by the MEXT/JSPS KAKENHI Grant Number 17K05398. K.S. is financially supported by the MEXT/JSPS KAKENHI Grant Number 19K03921. Data analysis were in part carried out on common use data analysis computer system at the Astronomy Data Center, ADC, of NAOJ.

Appendix A Proper motion measurements and fitting results

First, we measured the absolute position and proper motion of one of the selected maser spots in the feature ID19 at ~14 km s−1 (Table 2), which had stable spatial and velocity structures during our observation period. We could obtain phase-referenced images for five epochs from seven observing sessions (February 23, February 25, April 19, May 14, and June 09, 2017). The derived absolute coordinates of this maser spot are

for the first epoch (February 23, 2017), in which the error bars are formal uncertainties in the Gaussian fitting, and hence do not include systematics due to calibration errors. We compared the absolute positions of this maser spot in the first and second epochs, February 23 and 25, 2017, which are separated by only a two-day interval, and we found that the position difference between these two epochs is 0.03 mas. We regard this position difference as the typical uncertainty in the absolute positional measurements, which is taken into account in the absolute proper motion fitting as a systematic error.

We used the results from all five epochs in which we succeeded in phase-referencing astrometry to obtain the absolute proper motion of this spot by fitting the position movement as a function of time with the fixed annual parallax value of 0.563 mas (1.78 kpc; Burns et al. 2016). The best-fit absolute proper motion of this spot in the feature ID19 is − 1.08 ± 0.13 mas yr−1 in right ascension and − 0.42 ± 0.13 mas yr−1 in declination, with a post-fit residual of 0.04 mas in both coordinates.

Next, we performed fringe fitting and self-calibration for the strongest maser spot in the feature ID13 at ~6 km s−1 (see Table 2) to achieve a better phase calibration than that obtained in phase-referencing mode. In this procedure, the absolute position information was lost and the reference maser spot in the feature ID13 was placed at the (0,0) position, which is the phase-tracking center. Thus, the absolute position of the ID13 feature was registered by measuring its position offset from the spot in feature ID19 with absolute position known from the phase-referencing astrometry measurements.

We only measured the proper motions of those features that were detected in more than three epochs. The position of a feature at each epoch was determined by the intensity-weighted average of all maser spots in the feature. As a result, we determined the proper motions for 15 features among the 25 identified, as shown in Table 2. The average post-fit residual of the proper motion fitting for each feature is 0.12 mas (0.14 mas in right ascension and 0.08 mas in declination). The residuals are larger than those of the absolute position measurement for the maser spot in the ID19 feature (0.04 mas). This is most likely due to structural changes in the maser features, which are sometimes seen in VLBI astrometry (e.g., VERA collaboration 2020).

In order to estimate the systemic motion of the outflow with respect to the barycenter, we calculated the average proper motions of the maser features associated with the NE and SW lobe of the outflow. As done in Burns et al. (2016), we excluded the redshifted maser features located between the NE and SW lobes because they could not trace the NE-SW outflow. If the outflow motion is symmetric with respect to the barycenter, these motions should point toward the opposite direction at the same velocity. Finally, by subtracting the proper motion of the barycenter, we were able to obtain the proper motion of each maser feature with respect to the barycenter of the outflow ejected from S255 NIRS 3, as listed in Table 2. Due to the short monitoring period of only three months, some of the features have large uncertainties. However, we could see systematic proper motions in the outflow, as shown in Fig. 5. The less ordered proper motions in the NE lobe could be due to interaction with the burst-regenerated radio jet (see Fig. 4a).

The absolute proper motion of the barycenter was calculated from the absolute proper motion of the maser spot in the ID19 feature, and that relative to the barycenter and is equal to − 0.03 ± 0.40 mas yr−1 in right ascension and − 0.29 ± 0.30 mas yr−1 in declination. This means that the absolute proper motion of the barycenter is negligible. The motion is consistent with those obtained by Rygl et al. (2010) ( − 0.14 ± 0.54, − 0.84 ± 1.76 mas yr−1) and Burns et al. (2016) ( − 0.13 ± 0.20 , − 0.06 ± 0.27 mas yr−1). Therefore, we do not correct for the absolute positional offset due to the systemic motion with respect to the absolute coordinate frame when we compare our present results to those of previous works (Goddi et al. 2007; Burns et al. 2016).

The proper motion fitting of all maser features are presented in Figs. A.1A.15. The positions are measured with respect to that of the ID13 feature (not with respect to the barycenter or the absolute frame). The error bars are estimated from the Gaussian fitting errors, and the systematic errors of 0.1 mas summed in quadrature. In the same figures, we also plot the flux densities of the peak channels for each feature. In Fig. A.16, we plot the flux densities of those features that were detected in fewer than three epochs, which are insufficient for a proper motion measurement.

thumbnail Fig. A.1

Positions in right ascension (left) and declination (middle), and peak flux density (right) of maser feature 1 as a function of time.

thumbnail Fig. A.2

Same as Fig. A.1, but for feature 2.

thumbnail Fig. A.3

Same as Fig. A.1, but for feature 3. Triangles in the right panel represent the upper limit of the flux densities.

thumbnail Fig. A.4

Same as Fig. A.1, but for feature 4.

thumbnail Fig. A.5

Same as Fig. A.1, but for feature 5. Triangles in the right panel represent the upper limit of the flux densities.

thumbnail Fig. A.6

Same as Fig. A.1, but for feature 6.

thumbnail Fig. A.7

Same as Fig. A.1, but for feature 9. Triangles in the right panel represent the upper limit of the flux densities.

thumbnail Fig. A.8

Same as Fig. A.1, but for feature 10. Triangles in the right panel represent the upper limit of the flux densities.

thumbnail Fig. A.9

Same as Fig. A.1, but for feature 13.

thumbnail Fig. A.10

Same as Fig. A.1, but for feature 15.

thumbnail Fig. A.11

Same as Fig. A.1, but for feature 17. Triangles in the right panel represent the upper limit of the flux densities.

thumbnail Fig. A.12

Same as Fig. A.1, but for feature 18. Triangles in the right panel represent the upper limit of the flux densities.

thumbnail Fig. A.13

Same as Fig. A.1, but for feature 19.

thumbnail Fig. A.14

Same as Fig. A.1, but for feature 20. Triangles in the right panel represent the upper limit of the flux densities.

thumbnail Fig. A.15

Same as Fig. A.1, but for feature 22.

thumbnail Fig. A.16

Flux density variation for the features detected only in 1 or 2 VLBI epochs. Triangles in each panel represent the upper limit of the flux densities.

References

  1. Acosta-Pulido, J. A., Kun, M., Ábrahám, P., et al. 2007, AJ, 133, 2020 [NASA ADS] [CrossRef] [Google Scholar]
  2. Breen, S. L., Sobolev, A. M., Kaczmarek, J. F., et al. 2019, ApJ, 876, L25 [CrossRef] [Google Scholar]
  3. Brogan, C. L., Hunter, T. R., Cyganowski, C. J., et al. 2018, ApJ, 866, 87 [NASA ADS] [CrossRef] [Google Scholar]
  4. Brogan, C. L., Hunter, T. R., Towner, A. P. M., et al. 2019, ApJ, 881, L39 [CrossRef] [Google Scholar]
  5. Burns, R. A., Handa, T., Nagayama, T., Sunada, K., & Omodaka, T. 2016, MNRAS, 460, 283 [NASA ADS] [CrossRef] [Google Scholar]
  6. Burns, R. A., Sugiyama, K., Hirota, T., et al. 2020, Nat. Astron., 10, 506 [Google Scholar]
  7. Carattio Garatti, A., Garcia Lopez, R., Scholz, A., et al. 2011, A&A, 526, L1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  8. Carattio Garatti, A., Stecklum, B., Garcia Lopez, R., et al. 2017, Nat. Phys., 13, 276 [NASA ADS] [CrossRef] [Google Scholar]
  9. Cesaroni, R., Moscadelli, L., Neri, et al. 2018, A&A, 612, A103 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  10. Chen, P., Pearson, J. C., Pickett, H. M., Matsuura, S., & Blake, G. A. 2000, ApJS, 128, 371 [NASA ADS] [CrossRef] [Google Scholar]
  11. Chen, X., Sobolev, A. M., Breen, S. L., et al. 2020, ApJ, 890, L22 [CrossRef] [Google Scholar]
  12. Contreras Peña, C., Lucas, P. W., Kurtev, R., et al. 2017, MNRAS, 465, 3039 [NASA ADS] [CrossRef] [Google Scholar]
  13. Felli, M., Brand, J., Cesaroni, R., et al. 2007, A&A, 476, 373 [Google Scholar]
  14. Fujisawa, K., Yonekura, Y., Sugiyama, K., et al. 2015, ATel 8286, 1 [NASA ADS] [Google Scholar]
  15. Goddi, C., Moscadelli, L., Sanna, A., Cesaroni, R., & Minier, V. 2007, A&A, 461, 1027 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  16. Gray, M. D., Baudry, A., Richards, A. M. S., et al. 2016, MNRAS, 456, 374 [NASA ADS] [CrossRef] [Google Scholar]
  17. Greisen, E. W. 2003, Astrophys. Space Sci. Lib., 285, 109 [NASA ADS] [Google Scholar]
  18. Hirota, T., Kim, M. K., Kurono, Y., & Honma, M. 2014a, ApJ, 782, L28 [NASA ADS] [CrossRef] [Google Scholar]
  19. Hirota, T., Tsuboi, M., Kurono, Y., et al. 2014b, PASJ, 66, 106 [NASA ADS] [CrossRef] [Google Scholar]
  20. Honma, M., Kijima, M., Suda, H., et al. 2008a, PASJ, 60, 935 [NASA ADS] [CrossRef] [Google Scholar]
  21. Honma, M., Tamura, Y., & Reid, M. J. 2008b, PASJ, 60, 951 [NASA ADS] [Google Scholar]
  22. Hunter, T. R., Brogan, C. L., MacLeod, G., et al. 2017, ApJ, 837, L29 [NASA ADS] [CrossRef] [Google Scholar]
  23. Hunter, T. R., Brogan, C. L., MacLeod, G. C. et al. 2018, ApJ, 854, 170 [NASA ADS] [CrossRef] [Google Scholar]
  24. Kramer, B. H., Menten, K. M., Kamiński, T., et al. 2017, IAU Symp., 316, 155 [NASA ADS] [Google Scholar]
  25. Krumholz, M. R., Klein, R. I., McKee, C. F., Offner, S. S. R., & Cunningham, A. J. 2009, Science 323, 754 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  26. Kuiper, R., Klahr, H., Beuther, H., & Henning, T. 2011, ApJ, 732, 20 [NASA ADS] [CrossRef] [Google Scholar]
  27. Liu, S. -Y., Su, Y. -N., Zinchenko, I., Wang, K. -S., & Wang, Y. 2018, ApJ, 863, L12 [CrossRef] [Google Scholar]
  28. MacLeod, G. C., Smits, D. P., Goedhart, S., et al. 2018, MNRAS, 478, 1077 [CrossRef] [Google Scholar]
  29. MacLeod, G. C., Sugiyama, K., Hunter, T. R., et al. 2019, MNRAS, 489, 3981 [CrossRef] [Google Scholar]
  30. McMullin, J. P., Waters, B., Schiebel, D., Young, W., & Golap, K. 2007, ASP Conf. Ser., 376, 127 [Google Scholar]
  31. Menten, K. M., Melnick, G. J., & Phillips, T. G. 1990, ApJ, 350, L41 [NASA ADS] [CrossRef] [Google Scholar]
  32. Moscadelli, L., Cesaroni, R., & Rioja, M. J. 2005, A&A, 438, 889 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  33. Moscadelli, L., Sanna, A., Goddi, C., et al. 2017, A&A, 600, L8 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  34. Neufeld, D. A., & Melnick, G. J. 1990, ApJ, 352, L9 [Google Scholar]
  35. Patel, N. A., Curiel, S., Zhang, Q., et al. 2007, ApJ, 658, L55 [Google Scholar]
  36. Reid, M. J., & Honma, M. 2014, ARA&A, 52, 339 [Google Scholar]
  37. Reynolds, S. P. 1986, ApJ, 304, 713 [Google Scholar]
  38. Richards, A. M. S., Impellizzeri, C. M. V., Humphreys, E. M., et al. 2014, A&A, 572, L9 [Google Scholar]
  39. Rygl, K. L. J., Brunthaler, A., Reid, M. J., et al. 2010, A&A, 511, A2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  40. Sugiyama, K., Saito, Y., Yonekura, Y., & Momose, M. 2019, ATel 12446, 1 [Google Scholar]
  41. Szymczak, M., Olech, M., Wolak, P., Bartkiewicz, A., & Gawroński, M. 2016, MNRAS, 459, L56 [Google Scholar]
  42. Szymczak, M., Olech, M., Wolak, P., Gérard, E., Bartkiewicz, A. 2018, A&A, 617, A80 [CrossRef] [EDP Sciences] [Google Scholar]
  43. Uchiyama, M., Yamashita, T., Sugiyama, K., et al. 2020, PASJ, 72, 4 [Google Scholar]
  44. VERA Collaboration (Hirota, T., et al.) 2020, PASJ, 72, 50 [Google Scholar]
  45. Wang, Y., Beuther, H., Bik, A., et al. 2011, A&A, 527, A32 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  46. Yonekura, Y., Saito, Y., Sugiyama, K., et al. 2016, PASJ, 68, 74 [Google Scholar]
  47. Zinchenko, I., Liu, S.-Y., Su, Y.-N., et al. 2015, ApJ, 810, 10 [Google Scholar]

1

Throughout the paper, we define a spot as an individual maser emission at a single spectral channel, while a feature refers to a group of spots that are detected in more than three consecutive spectral channels at positions coincident within the beam size. Thus, a feature is regarded as a physical gas clump emitting the 22 GHz maser line.

All Tables

Table 1

Summary of VERA observations.

Table 2

Positions and proper motions of all features.

Table 3

Summary of ALMA Band 7 continuum data.

Table 4

Summary of the 321 GHz and 22 GHz H2O maser features.

All Figures

thumbnail Fig. 1

Total-power spectra from the VERA 20-m single-dish monitoring observations. A vertical dashed line indicates the LSR velocityof 6.4 km s−1, corresponding to the brightest components during the VERA monitoring (Fig. 2).

In the text
thumbnail Fig. 2

Auto-correlation total-power (solid red lines, AC) and scalar-averaged cross-power (blue solid lines, XC) spectra of the H2O masers. Thespectra were computed by averaging all the antennas/baselines and time ranges. Baselines were subtracted by fitting the line-free channels with a polynomial function. A vertical dashed line indicates the LSR velocity of 6.4 km s−1, corresponding to the brightest components.

In the text
thumbnail Fig. 3

Velocity-integrated flux densities of the 22 GHz H2O masers observed with the VLA. A vertical dashed line indicates the LSR velocity of 6.4 km s−1, corresponding to the brightest components during the VERA monitoring (Fig. 2). The red and blue lines correspond to the spectrum integrated around the central source NIRS 3 and that of the newly found feature offset by 5′′ northwest (NW) of NIRS 3, respectively.

In the text
thumbnail Fig. 4

(a) Distribution of the 22 GHz H2O maser features (colored symbols; Table 2) superposed on the VLA K-band continuum observed on December 27, 2016 (purple contours; Cesaroni et al. 2018) and the moment 0 map of the 321 GHz H2O maser emission (gray scale). Contour levels are 4, 8, 16, ... times the rms noise level of the 0.13 mJy beam−1 and the 32 mJy beam−1 km s−1 for the VLA K-band continuum and the moment 0 map of the 321 GHz H2O, respectively. The (0, 0) position is the phase tracking center, RA(J2000) = 06h12m54.006s and Dec = +17°59′22.96′′. The synthesized beam size of the ALMA Band 7 observation is shown in the bottom left corner of each panel. (b) Distribution of peak centroids of the 321 GHz H2O maser features (colored symbols) superposed on the ALMA Band 7 continuum (magenta contours) and the moment 0 map of the 321 GHz H2O maser emission (gray scale). Contour levels for the ALMA Band 7 continuum are 32, 64, and 128 times the rms noise level of 0.33 mJy beam−1, showing only the brightest part around the continuum peak.

In the text
thumbnail Fig. 5

Distributions of the 22 GHz H2O masers indicated by the plus symbols with the proper motion vectors measured with respect to the barycenter. The color code indicates the LSR velocities of the maser features. The 22 GHz H2O maser features observed in previous observations by Goddi et al. (2007) in 2005 and Burns et al. (2016) in 2010 are also plottedwith open circles and filled squares, respectively. The solid line indicates the outflow axis derived from a linear fit to the positions of all maser features, including previous results.

In the text
thumbnail Fig. 6

Cross-power spectrum of the 321 GHz H2O maser integrated over the 0.5′′ × 0.5′′ region.

In the text
thumbnail Fig. 7

(a) Distributions of the 22 GHz H2O maser features obtained after rotating the plot in Fig. 5 by 43.5 degrees. The abscissa is parallel to the outflow as indicated by a solid line in Fig. 5. Horizontal bars at the bottom right corner indicate the minimum, maximum, and mean positions of the SW lobealong the outflow axis for each epoch (2005, 2010, and 2017 from bottom to top). (b) Same as (a), where the displacement along the Y axis has been replaced by the observing epoch. The black solid line represents the movement of the head of the bow shock along the outflow axis.

In the text
thumbnail Fig. 8

Light curves of maser lines and continuum emissions. A red vertical line indicates the start of the 6.7 GHz methanol maser flare (July 07, 2015) reported by Fujisawa et al. (2015). Dashed lines indicate the dates of January 01, from 2015 to 2019. (a) Flux densities of the H2O masers integrated between − 15 and + 25 km s−1, from the VLA and VERA observations. The VERA results include the single-dish monitoring and the integrated spectra extracted from the VLBI image cubes. The VLA data were extracted from the image cubes, excluding the NW component. (b) Flux densities of representative 6.7 GHz methanol maser features observed with the Hitachi 32-m radio telescope (Yonekura et al. 2016), showing different behaviors among the light curves. (c) IR light curve at Ks band observed with the 1.5-m Kanata telescope in Hiroshima, Japan (Uchiyama et al. 2020). (d) Variation of the radio continuum emission at 22 GHz observed with the VLA (Cesaroni et al. 2018) and ALMA Band 7 at ~330 GHz (Liu et al. 2018, present paper). For the 330 GHz continuum emission on December 02, 2017 measured in the present study, we used the value obtained by integrating over a 0.5′′ aperture (420 mJy) rather than the result of the Gaussian fit (see discussion in Sect. 3.3).

In the text
thumbnail Fig. 9

Flux densities of the H2O maser features in the whole region around NIRS 3, the NE lobe, central (continuum) position, and the SW lobe observed with VERA (the corresponding spectra are shown in Fig. 10). The flux densities were computed by integrating the emission over the velocity range between − 15 and + 25 km s−1. The error bars represent the 30% flux uncertainty. A dashed vertical line indicates January 01, 2017.

In the text
thumbnail Fig. 10

Spectra of the H2O maser features integrated over different regions in the image cubes from the VERA monitoring data: (a) the whole emission region around NIRS 3; (b) the NE lobe; (c) the central (continuum) position; and (d) the SW lobe.A vertical dashed line indicates the LSR velocity of 6.4 km s−1, corresponding to the brightest components. Panel b: one can see the profiles of the weakest components at ~15 km s−1, while the brightest features at 6.4 km s−1 are better seen in panel a.

In the text
thumbnail Fig. 11

Comparison between VERA single-dish and VLA spectra observed at almost the same time in October 2016 and November (VLA)/December(VERA) 2016. Black, red, and blue solid lines show, respectively, the total-power spectrum from the VERA single-dish monitoring, the integrated VLA spectrum of NIRS 3, and that of the newly found maser source NW of NIRS 3.

In the text
thumbnail Fig. 12

Difference between auto-correlation spectra and those from VLA/VERA images. Left and right panels: flux densities integrated over all velocity components and differential spectra. Left panel: error bars correspond to 30% uncertainty, and the dotted vertical line indicates January 01, 2017. Right panel: dotted vertical line marks the systemic velocity.

In the text
thumbnail Fig. A.1

Positions in right ascension (left) and declination (middle), and peak flux density (right) of maser feature 1 as a function of time.

In the text
thumbnail Fig. A.2

Same as Fig. A.1, but for feature 2.

In the text
thumbnail Fig. A.3

Same as Fig. A.1, but for feature 3. Triangles in the right panel represent the upper limit of the flux densities.

In the text
thumbnail Fig. A.4

Same as Fig. A.1, but for feature 4.

In the text
thumbnail Fig. A.5

Same as Fig. A.1, but for feature 5. Triangles in the right panel represent the upper limit of the flux densities.

In the text
thumbnail Fig. A.6

Same as Fig. A.1, but for feature 6.

In the text
thumbnail Fig. A.7

Same as Fig. A.1, but for feature 9. Triangles in the right panel represent the upper limit of the flux densities.

In the text
thumbnail Fig. A.8

Same as Fig. A.1, but for feature 10. Triangles in the right panel represent the upper limit of the flux densities.

In the text
thumbnail Fig. A.9

Same as Fig. A.1, but for feature 13.

In the text
thumbnail Fig. A.10

Same as Fig. A.1, but for feature 15.

In the text
thumbnail Fig. A.11

Same as Fig. A.1, but for feature 17. Triangles in the right panel represent the upper limit of the flux densities.

In the text
thumbnail Fig. A.12

Same as Fig. A.1, but for feature 18. Triangles in the right panel represent the upper limit of the flux densities.

In the text
thumbnail Fig. A.13

Same as Fig. A.1, but for feature 19.

In the text
thumbnail Fig. A.14

Same as Fig. A.1, but for feature 20. Triangles in the right panel represent the upper limit of the flux densities.

In the text
thumbnail Fig. A.15

Same as Fig. A.1, but for feature 22.

In the text
thumbnail Fig. A.16

Flux density variation for the features detected only in 1 or 2 VLBI epochs. Triangles in each panel represent the upper limit of the flux densities.

In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.