© The Authors 2023

1. Introduction

All the stars in the Universe are generally divided into two main groups according to their mass: low-mass stars (M < 8 M_⊙) and high-mass stars (M > 8 M_⊙). For the first group, the formation process has been quite well established thus far, yet for the high-mass stars there are still several open questions that need to be addressed (e.g., Tan et al. 2014; Motte et al. 2018). Based on observational constraints, an evolutionary scenario of high-mass star formation (HMSF) was proposed by Motte et al. (2018), summarized in the following. High-mass stars form in molecular complexes, in particular in parsec-scale massive clumps or clouds that first undergo a global controlled collapse, ultimately forming low-mass prestellar cores. This first phase is known as starless massive dense cores (MDCs). After about 10⁴ yr, low-mass prestellar cores become protostars with growing mass. This phase is called protostellar MDC phase. It is only after ∼10⁵ yr more that we have the high-mass protostellar phase whereby, thanks to the gas flow streams generated by the global collapse, the protostars become high-mass protostars, even though they still harbor low-mass stellar embryos. At this stage of evolution, the high-mass protostars are still quiet in the infrared (IR) and if their accretion rates are efficient and strong, then they drive outflows. In this phase, the accretion disks are already formed. As soon as the stellar embryos reach more than 8 M_⊙ the high-mass protostars become IR-bright and they develop ultra-compact H II region (UCH II) that are quenched by infalling gas or confined to the photoevaporating disks. Finally, we have the H II region phase that lasts about 10⁵ − 10⁶ yr. In this phase, the ultraviolet radiation from the stellar embryos produces the H II region and the gas accretion toward the newborn star first slows down and then stops by terminating the main accretion phase. Thus, the high-mass star is formed.

The great importance that the magnetic field holds for several phases of HMSF has been demonstrated thanks to magnetohydrodynamical (MHD) simulations (e.g., Myers et al. 2014; Kuiper et al. 2016; Matsushita et al. 2018; Machida & Hosokawa 2020; Rosen & Krumholz 2020; Oliva & Kuiper 2023). Despite the observational difficulties, among which the low number of high-mass protostars that are usually found to be densely clustered in molecular clouds and their long distances to the Sun are of the greatest importance, measurements of the morphology and strength of magnetic fields close to high-mass protostars are still possible. These parameters can be obtained by observing the polarized emission of dust and molecular lines with the Atacama Large Millimeter Array (ALMA; e.g., Dall’Olio et al. 2019; Sanhueza et al. 2021) and by observing the polarized maser emission with the very long baseline interferometry (VLBI) technique (e.g., Surcis et al. 2011a,b, 2014; Surcis et al. 2022). A relevant star-forming region where it is possible to measure the magnetic field close to high-mass young stellar objects (YSOs) in different evolutionary phases is W75N(B).

The active high-mass star-forming region (HMSFR) W75N is part of the Cygnus X complex (Westerhout 1958; Harvey et al. 1977; Habing & Israel 1979) at a distance of 1.30 ± 0.07 kpc (Rygl et al. 2012). Haschick et al. (1981) identified three compact radio regions within W75N at a spatial resolution of ∼1.​​″5: W75N(A), W75N(B), and W75N(C). Hunter et al. (1994) further mapped the continuum emission of W75N(B) at a resolution of ∼0.​​″5 revealing the presence of three very compact subregions (Ba, Bb, and Bc). For the first time Hunter et al. (1994) also underlined the great importance of the region for understanding the star formation process thanks to its significant activity. The subregions Ba and Bb were renamed as VLA 1 and VLA 3 in 1997 when Torrelles et al. (1997) imaged them at a resolution of ∼0″​​.1 with the Very Large Array (VLA). Torrelles et al. (1997) also imaged, for the first time, another weaker and more compact radio source between VLA 1 and VLA 3, known as VLA 2. While VLA 1 and VLA 3 show a fully elongated radio continuum emission along the northeast-southwesterly and northwest-southeasterly directions, respectively, that are consistent with the morphology of thermal radio jets, VLA 2 shows a quasi-circular morphology resembling an UCH II. Despite their small separation (VLA 2 is only ∼0″​​.8, ∼1000 au at 1.3 kpc, from VLA 1), the three radio sources are thought to be massive YSOs (Shepherd et al. 2003) at three different evolutionary stages, with VLA 1 as the most evolved and VLA 2 as the least evolved (Torrelles et al. 1997). In addition, Shepherd et al. (2003) reported a systemic velocity for W75N(B) of V_lsr = +10.0 km s⁻¹, which can also be considered as the systemic velocity for VLA 1, VLA 2, and VLA 3. Five more new radio sources were identified in W75N(B) thanks to highly sensitive observations made with the upgraded VLA, namely: VLA 4, which is about ∼1″​​.5 south of Bc (Carrasco-González et al. 2010), VLA [NE] (∼8″ northeast of VLA 2), VLA [SW] (∼6″ southwest of VLA 2), and Bd (∼0.5″ northeast of VLA 4; Rodríguez-Kamenetzky et al. 2020). Thanks to ALMA observations at 1.3 mm (spatial resolution of ∼1.2″), Rodríguez-Kamenetzky et al. (2020) were able to associate VLA [NE] and VLA [SW] with the millimeter cores MM3 and MM2, respectively, indicating that they are also embedded YSOs. Furthermore, VLA 1, VLA 2, and VLA 3 are also associated with a millimeter core (MM1) but VLA 4, Bc, and Bd are not, which suggests that they are not embedded YSOs but, rather, shock-ionized gas (Rodríguez-Kamenetzky et al. 2020). In addition, multiwavelength VLA observations carried out by Carrasco-González et al. (2015) showed that VLA 2 had changed its morphology between 1996 and 2014 from a compact roundish source (≤160 au; Torrelles et al. 1997) to an extended source that is elongated in the northeast-southwest direction (220 × ≤160 au, with a position angle of PA = 65°), while VLA 1 still has continued to show an unchanged morphology since 1996 (Rodríguez-Kamenetzky et al. 2020). The spectral index analysis indicates that VLA 2 is a thermal, collimated ionized wind surrounded by a dusty disk or envelope (Carrasco-González et al. 2015), VLA 1 is at the early stage of the photoionization and it is driving a thermal radio jet, and VLA 3 is also driving a thermal radio jet whose shocks are traced by the obscured Herbig-Haro (HH) objects Bc and VLA 4 (Rodríguez-Kamenetzky et al. 2020).

A large-scale high-velocity CO-outflow, with an extension greater than 3 pc (PA = 66°) and with a total molecular mass greater than 255 M_⊙, has been detected from W75N(B) (e.g., Hunter et al. 1994; Davis et al. 1998; Shepherd et al. 2003; Makin & Froebrich 2018). Shepherd et al. (2003) found that the entire CO emission of W75N uncovers a complex morphology of multiple, overlapping outflows. In particular, they have suggested that VLA 2 may be driving the large-scale CO-outflow, while VLA 1 and VLA 3 are, instead, the centers of two additional, more compact outflows (extension of ∼0.21 pc and ∼0.15 pc, respectively). However, this scenario is not yet entirely clear and the main powering source of the large-scale CO-outflow remains unknown (e.g., Qiu et al. 2008). In addition, Shepherd et al. (2003) also determined that more than 10% of the molecular gas in W75N is outflowing material, and the combined outflow energy is roughly half the gravitational binding energy of the cloud, thus preventing its further collapse.

The great interest that W75N(B) has aroused in the past is also due to the presence of several maser species (OH, CH₃OH, and H₂O) around the two sources VLA 1 and VLA 2 (e.g., Haschick et al. 1981; Hunter et al. 1994; Torrelles et al. 1997; Lekht & Krasnov 2000; Minier et al. 2000; Surcis et al. 2009; Fish et al. 2011; Kang et al. 2016; Colom et al. 2018, 2021). In particular, OH masers are distributed throughout the region, with the majority of the maser spots associated with VLA 1 (Hutawarakorn et al. 2002; Fish et al. 2005). Nevertheless, VLA 2 is the site of the most intensive OH flare ever registered in a star-forming region (1000 Jy, Alakoz et al. 2005; Slysh et al. 2010), while other OH maser emission sites are situated on a ring structure around VLA 1, VLA 2, and VLA 3 (Hutawarakorn et al. 2002). Fish et al. (2011) measured the proper motions of the OH masers showing that most of those near VLA 1 (located ∼1″ northwest) are moving northward (V < +5 km s⁻¹) and those associated with VLA 2, and located southwest, are moving both toward the southwest and southeast (V < +5 km s⁻¹). The detection of OH maser Zeeman-pairs provided measurements of magnetic field strength between 6 and 8 mG close to VLA 1, where the 22 GHz H₂O maser are detected, and up to about 17 mG around VLA 2 (Fish et al. 2011). However, higher values (40–70 mG; Slysh & Migenes 2006; Slysh et al. 2010) were measured during the strong OH maser flare. The 6.7 GHz CH₃OH masers are only associated with VLA 1 and they are distributed parallel to the thermal radio jet of VLA 1 (Minier et al. 2000, 2001; Surcis et al. 2009). No CH₃OH masers had been detected around VLA 2 (Surcis et al. 2009; Rygl et al. 2012) until 2014, when three maser spots were detected in the southwest of the source (Carrasco-González et al. 2015). Thanks to the polarized emission of the CH₃OH masers, Surcis et al. (2009) were able to measure a magnetic field oriented southwest-northeast that is perfectly aligned with VLA 1, and whose strength on the line of the sight was found to be |B_|||> 10 mG (Surcis et al. 2019).

The 22 GHz H₂O masers have been widely studied and monitored both with single dishes and interferometers, revealing high-intensity variations, with extreme maser flares (up to 10³ Jy) and important variations in the maser distribution (e.g., Lekht & Sorochenko 1984; Lekht 1994; Lekht & Krasnov 2000; Torrelles et al. 2003; Surcis et al. 2011a, 2014; Kim et al. 2013; Krasnov et al. 2015; Kim & Kim 2018). The H₂O masers are associated with VLA 1 and VLA 2, and only one maser was associated with VLA 3 in 1996 (Torrelles et al. 1997), but it was never detected again (e.g., Torrelles et al. 2003; Surcis et al. 2014). Over a period of 16 yr, VLBI observations have shown that the evolution of the H₂O masers around VLA 1 and VLA 2, despite their close separation, is completely different. While the H₂O masers around VLA 1 are always linearly distributed (PA ≈ 43°) along the thermal radio jet, those detected around VLA 2 are instead tracing an expanding shell (expanding velocity of ∼30 km s⁻¹, Surcis et al. 2014) that evolved from a quasi-circular (Torrelles et al. 2003; Surcis et al. 2011a) to an elliptical structure (Kim et al. 2013; Surcis et al. 2014), following the morphology change in the continuum emission observed by Carrasco-González et al. (2015). Therefore, in VLA 2, the H₂O masers might be tracing the evolution from a non-collimated to a collimated outflow (Surcis et al. 2014). Furthermore, Surcis et al. (2014) also showed that the magnetic field around VLA 1 has not changed from 2005 to 2012 and it is always oriented along the direction of the thermal radio jet. On the other hand, the orientation of the magnetic field around VLA 2 has changed in a way that is consistent with the new direction of the major-axis of the shell-like structure that is now aligned with the thermal radio jet of VLA 1.

The peculiarity of the H₂O maser shell expansion with the contemporary variation of the magnetic field around VLA 2, together with the presence of a nearby immutable VLA 1 source, has made W75N(B) one of the most interesting cases with regard to investigating the evolution of early massive YSOs. For this reason, we performed VLBI monitoring observations of 22 GHz H₂O maser emission in full polarization mode every two years from 2014 to 2020 for a total of four epochs. In Sect. 3, we report the results of the monitoring observations from Sect. 2. We discuss the magnetic fields around VLA 1 and VLA 2 in Sect. 4. Finally, we present a full picture of the two massive YSOs in Sect. 5. Our conclusions are given in Sect. 6.

2. Observations and analysis

W75N(B) was observed at 22 GHz in full polarization spectral mode with several European VLBI Network¹ (EVN) antennas on four epochs separated by two years (see Table 1). The observations were carried out in June 2014 (epoch 2014.46), 2016 (epoch 2016.45), and 2018 (epoch 2018.44) and in October 2020 (epoch 2020.82), for a total observation time per epoch of 12 h. The bandwidth was 4 MHz in epochs 2014.46 and 2016.45, providing a local standard of rest velocity (V_lsr) range of ∼55 km s⁻¹ (after calibration ranging from −15.7 km s⁻¹ to +31.7 km s⁻¹), and 8 MHz in epochs 2018.44 and 2020.82 (after calibration ranging from −30.3 km s⁻¹ to +48.1 km s⁻¹). We carried out observations using a bandwidth that is two times wider in the last two epochs to search for maser emission at velocities <  − 15.7 km s⁻¹, as indicated by the results obtained in epoch 2016.45 (see Sect. A.2 and Table A.4). To measure the absolute positions of the H₂O masers, observations were conducted in phase-reference mode (with cycles phase-calibrator – target of 45 s–45 s). The phase-reference calibrator was J2040+4527 (separation = 2° ​​.856). The data were correlated with the EVN software correlator (SFXC; Keimpema et al. 2015) at the Joint Institute for VLBI ERIC (JIVE), using 2048 channels in epochs 2014.46 and 2016.45, and 4096 channels in epochs 2018.44 and 2020.82, generating all four polarization combinations (RR, LL, RL, and LR), with a spectral resolution in all epochs of ∼2 kHz (∼0.03 km s⁻¹).

