© The Authors 2024

1 Introduction

Outflows in the form of either slow and poorly collimated winds or fast jets are commonly observed toward forming stars, which suggests the existence of a tight connection between mass accretion and ejection in the process of star formation. The ability of protostellar winds to extract angular momentum from the circumstellar disk and enable mass accretion, driving disk evolution and dispersal and ultimately setting the condition for planet formation, critically depends on their launching mechanism. Several mechanisms have been proposed to link accretion and ejection, and to account for the high degree of jet collimation, the magnetic field plays an important role in most of them. The proposed mechanisms include magneto-centrifugal acceleration either at the interface of the stellar magnetosphere with the disk (X-winds, Shu et al. 1994) or from the centrifugally supported part of the disk (disk wind, DW; Blandford & Payne 1982), magnetic pressure alone (Lynden-Bell 2003), and thermal pressure following the disk photo-evaporation by intense stellar UV radiation (Hollenbach et al. 2000). Recent results on low-mass (~ 1 M_⊙) star formation favor models of magnetohydrodynamic (MHD) DWs since they can supersede magnetic turbulence and extract disk angular momentum in regions of low ionization and account for the emission of both winds and jets (Pascucci et al. 2023). Protostellar jets have also been investigated in intermediate (a few solar masses) and high-mass (≥10 M_⊙) pro-tostars through high-angular resolution observations of the radio continuum and the measurement of proper motions of the radio knots (Carrasco-González et al. 2021; Rodríguez-Kamenetzky et al. 2022), with indications being found in favor of X-winds and DWs or jet collimation resulting from ambient gas pressure.

In the past ten years, high-angular resolution (≈20–50 mas) Atacama Large Millimeter Array (ALMA) observations have identified several low- and high-mass young stellar objects (YSOs) where an LSR velocity (V_LSR) gradient transversal to the jet axis has clearly been detected. The latter is generally interpreted as jet rotation due to the magneto-centrifugal launching, and through various methods (Tabone et al. 2020), the launch radius of the wind has been derived. The launch radius, R_K, is the critical parameter to distinguish between the quasi-polar X-winds (Shu et al. 1994), R_K ≤ 0.1 au, from the radially extended DWs (Pudritz et al. 2007), R_K ~1–50 au, and this distinction is essential to understanding the wind’s role in driving mass accretion and disk evolution. By targeting thermal lines, these ALMA observations reach linear resolutions of > 10–50 au at best and can measure only line-of-sight velocities. Owing to the insufficient linear resolution and confusion between axial and rotation velocities, the ALMA measurements of the launch radii suffer from severe systematic errors and uncertainties. An emblematic case is that of the low-mass protostar HH212 (Lee et al. 2017), which is probably the clearest example of jet rotation, where the observers estimated a launch radius for the SO wind of 1 au and an ad-hoc model reproducing the ALMA observations indicated a much larger value of 40 au (Tabone et al. 2017, 2020).

By achieving linear resolutions of ~ 1 au and measuring 3D velocities, very long baseline interferometry (VLBI) observations of the 22 GHz water masers can directly trace the velocity field of protostellar winds. Our previous studies (Moscadelli et al. 2007, 2011; Sanna et al. 2010b; Goddi et al. 2011) have demonstrated the ability of 22 GHz water masers in tracing the velocity structure of the outflows near YSOs. Recently, we carried out the Protostellar Outflows at the EarliesT Stages (POETS) survey (Moscadelli et al. 2016; Sanna et al. 2018) to image the disk-outflow interface on scales of 10–100 au in a statistically significant sample (37) of luminous YSOs. We employed multi-frequency Jansky Very Large Array (JVLA) observations to determine the spatial structure of the ionized emission and multi-epoch Very Long Baseline Array (VLBA) observations (from the BeSSeL¹ survey) to derive the 3D velocity distribution of the 22 GHz water masers. Among the main results of POETS, we find that the water masers are always associated with weak, slightly resolved continuum emission, tracing the thermal jet emitted from the YSO, and the 3D velocity distribution of the water masers in most of the sources of the sample can be interpreted in terms of an MHD DW.

IRAS 21078+5211 is the most promising MHD DW candidate from the POETS survey. The star-forming region has a bolometric luminosity of 5 × 10³ L_⊙ (Moscadelli et al. 2016) at a distance of 1.63±0.05 kpc (Xu et al. 2013) and harbors a cluster of forming massive stars. Around the most massive YSO (5.6 + 2 M_⊙) of the cluster, on scales of a few 100 au, a rotating disk has been observed in high-density molecular tracers by employing the NOrthern Extended Millimeter Array (NOEMA; Moscadelli et al. 2021). The POETS JVLA 5 cm continuum observations have revealed a thermal jet (Moscadelli et al. 2016) directed toward the NE-SW (PA ≈ 44°) and emerging from the YSO, whose position at the center of the disk is pinpointed by compact JVLA continuum at 1.3 cm. The BeSSeL VLBA observations of the 22 GHz water masers, performed in 2010–2011, have discovered a cluster of masers placed ≈ 100 au NE from the YSO and whose proper motions are collimated NE-SW (PA = 49°) and trace the base of the jet from the YSO. In October 2020, we re-observed the water masers in IRAS 21078+5211 with sensitive Global VLBI observations and found that the maser spatial distribution has significantly changed with respect to the previous BeSSeL epochs (Moscadelli et al. 2022). Remarkably, close to the YSO, the water masers trace sinusoids in the plane of the sky, which can be interpreted as the sky-projection of 3D helixes seen almost edge-on. The water masers could then be tracing the streamlines of an MHD DW, where the gas is predicted to flow along helical magnetic field lines anchored to the rotating disk. Indeed, Moscadelli et al. (2023) have shown that resistive-radiative-gravito-MHD simulations of a DW can model the maser streamlines with great accuracy (rms deviation of ~1 au; see Fig. 1).

This article reports on new multi-epoch VLBA observations of the 22 GHz water masers in IRAS 21078+5211 carried out from March to October 2023. These new observations allow us to determine the 3D velocities of the water masers in the DW streamlines, which are key to accurately measuring the launch radii of the streamlines and reliably identifying the launching mechanism of the protostellar wind. Sect. 2 describes the new VLBA observations, and the distribution of the water maser 3D velocities is presented in Sect. 3. The discussion of the new results is provided in Sect. 4 and the conclusions are drawn in Sect. 5.

