Proper Motions of Water Masers in W49 N Measured by KaVA

We report the proper motions of 22 GHz water masers toward W49 N that were observed by the KVN and VERA Array (KaVA) during 2017 February–May. We found 263, 268, and 310 features in three successive epochs; they were distributed in a region of size 4 × 4 arcsec2. The strongest flux density was in the third epoch, and its averaged value was 18,090 Jy at V LSR +0.47 km s−1. For 102 H2O maser features, proper motion was detected across all three epochs. The average proper motions in R.A. and decl. offset were −0.352 and +0.890 mas yr−1, respectively. The morphology of the distribution of the H2O maser features was found to be a bipolar outflow structure with an inclination angle of 37° ± 13° to the line of sight, and the features were expanding from a well-defined outflow center. A model of the source combining expansion and rotation yielded a distance to W49 of 11.12 ± 0.96 kpc that is consistent with the results from trigonometric parallax. A redshifted lobe was situated in the northeast direction and a blueshifted lobe in the southwest direction. We also discussed the location of the powerful flaring H2O maser feature at V LSR = + 6 km s−1 and its possible mechanisms on the basis of spatial structures for the maser feature in VLBI maps observed with the KaVA, timed just before and during the rebrightening phase.


Introduction
The star formation process is an important problem that astronomers have been addressing for decades. The process is especially problematic in the case of high-mass stars (M *  8 M e ; Crowther et al. 2010) owing to their short lifetime and consequent scarcity of nearby objects and their frequent shrouding by interstellar material-so they are not easy to observe directly at short wavelengths, i.e., visible light-and because core nuclear fusion begins before the mass accretion process is complete. In order to observe inside the clouds, astronomers use a special technique of radio astronomy, very long baseline interferometry (VLBI; Rogers 1976; Thompson et al. 2017), to achieve very high angular resolution. In this technique, data from two or more widely separated radio telescopes are recorded on memory units for subsequent correlation.
Radio astronomers use a unique tool, cosmic masers, to penetrate the optically obscured star-forming regions (SFRs). The word "maser" was originally an acronym (Microwave Amplification by Stimulated Emission of Radiation) for the microwave analog of the laser. Studying the physical properties of masers toward the SFRs helps us to understand the properties of their environment, such as morphology, molecular abundance, density, dust and gas temperatures, and velocity fields.
The 22,235 MHz H 2 O maser is a common type of cosmic maser and easily detected in two types of Galactic objects: the envelopes of red supergiant and asymptotic giant branch stars and the disks and outflows associated with young stellar objects (e.g., Gray 2012). One excellent example of an SFR of this type is W49 North (W49 N, or G043.16+0.01), which is a region inside the high-mass SFR W49 A. It is the most luminous and complex SFR with H 2 O masers in our galaxy (Cheung et al. 1969;Burke et al. 1970). An H 2 O maser parallax distance of -+ 11.11 0.69 0.79 kpc has been measured (Zhang et al. 2013). The continuum sources were first identified by Dreher et al. (1984), i.e., seven main distinct components in W49 N and two subcomponents, which were observed at 1.5 and 2 cm using the Karl G. Jansky Very Large Array (VLA). 9 Those components (hereafter called regions) were A, B, C, D, F, and G. Dreher et al. (1984) defined two edge-brightened shells on the west side of the G region: the G1 and G2 sources. Using the VLA, De Pree et al. (2000) observed 7 mm radio continuum emission toward W49 N at 43.3 GHz with 0 3 resolution and found that the G2 source can be divided into three subsources, i.e., G2a, G2b, and G2c (see Figure 2 of De Pree et al. 2000). Extensive high-frequency and high-resolution imaging with the VLA and the Berkeley-Illinois-Maryland Association array was more recently completed and newly classified by De Pree et al. (2004), including a more refined set of definitions that subdivide the earlier designations of regions,with a further subdivision indicated by a lower-case letter following each subsource, for example G1s, G2a, G2b, and G2c and similarly for other regions. The 22 GHz H 2 O maser emission in W49 N has been observed and discussed in many previous works (Gwinn et al. 1992;Mac Low & Elitzur 1992, and references therein). The masers are found closely associated with the G1/ G2 sources and exhibit their strongest flux from these regions. Almost all works have proposed that the maser emission arises from an outflow centered near the G region. However, there were numerous H 2 O maser flares toward this region. For example, Honma et al. (2004) used a 20 m single dish of the VLBI Exploration of Radio Astrometry (VERA 10 ) to monitor the maximum phase of the outburst of H 2 O masers during 2003 October and also VLBI observations with the VERA array near its maximum phase. They found that the strongest outburst was at the local standard of rest (LSR) velocity of −30.7 km s −1 and located on the arc-like structure that was found in a region between G2a and G1, indicating a possible connection of the outburst to the shock phenomenon. Besides lots of detections of flux variability on the H 2 O masers in W49N (e.g., Boboltz et al. 1998;Liljeström & Gwinn 2000;Zhou et al. 2002), other powerful flaring events were reported for three maser spectral components at LSR velocities of −81, −60, and +6 km s −1 in 2017-2020 (Volvach et al. 2019a(Volvach et al. , 2020(Volvach et al. , 2019b. However, accurate positions for these flare components have not been determined using VLBI. There are a few intensive studies on the proper motions of H 2 O maser features, starting with Gwinn et al. (1992). They used five epochs of five-station VLBI observations during 1980-1982. Those five telescopes were the 37 m Haystack antenna in Westford, Massachusetts; the 43 m telescope of the National Radio Astronomy Observatory (NRAO) at Green Bank, West Virginia; one of the 25 m antennas of the VLA near Socorro, New Mexico; the Caltech 40 m telescope at Owens Valley; and either the 100 m telescope of the Max-Planck-Institut für Radioastronomie at Effelsberg, (then) West Germany, or the 26 m telescope at Onsala, Sweden. They found 105 H 2 O maser features around the G region and also computed expansion and rotation models of the H 2 O maser bipolar outflow. Another previous work (Zhang et al. 2013) used very high precision data, which were taken by the Very Long Baseline Array (VLBA) 11 to study the proper motion of H 2 O masers toward W49 N over the period 2010 March to 2011 April. They found numerous maser spots that appeared at four or more epochs within 1 yr and then selected only 11 maser spots detected at all 12 epochs in the arc-like structure of the G2a source. The absolute proper motion of the reference maser spot was found to be μ x = −4.49 ± 0.13 and μ y = −6.42 ± 0.12 mas yr −1 . They also fitted the data to models of expansion with and without rotation, which were based on the previous study by Gwinn et al. (1992).
In this work, we will report the flux variability and propermotion studies of H 2 O masers toward the W49 N SFR by using data from the KaVA, 12 which observed this source at three epochs within the period 2017 February-May. We have computed expansion and rotation models of the H 2 O maser bipolar outflow with a method that is independent of that employed by Gwinn et al. (1992).

Observations and Constructing the Image Cubes
Observations of the H 2 O (6 1,6 -5 2,3 ) line at 22.235080 GHz were carried out with KaVA, which consists of four 20 m antennas of VERA in Japan and three 21 m antennas of the Korean VLBI Network (KVN). The maximum baseline length of KaVA is 2270 km, providing the highest angular resolution of 1.2 mas. The pointing and phase-tracking center position of the target source, W49 N, was taken to be We observed W49 N in three epochs, each of which consisted of a time from horizon to horizon of 4.5 hr on two consecutive days. In the day-of-the-year (doy) system, these epochs were 058-059, 102-103, and 130-131 in 2017. The net on-source time was 3.2 hr epoch -1 . Left-handed circular polarization data were recorded on the hard disk devices at 1 Gbps. The total bandwidth of each session was 8 × 32 MHz, and a spectral resolution for the 22 GHz maser line was set to be 31.25 kHz, which corresponded to a total velocity coverage and resolution of 430 and 0.42 km s −1 , respectively. We assigned only one of the 32 MHz spectral windows to the 22 GHz H 2 O masers. The systemic velocity of W49 N is equal to +10.0 km s −1 (Cesaroni et al. 1988). To cover a velocity range as wide as ±300 km s −1 , with respect to the systemic velocity, we shifted the central frequency by 26 MHz (350 km s −1 ) in the observations of two consecutive days. We observed a nearby calibrator, DA406, for bandpass and delay calibration once every hour.
Correlation processing was carried out at the Korea-Japan Correlation Center in KASI, Korea. Calibration and synthesis imaging were performed using the NRAO Astronomical Image Processing System (AIPS) software package (Greisen 2003) in the standard manner. We calibrated bandpass and delay by using a nearby (in angular separation from the source) calibrator, DA406. Phase calibration was done using the strongest H 2 O maser spot in the first epoch, and it was the same position for every epoch 13 in W49 N by fringe fitting and selfcalibration. We adopted the V LSR of the strongest channel in every epoch for doing self-calibration; these had values in the range [+0.81:+1.20] km s −1 .
In order to search for maser emission in W49 N that was extended both in space (i.e., 2″ × 2″ square) and in velocity (∼600 km s −1 ), we made VLBI images using the following procedure. First, we made coarse images of W49 N for the first sessions (doy of 058 and 059 in 2017). All of the spectral channels in the 32 MHz bandwidth were analyzed to search for maser emission with a field of view of 800 mas and a 0.2 mas pixel size (4096 × 4096 pixels). To cover the 2″ × 2″ region, the center position of each image cube was shifted by 800 mas. Because of the extremely strong maser emission of >10 4 Jy, the image dynamic range was limited by the strong side lobes. In order to exclude these side lobes, we set a detection threshold of five times the rms noise in each channel to identify the real maser spots. Next, we made full images for all of the sessions in all of the channels but only around the selected positions of identified maser features in the first session (doy of 058 and 059). We regard spots as real for the purposes of further data analysis if they were detected at least in two consecutive sessions and were spatially coincident with each other within the beam size and spectrally coincident within the velocity width. Because there were plenty of H 2 O maser spots and many image fields, we developed program-specific software in the form of simple Fortran and bash shell scripts (see Phetra et al. 2019). After that, we applied Python scripts to classify the maser features and record the proper motion obtained from those features. Feature classification and propermotion calculations are discussed in more detail in Sections 2.2.1 and 2.2.2.

How to Define Maser Features
We measured the peak flux of H 2 O maser spots by using the Gaussian fitting technique from the SAD task in AIPS. The results will be reported in terms of the position offset in R.A. and decl. and the LSR velocities in each channel. The position of a feature is defined as a weighted mean of the positions over the channels in which its constituent spots are detected, with the weighting being equal to the corresponding channel flux density of that spot,¯( where x i is the measured value in each channel, either the offset coordinate position or the v LSR ; S i is the channel flux density; and n is the number of data (spots) in the same feature. The standard deviation on the mean is used for the error value in each parameter, which is defined by

The Proper-motion Calculation
In order to define the movement of a feature over all three epochs, we used the following process.
1. The reference position of a stable feature should have the same position offset in every epoch. 2. All features will be shifted with respect to the stable reference feature. 3. As explained in Section 2, all weighted LSR velocities for each moving feature should have the same range, and the standard deviation of velocities should be within 2.0 km s −1 (i.e., 0.4 km s −1 × 5 consecutive channels).
Even though we had only three epochs of data, we compared two methods for calculating the proper-motion results, following linear and quadratic equations (see Appendix A for the explanation of the equations), which were obtained from the KaVA observations. However, the results calculated from the linear equation were realistic, while just some features were found to conform with the quadratic equation. Therefore, the linear relationship was applied and fitted with the angular displacement of features in both R.A. and decl. offset. In order to avoid systematic errors, all of the proper motions in each direction were averaged, and the result was subtracted from each feature again. After that, square-root values of the proper motion in both R.A. and decl. offset (μ abs ) can be summarized from the equation ( ) m m m = + mas x y abs 2 2 . The transverse velocities were the vector quantity that was calculated from the basic equation, v t (km s −1 ) = 4.7μd, where μ is the proper motion (in units of arcseconds per year), and d is the distance to W49 N (we adopt the value of ∼11.1 kpc; see Zhang et al. 2013). Radial velocities were calculated from the "radio definition" nonrelativistic Doppler shift equation, where v obs is the observation velocity in kilometers per second; f obs and f rest are the observation and rest frequency, respectively; and c is the velocity of light in vacuum, 299,792 km s −1 . The transverse and radial velocities will be used in the modeling of the H 2 O maser bipolar outflow.

Comparison of the Maser Positions of KaVA with the High-accuracy VLBA Data
Since this observation lacked external phase-referencing astrometry, we could not directly obtain the absolute positions of masers. In order to compare the position offsets of maser features from the KaVA with the high-accuracy data of the VLBA (particularly the previous work by Zhang et al. 2013), we have applied some basic techniques to overlay and compare both data sets. The main processes can be explained as follows.
1. From a total of 102 moving features from the KaVA data (see Tables 3 and 4 in Appendix C), we select from the same range of data sets as those from Zhang et al. (2013, private communication) a small subset of at least three or four features in various well-separated positions, e.g., in the north, east, west, and south directions. 2. We make sure that they are the same features in both data sets; that is, the V LSR should be within ±0.4 km s −1 or equal to one channel of the spectral resolution of the KaVA. 3. We calculate the difference of the position offsets from both data sets. Note that the values should not be dramatically different from each other in R.A. and decl. offset. After that, we average them together to get the finalized value of the offset.
The position offsets will be applied and adopted as KaVA absolute positions (i.e., ∼10 mas accuracy) before overlaying with another facility, e.g., VLBA.

Total-power Spectrum
The 22 GHz H 2 O masers in W49 N were observed with the KaVA at three epochs, i.e., 2017 February 27-28, 2017 March 12-13, and 2017 May 10-11 (hereafter called the first, second, and third epoch, respectively). The LSR velocities were between −300 and +300 km s −1 . The total-power spectrum in each epoch was measured using the POSSM task in AIPS and plotted as shown in Figure 1. Blue and red lines represent two data sets that covered the LSR velocities as mentioned in Section 2. The dotted vertical line refers to the systemic velocity of this source, which is +10 km s −1 . We found that the strongest flux was at +0.47 km s −1 , and the highest flux density was 18,090 Jy in 2017 May. Other high-flux spectral regions were within the highly blueshifted LSR velocity range [−50:−100] km s −1 in every epoch.

Number of Maser Features and Their Distribution
We determined the number of maser features using the criteria that are described in Section 2.2.1, and the number of features found in epochs 1-3 were 263, 268, and 310, respectively. The number of features tended to increase, and we note especially that the third epoch has the highest number. These results are in good agreement with the previous work of Volvach et al. (2019a), in which they found significant flares of the 22 GHz H 2 O maser emission during the period 2017 July to 2018 November. We plotted only the number of features that were found in all three epochs, as shown in panel (a) of Figure 2. Their distribution offsets are in a region of size 4 × 4  arcsec 2 . One maser feature is separated westward by (−1.8, 0.0) arcsec. The arc-like structure and the clumpy features that showed the strongest flux are marked inside the dashed box of panel (a) and magnified as shown in panel (b). Mostly, the redshifted and blueshifted LSR velocities were found in the eastern and western parts of the source, respectively. Figure 3 shows the proper motions of the subset of 22 GHz H 2 O maser features that were detected by the KaVA at all three epochs. The length of each vector represents the magnitude of the proper motion compared to the standard (black arrow). The direction of each vector represents the direction of proper motion in R.A. and decl. offset coordinates. The colors of the arrows show the average radial velocity of each moving feature according to the color bar to the right of the figure. Only 102 features were classified as moving features over all three epochs (more details can be found in the Appendix C), and their weighted positions appear in Figure 2. The overall feature distribution is that of a red-blue divided shape, similar to a bipolar outflow, which is strong support to previous works, e.g., Gwinn et al. (1992) and Zhang et al. (2013). We mostly found that the features are moving outward (expansion) from their weighted mean positions. Predominantly blueshifted and redshifted velocities are found in the western and eastern parts of the source, respectively.

Discussion
In this discussion, we will focus on the following topics: (1) the location of the flux variable H 2 O maser feature and (2) H 2 O outflow modeling.

Location of the Flux Variable H 2 O Maser Feature
The H 2 O maser features at V LSR = −81,−60, and +6 km s −1 have been reported as exhibiting powerful flaring phenomena since 2017 (Volvach et al. 2019a(Volvach et al. , 2020(Volvach et al. , 2019b. These flares began in ∼2017 September, ∼2017 May, and ∼2017 April, respectively. We conducted the KaVA monitoring observations in this paper on 2017 February 27-28, April 12-13, and May 10-11. The KaVA monitoring data thus potentially contribute to establishing the location of the feature responsible for the powerful H 2 O maser flare at V LSR = +6 km s −1 that began in ∼2017 April. Indeed, the integrated fluxes of the KaVA total-power spectra, recorded in the LSR velocity range from −300 to +300 km s −1 , are among the weakest found during the period from 2016 August to 2017 September, as shown in  Figure 5 shows a close-up of spatial distributions for the brightest channel at V LSR = +5.90 km s −1 in H 2 O maser feature 27 located at (x-offset, y-offset) ∼ (+507, −248) mas (see Tables 3 and 4) in epochs 1, 2, and 3. A peak maser spot in each epoch presented an intensity of 1920, 1973, and 2047 Jy beam −1 , respectively. Comparison of these intensities, taking into account an uncertainty in each absolute intensity value of 20%, suggests no clear flux variation in the spot based on VLBI maps in the KaVA observation period. We thus assume that this maser spot at V LSR = +5.90 km s −1 could be identified as the one corresponding to the powerful flare feature at V LSR = +6 km s −1 in Volvach et al. (2019b) on the basis of their coincidence of LSR velocities within the velocity resolution in the KaVA observations of 0.42 km s −1 . We cannot rule out the possibility of a maser-maser interaction   Kramer et al. 2022 in preparation, and private communication) and the KaVA as total-power spectra in this paper is presented by white and black circles, respectively. These integrations were calculated over the LSR velocity range from −300 to +300 km s −1 . Error bars for integrated fluxes correspond to a typical value for scaling absolute flux of 10% and 20% in the Effelsberg and KaVA data, respectively. scheme causing the powerful flares at V LSR = +6 km s −1 because of the pointlike structure of the spots in all three epochs. Sizes are 1.17 × 1.05, 1.29 × 1.19, and 1.29 × 1.14 mas 2 , comparable to the synthesized beam sizes of 1.14 × 1.01, 1.30 × 1.14, and 1.23 × 1.07 mas 2 , and prevent us from understanding whether multiple maser spots overlap along the line of sight. Furthermore, epochs 2 and 3 of the KaVA observations coincide with the earliest period of the initiation for the powerful flare. Nevertheless, we think it reasonable to claim that the powerful flare at V LSR = +6 km s −1 could be excited by the shock front in the arc-like structured clump where the outburst at V LSR = −30.7 km s −1 was detected in Honma et al. (2004), because this maser spot in the powerful flaring feature and the outbursting maser feature in Honma et al. (2004) coexist in a region of size ∼72 mas in the arc-like structure, corresponding to ∼800 au at the source distance.

Model Fits to the Data
We fitted a small number of models of the outflow structure to the data. These were, in increasing order of sophistication, a model in which the unknown z-position is computed from its 3D velocity and 2D sky position, models where the expansion velocity is described by single and double power laws, and a model where rotation of the bipolar structure is allowed in addition to a two-power expansion law. The theory for these models largely follows Equations (2)-(13) of Gwinn et al. (1992), but our analysis differs in a number of ways that are explained below. The first difference is in the method of solution that was used to determine the optimized values of the free parameters for each model. While Gwinn et al. (1992) used a Levenberg-Marquardt method on N p partial derivatives of the χ 2 statistic, where N p is the number of free parameters, we used a combination of the simulated annealing and downhill simplex methods, implemented through the routines AMOTSA and AMOEBA (Press et al. 1997). This method has the advantage that only the χ 2 function itself needs to be computed, rather than the N p derivatives. Another difference is that the z-offset velocity is treated as a free parameter in the most basic model by Gwinn et al. (1992), but we note that χ 2 is formally independent of this parameter in this model (see Appendix B). Notes on the individual models are listed below, and the results are summarized in Table 1. 1. Basic model. The χ 2 statistic is given by Equation (2) of Gwinn et al. (1992), which we show here as where N p is the number of free parameters to be fitted for by comparing observational velocities u i and model velocities w i for i = 1...N m maser features, each supplying three components of the observed velocity vector, u i . The model velocities take the form where v 0 is a global velocity offset, and v i is the model-dependent flow velocity computed for the ith maser feature. The sole assumption about the flow velocity in this model, as in Equation (7) of Gwinn et al. (1992), is that it is radial in direction with respect to an offset position. Each feature can have an independent magnitude of the velocity. A consequence is that the source distance is indeterminate. A distance of 12.0 kpc was assumed to convert angular to linear distances in the sky plane. The set of parameters, in principle, includes the three components of v 0 ; the two sky components, x 0 and y 0 , of the offset position; and the z-position of every maser feature. However, by taking a partial derivative of Equation (3) with respect to z i , the z-position of feature i, it is possible to show that the unknown z-positions may be determined analytically from the five remaining independent parameters from the equation This is a corrected version of Equation (8) of Gwinn et al., noting that the original is dimensionally incorrect. We go on to show in Appendix B that the final model velocity w i of each feature is also independent of v 0z , so this model ends up with four independent parameters: the x-and y-components of the position and velocity offsets from the reference position and velocity, respectively. 2. Single power law. In this and all later models, χ 2 follows from Equation (10)   measurement error σ i . Comparing Equation (5) with Equation (3), there is now a distinction between the distance-independent velocity measurements in the zcomponent and the distance-dependent motions in the sky plane, making the source distance a free parameter. This is written in the dimensionless form d/d 0 , where the scaling distance d 0 = 12.0 kpc. Flow velocities now have an assumed functional form for each model (N) with N > 1. In this case, model 2, the flow velocity is given bŷ and r i is the magnitude of this vector. The expressionr i is the unit form of r i . Equation (6) adds four fitting parameters: the sky position offsets x 0 and y 0 , the bulk flow constant V 0 , and the expansion power law α 0 . The dimensions of V 0 are fixed by the requirement that the expression a V r i 0 0 must be a velocity. The single power-law case was not explicitly investigated by Gwinn et al. (1992) but follows from their Equation (11) by the simple expedient of setting V 1 = 0. A solution for this model is computed in Zhang et al. (2013). There are eight free parameters after elimination of the z-positions: the two sky position offsets, three velocity offsets, the distance, V 0 , and α 0 .
As in model 1, we may take a partial derivative of the formula for χ 2 , in this case, Equation (5), with respect to a z-position, z j . The resulting expression still allows the elimination of the z-positions in terms of the eight free parameters listed above. However, unlike Equation (4), the equation is nonlinear and must be solved numerically. A similar scheme applies to models 3 and 4 below. 3. Double power law. This follows model 2 of Gwinn et al.
The flow velocities are now computed from which is designed to better fit an expansion with a fast and a slow component. Equation (7) adds two additional fitting parameters compared with Equation (6), so there are now 10 free parameters after elimination of the zpositions, and this model offers the best possibility for a direct comparison of the solution methods by substituting the KaVA data with the data from Table 4 of Gwinn et al. (1992).  Table 5). The flow velocity is now represented aŝ where the group w g r i has the dimensions of an angular velocity and represents rotation of the outflow axis. The rotational term introduces a radial dependence, controlled by the index γ and a direction, represented by ω, so there Note. a In (km s −1 ) 2 , χ 2 is dimensionless in all other models. Model 1 has an assumed distance of 12.0 kpc. are four additional free parameters compared with model 3, bringing the total to 14 after elimination of the zpositions.
The uncertainty estimates on the model parameters in Table 1 were determined by running the model from different starting positions; the uncertainty is the sample standard deviation of 10 trials for each model. The results for our models are not completely in accord with those of either Gwinn et al. (1992) or Zhang et al. (2013). The most important difference is in the velocity offset. However, the most important parameter, the source distance, is reasonably robust in the more sophisticated models that include an adequate number of fitting parameters. The value 11.12 ± 0.96 kpc from our model 4 is consistent with the value of -+ 11.11 0.69 0.79 kpc derived from parallax by Zhang et al. (2013). The distance estimate from our model 2, though less accurate, is also consistent with the parallax distance (see comparison of the distance measurements toward W49 N from previous works and the current study in Table 2). Our model with rotation is consistent with a radial power of γ = 0.0, which makes our model closest to model 3 of Gwinn et al. in this respect. Our fits to the components of the rotation vector are consistent with those of model 3 of Gwinn et al., so its orientation with respect to the axis of the maser outflow looks the same within the errors as in Figure 8 of Gwinn et al. (1992).

Tests of the Algorithm
The optimization algorithm containing the AMOTSA and AMOEBA codes from Press et al. (1997) was tested on a generic function not related to the present work. The function is the Ddimensional multispike function, which is used because it contains a large number of local minima controlled by the parameter n and a global minimum at x i = 0.5 for all 1 i D. It is therefore a challenging problem to identify the global minimum among all of its close rivals. For the case where n = 21 and D = 4, the same as in model 1, the algorithm was able to detect the global minimum to an accuracy of 0.01 after 10 attempts, starting from different random points in the 4D unit volume. The algorithm was allowed a total of 10 8 moves, with 10 5 at any one of 1000 "temperatures." The stiffness factor that controls the number of high temperatures relative to the number of low temperatures was set to 1.2.

Conclusion
We have studied the proper motions of 22 GHz H 2 O masers toward the W49 N region using the KaVA over three epochs during 2017 February to May. The main results of this study are as follows.
1. The H 2 O maser features were distributed in a region of 4 × 4 arcsec 2 or 0.30 × 0.30 pc 2 , and the highest flux density was 18,090 Jy. They also show a bipolar outflow structure with wide-angle wind morphology in the northeast direction and an inclination angle of 37°± 13°to the line of sight.
2. The averaged internal proper motions in the R.A. and decl. offsets were −0.352 and +0.890 mas yr -1 , respectively. 3. On the basis of VLBI maps with the KaVA catching the timing at just before and during the rebrightening phase, the powerful flaring H 2 O maser feature at +6 km s −1 is possibly located at (x-offset, y-offset), i.e., ∼ (+507, −248) mas (seeTables3 and 4in Appendix C) from the LSR velocity correspondence basis and could be excited by the shock front in the arc-like structured clump, where the outburst was detected by Honma et al. (2004), i.e., VLSR = −30.7 km s −1 , although the possibility of a maser-maser interaction scheme could not be ruled out due to the pointlike spatial structures of the maser spots in all three VLBI epochs. 4. Modeling of the H 2 O masers as an expanding region with rotation yields a distance to W49 N of 11.12 ± 0.96 kpc. This value is consistent with the value recovered from trigonometric parallax observations of masers in Zhang et al. (2013) and earlier modeling results from Gwinn et al. (1992).

Appendix A Analysis of the Proper Motion
To analyze the proper motion, we need to consider a suitable equation for calculating the displacement angle in both R.A. and decl. offset. Even though the time range of observations was rather short, a total of 72 days, we have still considered two types of evolution equations: quadratic and linear.

A.1. Quadratic Equation
The basic evolution equation of the quadratic system is We used Cramer's rule to calculate the unknown parameters. The solutions are   Then, the unknown parameters are Finally, the proper motion at epoch i is the slope from Equation (A1) as follows: and the units of m i are milliarcseconds per day.

A.2. Linear Equation
The main evolution equation of the linear system is

Appendix B Proof that v 0z is Redundant in Model 1
We begin with Equation (4) of the main text, a formula for the unknown z-position of the ith maser feature, where each feature has a measured velocity vector (u ix , u iy , u iz ) and a measured sky position (x i , y i ), and the model has a velocity offset (u 0x , u 0y , u 0z ) and sky position offset (x 0 , y 0 ). Equation (4) is essentially Equation (8) of Gwinn et al. (1992) but with the position offsets included and the dimensions corrected. We simplify Equation (4) by defining the 2D radius ρ = ((x i − x 0 ), (y i − y 0 )), and the 2D dot product, noting that q i is the sky plane velocity relative to its component offsets. The simplified form is therefore We ultimately seek to show that components of the model velocity Gwinn et al. 1992) are all independent of v 0z , since Equation (2) of Gwinn et al. contains only the w i and measured velocities. It is useful to take Equation (5) of Gwinn et al. (1992) with p j = V 0 ( j), the flow velocity of a particular maser feature. As these velocities are individual to each feature, the sum over the features collapses, and the problem reduces to evaluating partial differentials of the form ∂v i /∂V 0 (i). From the definition of the flow velocity (Equation (7) of Gwinn et al. 1992), these differentials reduce to r i /r i , and Equation (5) of that work becomes the dot product, This equation does not appear explicitly in Gwinn et al. (1992) but is alluded to in the text following their Equation (9). A more useful version of Equation (B3) is obtained by expanding w i in terms of v 0i and v i and further expanding v i as a flow velocity from Equation (5) of Gwinn et al. (1992). With the flow speed as the subject, It is useful to divide Equation (B4) by r i before expanding the dot product and eliminating all occurrences of z i with the aid of Equation (B2), yielding 2 for the squared magnitude of the 3D position. After identification of the 2D dot product in the numerator and r i 2 in the denominator, a cancellation is possible, and Equation (B5) reduces to which is obviously independent of v 0z because the vectors q i and ρ i are both 2D sky plane quantities.
To complete the proof, we note that the x-component of the flow velocity is v ix = V 0 (i)(x i − x 0 )/r i from Equations (3), (4), and (7) of Gwinn et al. (1992). The right-hand side of this expression is the right-hand side of Equation (B6) multiplied by a relative x-position, neither of which contain v 0z . To form the model x-velocity, w ix , we add only v 0x , so this model velocity component is also independent of v 0z . A very similar argument holds for the model y-velocity. The z-component is slightly different; we begin in the same way with but it is then best to eliminate the dot product and the factor of r 1 i 2 with the aid of Equation (B2), leaving, simply, v iz = u iz − v 0z . In this case, the flow velocity does depend on v 0z , but upon adding v 0z to generate w iz , this dependence is lost.
We have now demonstrated that all three components of w i are independent of v 0z and so, therefore, is the χ 2 , or S 2 , statistic, given by Equation (2) of Gwinn et al. (1992).

Appendix C Results from the Proper-motion Analysis
These tables show the results obtained from the propermotion analysis according to the quadratic and linear fitting (Tables 3 and 4). Fea.
x-offset a δx b y-offset c δy μ No.
x-offset a δx b y-offset c δy μ  To compare the maser position offsets of KaVA with the high-accuracy VLBA data obtained from Zhang et al. (2013), however, each offset will be shifted with the value of (−516.464, +111.798) mas, i.e., R.A. and decl. offset, respectively.