The data were calibrated using the Astronomical Image Processing Software package (AIPS) by following the standard calibration procedure (e.g., Surcis et al. 2011a). Specifically, the bandpass, the delay, the phase, and the polarization calibration were performed in all epochs on the calibrator J2202+4216. Then we performed the fringe-fitting and the self-calibration on the brightest maser feature of each epoch. We note that the reference maser features VLA1.1.15, VLA1.2.07, VLA1.3.04, and VLA1.4.05 are given in Tables A.1, A.3, A.5, and A.7, respectively (for the notation definition see Sect. 3). The I, Q, U, and V Stokes cubes were then imaged using the AIPS task IMAGR. Afterward, the Q and U cubes were combined to produce cubes of linearly polarized intensity ($POLC=\sqrt{Q^2+U^2-{\sigma }_P^2}$) and polarization angle (POLA = 1/2 ×  atan(U/Q)). The polarized intensity cubes were corrected according to the noise ${\sigma }_P=\sqrt{[{(Q\times {\sigma }_Q)}^2+{(U\times {\sigma }_U)}^2]/(Q^2+U^2)}$, where σ_Q and σ_U are the noise of the Q and U Stokes cubes, respectively. The formal error on POLA due to the thermal noise is given by σ_POLA = 0.5 (σ_P/POLC)×(180° /π) (Wardle & Kronberg 1974).

To measure the absolute positions of the H₂O maser features, in all four epochs, we self-calibrated the phase-reference source J2040+4527 and the amplitude and phase solutions were applied only to the uncalibrated peak channel of the brightest maser feature of each epoch. We show the contours maps of J2040+4527 in Fig. 1 and the absolute positions (along with their uncertainties) of the reference maser feature of each epoch in Table 2. The uncertainties of the absolute positions of the reference maser features are estimated by adding quadratically the systematic errors (ε_sys = 0.028 mas) due to the source elevation limit and the separation between the calibrator and the target (Reid & Honma 2014, maximum and minimum elevation of W75N(B) at each station and in all epochs were ∼80° and ∼35°, respectively), the errors due to the thermal noise (ε_rms), the errors of the Gaussian fit of the peak spot of the reference maser feature (ε_gauss), and the errors of α₂₀₀₀ and δ₂₀₀₀ of J2040+4527 (0.68 mas and 0.80 mas, respectively, Petrov et al. 2011). The latter is usually unnecessary in relative astrometry studies (as in this one), however, due to the broad time interval between the EVN epochs we include it in our analysis to account for possible variation of the calibrator position from one epoch to the others. Another source of uncertainty that is taken into consideration here is due to the identification of the maser features that is meant to account for the maser spots scatter of each maser features (see Surcis et al. 2011a, hereafter S11a) and it is equal to half of the restoring beam of W75N(B).

Fig. 1. Contour maps of the phase-reference calibrator J2040+4527 in the EVN epochs 2014.46, 2016.45, 2018.44, and 2020.82 from left to right, respectively. Contours are 4, 8, 16, and 32 times the 0.6 mJy beam−1. The restoring beam is overplotted on the bottom left corner of each panel (see Table 1).

We had to use different approaches than we did in the past to calibrate the linear polarization angles of the H₂O maser features. This is because the last National Radio Astronomy Observatory (NRAO) POLCAL observations² of J2202+4216 were made in May-June 2012. For epoch 2014.46, we assumed that the magnetic field orientation on the plane of the sky around VLA 1 has not changed from the last VLBI epoch (i.e., 2012.54, Surcis et al. 2014. hereafter S14), as was the case between epochs 2005.89 and 2012.54 (S14). We aligned the H₂O maser features detected toward W75N(B) in epoch 2012.54, for which the absolute positions were unknown, with those in epoch 2014.46 by associating the maser feature VLA1.1.21 (V_lsr = 11.37 km s⁻¹; Table A.1) with the maser feature VLA1.28 (V_lsr = 11.35 km s⁻¹; S14), because both are spatially coincident and have similar radial velocities. As a consequence we found that the maser feature VLA1.1.05 (Table A.1) can be considered to be part of the same maser clump gas of VLA1.07 (S14), although they are not exactly the same maser feature. Therefore, we can assume that VLA1.1.05 has a mean linear polarization angle (χ) equal to −25°, which is the linear polarization angle measured in the maser clump gas of VLA1.07 by S14. We were thus able to estimate the polarization angles of the H₂O maser features in epoch 2014.46, with a systemic error of no more than ∼8°. For epoch 2016.46, we calibrated the linear polarization angles of the H₂O maser features by rotating the linear polarization angle measured for J2202+4216 from our EVN data to the one measured on 1st June 2016 by one of the Korean VLBI Network (KVN) antennas (χ_{J2202 + 4216, KVN} = −61° ±12° calibrated by using 3C286, priv. comm.). In this case, the systemic error was of no more than ∼12°. In this epoch, we also tried to calibrate the linear polarization angles by observing the well known polarization calibrator 3C48, but the rms was not sufficient to detect the linear polarization intensity above 3σ (see Table 1). Thanks to the wider bandwidth in the last two epochs (2018.44 and 2020.82), we were instead able to calibrate the linear polarization angles by using the polarization calibrator 3C48. For 3C48, we assumed a polarization angle at K-band equal to −70° (Perley & Butler 2013).

Similarly to S14, we analyzed the polarimetric data following the procedure reported in S11a. Therefore, we first identified the H₂O maser features and then we measured the mean linear polarization fraction (P_l) as well as the mean linear polarization angle (χ) for each identified H₂O maser feature, considering only the consecutive channels (more than two) across the total intensity spectrum for which the polarized intensity is ≥10 rms. Afterward, by using the full radiative transfer method (FRTM) code for 22 GHz H₂O masers (Vlemmings et al. 2006), which is based on the model for unsaturated 22 GHz H₂O masers of Nedoluha & Watson (1992), we modeled the observed total intensity and linear polarization spectra of the linearly polarized maser features by gridding the intrinsic maser linewidth (ΔV_i) between 0.4 and 4.5 km s⁻¹, in steps of 0.025 km s⁻¹. We did this using a least square fitting routine (χ²-model) with an upper limit of the emerging brightness temperature (T_bΔΩ, where ΔΩ is the maser beaming) of 10¹¹ K sr (for more details see Appendix A of Vlemmings et al. 2006). Thus, we were able to obtain as outputs of the FRTM code the values of T_bΔΩ and ΔV_i that produce the best fit models for our linearly polarized maser features. Because the FRTM code is based on a model for unsaturated H₂O maser, it cannot properly disentangle the values of T_bΔΩ and ΔV_i in the case of saturated maser features and therefore it provides only a lower limit for T_bΔΩ and an upper limit for ΔV_i. An upper limit for T_bΔΩ below which the maser features can be considered unsaturated is T_bΔΩ < 6.7 × 10⁹ K sr (Surcis et al. 2011b). Then, from T_bΔΩ and P_l we could estimate the angle between the maser propagation direction and the magnetic field (θ) from which the 90° ambiguity of the magnetic field orientation with respect to the linear polarization vectors can be solved.

Indeed, if θ > θ_crit = 55° the magnetic field appears to be perpendicular to the linear polarization vectors; otherwise, it is parallel (Goldreich et al. 1973). Nedoluha & Watson (1992) found that T_bΔΩ scaled linearly with Γ + Γ_ν, which are the maser decay rate and cross-relaxation rate, respectively. As explained in Surcis et al. (2011a), Γ_ν varies with the temperature and Γ = 1 s⁻¹ for the H₂O maser emission, therefore, in our fit, we consider a value of Γ + Γ_ν = 1 s⁻¹ that allows us to adjust the fitted T_bΔΩ values by simply scaling it according to the real Γ + Γ_ν value, as described in Anderson & Watson (1993). We note that ΔV_i and θ do not need to be adjusted. The errors of T_bΔΩ, ΔV_i, and θ were determined by analyzing the probability distribution function of the full radiative transfer χ²-model fits. In cases where a maser feature is also circularly polarized, we can use the best estimates of T_bΔΩ and ΔV_i in the FRTM code to produce I and V models that can be used to fit the V spectra of the circularly polarized maser features from which we can measure the Zeeman splitting. Due to the typical weak circularly polarized emission of H₂O masers (< 1%), it is important to consider the self-noise (σ_s.−n.) produced by the maser features in their channels to determine whether the circularly polarized emission is real. The self-noise becomes important when the power contributed by the astronomical maser is a significant portion of the total received power (Sault 2012). Therefore, a detection of circularly polarized emission can be considered real only when the V peak intensity of a maser feature is both > 3rms and > 3σ_s.−n. (see Table 1).

3. Results

We report in Table 3 the number of H₂O maser features detected around VLA 1 and VLA 2 in the four EVN epochs with their corresponding local standard of rest velocity (V_lsr) and peak intensity (I) ranges, the mean linewidth of the maser features (⟨Δv_L⟩), the ranges of P_l and of the circular polarization fraction (P_V), the ranges of χ, and all the outputs of the FRTM code (T_bΔΩ, ΔV_i, and θ) with the derived magnetic field parameters. These are the ranges of the estimated orientation of the magnetic field on the plane of the sky (Φ_B) and its error-weighted value (⟨Φ_B⟩), the range of the magnetic field strength along the line of sight in absolute values (|B_|||) and their error-weighted values (⟨|B_|||⟩), and the error-weighted values of the estimated 3D magnetic field strength (⟨|B|⟩). In addition, we report the detailed results obtained from the four EVN epochs with their plots and tables in Appendix A. We should mention here that the detected H₂O maser features are called throughout the paper as VLA1.x.yy and VLA2.x.yy, where x is a number from 1 to 4 indicating the EVN epoch from the first (2014.46) to the fourth (2020.82) and yy is the number of the H₂O maser feature counted from west to east in each epoch. Maser features with the same yy value but with different x value are not necessarily related to each others, that is they are not necessarily the same maser feature detected in different epochs.

The main objective of our monitoring project is to determine whether the 22 GHz H₂O maser distributions (Sect. 3.1) and the maser features characteristics such intensity (Sect. 3.2) and polarization (Sect. 3.3) around VLA 1 and VLA 2 show any variation over time and how the magnetic field behaves accordingly. To do this, we have to correct the positions of the detected H₂O maser features from epoch 2016.45 to epoch 2020.82 by considering the proper motion of the entire region with respect to the Earth and assuming epoch 2014.46 as the reference epoch. Rygl et al. (2012) measured the median proper motion of the 6.7 GHz CH₃OH maser features associated with VLA 1 and VLA 2. They assumed that this motion represents the proper motion of the entire region W75N(B). The components of this proper motion along the right ascension and declination are ⟨μ_α⟩ = (−1.97 ± 0.10) mas yr⁻¹ and ⟨μ_δ⟩ = (−4.16 ± 0.15) mas yr⁻¹. We therefore corrected the positions of the H₂O maser features of the last three EVN epochs by assuming that both VLA 1 and VLA 2 moved from epoch 2014.45 with a proper motion equal to that measured by Rygl et al. (2012). Furthermore, a comparison of the absolute positions of the continuum emission of VLA 1 and VLA 2 at K-band, as measured by Torrelles et al. (1997, epoch 1996.96), Carrasco-González et al. (2015) and Rodríguez-Kamenetzky et al. (2020, epoch 2014.20) with the VLA, has shown that, while VLA 2 does not show any further motion within the region W75N(B), VLA 1 does actually move. The proper motion of VLA 1 within W75N(B) is ⟨μ_α⟩^VLA 1 = ( − 6.3 ± 0.4) mas yr⁻¹ and ⟨μ_δ⟩^VLA 1 = ( + 5.0 ± 0.4) mas yr⁻¹ (see Appendix B). Therefore, before comparing the maser features around VLA 1, we must first apply a further correction to their positions.

Unfortunately, we cannot compare the maser distributions of the EVN epochs with those observed previously with the VLBA in epochs 2005.89 and 2012.54, because we do not have any information on the absolute positions of the maser features in those epochs (S11a; S14). Nevertheless, when necessary refer to Table A.3 of S14 for the parameters of the VLBA epochs 2005.89 and 2012.54.

3.1. Spatial and velocity distribution of the H₂O masers

3.1.1. VLA 1