Fig. 1 Three-dimensional view of the MHD DW in IRAS 21078+5211 from simulations in Moscadelli et al. (2023). The observed maser positions (red and orange dots) are overlaid on streamlines (blue lines) computed from resistive-radiative-gravito-MHD simulations of a jet around a forming massive star.

2 Very Long Baseline Array observations

We observed the 6_ı6 − 5₂₃ H₂O maser transition (rest frequency 22.235079 GHz) toward IRAS 21078+5211 (tracking center: RA(J2000) = 21^h9^m21.720 and Dec(J2000) = +52°22′37″.08) with the ten antennae of the VLBA of the National Radio Astronomy Observatory (NRAO²) at four epochs: March 6, June 5, August 5, October 19, 2023. Each epoch lasted 14 hr. Because of technical problems, the antenna of Mauna Kea, Owens Valley, and Saint Croix respectively did not take part in the observations in the first, second, and fourth epochs. To determine the absolute positions and velocities, we performed phase-referencing observations, alternating scans on the target and the phase-reference calibrators every 45 s. The phase-reference calibrators were the quasars 2116+543 and 2051+528, both within 2.5° from the target and with a correlated flux of ~0.1 Jy beam⁻¹ at 22 GHz. Calibration scans of ≈8 min were observed every hour on the fringe-finder and bandpass calibrators J2202+4216, 2007+777, 3C84, and 3C48.

We recorded dual circular polarization through four adjacent bandwidths of 16 MHz, one of them centered at the YSO’s V_LSR of −6.4 km s⁻¹. The four 16 MHz bandwidths were used to increase the signal to noise ratio of the weak continuum (phase-reference) calibrators. The data were correlated with the VLBA-DiFX software correlator at Socorro (New Mexico, USA) in two correlation passes using 1024 and 128 spectral channels to correlate the maser 16 MHz bandwidth and the whole set of four 16 MHz bandwidths, respectively. The spectral resolution attained across the maser 16 MHz band was 0.21 km s⁻¹. The correlator averaging time was 1 s.

Data were reduced with the Astronomical Image Processing System (AIPS, Greisen 2003) package following the VLBI spectral line procedures in the AIPS COOKBOOK³. At each epoch, the emission of an intense and compact maser channel was self-calibrated, and the derived (amplitude and phase) corrections were applied to all the maser channels before imaging. To cover the whole maser emission, we produced images extending 1”.23 in both RA cos δ and Dec, and 84 km s⁻¹ in V_LSR. Reflecting the changes in the antennae taking part in the observations, using natural weighting, the FWHM major and minor sizes of the beam vary in the ranges 0.6–1.2 mas and 0.4–0.8 mas, respectively, and the beam PA in the range −3°–75°. In channel maps with a (relatively) weak signal, the 1σ rms noise is ≈4 mJy beam⁻¹, which is close to the expected thermal noise.

For a description of the criteria used to identify individual maser features, derive their parameters (position, intensity, flux, and size), and measure their (relative and absolute) proper motions, we refer to Sanna et al. (2010a). Individual maser features are a collection of quasi-compact spots observed on contiguous channel maps and spatially overlapping (within their FWHM size). The spot positions are determined by fitting a two-dimensional elliptical Gaussian to their spatial emissions. The uncertainty of the spot position relative to the reference maser channel is the contribution of two terms: $\text{Δ}{\theta }_{\text{spot }}=\sqrt{\text{Δ}{\theta }_{\text{fit }}^2+\text{Δ}{\theta }_{\text{bandpass }}^2}$. The first term depends on the S/N of the data, following Reid et al. (1988): Δθ_fit = θ_beam/(2 S/N), where θ_beam is the resolution beam size, conservatively taken equal to the maximum value of the FWHM major beam size of 1.2 mas. The second term depends on the accuracy of the bandpass calibration through the expression (Zhang et al. 2017): Δθ_bandpass = θ_beam(ΔΨ/360°), where ΔΨ in degrees is the phase stability across the observing band. In our case, ΔΨ ≲ 5° and Δθ_bandpass ≲ 0.02 mas becomes the dominant error term for spot intensity ≥200 mJy beam⁻¹. The maser feature position (and corresponding error) is estimated from the error-weighted mean of the spot positions (and corresponding errors), and the feature V_LSR from the intensity-weighted mean of the spots’ V_LSR. The error on the line-of-sight velocities, defined as the difference between the maser V_LSR and the YSO’s V_LSR, is taken equal to the uncertainty of 0.4 km s⁻¹ on the YSO’s V_LSR, exceeding the maser V_LSR error of ≲ 0.2 km s⁻¹.

Table 1 reports the parameters (intensity, V_LSR, position, and relative proper motion) of the 22 GHz water masers in IRAS 21078+5211. As described in Sanna et al. (2010a), at least three observing epochs are necessary to assess the time persistence of a maser feature and accurately derive its proper motion. Relative proper motions are calculated with respect to the geometric center (hereafter “center of motion”, identified with label #0 in Table 1) of the features with a stable spatial and spectral structure, persisting over the four observing epochs. Errors of the relative proper motions are typically of ≈0.6 km s⁻¹. At each epoch, inverse phase-referencing (Reid et al. 2009) produced good S/N (≥10) images of the two phase-reference calibrators. Taking into account that the calibrators are relatively compact, with size ≲1 mas, and that the absolute position of the calibrators is known within a few 0.1 mas, we estimate that the error on the absolute position of the masers is ≲0.5 mas. Without correction for the uncompensated delay introduced by the Earth’s atmosphere (Reid et al. 2009), our VLBA observations reach an astrometric accuracy of ≲0.1 mas, and the absolute proper motions have typical errors of 2–3 km s⁻¹. The derived absolute proper motions have been corrected for the apparent proper motion due to the Earth’s orbit around the Sun (parallax), the solar motion, and the differential Galactic rotation between our LSR and that of the maser source. We adopted a flat Galaxy rotation curve (R₀ = 8.33 ± 0.16 kpc, Θ₀ = 243 ± 6 km s⁻¹ (Reid et al. 2014)) and the solar motion $U={11.1}_{-0.75}^{+0.69},V={12.24}_{-0.47}^{+0.47},$ and $W={7.25}_{-0.36}^{+0.37}$ km s⁻¹ by Schönrich et al. (2010), who revised the Hipparcos satellite results.

3 Results

Figure 2 illustrates the reference frames employed to describe the maser positions and velocities, and magnetic fields. The jet axis z is inclined 15° with the plane of the sky, and its projection on the sky xl is directed at PA = 50° (see Appendix A). The axis R is orthogonal to xl in the plane of the sky, and the axis oɀ is orthogonal to ɀ in the plane containing ɀ and the line-of-sight los. The maser positions in the reference frame (R, oɀ, ɀ) can also be described using cylindrical coordinates: the inclination angle θ between the position vector r and the axis ɀ; ρ, the projection of r into the plane (R, oɀ), oriented at the azimuthal angle ϕ with respect to the axis oɀ. We note that ϕ increases clockwise in agreement with the disk rotation. Finally, velocity and magnetic field vectors can be decomposed in the reference frame co-moving with the masers along the axial (ê_ɀ), radial (ê_ρ), and azimuthal (ê_ϕ) directions.

Figure 3 shows the maser positions and proper motions from the new 2023 VLBA observations. We note that the observing parameters of the 2023 VLBA observations, including beam size, spectral resolution, and sensitivity, are comparable to those of the previous October 2020 Global VLBI observations. The new observations continue to show a bipolar distribution of masers, both spatially and kinematically (in terms of proper motions and radial velocities), consistent with the findings of Moscadelli et al. (2022) from the 2020 observations. Maser emission is concentrated in the same three regions as in the previous Global VLBI observations, namely, to the NE, SW, and N of the YSO. However, inside each of these regions, the spatial distribution of the masers has changed. Figure 4 compares the positions of the water masers in the NE and SW regions between the 2020 Global VLBI and 2023 VLBA observations. In the NE region at 20 au ≤ xl ≤ 90 au, the masers from 2020 traced a continuous curved pattern, while two separated streams at different distances from the jet axis are visible in the observations of 2023. In the SW region, the 2020 masers drew a single arc extending over −100 au ≤ xl ≤ −20 au; in 2023 there are hints of multiple arcs spreading over a larger distance of −170 au ≤ xl ≤ −20 au.

As in the 2020 observations, the maser patterns of the new 2023 VLBA observations still present sinusoidal shapes, which can be interpreted in terms of sky-projection of 3D helical streamlines. Figure 5 shows a detailed view of these streamlines and the maser velocities along them. In the 2023 observations, inside the NE region over the xl range 20–90 au, the masers trace two sinusoids of different amplitudes, named B1 and B2. In the SW region, only a single sinusoid, named R1, can be reliably identified with the masers. In Appendix B we discuss the procedure to fit the sinusoids to the maser positions, and the parameters of the sinusoidal fits (ℜ, the amplitude of the sinusoid; f_ɀ, the spatial frequency; and xl₀, the position of zero phase) are reported in Table 2. A comparison of this table with Table 1 in Moscadelli et al. (2022) confirms that the DW streamlines in IRAS 21078+5211 have changed over a timescale ≲1 yr.

As discussed in Appendix C, the proper motions presented in Fig. 3 are calculated relative to the YSO. The uncertainty is of a few kilometers per second and velocity amplitudes range from a few kilometers per second up to ≈60 km s⁻¹, with the fastest motions found close to the jet axis at high xl ≥120 au. Some of the masers with the highest velocities belong to the B3 stream, an elongated distribution of masers that extend from the NE region northward. We note that the spatial distribution of the masers of the B3 stream inside the NE region (over 105 au ≤ xl ≤ 150 au; see Fig. 5, left panel) is not actually linear but slightly arched. However, as discussed in Appendix B, the positional scatter of these masers is too large to permit a reliable determination of a streamline. The masers in the N region trace an extended (≈50 au) arc (see Fig. 3), with proper motions directed about perpendicular to the arc, suggesting that it is expanding away.

Fig. 2 Reference frames used for maser positions and velocities, and magnetic fields. (a) Transformation of coordinates by rotation around the line-of-sight axis los over the position angle of the jet to bring Dec onto the sky-projected jet axis xl. (b) Transformation of coordinates by rotation around the axis R over the inclination angle of the jet to bring xl onto the jet axis ɀ. (c) Transformation of coordinates by rotation around the axis ɀ over the rotation angle ϕ of the maser position to bring oɀ along the direction êρ. (d) Decomposition of the maser velocity into the axial (vɀ), radial (vρ), and azimuthal (vϕ) components.

4 Discussion

4.1 The magneto-centrifugal disk wind

Figure 6 presents a sketch of an MHD DW streamline where the gas flows along a helical magnetic field line anchored to the rotating disk. According to the MHD DW theory, the gas is magneto-centrifugally launched along a quasi-poloidal magnetic field line departing from the disk. The field line behaves like a rigid flow tube co-rotating with the disk until reaching the Alfvén point where the increasing kinetic energy of the outflowing gas equals the magnetic energy. Beyond that point, the inertia of the flow becomes important, the magnetic field line starts lagging behind the disk rotation, and the azimuthal component of the magnetic field grows in the rotation verse opposite to the disk. Thus, the geometry of the magnetic field line approaches that of a helix winding around the disk axis in counter-rotation with the disk.

The shape of the helical trajectories of the masers is characterized by the magnetic field helix angle α_B = arctan(B_ϕ/B_ɀ), which is the angle with which the helical magnetic field winds around the jet axis (see Fig. 6). The terms B_ɀ and B_ϕ are the axial and azimuthal components of the magnetic field, respectively. In particular, by performing a sinusoidal fit of the sky-projected maser streamlines, one can derive ℜ, the amplitude of the sinusoid, corresponding to the radius of the helix and f_ɀ, the spatial frequency of the sinusoid (see Sect. 3 and Table 2). Following Eq. (12) in Moscadelli et al. (2022), ‖B_ϕ‖/‖B_ɀ|| = ‖ tan(α_B)‖ = f_ɀℜ, that is, the product f_ɀℜ defines the shape of the helix. In the above equations, we have accounted for the negative sign of B_ϕ, as it is directed opposite to the verse of increasing ϕ.

An observer co-rotating with the launch point of the streamline (at the angular velocity ω_K) sees the gas moving along the helical field line. That is, in such a reference frame, the flow velocity is parallel to the magnetic field. Therefore, we can write $$\frac{\Vert v_{\varphi }^{\prime }\Vert }{\Vert v_z\Vert }=\frac{\Vert B_{\varphi }\Vert }{\Vert B_Z\Vert }=f_z\text{\hspace{0.17em}ℜ},$$(1)

where we use the apical prime to indicate that the azimuthal velocity is measured in the co-rotating reference frame. On the other hand, the maser motion we observe is the composition of the motion along the streamline and the rotation of the streamline anchored to the disk. While the axial component of the velocity is not affected by such a rotation, the azimuthal velocity component we measure is given by $$v_{\varphi }=v_{\varphi }^{\prime }+{\omega }_{\text{K}}\text{ℜ}.$$(2)

Combining Eqs. (1) and (2), we derive $${\omega }_{\text{K}}=\frac{v_{\varphi }}{\text{ℜ}}+f_z\Vert v_z\Vert ,$$(3)

which shows that the knowledge of the streamline shape and maser 3D velocity allows us a direct determination of the angular velocity of the launch radius. Our previous Global VLBI observations (Moscadelli et al. 2022) of IRAS 21078+5211 measured just the maser line-of-sight velocity, and we could only set lower limits for the axial velocities of the masers and the angular velocities of the launch radii of the streamlines and only derive upper limits for the launch radii.

Knowing the jet PA = 50° and its inclination i_jet = 15° with the plane of the sky (see Appendix A), we can convert the velocity components from the reference system (RA, Dec, los) in which they are originally measured to the reference system (R, oɀ, ɀ; see Figs. 2a and 2b). The rotation (or azimuthal) angle, ϕ, of the masers belonging to a streamline is calculated with the formula: ϕ = arcsin(R/ℜ), where R is the sky-projected radius of the maser and R is the rotation radius of the streamline. By fitting sinusoids to the sky-projection of helical motions and varying the inclination angle of the helix axis, we verified that this formula is accurate to within a few degrees for small inclinations ≤20°. As the streamline winds counterclockwise around the disk axis, the rotation angle goes from 180° to 0° for increasing xl along the two NE streamlines, and from 360° to 180° for decreasing xl along the SW streamline. The knowledge of the rotation angle of a given maser allows us to determine the radial, v_ρ, and azimuthal, v_ϕ, components of its velocity (see Figs. 2c and 2d). Table 3 reports the mean and standard deviation of the axial, radial, and azimuthal velocity components for each streamline. While v_ɀ dominates over v_ρ and V_ϕ in the B1 and B2 streamlines, v_ρ and v_ɀ are comparable in the R1 streamline.

Figure 7a shows the distribution of the Keplerian angular velocities computed from all masers using Eq. (3). In agreement with our expectation, the large majority (75–80%) of the masers belonging to the same streamline provide consistent values of the angular velocity. We emphasize that the found consistency provides a direct test of the MHD DW model employed to interpret the data. A low fraction of outliers⁴ can be expected considering both the high time variability and the possibility of spatial overlap of the streamlines (see Fig. 4 and Appendix B). Our method of filtering out the outliers is mainly based on the scattering of the maser positions from the sinusoidal fit and can be ineffective at the loci of intersection of two different streamlines.

Table 2 lists the value (and corresponding error) of the angular velocity of the launch radius for each streamline, calculated by taking the mean (and standard deviation) of the consistent values of angular velocities in Fig. 7a. By reporting the derived angular velocities on the Keplerian profile of the YSO, Figure 7b shows that the launch radii of the three streamlines can be determined with high accuracy: 10.8±0.8 au, 19.2±2.8 au, and 49±5 au for B1, B2, and R1, respectively. If the uncertainty of 2 M_⊙ on the YSO mass is taken into account, the values of the launch radii change by ≈12%. We stress that our findings are crucial to favor the DW model over the X-wind model, as the latter predicts launch radii <0.1 au. Besides, the high velocities (up to 60 km s⁻¹) observed close to the jet axis (see Fig. 3) are inconsistent with acceleration due to (only) magnetic or thermal pressure, which produce typical terminal flow speeds <20 km s⁻¹ (Oliva & Kuiper 2023b). In conclusion, magneto-centrifugal launching in a radially extended MHD DW appears to be the only viable process to explain the maser 3D velocity field in this source.

Fig. 3 Water maser 3D velocities in IRAS 21078+5211. Colored dots and arrows give water maser positions (relative to the reference feature # 0) and proper motions (relative to the YSO), respectively, with colors denoting the maser VLSR. We show only the masers for which the position at the first observing epoch can be derived. The amplitude scale for the velocity is reported in the bottom-left corner of the panel. The black dotted rectangles encompass the three regions to the N, NE, and SW of the YSO, where maser emission concentrates. The sky-projected jet and disk axes are in red and black, respectively, with distances labeled in astronomical units. An elongated maser stream directed approximately S-N is traced with a dashed-line segment labeled “B3.”

Fig. 4 Time variability of the maser streamlines. Comparison of the positions of water masers observed in October 2020 (green triangles) and March-October 2023 (red dots) for the NE (left) and SW (right) regions.

Fig. 5 Three-dimensional velocities along the maser streamlines. The colored dots and arrows and the red axes have the same meaning as in Fig. 3. (Left panel) Expanding view of maser positions and proper motions in the NE region. The green and blue dashed curves are the fitted sinusoids to the B1 and B2 streamlines, respectively, whose parameters are reported in Table 2. (Right panel) Expanding view of maser positions and proper motions in the SW region. The red dashed curve is the fitted sinusoid to the R1 streamline, whose parameters are reported in Table 2.

Fig. 6 Co-rotating helical magnetic field line. The water masers (red balls) move along a helical magnetic field line (blue) anchored to the rotating disk. The rotation radius and angle, ℜ and ϕ, and the axial, Bɀ and vɀ, and azimuthal, Bϕ and v′ϕ, components of the magnetic field and velocity in the reference frame co-rotating with the disk are labeled. The tangent to the helix is inclined by an angle αB with respect to the ɀ axis.

Fig. 7 Determination of the launch radii of the maser streamlines. (a) The streamline to which each panel refers is indicated in the top-left corner of the panel. Maser labels (see Table 1) and error bars indicate the values (and corresponding errors) of the angular velocity of the launch radius inferred with the masers belonging to the streamline, which are plotted versus the position along the sky-projected jet axis xl. The dashed blue line and error bar give the mean and standard deviation of the consistent values, excluding the outliers denoted with red labels. (b) Vertical and horizontal error bars plot the values (and associated errors) of the angular velocity and rotation radius of the streamlines. The heavy and thin curves are the Keplerian profiles evaluated at three different YSO masses, as labeled. The launch radius of each streamline for the nominal value of the YSO mass of 5.6 M⊙ is identified with horizontal and vertical dotted lines, and the associated measurement error is marked by the horizontal error bar in the bottom.

4.2 Magneto-centrifugal launching versus magnetic pressure

Figure 8 compares the velocities of the masers in the B3 stream and N region. In the B3 stream, for increasing z, the velocities of the masers grow and collimate toward the jet axis, from v_ɀ ≤10 km s⁻¹ and θ_v ≥ 60° at ɀ ≈ 100 au to v_ɀ ≈60 km s⁻¹ and θ_V ≈10° at ɀ ≥ 200 au. The large change of ≈100° in the azimuthal angle ϕ_V can be interpreted in terms of rotation of the velocity vector around the jet axis. The acceleration to large (≳100 km s⁻¹) speeds together with the collimation of the velocities close to the jet axis represent a distinct characteristic of the magneto-centrifugal launching of an MHD DW (Pudritz et al. 2007). Therefore, the analysis of the maser kinematics in the B3 stream and in the B1, B2, and R1 streamlines leads us to conclude that all the outflowing gas close (R ≤ 50 au) to the jet axis is magneto-centrifugally launched by an MHD DW. As previously investigated in simulations (Oliva & Kuiper 2023b), the acceleration of the flow occurs when the magnetic field guides the flow outward, where the gravitational energy is less negative. The flow then becomes super-Keplerian and runs away: the kinetic energy keeps increasing along the magnetic field line, carrying angular momentum outward in the process. The kinematics of the masers in the N region, at R ≥ 90 au, is remarkably different from that in the B3 stream. In the N region, masers are found at relatively low heights, ɀ ≤60 au, and move at small speeds, ≤20 km s⁻¹, along directions about parallel to the disk (70° ≲ θ_V ≲ 110°). The arc-like spatial distribution (see Fig. 3) together with the small spread in the azimuthal angle of the velocities, ϕ_V ≈ 310°±10° (see Fig. 8), suggest that these masers trace a radially expanding narrow shock-front placed at relatively large radii, 90 ≤ R ≤ 130 au (see Fig. 3), and extending from the disk plane up to ɀ ≈ 60 au. This kinematic pattern is clearly inconsistent with magneto-centrifugal launching, and instead it could be explained in terms of acceleration by magnetic pressure, which, according to our MHD DW model (Moscadelli et al. 2023), dominates at R ≥ 50 au.

We previously presented (Moscadelli et al. 2023) a resistivemagnetohydrodynamical simulation (including radiation transport and self-gravity) that closely (rms deviation of ~1 au) reproduces the maser spatial patterns of the 2020 Global VLBI observations (see Fig. 1). We looked for simulations that simultaneously reproduce the mass of the protostar, the disk radius, the presence of a high velocity jet, and the accretion rate observed in IRAS 21078+5211. We then explored several launch radii for the velocity streamlines (in the frame of reference co-rotating with the launch point) and compared against the maser spatial patterns. In the simulations, we previously found an Alfvén surface ≲50 au wide that becomes more conical toward the protostar because the outflow is constrained by the more dense material from the disk and infalling envelope. The arc-like maser distribution in the N region is located at R ≳ 100 au away from the protostar. Even though we cannot rule out that the magneto-centrifugal mechanism can operate at such distances (see discussion on the caveats of the simulations in our previous articles (Moscadelli et al. 2023; Oliva & Kuiper 2023a)), the condition that the Alfvén speed in the launching region be sufficiently large must be satisfied. This means that the density must be sufficiently low and the magnetic field sufficiently strong. Both conditions are less likely to be satisfied at large distances from the protostar as the infalling material is encountered.

On top of the accretion disk, the rotational dragging of the magnetic field lines creates a magnetic pressure gradient that can be stronger than gravity (see our previous discussion in Oliva & Kuiper (2023b) on magnetic tower flows). Given that the maser kinematics in the N region do not exhibit the signatures of magneto-centrifugal motion, we examined whether magnetic pressure could be the driving force behind their kinematics. Figure 9 shows a snapshot of the simulation including an analysis of forces (right panel). The label a in the figure indicates the magneto-centrifugal jet, which is launched and accelerated in the region where the flow is simultaneously sub-Alfvénic (shadowed yellow contour) and centrifugal force dominates over gravity (blue shadowed contour). The region where the magnetic pressure gradient dominates over gravity is bigger than the Alfvén surface (cf. unshadowed part of the red contour). As a consequence, regions b and c are outflowing but not magneto-centrifugally launched. Region b is an episodic phenomenon we previously described as a cavity wall ejection (Oliva & Kuiper 2023b). At this stage, the radially outward motion in b is driven by magnetic pressure, which is accumulated on top of the accretion disk (d) at a similar height (but not radial distance) as the masers in the N region. Additionally, most of the material in b is super-Alfvénic, meaning that we expect no correlation of the magnetic field direction and flow, just as observed in the polarization data from the N region (Moscadelli et al. 2023). A magnetic tower flow (c) should have an important vertical velocity component. Because the proper motions in the N region exhibit a rather weak v_ɀ, we conclude that a plausible production mechanism for the kinematical features of the masers in the N region is an episodic ejection event. In the simulations, such episodic ejection phenomena have lifetimes of 10–100 years. Future observations will provide a better understanding of the nature of the flow in the N region.

Fig. 8 Change in the maser velocity direction and amplitude moving along the jet axis. The left and right panels refer to the B3 stream and N region, respectively. The inclination, θ, and azimuthal, ϕ, angles of the direction of the velocity (see Fig. 2c) and the velocity component along the jet, vɀ, are plotted versus the position along the jet, ɀ, in the top, middle, and bottom panels, respectively. Values and corresponding errors are denoted with dots and error bars, respectively, where red marks velocities ≲6 km s−1 in the B3 stream.

Fig. 9 Snapshot of the simulation when the mass of the protostar is 5 M⊙. The left panel shows the density, velocity, and magnetic field lines in the circumstellar gas. The right panel contains contours (solid lines) accompanied by shadows (semi-transparent thick lines of the same color as the contour). The shadow indicates in which direction the inequality listed in the legend holds true. A smoothing algorithm has been applied to the contours to improve clarity. For more details, we refer to Oliva & Kuiper (2023b).

5 Conclusions

We have recently carried out the POETS survey, where the 3D velocities of the protostellar winds in a relatively large (37) sample of luminous YSOs have been mapped at linear resolutions of ~ 1 au through VLBI observations of the 22 GHz water masers. This technique overcomes the limitations of interfero-metric observations of thermal tracers, which, being significantly more extended than the maser cloudlets, are less suitable to act as kinematic probes of the small regions where protostellar winds are launched and collimated. The POETS results suggest that MHD DWs can be the common launching mechanism of protostellar winds in high-mass YSOs. In this article, we have tested the MHD DW scenario in the YSO IRAS 21078+5211, the best MHD DW candidate from the POETS survey. We have performed multi-epoch VLBA observations of the 22 GHz water masers in IRAS 21078+5211 to measure the 3D velocities of the maser cloudlets tracing the DW streamlines.

In IRAS 21078+5211, the fast (10–60 km s⁻¹) outflowing gas moving at small angles (≲30°) and separation (R ≲40 au) from the jet axis is magneto-centrifugally launched in a radially extended MHD DW. The water masers trace three DW streamlines whose launch radii are measured with an accuracy of ~1 au: 10.8±0.8 au, 19.2±2.8 au, and 49±5 au. We have directly tested the model of an MHD DW by verifying that the water masers move along individual DW streamlines co-rotating with the disk at the angular velocity of their launch radius. At larger radii (R ≳100 au), the relatively slow (≲20 km s⁻¹) flow traced with the masers propagates outward mainly along the radial direction, and its properties are clearly inconsistent with magneto-centrifugal acceleration. The comparison with our resistive-magnetohydrodynamical simulations suggests that these masers can trace an episodic ejection event driven by the magnetic pressure accumulated on top of the accretion disk. This work proves that water maser VLBI observations can determine the launching mechanism of protostellar winds in a trustworthy manner, as the model assumptions can also be tested, and it calls for applying this technique across the range of protostellar masses and ages traceable with water masers. Our findings are in line with recent results of low-mass star formation and suggest that MHD DWs could be the common launching mechanism of protostellar winds from low to high-mass protostars.

Data availability

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

Acknowledgements

A.O. acknowledges financial support from the Federal Commission for Scholarships for Foreign Students for the Swiss Government Excellence Scholarship (ESKAS No. 2023.0405) for the academic year 2023–2024, the Office of International Affairs and External Cooperation of the University of Costa Rica, and the European Research Council under the European Union’s Horizon 2020 research and innovation program (project ID 833925, project STAREX). O.B. acknowledges financial support from the Italian Ministry of University and Research – Project Proposal CIR01_00010. Scientific results from data presented in this publication are derived from the following VLBA (NRAO) project code: BM534. NRAO is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc.

Appendix A The orientation of the jet axis

In Moscadelli et al. (2022) we have fixed the PA of the jet, PA_jet, to be ≈50° on the basis of several evidences: 1) the elongation at PA = 56°±12° of the slightly resolved (size ≈ 150 au) JVLA continuum at 1.3 cm; 2) the collimated, PA = 49° ± 18°, proper motions of the water masers near (100–200 au) the YSO from the previous VLBA (BeSSeL) observations; 3) the direction at PA ≈ 44° of the non-thermal radio jet traced by the extended (size ≈ 1000 au) JVLA continuum at 5 cm. On top of this evidence, we note that a change of ≳5° in PA_jet would compromise the sinusoidal fit of the maser streamlines. In fact, if the jet axis was oriented at PA ≳55° it would intersect the streamlines, while in correspondence of a jet PA ≲45° the lower and upper extremes of the streamlines would be placed at significant (and remarkably different) distances from the axis. Therefore, we can set PA_jet = 50°±5°.