The number of 22 GHz H₂O maser features detected around VLA 1 has decreased from the EVN epoch 2014.46 to the EVN epoch 2020.82 (see Table 3). If we plot all the H₂O maser features detected in all the four EVN epochs, after correcting their positions as reported above, we see that these are always distributed along the radio continuum emission of VLA 1 presented by Rodríguez-Kamenetzky et al. (2020). This continuum emission was obtained at a resolution of tens of milliarcseconds by observing with the VLA a wide range of frequencies (4 − 48 GHz) in 2014. Rodríguez-Kamenetzky et al. (2020) concluded that VLA 1 is at the early stage of photoionization and it is driving a thermal radio jet at scale of about 0.1 arcsec (≈130 au). In Fig. 2 we overplot the H₂O maser features detected in the four EVN epochs to the continuum emission measured at Q-band by Rodríguez-Kamenetzky et al. (2020). The velocities of all the H₂O maser features are not shown in this figure, however, they are shown in the left panels of Fig. A.1, where every EVN epoch is reported in a different panel. In order to better visualize the accordance of the maser features distribution with the position angle of the thermal radio jet, we made a linear fit of the maser features in each epoch. The results of these fits are shown in the right panel of Fig. 2 and their parameters are reported in Table 4. For the linear fit of epoch 2014.46, we did not consider the group of five maser features located south because their positions would strongly affect the position angle of the fitting line. This is not the case for the northeast and southeast maser features of epochs 2016.45 and 2018.44, respectively, which we did indeed include in our linear fits. We see that the position angle of the fitted lines (see Fig. 2 and Table 4 for comparison) are consistent with the position angle of the thermal radio jet (solid gray line in Fig. 2, PA = +42° ± 5°; Rodríguez-Kamenetzky et al. 2020).

Fig. 2. Comparison of the H2O maser features detected toward VLA 1 (left panel) in the four EVN epochs (2014.46, 2016.45, 2018.44, and 2020.82) and superimposed to the uniform-weighted continuum map at Q-band (central frequency 44 GHz) of the thermal radio jet driven by VLA 1 obtained with the VLA by Rodríguez-Kamenetzky et al. (2020). The light gray contours are 5, 6, 7, 8, 10, 12, 14 times σ = 100 μJy beam−1 and the VLA beam is shown on the bottom right corner. The positions of the maser features are corrected assuming the proper motion of the region W75N(B) equal to the median proper motion measured for the 6.7 GHz CH3OH maser features by Rygl et al. (2012), ⟨μα⟩ = (−1.97 ± 0.10) mas yr−1 and ⟨μδ⟩ = (−4.16 ± 0.15) mas yr−1. In addition, the positions of the maser features associated with VLA 1 and detected in the last three epochs have been further corrected considering the proper motion of VLA 1 within W75N(B). This proper motion is ⟨μα⟩VLA 1 = ( − 6.3 ± 0.4) mas yr−1 and ⟨μδ⟩VLA 1 = ( + 5.0 ± 0.4) mas yr−1 (see Appendix B). The size of the symbols are ten times the uncertainties of the absolute positions of the maser features (see Table 2). For clarity we show the linear fit of the maser features for each epoch (right panel), along with the length of the segments correspond to the position of the maser features used in the fitting. The reference position is α2000 = 20h38m36s​.43399 and δ2000 = +42° 37′34″​​.8710. The parameters of the linear fit are reported in Table 4. The solid light-gray line indicates the orientation of the thermal radio jet (PA = +42°±5°; Rodríguez-Kamenetzky et al. 2020).

We note that the mean velocities along the line of sight of the maser features (⟨V⟩, as in Col. 7 of Table 4) in the first three epochs, which are all around +11.5 km s⁻¹, and their similar maximum line of sight velocities (V_max ≈ +16 km s⁻¹, Col. 8 of Table 4) are largely different than those of epoch 2020.82 (⟨V⟩= + 18.5 km s⁻¹ and V_max = +26.3 km s⁻¹). These differences might be explained either with an acceleration of the motion of VLA 1 along the line of sight and away from us, or with the variation of the masing conditions – even though a combination of the two would appear to be the realistic case (see Fig. 3). In addition, we also note that the maser velocities are spatially mixed on the plane of the sky in each epoch.

Fig. 3. Comparison of the total sum of the H2O maser spectra detected toward VLA 1 with the EVN in epochs 2014.46, 2016.45, 2018.44, and 2020.82. The vertical dashed dark-green line indicates the assumed systemic velocity of the region (Vlsr = +10.0 km s−1; Shepherd et al. 2003).

3.1.2. VLA 2

The number of 22 GHz H₂O maser features detected around VLA 2 in the four EVN epochs ranges between 37 and 44 (see Table 3), which is roughly half of the number of maser features detected in the two previous VLBA epochs (88 and 68 in 2005.89 and 2012.54, respectively; S14). These differences might not be because of the sensitivity of the EVN observations, which was better than that of the VLBA epochs (see Table 1 and S11a; S14), but of the maser activity of the region. The maser distribution in the four EVN epochs is consistent with the elliptical shell observed for the first time by Kim et al. (2013) and confirmed by S14. S14 also measured a mean expansion velocity of the maser shell on the plane of the sky, which begun to expand in 1999 when the shell was quasi-circular and continued by becoming elliptical in 2007 (Torrelles et al. 2003; S11a; Kim et al. 2013), of 4.9 mas yr⁻¹ (30 km s⁻¹ at a distance of 1.3 kpc). However, the previous observations could not be properly compared because only those presented in Kim et al. (2013) were made in phase-reference mode and therefore the absolute positions of most of the maser features in the other epochs were unknown. For estimating the expansion velocity, S14 assumed that the center of the shells coincides in the different epochs and consequently the measured expansion is radial. Since we were able to measure the absolute positions of the maser features in the four EVN epochs, we are now able to properly measure the expansion velocity. Nevertheless, we can compare our expansion velocities only in magnitude and not in direction with that measured by S14. We plot all the H₂O maser features detected with the EVN and overplotted to the continuum emission at K-band (Carrasco-González et al. 2015) in the left panel of Fig. 4. Here, we group the maser features in four different zones. As for VLA 1, the velocities of all the H₂O maser features detected in VLA 2 are shown in the right panels of Fig. A.1.

Fig. 4. Comparison of the H2O maser features detected toward VLA 2 (left panel) in the four EVN epochs (2014.46, 2016.45, 2018.44, and 2020.82) and superimposed to the natural-weighted continuum map at K-band (central frequency 22 GHz) of the thermal, collimated ionized wind emitted by VLA 2 obtained with the VLA by Carrasco-González et al. (2015). The light gray contours are 5, 10, 15, 20, 25, 50, and 75 times σ = 10 μJy beam−1 and the VLA beam is shown on the bottom right corner. The positions of the maser features are corrected assuming the proper motion of the subregion W75N(B) equal to the median proper motion measured for the 6.7 GHz CH3OH maser features by Rygl et al. (2012), ⟨μα⟩ = (−1.97 ± 0.10) mas yr−1 and ⟨μδ⟩ = (−4.16 ± 0.15) mas yr−1. The size of the symbols are ten times the uncertainties of the absolute positions of the maser features (see Table 2). For clarity, we show the parabolic fit of the maser features of zone 1 and the linear fit of the maser features of zone 2 (north features) and zone 4 for each epoch (right panel). The reference position is α2000 = 20h38m36s​.43399 and δ2000 = +42° 37′34.​​″8710. The parameters of the linear and parabolic fits are reported in Tables 4 and 5, respectively.

Zone 1. This is located northeast. Interestingly, while S14 found expanding motions, our data now shows the opposite trend, with apparent motions toward the central source. Indeed, the maser features seem to not trace an expansion, but they actually seem to “bounce”. The maser features of epoch 2014.46 are generally slightly at northeast of those of epochs 2016.45 and 2018.44, which are more northeast than those of epoch 2020.82. This apparent “bouncing” could be interpreted as evidence that the outflowing gas, where the masers arise, encounters an obstacle, as already proposed by Kim & Kim (2018). This obstacle can be either a much denser medium that stops the expansion of the gas (case A) or the absence of physical conditions, such as density and temperature, for producing the H₂O maser emission (case B). In case A, the impact of the gas with a denser medium might have produced an additional slow inward shock that pumped the maser features in the four EVN epochs. To estimate the proper motion of the gas on the plane of the sky due to this inward shock we fit the maser features of zone 1 with polynomials of second order, one per EVN epoch, and the results are shown in the right panel of Fig. 4 and the parameters are listed in Table 5. In addition in Fig. 5 we show the areas where the maser features considered in the polynomial fit are located. We note that only the areas of epochs 2014.46 and 2020.82 do not overlap and, therefore, their polynomial fits can be used to estimate the proper motion.

Fig. 5. Visual comparison of the polynomial fit areas of zone 1 of VLA2. Each area represents the area in each epoch where the maser features considered in the polynomial fit (see Fig. 4 and Table 5) are located. The colors correspond to those given in Fig. 4.

Hence, we estimated the proper motion by measuring a mean distance between the curves of epochs 2014.46 and 2020.82 and the resulting velocity is equal to −2.5 mas yr⁻¹ that corresponds to ∼ − 15.4 km s⁻¹ at a distance of 1.3 kpc, the minus sign indicates that the motion is opposite to the expansion measured by S14. Although case A is plausible, we should note that the curvature of the maser distribution in zone 1 is contrary to what would be expected (see right panel of Fig. 4). In case B, the physical conditions of the gas in the northeast are not suitable for the maser emission and what we observe are maser features pumped by different outward shocks in each epoch. However, both cases are undermined by the presence of the broad maser feature VLA2.4.39 (I = 1.94 Jy beam⁻¹, Δv_L = 1.65 km s⁻¹; see Table A.8) that is located farther northeast of all the other maser features of zone 1. This maser feature shows physical parameters, such as the velocity (V_lsr = +16.35 km s⁻¹), consistent with those of the other maser features of the zone. However, the presence of this isolated maser feature can also be justified by the presence of a belt of gas, between VLA2.4.39 and the other maser features, where the maser conditions are not met. Therefore, if the shock that pumped VLA2.4.39 is the same that pumped the maser features in the EVN epoch 2014.46, we can estimate the proper motion due to this shock. VLA2.4.39 is at about 53 mas from the front of the maser features detected in epoch 2014.46, therefore the proper motion is 7.9 mas yr⁻¹ that corresponds to ∼49 km s⁻¹ at a distance of 1.3 kpc. This is higher than the expansion velocity measured by S14 and the proper motions measured by Kim et al. (2013), ∼20 − 30 km s⁻¹.

Carrasco-González et al. (2015) presented the most recent continuum maps of VLA 2, which were obtained by observing four different frequency bands (C, U, K, and Q) with the VLA in 2014. From these new maps, it was possible to verify the variation of collimation of the outflow emitted from VLA 2 as traced by the H₂O maser features from 1999 to 2012 (S14). We can therefore compare our maser distributions with the VLA continuum emission at K-band (see Fig. 4). The asymmetric morphology of the K-band continuum emission, which shows a weaker emission toward southwest and none toward northeast, suggests the presence of an obstacle toward the northeast, as we describe earlier in this paper. This obstacle might be an inhomogeneity within the dusty disk or envelope supposed by Carrasco-González et al. (2015). In particular, the distribution of the maser features of zone 1 seems to be the continuation toward northwest of the last external contour of the continuum emission at 50 μJy beam⁻¹ (see Fig. 4), even though no continuum emission is detected where these maser features arise.

We note that the velocity range of the maser features of zone 1 in the four EVN epochs (−11 $<V_{\mathrm{lsr}}^{2014.46}<+17$ km s⁻¹, −16 km s⁻¹ $<V_{\mathrm{lsr}}^{2016.45}<+17$ km s⁻¹, −16 km s⁻¹ $<V_{\mathrm{lsr}}^{2018.44}<+16$ km s⁻¹, and +6 km s⁻¹ $<V_{\mathrm{lsr}}^{2020.82}<+17$ km s⁻¹) is similar to that covered in the VLBA epoch 2012.54 (−10 km s⁻¹ $<V_{\mathrm{lsr}}^{2012.54}<+18$ km s⁻¹; S14), even though the range observed in the EVN epoch 2020.82 shows only redshifted velocities. This might suggest that the maser features of zone 1 traced the same shock, that moved outward, from 2012.54 to 2018.44, and then they quenched any time between 2018.44 to 2020.82. We also note that, in the first three EVN epochs, the blue- and red-shifted maser features do not follow any particular spatial distribution, that is blue- and redshifted maser features are spatially coincident and do not show any velocity gradient. This might further indicate that the gas expands along the walls of the denser medium when encountering it. Consequently, the maser features in the last EVN epoch (2020.82) might trace a different shock that still move outward rather than inward, according to the morphology of the maser distribution. This might exclude the possible inward shock supposed above in the discussion of case A.