Considering the small inclination with the plane of the sky of the maser 3D velocities from previous VLBA observations, in Moscadelli et al. (2022) we have set an upper limit for the jet inclination i_jet ≤ 30°. Now, we improve on the estimate of the jet inclination by arguing that a selection of (5) masers from previous VLBA observations should accurately trace the jet direction in light of the following evidences: 1) they are placed along the jet axis at the largest distance from the YSO; 2) they move at relatively high speed, between 23 and 37 km s⁻¹; 3) the PA of their proper motions is in the range 45°– 55°; 4) the inclination angle of their 3D velocities varies in the small range 12°– 19°, with an average inclination of 15°. Therefore, we can set i_jet = 15°±5°, the jet axis being tilted toward us at PA = PA_jet.

Appendix B Streamlines

In Figs. 3 and 5, the plotted maser positions correspond to the first observing epoch, since we plot either masers detected at the first epoch or those observed at least in two successive epochs for which the first epoch position can be recovered through their proper motion. Instead, Fig. B.1 shows all the masers detected across the four epochs. The comparison between Fig. B.1 and the left panel of Fig. 5 shows the importance of using simultaneous masers when tracing the DW streamlines. The positions of the masers observed only at either the second, third or fourth epoch fall in between the B1 and B2 streams, thus making it hard to recognize the two streamlines.