Zone 2. The low number of H₂O maser features of zone 2 (northwest, see Fig. 4) detected from the VLBA epoch 2012.54 to the EVN epoch 2018.44 and their total H₂O maser intensity (see Sect. 3.2.2) indicate a low maser activity in this zone since the appearance of the elliptical maser distribution. However, the maser features detected in the last EVN epoch 2020.82 are about three times greater in number and their total H₂O maser intensity is almost ten times higher than previously detected (see Sect. 3.2.2), suggesting a sudden increment of the maser activity as never observed before in the northeast part of the maser distribution, neither when the distribution was quasi-circular (Torrelles et al. 1997, 2003; S11a) nor afterward (Kim et al. 2013; S11a; S14). Comparing only the EVN epochs, we note that all the maser features in zone 2 north of declination 42°37′34.​​″2 in Fig. 4 have a velocity range between −1 km s⁻¹ and +9 km s⁻¹, that is, they are all blueshifted, with the exception of VLA2.4.12 that shows a redshifted velocity of V_lsr = +12.77 km s⁻¹. The rest of the maser features of zone 2 (south of declination 42°37′34.​​″2 in Fig. 4), which are the only ones detected on the continuum emission at K-band (see Fig. 4), show much higher redshifted velocities (+15 km s⁻¹< V_lsr < +23 km s⁻¹). In contrast to zone 1, we see that the maser features in zone 2 detected in one epoch are always detected outward with respect the previous epoch. To better estimate the apparent motion between different epochs, we made a linear fit of the maser features in each epoch. In particular, we were able to measure the expansion velocity of the north maser features in zone 2 by considering the distance of the median point of the line of one epoch from the line of the next epoch. The results are reported in Table 4 and in the right panel of Fig. 4. The expansion velocities on the plane of the sky measured between the epochs 2016.45 and 2018.44 (3.8 mas yr⁻¹) and between epochs 2018.44 and 2020.82 (4.3 mas yr⁻¹) are consistent with the expansion velocity of 30 km s⁻¹ (4.9 mas yr⁻¹) measured by S14. In zone 2, it is also possible to identify four maser features, each detected in a different EVN epoch, located below the linear fits and at the center of the dashed rectangle that highlights zone 2 in Fig. 4, apparently seeming to trace an outward motion. These maser features are VLA2.1.04, VLA2.2.07, VLA2.3.06, and VLA2.4.13 (see Table 6). Although these maser features have different line of sight velocities (see Col. 3 of Table 6), we can estimate the expansion velocity between the four EVN epochs by assuming that the shock that pumped them is the same. We find a proper motion of 4.8 mas yr⁻¹ (∼30 km s⁻¹ at 1.3 kpc) between the epochs 2014.46 and 2016.45, 4.5 mas yr⁻¹ (∼28 km s⁻¹) between epochs 2016.45 and 2018.44, and 4.7 mas yr⁻¹ (∼29 km s⁻¹) between the epochs 2018.44 and 2020.82 (see Table 6). The mean proper motion between epochs 2014.46 and 2020.82 is equal to 4.6 mas yr⁻¹ (∼28 km s⁻¹), which is again consistent with the expansion velocity measured by S14 on the plane of the sky, that is 30 km s⁻¹. The consistency of the expansion velocities measured from the maser features of zone 2 with that measured previously by S14 is a further clue of the presence, since epoch 2014.46, of an obstacle in front of the expanding gas in zone 1 that is absent in front of the gas in zone 2 that can freely expand.

Zone 3. The H₂O maser features of zone 3 are all located toward the bright core of the continuum emission observed by Carrasco-González et al. (2015; see Fig. 4) and show velocities in the range between +1 km s⁻¹ and +22 km s⁻¹. We note that the blueshifted maser features are always located outward than the redshifted in all the epochs. Looking at Fig. 4 we can divide the maser features in zone 3 into four groups: one in the east (group 3A), one in the center (group 3B), one in the south (group 3C), and the fourth in the southwest (group 3D). Group 3B coincides with the peak of the continuum emission at K-band (∼800 μJy beam⁻¹), groups 3A and 3C are located between the contours at 25σ (250 μJy beam⁻¹) and 50σ (500 μJy beam⁻¹) of Fig. 4, while group 3D between the contours at 20σ (200 μJy beam⁻¹) and 25σ (250 μJy beam⁻¹). For the first three groups (3A-3C), we can identify maser features with similar velocity in at least two epochs and we are therefore able to estimate their proper motions. In particular, in group 3A, we identified a maser feature in three consecutive epochs (from 2014.46 to 2018.44), corresponding to VLA2.1.22, VLA2.2.18, and VLA2.3.32 (see Table 6). We measured a constant velocity of ∼12 km s⁻¹ (2.0 mas yr⁻¹), pointing slightly toward southwest between epochs 2014.46 and 2016.45 and between epochs 2016.45 and 2018.44 (see Table 6). The maser features VLA2.1.03 and VLA2.3.09 of group 3C can be considered to be tracing the same gas, therefore, the estimated proper motion is ∼8 km s⁻¹ (1.3 mas yr⁻¹; see Table 6). Both proper motions of groups 3A and 3C are much lower than that measured by S14. The identification of common maser features with similar V_lsr in group 3B is very difficult. However, we can identify three maser features by considering their relative position in Fig. 4. These are VLA2.1.09, VLA2.2.12, and VLA2.3.13 (see Table 6). The proper motion is then ∼31 km s⁻¹ (5.0 mas yr⁻¹) between epochs 2014.46 and 2016.45, and ∼46 km s⁻¹ (7.5 mas yr⁻¹) between epochs 2016.45 and 2018.44 (see Table 6). The mean proper motion between epochs 2014.46 and 2018.44 is about 38 km s⁻¹ (6.2 mas yr⁻¹) that is slightly higher than the expansion velocity of 30 km s⁻¹ measured by S14.

Zone 4. This zone is located southwest and the H₂O maser features distribution is roughly aligned north-south as observed previously (Torrelles et al. 2003; Kim et al. 2013; S14). No H₂O maser emission was detected toward this zone in the EVN epoch 2014.46. A comparison with the continuum emission at K-band reported by Carrasco-González et al. (2015) reveals that the maser features of zone 4, which are all redshifted (+16 km s⁻¹< V_lsr < +29 km s⁻¹) and are all associated with the weak continuum emission at around 200-300 μJy beam⁻¹ (see Fig. 4). The results of the linear fit of the maser features are reported in Table 4 and displayed in the right panel of Fig. 4. Because the linear fit of epoch 2016.45 crosses that of epoch 2018.44 (see right panel of Fig. 4), we can estimate an expansion velocity on the plane of the sky only between epochs 2018.44 and 2020.82. This is the largest ever value measured toward the maser features around VLA 2, equal to 12.7 mas yr⁻¹ (∼78 km s⁻¹). This high velocity might suggest that the gas does not encounter any dense matter toward southwest that could slow it down.

3.2. Intensity variability

3.2.1. VLA 1

The total H₂O maser intensity (∑I) of the maser features first had an increment between the VLBA epoch 2012.54 (S14) and the EVN epoch 2016.45 and then it decreased again till the last EVN epoch 2020.82. This variation of the intensity can be seen in Fig. 6, where we also show the intensities measured in epochs 2005.89 and 2012.54 with the VLBA (S11a; S14). The lowest intensity was measured in 2012.54 (∑I = 206.31 Jy beam⁻¹) by S14 and it coincides with the minimum of activity registered toward W75N(B) by the 22-m Pushchino telescope in May through July 2012 (Krasnov et al. 2015). The H₂O maser features cover a broad line-of-sight velocity range from +7.8 km s⁻¹ to +26.4 km s⁻¹ over the four EVN epochs and, as previously mentioned, with maser features with line-of-sight velocities of ≥ + 20 km s⁻¹ detected only in epoch 2020.82 (see Table 3 and Fig. 3). We therefore compared the intensities for five different ranges of velocities along the line of sight (ranges I–V) in Fig. 7 and we plot the corresponding maser features superimposed to the continuum emission at Q-band (Rodríguez-Kamenetzky et al. 2020) in the five panels of Fig. 8. We note that each range of velocities has a maximum for the total intensity at different epochs (the corresponding epoch is colored in light gray in Fig. 8). In particular, the maximum is reached earlier for the ranges with the highest velocities, except for the maser features with velocities ≥ + 20 km s⁻¹ (range V) that are detected only in epoch 2020.82. Indeed, we see that the ∑I values of the maser features with +7.8 km s⁻¹ < V_lsr < +10 km s⁻¹ (range I) – which are the only blueshifted features with respect to the systemic velocity of the region (V_lsr = +10.0 km s⁻¹; Shepherd et al. 2003) – and with +10 km s⁻¹ < V_lsr < +13 km s⁻¹ (range II) show the maximum in epochs 2018.44 and 2016.45, respectively. Whereas the ranges +13 km s⁻¹ < V_lsr < +16 km s⁻¹ (range III) and +16 km s⁻¹ < V_lsr < +20 km s⁻¹ (range IV) show their maximum in epoch 2014.46. Furthermore, we note that the maser features of ranges I and II, which are, respectively, within ∓3 km s⁻¹ from the systemic velocity, and those of IV are mainly located along the southwestern tail of the thermal radio jet at Q-band; in particular, the intensity of ranges I and II reaches its maximum when the maser features are along this tail (see Fig. 8). A few maser features of ranges II and IV are close to the position of the strong core of the continuum emission at Q-band and these all but one are detected in the last two epochs, the one is detected in epoch 2014.46 (range II). Most of the maser features of range III are instead aligned east-west below the core of the continuum emission at Q-band and its intensity reaches the maximum in epoch 2014.46. The most redshifted maser features (range V) are the weakest ones among all of those detected in the VLBA and EVN epochs, which can be due to the fact that they have arisen only recently. These maser features are located in the very west edge of the southern tail and slightly north of the core. The spatial coincidence of the blue- and redshifted maser features of ranges I and II may suggest that the maser features are tracing the central part of one of the lobes of the thermal jet and that this has such an inclination that the maser features located on the surface closer to us appear slightly blueshifted, while those located on the opposite surface appear slightly redshifted.

Fig. 6. Comparison of the total sum of the H2O maser intensity (I) measured toward VLA 1 and VLA 2 with the VLBI in epochs 2005.89 (VLBA, S11a), 2012.54 (VLBA, S14), 2014.46 (EVN), 2016.45 (EVN), 2018.44 (EVN), and 2020.82 (EVN). The lowest intensities coincide with the minimum of activity registered toward W75N(B) by the 22-m Pushchino telescope in May-July 2012 (Krasnov et al. 2015).

Fig. 7. Comparison of the total sum of the H2O maser intensity (I) for five different ranges of velocities in VLA 1 as measured in epochs 2012.54 (VLBA, S14), 2014.46 (EVN), 2016.45 (EVN), 2018.44 (EVN), and 2020.82 (EVN).

Fig. 8. Comparison of the H2O maser features detected for five different ranges of velocities in VLA 1 in the four EVN epochs 2014.46, 2016.45, 2018.44, and 2020.82. The size of the symbols are scaled logarithmically according to their peak intensities. The epoch for which the maximum of the intensities is registered for each velocities range is colored in light gray. The assumed velocity of the region is Vlsr = +10.0 km s−1 (Shepherd et al. 2003).

3.2.2. VLA 2

The total intensity of the H₂O maser emission almost constantly increased between the VLBA epoch 2012.54 and the EVN epoch 2020.82 (see Fig. 6). This reflects the fact that the maser features in the four EVN epochs, although there are only half as many of them, are actually brighter than those detected in the two previous VLBA epochs (see Table 3 and S11a; S14). We show the total sum of the H₂O maser spectra detected towards VLA 2 in the four EVN epochs in Fig. 9, whereby the complexity of the H₂O maser emission around VLA 2 is further confirmed.

Fig. 9. Comparison of the total sum of the H2O maser spectra detected toward VLA 2 with the EVN in epochs 2014.46, 2016.45, 2018.44, and 2020.82. The vertical dashed dark-green line indicates the assumed systemic velocity of the region (Vlsr = +10.0 km s−1; Shepherd et al. 2003).

We compare the ∑I for the four zones (see Sect. 3.2.2) in Fig. 10. Here, we see that zone 1 in the four EVN epochs shows increments of intensity between epochs 2014.46 and 2016.45, as well as between epochs 2018.44 and 2020.82, with a slight decrement between epochs 2016.45 and 2018.44. The increments of total intensity are due to the high intensity of individual maser features rather than to their number (see Tables A.2, A.4, A.6, and A.8), indicating that the masing conditions are more favorable inward than outward. Also, the maser features of zone 2 show an increment of ∑I between epochs 2018.44 and 2020.82, similarly to what we see in zone 1 between the same two epochs, this is 127 Jy beam⁻¹ for zone 1 and 137 Jy beam⁻¹ for zone 2.

Fig. 10. Comparison of the total H2O maser intensity in the four zones of VLA 2 (see Fig. 4) measured in epochs 2012.54 (VLBA; S14), 2014.46 (EVN), 2016.45 (EVN), 2018.44 (EVN), and 2020.82 (EVN).

In zone 3, we instead observe a decrement of ∑I between the VLBA epoch 2012.54 and the four EVN epochs. We detected the brightest maser features in the EVN epoch 2014.46 toward south (VLA2.1.01 and VLA2.1.03, each with I ≈ 21 Jy beam⁻¹), while the rest of the maser features of zone 3 show a maser intensity I < 13 Jy beam⁻¹ in all four EVN epochs, with the brightest of them detected toward the center of the continuum emission. The trend of ∑I in zone 4 between epochs 2016.45 and 2020.82 is identical to that observed for the maser features of zone 3.

3.3. Maser polarization

The analysis of the polarized emission from the H₂O maser features detected in the four EVN epochs allow us to compare, in addition to the magnetic field described in Sect. 4, some physical parameters of the maser features and of the gas where they arise.

3.3.1. VLA 1

We note that the highest values of P_l were measured in epochs 2014.46 (15.6%) and 2018.44 (10.6%) and the corresponding maser features are located far southwest of the continuum emission at Q-band (2014.46) and to the west of the bright core of this continuum emission (2018.44). Meanwhile, the lowest values of P_l are all measured in epoch 2020.82 (see Table 3).

As described in Sect. 2, we can estimate some intrinsic characteristics of the maser features by modeling the linearly polarized emission with the FRTM code. The averaged intrinsic linewidth is larger in epoch 2014.46 (⟨ΔV_i⟩ $=3.0_{-0.3}^{+0.2}$ km s⁻¹) and then consistently decreases, to reach a minimum value of $2.6_{-0.3}^{+0.1}$ km s⁻¹ in epoch 2020.82. This implies that the gas temperature of the region, T, has also decreased from epoch 2014.46 to epoch 2020.82. Indeed, from the equation T ≈ 100 × (⟨ΔV_i⟩/0.5)² (Nedoluha & Watson 1992), this is valid if the broadening of the maser line due to the turbulence is negligible, we have: $T_{2014.46}={3600}_{-684}^{+496}$ K and $T_{2020.82}={2704}_{-588}^{+212}$ K. If the H₂O masers are pumped by non-dissociative shocks, the gas can reach temperatures around 4000 K and above this threshold the H₂O molecule is dissociated (Kaufman & Neufeld 1996). However, the estimated high temperatures indicate that the contribution of the turbulence to ΔV_i is not actually negligible. However, if we assume that this contribution is constant in time and everywhere in the source we can still qualitatively compare the estimated temperatures between the different epochs. The averaged emerging brightness temperature (⟨T_bΔΩ⟩) instead increases from epoch 2014.46 to epoch 2020.82 (see Table 3), this might indicate that the saturation level of the maser features increases from one epoch to the other. This is related to the efficiency of the pumping mechanism. Indeed the saturation regime is reached when the stimulated emission rate becomes larger than the decay rate to the upper level of the maser transition, in other words the pumping mechanism is no longer able to provide enough population inversion between the maser levels. In the case of the H₂O maser, the pumping mechanism is due to the shocks produced by the outflow hitting the surrounding matter. Therefore, the shocks should have lost part of their energy from one epoch to the other, which might also be suggested by the estimated gas temperature.

Circular polarization was measured only in epochs 2016.45 and 2018.44 (see Appendix A) when the total H₂O maser intensity reached its highest values. From the three circularly polarized maser features, we measured a magnetic field strength along the line of sight between −676 ± 102 mG (epoch 2016.45; see Table A.3) and −764 ± 22 mG (epoch 2018.44; see Table A.5), which are consistent with what S11a and S14 measured in the previous two VLBA epochs ($-400\mathrm{mG}<B_{||}^{2005.89}<+810$ mG and $-540\mathrm{mG}<B_{||}^{2012.54}<+18$ mG). We note that the positive and negative signs of B_|| indicate that the magnetic field is either pointing away or toward the observer, respectively.

3.3.2. VLA 2

We note that all the linearly polarized maser features in the four EVN epochs, except for VLA2.1.24 (P_l = 4.6 %), show a P_l value in the range between 0.2% and 2.7% (see Table 3). The outputs of the FRTM code provides consistent T_bΔΩ values in all EVN epochs, only four and one maser features show T_bΔΩ on the order of 10⁶ and 10¹⁰ K sr, respectively, and consistent ΔV_i in three of the four EVN epochs. In epoch 2014.46, we have ΔV_i< 2 km s⁻¹ while in the other three EVN epochs, ΔV_i> 2 km s⁻¹ with a few exceptions (VLA2.2.22, VLA2.3.21, and VLA2.3.23). As we already made previously for VLA 1 (Sect. 3.3.1), and keeping in mind that the values do not actually indicate the actual temperature of the gas, we can estimate the mean temperature from the ΔV_i values and we then have $T_{2014.46}={1024}_{-240}^{+272}$ K, $T_{2016.45}={3364}_{-448}^{+480}$ K, $T_{2018.44}={1764}_{-608}^{+1600}$ K, and $T_{2020.82}={3364}_{-864}^{+236}$ K. We note that all the linearly polarized maser features in epoch 2014.46 are located northeast, where the expanding gas is supposed to encounter a denser medium.

We were also able to measure B_|| from the circularly polarized maser emission of a total of five maser features in the three EVN epochs 2016.45 (−1498 ± 225 mG and −2426 ± 364 mG, Table A.4), 2018.44 (−355 ± 69 mG and +439 ± 66 mG, Table A.6), and 2020.82 (−452 ± 68 mG, Table A.8), all of them located north to northeast. These values are larger than those measured in epochs 2005.89 ($-186\mathrm{mG}<B_{||}^{2005.89}<+957$ mG; S11a) and 2012.54 ($-152\mathrm{mG}<B_{||}^{2005.89}<-103$ mG; S14) with the VLBA.

4. Magnetic field

We discuss the magnetic field around VLA 1 and VLA 2 in Sects. 4.1 and 4.2, respectively. Here, we are able to estimate and compare the orientation of the magnetic field from the linear polarization vectors measured in the four EVN epochs. For each epoch, we determined the error-weighted orientation of the magnetic field (⟨Φ_B⟩), listed in Table 3. We note that the position angle of a magnetic field vector on the plane of the sky has three values contemporary: Φ_B and Φ_B ± 180°. Therefore, we can state that the magnetic field vectors have a sort of periodicity of 180° that exists only as a consequence of how Φ_B is defined on the plane of the sky: positive if measured counterclockwise from north and negative if measured clockwise from the north.

4.1. VLA 1

We plot ⟨Φ_B⟩ as function of time in the left panel of Fig. 11, as measured in the two previous VLBA epochs (2005.89 and 2012.54) and in the four EVN epochs. Here, we also show a 180° periodicity linear fit (dashed gray lines) that has a slope of $m_{\mathrm{VLA}1}=(-15°.38\pm 1°.80){\mathrm{yr}}^{-1}$. However, the apparent rotation is not due to calibration uncertainties (see Table 1) or to maser polarization variability, which can influence the relation between the magnetic field orientation and the polarization vectors only if the H₂O maser features are highly saturated (this is not our case), but it is simply the consequence of averaging the Φ_B angles that are estimated from linear polarization vectors measured in different locations along VLA 1, where the magnetic field is actually differently oriented. Therefore, a punctual comparison between the magnetic field vectors measured in the different epochs is necessary, but this cannot be done with common polarized maser features detected in consecutive epochs because of the maser variability on short timescales. Instead, what we can do is a comparison of all the magnetic field vectors estimated in all the epochs. Indeed, we plot each single magnetic field vector estimated from all the linearly polarized maser features detected in the four EVN epochs in the right panel of Fig. 11. Here, we can see that the magnetic field vectors estimated in similar location, but in different epochs, seem to represent a quasi-static magnetic field. We can therefore consider the magnetic field vectors estimated from the H₂O maser features in one epoch as representative of the magnetic field in those locations rather than in that time. Consequently, we can gather the magnetic field vectors of all the EVN epochs and consider them as measurements carried out at the same time. This allows us to compare the magnetic field with the continuum emission at Q-band (see Fig. 11). The magnetic field vectors seem to follow the morphology of the continuum emission, in particular, the internal vectors are in accordance with the radio continuum contours (right panel of Fig. 11). This suggests that the magnetic field is along the thermal jet and it bends toward south at the southwest end of the thermal jet and toward the north at the northeastern end. In addition, we note that the magnetic field orientation in the northeast and in the far southwest coincides with that of the large-scale magnetic field vectors reported by Palau et al. (2021, their Fig. 2), showing the 1.3 mm polarized continuum emission of W75N(B) as observed by Alves et al. (in prep.).

Fig. 11. Multi-epoch comparison of ⟨ΦB⟩ angles for VLA 1 (left panel). The epochs are 2005.89 (VLBA, S11a), 2012.54 (VLBA, S14), 2014.46 (EVN), 2016.45 (EVN), 2018.44 (EVN), and 2020.82 (EVN). The two horizontal solid black lines indicate the ±90° angles. The dashed gray lines represents the best linear fit results of the data considering that the position angle of a magnetic field vector has contemporary multiple values equal to ⟨ΦB⟩=⟨ΦB⟩±180°. The slopes of these lines are mVLA 1 = ( − 15° ​.38 ± 1° ​.80) yr−1. Magnetic field vectors and strength along the line of sight (next to the corresponding maser feature) as estimated from all the linearly and circularly polarized maser features detected around VLA 1 (right panel) during the four EVN epochs (2014.46, 2016.45, 2018.44, and 2020.82). The vectors are superimposed to the uniform-weighted continuum map at Q-band (central frequency 44 GHz) of the thermal jet driven by VLA 1 (Rodríguez-Kamenetzky et al. 2020). For more details, see the caption for Fig. 2.

The three circularly polarized maser features are associated with the ends of the thermal jet as observed at Q-band (see Fig. 11) and the estimated magnetic field points toward us (see Sect. 3.3.1). The 3D magnetic field strength can be estimated from |B| = ⟨|B_|||⟩/cos⟨θ⟩ (if θ ≠ ±90° ), where ⟨|B_|||⟩ and ⟨θ⟩ are the error-weighted mean values of |B_||| and θ, respectively. We must note that the magnetic field derived from the circularly polarized emission of the H₂O maser features can be very high because these masers probe shocked gas, so it is not representative of the whole region. We then have in the two EVN epochs |B|^2016.45 > 2.0 G, which is a lower limit obtained by considering θ = 71°, and ⟨|B|^2018.44⟩ = 3.6 G. Both these values are much larger than what was measured on average in the VLBA epochs (⟨|B|^2005.89⟩≈0.7 G and ⟨|B|^2012.54⟩> 0.1 G; S11a; S14). The differences might be due to the increment of the magnetic field strength or to the hydrogen number density of the 22 GHz H₂O maser. In the latter case, by considering the relation $\left|B\right|\propto n_{{\text{H}}_2}^{0.65}$ (Crutcher et al. 2010), we can estimate an increment of $n_{{\text{H}}_2}^{{\text{H}}_2\text{O}}$ of one order of magnitude between the VLBA epoch 2005.89 and the EVN epoch 2018.44.

4.2. VLA 2

Similarly to VLA 1, we observe an apparent clockwise rotation of ⟨Φ_B⟩ from the VLBA epoch 2015.89 to the EVN epoch 2020.82 (see left panel of Fig. 12). In this case we notice that ⟨Φ_B⟩ measured in the EVN epoch 2020.82 is on a different 180° periodicity line with respect to those of the previous EVN epochs. The slope of the 180° periodicity lines is $m_{\mathrm{VLA}2}=(-19°.73\pm 1°.89){\mathrm{yr}}^{-1}$, which might indicate that the rotation in VLA 2 is faster than in VLA 1. As explained in the case of VLA 1, the apparent rotation of ⟨Φ_B⟩ can be affected by the location of each maser feature for which we are able to estimate the magnetic field vectors and it cannot therefore be taken at face value. Indeed, despite our earlier claims (S14), we now think there is no evidence for such a rotation. A comparison of the magnetic field vectors estimated in the four EVN epochs is reported on the right panel of Fig. 12. Two main aspects can be observed here, firstly almost all the magnetic field vectors are located northeast and secondly the magnetic field vectors estimated in one epoch can be considered representative of a quasi-static magnetic field in those locations rather than in that time. The magnetic field is generally (within the errors) perpendicular to the expansion velocities of the gas, which are represented with arrows in Fig. 12, all around VLA 2, and it is parallel only in the northeast where the expanding gas encounters the denser matter. The fact that the magnetic field vectors are along the shock front is not surprising. Indeed, H₂O masers arise behind fast C- or J-type shocks (e.g., Kaufman & Neufeld 1996; Hollenbach et al. 2013) that alter the initial magnetic field configuration, while propagating in the ambient medium. As explained in detail in Goddi et al. (2017), the magnetic field component perpendicular to the shock velocity is compressed and might dominate the parallel component, which remains unaffected; consequently, the resulting magnetic field probed by the H₂O maser features might be expected to be along the shock front. The sudden change of the orientation of the magnetic field vectors in the northeast cannot be explained with the typical 90°-flip because the estimated θ angles of those H₂O maser features are much greater than 55° (see Sect. 2). Therefore, this can be justified with a variation of the ratio between the perpendicular and the parallel components of the magnetic field with respect to the expanding velocity. In other words, even though the shock is still able to pump the H₂O molecules in order to produce the maser emission, it is not able anymore to compress enough the perpendicular component of the magnetic field in order to make it dominating over the parallel component. This can be considered a further clue in favor of the presence of a very dense gas in the northeast where the magnetic field is oriented in a northeast-southwestly direction.