Figure 4 shows that, in both the NE and SW regions, the older (2020) and newer (2023) maser patterns are spatially contiguous, intersect and partly overlap. As the maser streamlines change on small timescales ≲1 yr comparable or lower than the maser life time, the spatial and time overlap of different streamlines can make it difficult to discern the masers belonging to a given streamline. In the following, we briefly describe the procedure employed to fit the three sinusoids to the maser positions tracing the R1, B1 and B2 streamlines. Our main criterion to select the masers is the quality of the sinusoidal fit. Since the parameters to fit are three, a streamline can be identified only if it is adequately sampled by a minimum of four maser clusters.

As shown in the left panel of Fig. B.2, the R1 streamline is traced by four, spatially distinct maser clusters. The sinusoidal fit of the maser positions has a very low rms deviation of 0.35 au. Besides, the sinusoidal fit of the positions of the three maser clusters closer to the YSO (over −20 au ≥ xl ≥ −120 au) intercepts the fourth, more distant cluster (at xl ≈ −150 au), too. Any other combination of masers of the SW region results in a sinusoidal fit of much lower quality. For instance, replacing the first two clusters of the R1 streamline (over −20 au ≥ xl ≥ −60 au) with those nearby but at slightly smaller radius (see Fig. 5, right panel), degrades the rms deviation of the sinusoidal fit to 1 au.