Fig. 12. Multi-epoch comparison of ⟨ΦB⟩ angles for VLA 2 in the left panel. The epochs are 2005.89 (VLBA, S11a), 2012.54 (VLBA, S14), 2014.46 (EVN), 2016.45 (EVN), 2018.44 (EVN), and 2020.82 (EVN). The gray point indicate the ⟨ΦB⟩ angle measured in 2020.82 after a rotation of 180°, considering that ⟨ΦB⟩=⟨ΦB⟩±180°. The two horizontal solid black lines indicates the ±90° angles. The dashed gray lines represents the best linear fit results of the data considering that the position angle of a magnetic field vector has contemporary multiple values equal to ⟨ΦB⟩=⟨ΦB⟩±180°. The slopes of these lines are mVLA 2 = ( − 19° ​.73 ± 1° ​.89) yr−1. Right panel: magnetic field vectors and strength along the line of sight (next to the corresponding maser feature) as estimated from the linearly and circularly polarized maser features detected around VLA 2 during the four EVN epochs (2014.46, 2016.45, 2018.44, and 2020.82) and superimposed to the natural-weighted continuum maps at K-band (central frequency 22 GHz) of the thermal, collimated ionized wind emitted by VLA 2 (Carrasco-González et al. 2015). For more details see caption of Fig. 4. The arrows represent the direction of the proper motions on the plane of the sky reported in Tables 4–6, and the arrows length is proportional to their magnitudes. The double arrow in the northeast indicates the uncertainties of the motion direction, due to the presence of an inward motion (case A) or of an outward motion (case B) of the gas, as explained in Sect. 3.1.2.

We also compared the magnetic field vectors with the morphology of the continuum emission at K-band obtained by Carrasco-González et al. (2015) with the VLA in Fig. 12. Here, we see that some vectors seems to follow the morphology of the continuum emission suggesting a possible correlation between the magnetic field orientation and the shaping of the continuum emission. In addition, the orientation of the magnetic field vectors in the northeast of VLA 2 are in good agreement with the values reported by Palau et al. (2021) at large scales (∼1.3″), using submillimeter continuum data thus indicating that the magnetic field at small scales might connect with the magnetic field at large scales.

To properly compare the magnetic field strength measured in different epochs, we need to estimate |B| (see Sect. 4.1). We have |B_VLA2.2.17| = 9.5 G and |B_VLA2.2.27| = 27.8 G in epoch 2016.45, |B_VLA2.3.14|> 0.7 G and |B_VLA2.3.39|> 0.6 G in epoch 2018.44 (assuming a lower value of θ = 53°, see Table 3), and |B_VLA2.4.27| = 0.6 G in epoch 2020.82. As we mention in Sect. 4.1, a variation of the magnetic field strength might largely be due to a variation of n_H2 of the gas where the H₂O maser features arise. This variation follows the empirical equation of Crutcher et al. (2010) reported in Sect. 4.1. Considering the measurements of B made in the northeast of VLA 2 (zone 1), we can estimate n_H2 of the gas in the five epochs by assuming that the lowest value must be at least 10⁸ cm⁻³, which is the lowest limit for having H₂O maser emission at 22 GHz (Elitzur et al. 1989) and this can be assumed for the epoch with the weakest magnetic field strength that we measured. These values are reported in Table 7, where we name the n_H2 of the gas as n₁ − n₅ from epoch 2005.89 to epoch 2020.82, respectively. We note that n_H2 decreases from the VLBA epoch 2005.89 to the VLBA epoch 2012.54, when the morphology of the maser feature distribution changed from quasi-circular to elliptical (Kim et al. 2013; S14). This might be due to the dilution of the gas. Furthermore, in this transition we also observe a change of sign of the magnetic field from positive to negative and it therefore indicates that the magnetic field changes its pointing direction from away to toward the observer. From the VLBA epoch 2012.54 to the EVN epoch 2016.45, when the expansion of the gas in the northeast continues, n_H2 increases to the highest possible value for having H₂O maser emission and this might be due to the encounter of the expanding gas with the denser medium. This might be a further indication that the absence of H₂O maser emission farther in the northeast is due to the presence of an overly high density medium (case A in Sect. 3.1.2). Finally, n_H2 shows the lowest possible value in the EVN epochs 2018.44 and 2020.82 when the H₂O masers are probing the gas in the same region as they did in the VLBA epoch 2012.54. This suggests that in this zone not only the maser features in the EVN epoch 2020.82 are pumped by a different shock, as we suppose in Sect. 3.1.2, but that likely all the maser features in all the EVN epochs are pumped by different shocks in this zone.

5. Discussion

In this section, we put together all the pieces of information that we have presented in previous sections to provide a full picture of both VLA 1 and VLA 2.

5.1. VLA 1

We can determine that the 22 GHz H₂O maser features around VLA 1 are probing the passage throughout the gas where the masers arise of a nondissociative shock. This shock is produced by the expansion of a thermal jet from VLA 1. The magnetic field is along the axis of the thermal jet, pointing toward us, and it bends toward south and north at the southwest and northeast ends of the thermal jet, respectively, following the large-scale magnetic field morphology (Palau et al. 2021). From the different maser and gas characteristics observed between the VLBA and EVN epochs, we can also conclude that the compression shock probed by the maser features in the recent epochs is a different one than that probed in the past. It is also possible that the maser features in the six epochs were pumped by a series of shocks rather than a single one. The presence of more shocks might be justified by episodic variations in the velocity of the jet as observed for instance in the intermediate-mass protostar in the Serpens star-forming region (Rodríguez-Kamenetzky et al. 2022).

5.2. VLA 2

In VLA 2, we have an asymmetric expansion of the gas on the plane of the sky that depends on the particular direction of the motion. Indeed, the expansion of the gas is prevented by the presence of a dense gas in the northeast, in the northwest the gas is moving outward with a velocity around 26–28 km s⁻¹; while in the center, the motion is toward southeast with a velocity around 38 km s⁻¹; in the east, the motion of the gas is toward south with a velocity of 12 km s⁻¹; in the south, the gas is expanding southward with a velocity of 8 km s⁻¹; and in the southwest the gas is the fastest one with a westward velocity of 78 km s⁻¹. If we average the magnitude of all the velocities that we estimated we get |V|≈33 km s⁻¹ that is similar to the symmetric expanding velocity of 30 km s⁻¹ measured by S14, who considered the fitted ellipses centered to a common center. Furthermore, the maser features in the northeast show both blue- and redshifted V_lsr without an ordered spatial distribution suggesting that the gas is moving along the walls of the dense core that the gas encounters. In the northwest and southeast, the masers show blueshifted velocities, that is, the gas is moving toward us; while in the center and southwest the masers are all redshifted, that is, the gas is moving away from us. The comparison with the K-band continuum emission (Carrasco-González et al. 2015, see Figs. 4 and 12) suggests that the gas in the southwest can expand without encountering any obstacle, while in the north there might be a very dense medium that may be part of the envelope assumed by Carrasco-González et al. (2015), which is able to slow down the expanding gas in the northwest and to stop completely the expansion in the northeast. The presence of this denser medium is also responsible to alter the morphology and strength of the magnetic field. The magnetic field is generally perpendicular to the proper motions all around VLA 2, but in the northeast, it becomes parallel and stronger after encountering the denser medium. This might be the consequence of the inefficiency of the compression of the gas due to the passage of the shock, which makes it so that the perpendicular component of the magnetic field no longer dominates over the parallel component, with respect to the expansion direction. Furthermore, the magnetic field morphology in the northeast region of VLA 2 agrees with the large-scale magnetic field reported by Palau et al. (2021). Regarding the 3D magnetic field strength, we find that it is higher than previously measured only in epoch 2016.45, where the maser features are detected in the northeast high density region of VLA 2.

6. Summary

We observed the polarized emission of 22 GHz H₂O maser around the radio sources VLA 1 and VLA 2, located in the HMSFR W75N(B), with the EVN on four epochs, separated by two years from 2014 to 2020. We detected both linearly and circularly polarized emission in all epochs but one (epoch 2014.46) around both the radio sources. A comparison of the maser distributions and the magnetic fields among the four EVN epochs shows that:

• The 22 GHz H₂O maser emission does not probe the magnetic field as it is in a single observing epoch, but it probes a portion of a quasi-static magnetic field in the region where the maser arises. That is, the magnetic field vectors estimated in one epoch must be considered as a measurement of the magnetic field in those locations rather than in that time. Therefore, mapping the magnetic field in different epochs with the 22 GHz H₂O maser allows us to reconstruct the magnetic field in a larger area.

• In VLA 1, the 22 GHz H₂O maser features are probing the passage of a nondissociative shock produced by the expansion of the thermal radio jet of VLA 1. The magnetic field is along the axis of the thermal radio jet and bends toward south and north at the southwest and northeast ends of the jet. The magnetic field strength along the line of sight measured from the Zeeman splitting of the maser line is in the range $-764\mathrm{mG}<B_{||}^{\mathrm{EVN}}<-676$ mG and is consistent with previous VLBA measurements.

• The 22 GHz H₂O maser features around VLA 2 are tracing an asymmetric expansion of the gas that is actually halted in the northeast, where the gas likely encounters a very dense medium. The inferred magnetic field vectors are almost all perpendicular to the proper motion, but in the northeast, where the expanding gas encounters the supposed denser medium, the vectors become parallel to the expansion direction. The magnetic field strengths along the line of sight are all measured in the north-northeast and their values are in the range of $-355\mathrm{mG}<B_{||}^{\mathrm{EVN}}<-2426$ mG, which is greater than previous measurements with the VLBA.

In conclusion, W75N(B) is among the best HMSF laboratories and it is thanks to the presence of multiple YSOs at different evolutionary stages that it can play a crucial role in shedding light on the formation process of high-mass stars. In addition, we have demonstrated the importance of VLBI monitoring observations of polarized maser emission in reconstructing the morphology of magnetic field close to massive YSOs.

Acknowledgments

We wish to thank the anonymous referee for the useful suggestions that have improved the paper. The European VLBI Network is a joint facility of independent European, African, Asian, and North American radio astronomy institutes. Scientific results from data presented in this publication are derived from the following EVN project code(s): ES074. G.S. thanks Jorge Cantó for the useful discussion and suggestions. J.M.T. acknowledges partial support from the PID2020-117710GB-100 grant funded by MCIN/AEI/10.13039/501100011033, and by the program Unidad de Excelencia María de Maeztu CEX2020-001058-M. J.F.G. acknowledges support from grants PID2020-114461GB-I00 and CEX2021-001131-S, funded by MCIN/ AEI /10.13039/50110001103. S.C. acknowledges support from UNAM, and CONACyT, México. This work was partially supported by FAPESP (Fundação de Amparo à Pesquisa do Estado de São Paulo) under grant 2021/01183-8.

References

Appendix A: Detailed Results.

We report here the detailed results obtained for each EVN epoch, from epoch 2014.46 (Section A.1) to epoch 2020.82 (Section A.4). In Fig. A.1, we show the 22 GHz H₂O maser distributions around VLA 1 and VLA 2 for the four epochs. In Tables A.1–A.8, we list all the H₂O maser features detected towards VLA 1 and VLA 2 in the four EVN epochs. The tables are organized as follows. The name of the feature is reported in col. 1. The position offset with respect to the reference maser feature are reported in cols. 2 and 3, as well as in col. 4, when it is the case, we report the zone that the maser feature belongs to. The peak intensity, the LSR velocity, and the FWHM of the total intensity spectra of the maser features, that are obtained using a Gaussian fit, are reported in cols. 5, 6, and 7, respectively. The linear polarization fraction and the linear polarization angles are instead reported in cols. 8 and 9, respectively. The best-fitting results obtained by using the FRTM code are reported in cols. 10 (the intrinsic linewidth), 11 (the emerging brightness temperature), and 14 (the angle between the magnetic field and the maser propagation direction). The later is used to solve the 90° ambiguity of the magnetic field orientation with respect to the linear polarization vector (see Sect. 2). The circular polarization fraction and the estimated magnetic field strength along the line of sight are finally reported in cols. 12 and 13, respectively.

Fig. A.1. Close-up view of the 22 GHz H2O maser features detected around the radio sources VLA 1 (left panels) and VLA 2 (right panels) in epochs 2014.46 (top panels), 2016.45, 2018.44, and 2020.82 (bottom panels). The reference positions are reported in Table 2. The letters A, B, and C indicate the position of the three maser groups identified around VLA 1 by S11a and they are reported here only for reference. The octagonal symbols are the identified maser features scaled logarithmically according to their peak intensity (Tables A.1-A.8). Maser LSR radial velocities are indicated by color (the assumed systemic velocity of the region is Vlsr = +10.0 km s−1; Shepherd et al. 2003). The linear polarization vectors, scaled logarithmically according to the polarization fraction Pl, are overplotted. In the bottom-right corner of all panels the error-weighted orientation of the magnetic field (⟨ΦB⟩) is also reported; the two dashed segments indicate the uncertainties.

A.1. Epoch 2014.46

We detected 28 and 43 H₂O maser features around VLA 1 and VLA 2, respectively (see Tables A.1 and A.2). The maser features around VLA 1 have V_lsr ranging from +7 km s⁻¹ and +20 km s⁻¹, with peak intensities in the range 0.4 Jy beam⁻¹< I ≲ 1200 Jy beam⁻¹. The distribution is still linear from southwest to northeast as previously observed (S11a; S14), although we did not detect any maser of groups A and C as defined by S11a and indicated in Fig. A.1. The overall maser distribution around VLA 2, with H₂O maser features with −12 km s⁻¹< V_lsr < +22 km s⁻¹ and 0.2 Jy beam⁻¹< I < 32 Jy beam⁻¹, is identical to the distribution detected by Kim & Kim (2018) with the KVN & VERA Array (KaVA) and it is similar to the elliptical distribution observed in 2012.54 by S14, even though no maser feature was detected in the west-southwest. Actually, we note that most of the maser features (72%, 31/43) are detected along the arc structure in the north, which is less rounded than that observed in 2012.54 (S14).

We detected linearly polarized emission toward twelve and five maser features in VLA 1 (P_l = 0.4%−15.6%) and VLA 2 (P_l = 0.9%−4.6%), respectively. The maser feature VLA1.1.03 shows the highest mean linear polarization fraction measured among the four epochs, that is, P_l = 15.6%, but it is much lower than the highest P_l ever measured toward W75N(B) (P_l = 25.7%, S11a). The error-weighted linear polarization angles are ⟨χ⟩_VLA 1 = −47°±30° and ⟨χ⟩_VLA 2 = −40°±7°. The FRTM code was able to properly fit eight and four H₂O maser features in VLA 1 and VLA 2, respectively, and for all of them it provided θ angles much greater than 55° indicating that the magnetic field is perpendicular to the linear polarization vectors. No circular polarization was detected at 3σ_s.−n., which implies an upper limit for the brightest maser feature VLA1.1.15 of P_V < 0.01% (see Table 1).

A.2. Epoch 2016.45

We list in Tables A.3 and A.4 the 20 (0.2 Jy beam⁻¹< I < 1620 Jy beam⁻¹) and 37 (0.2 Jy beam⁻¹< I < 91 Jy beam⁻¹) H₂O maser features detected toward VLA 1 and VLA 2, respectively. All the maser features detected around VLA 1, but one, have LSR velocities between +7.8 km s⁻¹ and +16.4 km s⁻¹. The maser feature VLA1.2.16 (I = 0.24 Jy beam⁻¹) is the most blueshifted H₂O maser feature ever detected toward VLA 1, with V_lsr = +2.68 km s⁻¹. The maser distribution around VLA 1 is still linear and only VLA1.2.20 is located about 100 mas northeast of an area of 70 mas × 60 mas, where all the other maser features arose. The elliptical distribution of the H₂O maser features around VLA 2 (−16 km s⁻¹< V_lsr < +21 km s⁻¹) is also confirmed. In this epoch, we detected six maser features (VLA2.2.01–VLA2.2.06) in the southwest, these are all redshifted (V_lsr > 15 km s⁻¹, see Fig. A.1). One of the maser feature (VLA2.2.37) was actually partially detected because the 2-MHz bandwidth did not cover the entire line emission, and we are therefore able to provide only a lower limit of the peak intensity (I > 30 Jy beam⁻¹). This detection led us to double the observed bandwidth in the next two EVN epochs to recover the maser emission at V_lsr < −15.7 km s⁻¹ (see Sect. 2).

Seven and 13 H₂O maser features associated with VLA 1 (P_l = 0.7%−2.2%) and VLA 2 (P_l = 0.6%−2.7%) showed linearly polarized emission, respectively. The error-weighted linear polarization angles are ⟨χ⟩_VLA 1 = −60°±16° and ⟨χ⟩_VLA 2 = −74°±19°. We were able to properly fit with the FRTM code about 40% (3 out of 7) of the linearly polarized maser features in VLA 1 and ∼70% (9 out of 13) of those in VLA 2. We got θ > 55° for all of them, and the magnetic field is therefore perpendicular to the linear polarization vectors. In this epoch we were also able to measure circular polarization toward three maser features, one in VLA 1 (VLA1.2.01, P_V = 3.5%) and two in VLA 2 (VLA2.2.17 and VLA2.2.27, P_V = 4.9% and 7.8%, respectively). Their spectra are shown in Fig. A.2, where the best-fit model obtained from the FRTM code are overplotted as thick red lines. We note that P_V of VLA2.2.27 is the highest ever measured toward both sources VLA 1 and VLA 2, while the other two maser features have P_V consistent with previous detections (S11a; S14).

Fig. A.2. Total intensity spectra (I, upper panel) and circular polarization intensity spectra (V, lower panel) for the H2O masers VLA1.2.01, VLA2.2.17, VLA2.2.27, VLA1.3.04, VLA1.3.16, VLA2.3.14, VLA2.3.39, and VLA2.4.27 (see Tables A.3, A.4, A.5, A.6, and A.8). The thick red lines are the best-fit models of I and V emission obtained using the outputs of the FRTM code for 22 GHz H2O masers (see Sect. 2), in particular these models are used to fit the spectra (for more details see S11a). The maser features were centered on zero velocity.

A.3. Epoch 2018.44

Toward VLA 1 we detected 20 H₂O maser features in the velocity range +7.7 km s⁻¹≤V_lsr ≤ +15.0 km s⁻¹ and with peak intensities between 2 Jy beam⁻¹ and 450 Jy beam⁻¹ (see Table A.5). These maser features are linearly distributed from southwest to northeast (∼270 au) and are associated only with group B of S11a. We detected 44 H₂O maser features (0.04 Jy beam⁻¹≤I ≤ 128 Jy beam⁻¹) elliptically distributed around VLA 2 and with a wide range of velocities (−16.0 km s⁻¹< V_lsr < +28.3 km s⁻¹; see Table A.6).

The FRTM code was able to properly fit all the maser features that showed linearly polarized emission, these are 12 in VLA 1 (P_l = 0.7%−10.6%) and 10 in VLA 2 (P_l = 1.0%−2.3%). The error-weighted linear polarization angles are ⟨χ⟩_VLA 1 = −87°±5° and ⟨χ⟩_VLA 2 = +76°±12°. Differently from the previous two EVN epochs, the magnetic field is not always perpendicular to the linear polarization vectors. Indeed, for two maser features in VLA 1 (VLA1.3.06 and VLA1.3.07) and two in VLA 2 (VLA2.3.23 and VLA2.3.28) the FRTM code provided such θ angle values that the magnetic field is more likely parallel than perpendicular. To determine the relative orientation of the magnetic field, we consider the associated errors ε^± of θ, where the plus and minus signs indicate the positive and negative errors respectively, and in particular the quantities θ^± = θ + ε^±. If |θ⁺ − 55°| < |θ⁻ − 55°| the magnetic field is more likely parallel to the linear polarization vectors, while if |θ⁺ − 55°| > |θ⁻ − 55°| the magnetic field is more likely perpendicular (e.g., Surcis et al. 2015). We also detected circularly polarized emission toward two maser features in VLA 1 (VLA1.3.04 and VLA1.3.16) and two maser features in VLA 2 (VLA2.3.14 and VLA2.3.39). Whereas we were able to fit the V spectra of VLA1.3.04 and VLA1.3.16 with the model obtained by using the outputs of the FRTM code, we had to find the best models for VLA2.3.14 and VLA2.3.39 by considering the values of T_bΔΩ and ΔV_i that best fit their total intensity spectra. This was necessary because both VLA2.3.14 and VLA2.3.39 did not show any linearly polarized emission and therefore we could not use the FRTM code to determine them. In particular, we found that the best estimates of T_bΔΩ and ΔV_i are equal to 10⁶ K sr and 1.4 km s⁻¹ for VLA2.3.14, respectively, and T_bΔΩ  = 3.2 × 10⁹ K sr and ΔV_i  = 2.0 km s⁻¹ for VLA2.3.39. The results of the fit can be seen in Fig. A.2.

A.4. Epoch 2020.82

In the fourth and last EVN epoch, we detected 10 and 39 H₂O maser features toward VLA 1 and VLA 2 (see Tables A.7 and A.8), respectively. The maser features in VLA 1 (0.1 Jy beam⁻¹< I < 303 Jy beam⁻¹) are compactly located around the two brightest maser features VLA1.4.05 (V_lsr = 9.87 km s⁻¹) and VLA1.4.10 (V_lsr = 10.87 km s⁻¹), which are about 47 mas apart (∼60 au, see Fig. A.1) and are associated with group B of S11a. The 70% of the H₂O maser features detected along VLA 1 have velocities in the range +18.9 km s⁻¹< V_lsr < +26.4 km s⁻¹. The last time that H₂O maser emission with V_lsr > +19 km s⁻¹ was detected toward VLA 1 was in 2005 (S11a), but those maser features were associated only with group A of S11a and they were at ∼500 mas (∼650 au) northeast of group B. The H₂O maser features elliptically distributed around VLA 2 covers the smallest range of velocities (−1.0 km s⁻¹< V_lsr < +27.2 km s⁻¹) and the largest range of peak intensities (0.05 Jy beam⁻¹≤I < 286 Jy beam⁻¹) among the four EVN epochs.

The highest linear polarization fraction of epoch 2020.82 was measured in VLA 2 and its value is P_l = 1.9%. In particular, we detected linearly polarized emission toward a total of seven maser features, with two in VLA 1 and five in VLA 2 (see Tables A.7 and A.8). The FRTM code properly fit all of them, providing θ angles for VLA1.4.02 and VLA1.4.05 lower than 55°, which implies a magnetic field parallel to the linear polarization vectors. In the case of VLA2.4.31, we find that the magnetic field is more likely to be parallel and, indeed, for this feature, we have |θ⁺ − 55°| < |θ⁻ − 55°|, while for all the other four linearly polarized maser features in VLA 2, the magnetic field is perpendicular to the linear polarization vectors (see Table A.8). Only one maser feature shows circular polarization (VLA2.4.27, P_V = 1.3%), but because it does not show any linear polarization, we had to determine the best estimates of T_bΔΩ and ΔV_i for modeling its V spectra (see Fig. A.2). These values are 10⁸ K sr and 3.0 km s⁻¹, respectively.

Appendix B: Proper motion estimate of VLA 1

To properly compare the positions of the H₂O maser features detected in the four EVN epochs with the continuum emissions of VLA 1 and VLA 2 observed in 2014, we must verify whether the two massive YSOs have any proper motion within the region W75N(B). We report in Table B.1 the absolute positions of the three massive YSOs VLA 1, VLA 2, and VLA 3, as measured in the K-band with the VLA by Torrelles et al. (1997), Carrasco-González et al. (2015), and Rodríguez-Kamenetzky et al. (2020) in epochs 1996.96 and 2014.20. These positions were obtained by deconvolving the data of the two epochs with the same restoring circular beam of 0.​​″12, as done by Carrasco-González et al. (2015) and Rodríguez-Kamenetzky et al. (2020), and by fitting the emission with a Gaussian fit. Furthermore, as already done for the maser features of the last three EVN epochs, the positions of the three YSOs of epoch 2014.20 were also corrected, assuming the proper motion of the region W75N(B) equal to the median proper motion measured for the 6.7 GHz CH₃OH maser features by Rygl et al. (2012), ⟨μ_α⟩ = (−1.97 ± 0.10) mas yr⁻¹ and ⟨μ_δ⟩ = (−4.16 ± 0.15) mas yr⁻¹.

In the three panels of Fig. B.1, we plot the corrected positions of the three massive YSOs for the two epochs with the relative errors due to the thermal noise and to the Gaussian fit errors (Δα, Δδ), which are represented by squares and with the restoring circular beam centered on each position (circles). We note from Fig. B.1 that the positions of VLA 2 and VLA 3 in both epochs are consistent within the beam and from Table B.2 that the position shifts between the two epochs are similar suggesting a possible systematic uncertainty. We cannot verify if this uncertainty is real or not because the structure of the two sources varied between the two epochs, in particular, that of VLA 2. We can therefore assume that VLA 2 and VLA 3 did not move within W75N(B) over 17.24 years and that the positions of the H₂O maser features in VLA 2 in the last three EVN epochs do not need further corrections.