Fig. B.1 Non-simultaneous water masers. Colored dots and the red axis have the same meaning as in Fig. 3. At variance with the left panel of Fig. 5, all the detected masers are shown, including those observed only in a single epoch different from the first epoch. For these masers, missing the proper motion information, the first-epoch position cannot be derived.

In the case of the B1 stream (see Fig. B.2, right panel), after excluding four masers clearly separated in V_LSR (the yellow-green dots and arrows in Fig. 5, left panel), the remaining masers accurately sample (at more than 10 distinct positions) the sinusoid, whose fit has a low rms deviation of 0.65 au. Instead, in the case of the B2 stream, by including all the masers along the stream, the resulting sinusoidal fit has a significantly higher rms of 1.4 au. Differently from the R1 stream, the number of masers in the B2 stream is too high to select unambiguously the positions that provide the best sinusoidal fit. Therefore, we have considered the velocities to see whether it was possible to distinguish masers belonging to different streamlines on the basis of their motion.

Figure B.3 shows that a set of (8) masers (red labels) of the B2 stream have almost the same velocities in comparison with the rest of the masers presenting a much higher velocity scatter. Figure B.4 shows the proper motions (red arrows) of these eight masers, which are about parallel to each other. Assuming, conservatively, that this set of masers belongs to a different streamline, we have excluded them and repeated the sinusoidal fit of the B2 streamline. The result is presented in Fig. B.4, which shows that a few (blue labels) other masers can be also excluded on the basis of the very scattered positions. The sinusoidal fit to the remaining masers (black labels), defining the B2 streamline, has a rms deviation of 0.9 au.