Fig. B.1. Comparison of the positions of the continuum emission at K-band of VLA 1 (left panel), VLA 2 (central panel), and VLA 3 (right panel) as observed in epoch 1996.96 (solid black; Torrelles et al. 1997) and epoch 2014.20 (dashed red; Carrasco-González et al. 2015; Rodríguez-Kamenetzky et al. 2020). The squares represent the uncertainties of the peaks positions (see Table B.1) and the circles are the restoring circular beam of 0.​​″12 used to properly compare the two epochs. Each panel is centered at the corresponding position of the YSO in epoch 1996.96 (see Table B.1). The positions of epoch 2014.20 are corrected for the proper motion of the region W75N(B) as measured by Rygl et al. (2012), ⟨μα⟩ = (−1.97 ± 0.10) mas yr−1 and ⟨μδ⟩ = (−4.16 ± 0.15) mas yr−1.

It is evident from both Fig. B.1 and Table B.2 that VLA 1 actually moved within W75N(B) between the two epochs. From the shift reported in Table B.2 we measure a proper motion of VLA 1 within W75N(B) of ⟨μ_α⟩^VLA 1 = (−6.3 ± 0.4) mas yr⁻¹ (−38.9 ± 2.5 km s⁻¹) and ⟨μ_δ⟩^VLA 1 = (+5.0 ± 0.4) mas yr⁻¹ (+30.9 ± 2.5 km s⁻¹) along the right ascension and declination, respectively. This motion might suggest that VLA 1 is a runaway protostar, but to verify this possibility it is necessary to carry-on ad-hoc multi epoch observations. Nevertheless, we decide to correct the positions of the H₂O maser features detected in the last three EVN epochs and associated with VLA 1 accordingly to the estimated proper motion of VLA 1.

All Tables

All Figures

	Fig. 1. Contour maps of the phase-reference calibrator J2040+4527 in the EVN epochs 2014.46, 2016.45, 2018.44, and 2020.82 from left to right, respectively. Contours are 4, 8, 16, and 32 times the 0.6 mJy beam−1. The restoring beam is overplotted on the bottom left corner of each panel (see Table 1).
In the text

Fig. 2. Comparison of the H2O maser features detected toward VLA 1 (left panel) in the four EVN epochs (2014.46, 2016.45, 2018.44, and 2020.82) and superimposed to the uniform-weighted continuum map at Q-band (central frequency 44 GHz) of the thermal radio jet driven by VLA 1 obtained with the VLA by Rodríguez-Kamenetzky et al. (2020). The light gray contours are 5, 6, 7, 8, 10, 12, 14 times σ = 100 μJy beam−1 and the VLA beam is shown on the bottom right corner. The positions of the maser features are corrected assuming the proper motion of the region W75N(B) equal to the median proper motion measured for the 6.7 GHz CH3OH maser features by Rygl et al. (2012), ⟨μα⟩ = (−1.97 ± 0.10) mas yr−1 and ⟨μδ⟩ = (−4.16 ± 0.15) mas yr−1. In addition, the positions of the maser features associated with VLA 1 and detected in the last three epochs have been further corrected considering the proper motion of VLA 1 within W75N(B). This proper motion is ⟨μα⟩VLA 1 = ( − 6.3 ± 0.4) mas yr−1 and ⟨μδ⟩VLA 1 = ( + 5.0 ± 0.4) mas yr−1 (see Appendix B). The size of the symbols are ten times the uncertainties of the absolute positions of the maser features (see Table 2). For clarity we show the linear fit of the maser features for each epoch (right panel), along with the length of the segments correspond to the position of the maser features used in the fitting. The reference position is α2000 = 20h38m36s​.43399 and δ2000 = +42° 37′34″​​.8710. The parameters of the linear fit are reported in Table 4. The solid light-gray line indicates the orientation of the thermal radio jet (PA = +42°±5°; Rodríguez-Kamenetzky et al. 2020).

In the text

	Fig. 3. Comparison of the total sum of the H2O maser spectra detected toward VLA 1 with the EVN in epochs 2014.46, 2016.45, 2018.44, and 2020.82. The vertical dashed dark-green line indicates the assumed systemic velocity of the region (Vlsr = +10.0 km s−1; Shepherd et al. 2003).
In the text

Fig. 4. Comparison of the H2O maser features detected toward VLA 2 (left panel) in the four EVN epochs (2014.46, 2016.45, 2018.44, and 2020.82) and superimposed to the natural-weighted continuum map at K-band (central frequency 22 GHz) of the thermal, collimated ionized wind emitted by VLA 2 obtained with the VLA by Carrasco-González et al. (2015). The light gray contours are 5, 10, 15, 20, 25, 50, and 75 times σ = 10 μJy beam−1 and the VLA beam is shown on the bottom right corner. The positions of the maser features are corrected assuming the proper motion of the subregion W75N(B) equal to the median proper motion measured for the 6.7 GHz CH3OH maser features by Rygl et al. (2012), ⟨μα⟩ = (−1.97 ± 0.10) mas yr−1 and ⟨μδ⟩ = (−4.16 ± 0.15) mas yr−1. The size of the symbols are ten times the uncertainties of the absolute positions of the maser features (see Table 2). For clarity, we show the parabolic fit of the maser features of zone 1 and the linear fit of the maser features of zone 2 (north features) and zone 4 for each epoch (right panel). The reference position is α2000 = 20h38m36s​.43399 and δ2000 = +42° 37′34.​​″8710. The parameters of the linear and parabolic fits are reported in Tables 4 and 5, respectively.

In the text

	Fig. 5. Visual comparison of the polynomial fit areas of zone 1 of VLA2. Each area represents the area in each epoch where the maser features considered in the polynomial fit (see Fig. 4 and Table 5) are located. The colors correspond to those given in Fig. 4.
In the text

	Fig. 6. Comparison of the total sum of the H2O maser intensity (I) measured toward VLA 1 and VLA 2 with the VLBI in epochs 2005.89 (VLBA, S11a), 2012.54 (VLBA, S14), 2014.46 (EVN), 2016.45 (EVN), 2018.44 (EVN), and 2020.82 (EVN). The lowest intensities coincide with the minimum of activity registered toward W75N(B) by the 22-m Pushchino telescope in May-July 2012 (Krasnov et al. 2015).
In the text

	Fig. 7. Comparison of the total sum of the H2O maser intensity (I) for five different ranges of velocities in VLA 1 as measured in epochs 2012.54 (VLBA, S14), 2014.46 (EVN), 2016.45 (EVN), 2018.44 (EVN), and 2020.82 (EVN).
In the text

Fig. 8. Comparison of the H2O maser features detected for five different ranges of velocities in VLA 1 in the four EVN epochs 2014.46, 2016.45, 2018.44, and 2020.82. The size of the symbols are scaled logarithmically according to their peak intensities. The epoch for which the maximum of the intensities is registered for each velocities range is colored in light gray. The assumed velocity of the region is Vlsr = +10.0 km s−1 (Shepherd et al. 2003).

In the text

	Fig. 9. Comparison of the total sum of the H2O maser spectra detected toward VLA 2 with the EVN in epochs 2014.46, 2016.45, 2018.44, and 2020.82. The vertical dashed dark-green line indicates the assumed systemic velocity of the region (Vlsr = +10.0 km s−1; Shepherd et al. 2003).
In the text

	Fig. 10. Comparison of the total H2O maser intensity in the four zones of VLA 2 (see Fig. 4) measured in epochs 2012.54 (VLBA; S14), 2014.46 (EVN), 2016.45 (EVN), 2018.44 (EVN), and 2020.82 (EVN).
In the text

Fig. 11. Multi-epoch comparison of ⟨ΦB⟩ angles for VLA 1 (left panel). The epochs are 2005.89 (VLBA, S11a), 2012.54 (VLBA, S14), 2014.46 (EVN), 2016.45 (EVN), 2018.44 (EVN), and 2020.82 (EVN). The two horizontal solid black lines indicate the ±90° angles. The dashed gray lines represents the best linear fit results of the data considering that the position angle of a magnetic field vector has contemporary multiple values equal to ⟨ΦB⟩=⟨ΦB⟩±180°. The slopes of these lines are mVLA 1 = ( − 15° ​.38 ± 1° ​.80) yr−1. Magnetic field vectors and strength along the line of sight (next to the corresponding maser feature) as estimated from all the linearly and circularly polarized maser features detected around VLA 1 (right panel) during the four EVN epochs (2014.46, 2016.45, 2018.44, and 2020.82). The vectors are superimposed to the uniform-weighted continuum map at Q-band (central frequency 44 GHz) of the thermal jet driven by VLA 1 (Rodríguez-Kamenetzky et al. 2020). For more details, see the caption for Fig. 2.

In the text

Fig. 12. Multi-epoch comparison of ⟨ΦB⟩ angles for VLA 2 in the left panel. The epochs are 2005.89 (VLBA, S11a), 2012.54 (VLBA, S14), 2014.46 (EVN), 2016.45 (EVN), 2018.44 (EVN), and 2020.82 (EVN). The gray point indicate the ⟨ΦB⟩ angle measured in 2020.82 after a rotation of 180°, considering that ⟨ΦB⟩=⟨ΦB⟩±180°. The two horizontal solid black lines indicates the ±90° angles. The dashed gray lines represents the best linear fit results of the data considering that the position angle of a magnetic field vector has contemporary multiple values equal to ⟨ΦB⟩=⟨ΦB⟩±180°. The slopes of these lines are mVLA 2 = ( − 19° ​.73 ± 1° ​.89) yr−1. Right panel: magnetic field vectors and strength along the line of sight (next to the corresponding maser feature) as estimated from the linearly and circularly polarized maser features detected around VLA 2 during the four EVN epochs (2014.46, 2016.45, 2018.44, and 2020.82) and superimposed to the natural-weighted continuum maps at K-band (central frequency 22 GHz) of the thermal, collimated ionized wind emitted by VLA 2 (Carrasco-González et al. 2015). For more details see caption of Fig. 4. The arrows represent the direction of the proper motions on the plane of the sky reported in Tables 4–6, and the arrows length is proportional to their magnitudes. The double arrow in the northeast indicates the uncertainties of the motion direction, due to the presence of an inward motion (case A) or of an outward motion (case B) of the gas, as explained in Sect. 3.1.2.

In the text

Fig. A.1. Close-up view of the 22 GHz H2O maser features detected around the radio sources VLA 1 (left panels) and VLA 2 (right panels) in epochs 2014.46 (top panels), 2016.45, 2018.44, and 2020.82 (bottom panels). The reference positions are reported in Table 2. The letters A, B, and C indicate the position of the three maser groups identified around VLA 1 by S11a and they are reported here only for reference. The octagonal symbols are the identified maser features scaled logarithmically according to their peak intensity (Tables A.1-A.8). Maser LSR radial velocities are indicated by color (the assumed systemic velocity of the region is Vlsr = +10.0 km s−1; Shepherd et al. 2003). The linear polarization vectors, scaled logarithmically according to the polarization fraction Pl, are overplotted. In the bottom-right corner of all panels the error-weighted orientation of the magnetic field (⟨ΦB⟩) is also reported; the two dashed segments indicate the uncertainties.

In the text

Fig. A.2. Total intensity spectra (I, upper panel) and circular polarization intensity spectra (V, lower panel) for the H2O masers VLA1.2.01, VLA2.2.17, VLA2.2.27, VLA1.3.04, VLA1.3.16, VLA2.3.14, VLA2.3.39, and VLA2.4.27 (see Tables A.3, A.4, A.5, A.6, and A.8). The thick red lines are the best-fit models of I and V emission obtained using the outputs of the FRTM code for 22 GHz H2O masers (see Sect. 2), in particular these models are used to fit the spectra (for more details see S11a). The maser features were centered on zero velocity.

In the text

Fig. B.1. Comparison of the positions of the continuum emission at K-band of VLA 1 (left panel), VLA 2 (central panel), and VLA 3 (right panel) as observed in epoch 1996.96 (solid black; Torrelles et al. 1997) and epoch 2014.20 (dashed red; Carrasco-González et al. 2015; Rodríguez-Kamenetzky et al. 2020). The squares represent the uncertainties of the peaks positions (see Table B.1) and the circles are the restoring circular beam of 0.​​″12 used to properly compare the two epochs. Each panel is centered at the corresponding position of the YSO in epoch 1996.96 (see Table B.1). The positions of epoch 2014.20 are corrected for the proper motion of the region W75N(B) as measured by Rygl et al. (2012), ⟨μα⟩ = (−1.97 ± 0.10) mas yr−1 and ⟨μδ⟩ = (−4.16 ± 0.15) mas yr−1.

In the text