Finally, we have also tried to fit a sinusoid to the masers of the B3 stream inside the NE region (see Fig.5, left panel) whose spatial distribution is roughly arc-like. After excluding a few, most scattered masers, the rms deviation of the sinusoidal fit is 2 au, 2–6 times larger than that of the B2, B1 and R1 streamlines. Besides, the analysis of the maser velocities does not provide any clue to separate the masers into subsets belonging to different streamlines. Our conclusion is that we cannot employ these masers to identify a streamline with a confidence level as good as that of the previous three streamlines.

Fig. B.2 Sinusoidal fits of the R1 and B1 streamlines. Maser labels (see Table 1) mark the position of the masers in the R1 (left panel) and B1 (right) streamlines. We warn that the radial positions are stretched out to allow one to read the maser labels. Positional errors are not shown since they are smaller than the labels. In both panels, the black dashed line draws the sinusoid (see Table 2) fitted to the maser positions.

Fig. B.3 Distribution of velocities in the B2 stream. Maser labels (see Table 1) and error bars report values and corresponding errors for the velocity components along the R (upper panel), xl (middle panel) and los (lower panel) axes. The error bars of the los components are not shown, as they are smaller than the labels. Red labels denote a subset of masers with similar velocity components.

Fig. B.4 Sinusoidal fit of the B2 streamline. Maser labels (see Table 1) mark the position of the masers of the B2 stream. We warn that the radial positions are stretched out to allow one to read the maser labels. Positional errors are not shown since they are smaller than the labels. Red labels and arrows denote the positions and proper motions of a subset of masers showing similar velocities (see Fig. B.3). The black dashed line draws the sinusoid (see Table 2) fitted to the positions of the masers with black labels. The blue labels indicate a few masers scattered in position with respect to the sinusoid.

Appendix C The reference system co-moving with the young stellar object

As described in Sect. 2, the relative proper motions are significantly more accurate than the absolute proper motions. However, to use them to describe the motion with respect to the YSO, we need to determine the proper motion of the “center of motion” (feature #0) with respect to the YSO. A good reference for that is the absolute proper motion of feature #0: v_RA = 4.5 ± 1.7 km s⁻¹, v_DEC = 18.5 ± 1.7 km s⁻¹ (or, equivalently, v_xl = 15.3 ± 1.7 km s⁻¹, v_R = −11.3 ± 1.7 km s⁻¹). Although this proper motion is corrected for the parallax, solar motion and differential Galactic rotation, the virial motion of an individual YSO can still deviate from that of the harboring high-mass star-forming region up to ≈10 km s⁻¹ (Reid et al. 2019). In the following, we derive an independent estimate of the proper motion of feature #0 with respect to the YSO accurate within a few kilometers per second.

While the relative proper motions are calculated with respect to feature #0, the good knowledge of the YSO’s V_LSR allows us to determine accurate values of the velocity component along the line-of-sight, v_los. If for a set of masers v_los can be expected to be directly proportional to one of the two proper motion components, v_RA or v_DEC (or, equivalently, v_xl or v_R), the offset of that component can be determined.

First, we apply this method to the set of eight masers of the B2 stream whose proper motions are almost parallel (see Fig. B.4). The 3D velocities of these masers are also approximately parallel, as we have verified that they are inclined on average by only 13° with respect to a central direction. Let us indicate with M and δ the mean and standard deviation, respectively, of the values of one specific velocity component. We can expect that: M_xl/M_los = δ_xl/δ_los and M_R/M_los = δ_R/δ_los. Since the standard deviations are not affected by the presence of an offset, these relations allow us to calculate the corrected value of M_xl and M_R, and evaluate the proper motion of feature #0 with respect to the YSO. Using all or a subset of the eight masers, we always find consistent corrections resulting into velocity components of feature #0 varying within the small ranges: $v_{xl}^{f0}=17.5\pm 2.5$ km s⁻¹ and $v_R^{f0}=-12.5\pm 1.0$ km s⁻¹. We note that the derived correction agrees with the absolute proper motion of feature #0 reported above.

As discussed in Sect. 4.2, the masers belonging to the N region move radially away at a similar azimuthal angle of ≈304°. Since v_ρ dominates on the other velocity components (see Table 3), v_los and v_R are expected to be directly proportional. Using the relation: M_R/M_los = δ_R/δ_los, we derive $v_R^{f0}=-11.5$ km s⁻¹, consistent with the previous determination. In conclusion, the good agreement among these different methods make us confident that the corrected relative proper motions represent the velocities in the reference frame of the YSO with accuracies of a few kilometers per second.

References

All Tables

All Figures

	Fig. 1 Three-dimensional view of the MHD DW in IRAS 21078+5211 from simulations in Moscadelli et al. (2023). The observed maser positions (red and orange dots) are overlaid on streamlines (blue lines) computed from resistive-radiative-gravito-MHD simulations of a jet around a forming massive star.
In the text

Fig. 2 Reference frames used for maser positions and velocities, and magnetic fields. (a) Transformation of coordinates by rotation around the line-of-sight axis los over the position angle of the jet to bring Dec onto the sky-projected jet axis xl. (b) Transformation of coordinates by rotation around the axis R over the inclination angle of the jet to bring xl onto the jet axis ɀ. (c) Transformation of coordinates by rotation around the axis ɀ over the rotation angle ϕ of the maser position to bring oɀ along the direction êρ. (d) Decomposition of the maser velocity into the axial (vɀ), radial (vρ), and azimuthal (vϕ) components.

In the text

Fig. 3 Water maser 3D velocities in IRAS 21078+5211. Colored dots and arrows give water maser positions (relative to the reference feature # 0) and proper motions (relative to the YSO), respectively, with colors denoting the maser VLSR. We show only the masers for which the position at the first observing epoch can be derived. The amplitude scale for the velocity is reported in the bottom-left corner of the panel. The black dotted rectangles encompass the three regions to the N, NE, and SW of the YSO, where maser emission concentrates. The sky-projected jet and disk axes are in red and black, respectively, with distances labeled in astronomical units. An elongated maser stream directed approximately S-N is traced with a dashed-line segment labeled “B3.”

In the text

	Fig. 4 Time variability of the maser streamlines. Comparison of the positions of water masers observed in October 2020 (green triangles) and March-October 2023 (red dots) for the NE (left) and SW (right) regions.
In the text

Fig. 5 Three-dimensional velocities along the maser streamlines. The colored dots and arrows and the red axes have the same meaning as in Fig. 3. (Left panel) Expanding view of maser positions and proper motions in the NE region. The green and blue dashed curves are the fitted sinusoids to the B1 and B2 streamlines, respectively, whose parameters are reported in Table 2. (Right panel) Expanding view of maser positions and proper motions in the SW region. The red dashed curve is the fitted sinusoid to the R1 streamline, whose parameters are reported in Table 2.

In the text

Fig. 6 Co-rotating helical magnetic field line. The water masers (red balls) move along a helical magnetic field line (blue) anchored to the rotating disk. The rotation radius and angle, ℜ and ϕ, and the axial, Bɀ and vɀ, and azimuthal, Bϕ and v′ϕ, components of the magnetic field and velocity in the reference frame co-rotating with the disk are labeled. The tangent to the helix is inclined by an angle αB with respect to the ɀ axis.

In the text

Fig. 7 Determination of the launch radii of the maser streamlines. (a) The streamline to which each panel refers is indicated in the top-left corner of the panel. Maser labels (see Table 1) and error bars indicate the values (and corresponding errors) of the angular velocity of the launch radius inferred with the masers belonging to the streamline, which are plotted versus the position along the sky-projected jet axis xl. The dashed blue line and error bar give the mean and standard deviation of the consistent values, excluding the outliers denoted with red labels. (b) Vertical and horizontal error bars plot the values (and associated errors) of the angular velocity and rotation radius of the streamlines. The heavy and thin curves are the Keplerian profiles evaluated at three different YSO masses, as labeled. The launch radius of each streamline for the nominal value of the YSO mass of 5.6 M⊙ is identified with horizontal and vertical dotted lines, and the associated measurement error is marked by the horizontal error bar in the bottom.

In the text

Fig. 8 Change in the maser velocity direction and amplitude moving along the jet axis. The left and right panels refer to the B3 stream and N region, respectively. The inclination, θ, and azimuthal, ϕ, angles of the direction of the velocity (see Fig. 2c) and the velocity component along the jet, vɀ, are plotted versus the position along the jet, ɀ, in the top, middle, and bottom panels, respectively. Values and corresponding errors are denoted with dots and error bars, respectively, where red marks velocities ≲6 km s−1 in the B3 stream.

In the text

Fig. 9 Snapshot of the simulation when the mass of the protostar is 5 M⊙. The left panel shows the density, velocity, and magnetic field lines in the circumstellar gas. The right panel contains contours (solid lines) accompanied by shadows (semi-transparent thick lines of the same color as the contour). The shadow indicates in which direction the inequality listed in the legend holds true. A smoothing algorithm has been applied to the contours to improve clarity. For more details, we refer to Oliva & Kuiper (2023b).

In the text

	Fig. B.1 Non-simultaneous water masers. Colored dots and the red axis have the same meaning as in Fig. 3. At variance with the left panel of Fig. 5, all the detected masers are shown, including those observed only in a single epoch different from the first epoch. For these masers, missing the proper motion information, the first-epoch position cannot be derived.
In the text

Fig. B.2 Sinusoidal fits of the R1 and B1 streamlines. Maser labels (see Table 1) mark the position of the masers in the R1 (left panel) and B1 (right) streamlines. We warn that the radial positions are stretched out to allow one to read the maser labels. Positional errors are not shown since they are smaller than the labels. In both panels, the black dashed line draws the sinusoid (see Table 2) fitted to the maser positions.

In the text

	Fig. B.3 Distribution of velocities in the B2 stream. Maser labels (see Table 1) and error bars report values and corresponding errors for the velocity components along the R (upper panel), xl (middle panel) and los (lower panel) axes. The error bars of the los components are not shown, as they are smaller than the labels. Red labels denote a subset of masers with similar velocity components.
In the text

Fig. B.4 Sinusoidal fit of the B2 streamline. Maser labels (see Table 1) mark the position of the masers of the B2 stream. We warn that the radial positions are stretched out to allow one to read the maser labels. Positional errors are not shown since they are smaller than the labels. Red labels and arrows denote the positions and proper motions of a subset of masers showing similar velocities (see Fig. B.3). The black dashed line draws the sinusoid (see Table 2) fitted to the positions of the masers with black labels. The blue labels indicate a few masers scattered in position with respect to the sinusoid.

In the text