Snapshot of a magnetohydrodynamic disk wind traced by water maser observations

The formation of astrophysical objects of different nature and size, from black holes to gaseous giant planets, involves a disk-jet system, where the disk drives the mass accretion onto a central compact object and the jet is a fast collimated ejection along the disk rotation axis. Magnetohydrodynamic disk winds can provide the link between mass accretion and ejection, which is essential to ensure that the excess angular momentum is removed from the system and accretion onto the central object can proceed. However, up to now, we have been lacking direct observational proof of disk winds. This work presents a direct view of the velocity field of a disk wind around a forming massive star. Achieving a very high spatial resolution of ~0.05 au, our water maser observations trace the velocities of individual streamlines emerging from the disk orbiting the forming star. We find that, at low elevation above the disk midplane, the flow co-rotates with its launch point in the disk, in agreement with magneto-centrifugal acceleration where the gas is flung away along the magnetic field line anchored to the disk. Beyond the co-rotation point, the flow rises spiraling around the disk rotation axis along a helical magnetic field. We have performed (resistive-radiative-gravito-) magnetohydrodynamic simulations of the formation of a massive star and record the development of a magneto-centrifugally launched jet presenting many properties in agreement with our observations.

flow co-rotates with its launch point in the disk, in agreement with magnetocentrifugal acceleration where the gas is flung away along the magnetic field line anchored to the disk. Beyond the co-rotation point, the flow rises spiraling around the disk rotation axis along a helical magnetic field. We have performed (resistive-radiative-gravito-) magnetohydrodynamic simulations of the formation of a massive star starting from the gravitational collapse of a rotating cloud core threaded by a magnetic field. A magneto-centrifugally launched jet develops around the forming massive star which has properties matching many features of the maser and thermal (continuum and line) observations of our target. Our results are, presently, the clearest evidence for a magnetohydrodynamic disk wind, and show that water masers and forming massive stars provide a suitable combination of tracer and environment to allow us studying the disk-wind physics.
Magnetohydrodynamic (MHD) disk winds have been proposed to be the engines of the powerful jets observed at varying length scales in many diverse sources, from young stellar objects 32 (YSO) to black holes 4 . According to the classical model of an ideal MHD disk wind 4 , in the reference frame co-rotating with the launch point, the flow streams along the magnetic field line anchored to the accretion disk. An observer at rest sees magnetocentrifugal acceleration: the magnetic field keeps the flow in co-rotation with its launch point while its radial distance increases, till reaching the Alfvén point where the poloidal kinetic and magnetic energies are equal. Beyond the Alfvén point, the flow spirals outward Article number, page 2 of 36 along the rotation axis with a stably increasing ratio of the streaming onto the rotational velocity, until it gets eventually collimated into a fast jet 29,14 . So far, the best observational evidence for a MHD disk wind has been the finding of line of sight velocity gradients transversal to the jet axis, which are interpreted in terms of jet rotation and the imprint of the magneto-centrifugal acceleration 3,9,17,1 . However, that is an indirect evidence and the derivation of key parameters, as the launch radius and the magnetic lever arm, can be seriously affected by systematic biases 39 . On scales of ∼100 au, a few studies based on Very Long Baseline Interferometry (VLBI) maser (the laser equivalent at microwave band) observations have revealed rotating disk-like 19,21,37 , conical 22  In October 2020, we performed novel observations (see Fig. 2a) of the water maser emission in IRAS 21078+5211 by including all telescopes available in the VLBI network, with the aim to simulate next-generation radio interferometers which will improve current sensitivities by more than an order of magnitude (see Appendix A). In the following, we will show that these new observations prove that the water masers trace magnetized streams of gas emerging from the YSO's disk (see Fig. 2b and Appendix H). The maser emission concentrates in three regions to NE, N, and SW, inside the three dotted rectangles of  Fig. 2a), whose sky-projection is known from previous observations of the maser proper motions and radio jet (see Appendix E), we observe two elongated structures, blue-and red-shifted (with respect to the systemic V LSR of the YSO: V sys = −6.4 km s −1 ) to NE and SW, respectively. These structures are Article number, page 3 of 36 the opposite lobes of a collimated outflow from the YSO, located in between the two lobes; the disk axis (the black dashed line in Fig. 2a) is the intercept of the jet axis at the YSO position. From previous VLBA observations we know that the jet axis has to lay close to the plane of the sky, with an inclination ≤ 30 • . According to the maser V LSR , the jet is inclined towards us to NE, and away from us to SW.
The jet and disk axes provide a convenient coordinate system to refer the maser positions to. In the following, we present the interpretation of the maser kinematics, which is based on the analysis of the three independent observables: z, the elevation above the disk plane (or offset along the jet), R, the radial distance from the jet axis (or transversal offset), and the maser V LSR . As discussed in Appendix A, the accuracy of the maser positions is ≈ 0.05 au, and that of the maser V LSR ≈ 0.5 km s −1 . Without loss of generality, we can express the maser velocities as the sum of two terms, one associated with the toroidal component or rotation around the jet axis, V rot , and the other associated with the poloidal component including all the contributions owing to non-rotation, V off . Since the jet axis is close to the plane of the sky and we observe the rotation close to edge-on (see Fig 3), we can write: where φ is the angle between the rotation radius R and the line of sight, and ω and t are the angular velocity and the time, respectively. Fig. 4c shows the remarkable finding that the spatial coordinates z and R of the maser emission in the SW flow satisfy the relation: where C, the amplitude of the sinusoid, f z , the spatial frequency, and z 0 , the position of zero phase, are fitted constants (see Table 1). In Appendix B we demonstrate that the masers in the NE flow can be separated in three different streams, each of them satisfying the relation 4 (see Fig. 5c and Table 1). The comparison of Eqs. 2 and 4 leads to a straightforward interpretation of the sinusoidal relation between the coordinates by taking: The former equation indicates that the rotation radius is the same for all the masers, the latter shows that the motions of rotation around and streaming along the jet axis are locked together, which is the condition for a spiral motion.
Denoting with V z the streaming velocity along the jet axis, we can write |z − z 0 | = V z t and, comparing with Eq. 3, we derive the relation between the rotation and streaming Article number, page 4 of 36 velocities: According to Eq. 5, the observation of a well defined sinusoidal pattern requires that ω and V z are directly proportional, or constant. The constancy of V z implies that V off is also constant, because, if the rotation radius does not change, V off is the projection along the line of sight of V z . Following Eqs. 1 and 2, the constancy of ω and V off would result into a tight linear correlation between V LSR and transversal offsets R. While a good linear correlation between V LSR and R is observed for the SW and NE-1 spiral motions (see Figs. 4b and 5b, black symbols), the scatter in velocity is considerable for the NE-2 spiral motion (see Fig. 5b, red symbols). In Appendix D we investigate the physical reason for the observation of well defined sinusoidal patterns despite the presence of a significant velocity scatter. Applying the equations of motions for an axisymmetric MHD flow, we find that the magnetic field configuration has to be helical over the maser emission region, and the motion along such an helical field line, in the reference frame co-rotating with the launch point, leads to the sinusoidal pattern of maser positions.
We consider now the N region (see Fig. 2a) and show that, in this region as well, the maser kinematics is consistent with the predictions for a MHD disk wind. The N masers have a larger separation from the jet axis than the NE and SW masers. A few nearby masers show quite different V LSR , which could hint at distinct streams, as observed (see Fig. 5a) and discussed (see Appendix B) for the NE flow. In this case, however, only a single stream is reasonably well sampled in position and V LSR with the masers, and we focus our kinematical analysis on that. The spatial distribution of this stream presents an arc-like shape: a subset of maser features draws a line at small angle with the disk axis and another group extends at higher elevations about parallel to the jet axis (see Fig. 6a). Fig. 6c shows that the maser V LSR increases linearly with R in the range −95 au R −70 au up to an elevation z ≈ 60 au. The relatively large separation from the jet axis and radial extent, and the arc-like shape of the maser distribution lead us to think that the N emission is observed close to the plane of the sky. In this case, the maser V LSR should mainly trace rotation, which is also expected to be the dominant velocity component at low elevations above the disk. Then, the good linear correlation between V LSR and R indicates that the masers co-rotate at the same angular velocity, ω N = 0. Eqs. 1 and 2, the derivation of ω N does not depend on the maser geometry. Therefore, the finding that the masers lay along a line intercepting the disk at ≈ −40 au provides an "a posteriori" test of the assumption that the N emission is observed close to the plane of the sky.
The masers found at elevation z > 60 au appear to set aside of the arc-like distribution (see Fig. 6a) and their V LSR are significantly more negative and do not follow the linear correlation with the radius (see Fig. 6c). A natural interpretation is that the location at where ω K is the Keplerian angular velocity of the launch point. Eq. 6 shows that ω e is the angular velocity of the spiraling trajectory as observed in the reference frame corotating with the launch point. Based on axisymmetric MHD models 31,39 , the ratio ω K /ω Article number, page 6 of 36 increases stably from 1 up to a value ≈ 4 while the gas climbs from z A to 10 R K , where z A is the elevation of the Alfvén point and R K is the launch radius. Being ω ≤ ω K , the negative value of ω e indicates that the rotation angle of the maser positions decreases with z. Following the previous discussion, Eq. 5 has to be corrected by replacing ω with ω e : A good test of the above considerations comes directly from our data. Assuming that all the masers move along a single trajectory and using the fitted values of ω and f z in Eq. 5, we obtain implausibly small values for V z : 31, 37, 17 and 49 km s −1 , for the SW, NE-1, NE-2, and NE-3 spiral motions, respectively. There are two strong observational evidences that the derived streaming velocities are too small. First, comparing them with the values of V off (see Table 1 Since V off corresponds to the line of sight projection of V z , we can write: where V off is corrected for the systemic V LSR of the YSO, and i axi is the inclination angle of the jet axis with the plane of the sky. As we know that i axi ≤ 30 • , Eq. 8 allows us to derive a lower limit for V z , reported in Table 2. Using the derived lower limit of V z and the corresponding value of f z (see Table 1), by means of Eq. 7 we can calculate a lower limit for ω e . Finally, we use Eq. 6 and the fitted value of ω (see Table 1) to infer a lower limit for ω K = ω +|ω e | and, knowing the mass of the YSO, a corresponding upper limit for the launch radius R K (see Table 2). The NE-1 stream, which extends the most in elevation (from 20 to 130 au, see Fig. 5a), includes a group of maser features at elevation of ≈ 20 au, which should be located closer to the Alfvén point. In Appendix F, we study the change of V LSR versus R internal to this cluster and obtain a lower limit for ω K ≥ 2.4 km s −1 au −1 that agrees well with the corresponding value reported in Table 2 for the NE-1 stream.
The results shown in Table 2 are in general agreement with the theoretical predictions from MHD disk winds 34,33 : the inferred lower limits on the streaming velocity V z increase steeply with decreasing launch radius; for the SW and NE-3 streams, which are those with relatively lower uncertainty in the measurement of both ω and R (see Table 1), the ratio between V z and the rotation velocity ω R is ≥ 2-3.
In conclusion, our observations resolve, for the first time, the kinematics of a MHD disk wind on length scales of 1-100 au, allowing us to study the velocity pattern of individual streamlines launched from the disk. As represented in Fig. 2b, close to the disk rotation axis we observe flows spiraling outward along a helical magnetic field, launched from locations of the disk at radii ≤ 6-17 au. At larger separation from the rotation axis, we observe a stream of gas co-rotating with its launch point from the disk at radius of ≈ 40 au, in agreement with the predictions for magneto-centrifugal acceleration. Our interpretation is supported by (resistive-radiative-gravito-) MHD simulations of the formation of a massive star that lead to a magneto-centrifugally launched jet whose properties agree with our maser and thermal (continuum and line) observations of IRAS 21078+5211. These results provide the best evidence for a MHD disk wind to date. Since water maser emission is widespread in YSOs, sensitive VLBI observations of water masers can be a valuable tool to investigate the physics of disk winds.

Acknowledgments
We Notes. Column 1 denotes the maser stream; Cols. 2 and 3 provide the values of ω and V off from the linear fit of maser VLSR versus R; Cols. 4, 5 and 6 report the amplitude, the spatial frequency and the position of zero phase, respectively, of the sinusoidal fit of the maser coordinates R versus z.
a The determination of this error, smaller than the value, 1.7 km s −1 au −1 , from the linear fit, is discussed in Appendix C.
Notes. Column 1 indicates the maser stream; Col. 2 reports the estimated streaming velocity along the jet axis; Cols. 3 and 4 give the estimate of the angular velocity and the radius, respectively, at the launch point.   Table 1.  Table 1. We plot the linear transformation of the radii to reduce the overlap and improve the visibility of each of the three streams. The dashed curves are the fitted sinusoids. In Table 1, we report the parameters of the sinusoidal fits of the radii.
Article number, page 14 of 36 we estimate that the error on the absolute position of the masers is 0.5 mas. 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: ∆θ spot = ∆θ 2 fit + ∆θ 2 bandpass . The first term depends on the SNR of the data, following 35 :

Appendix B: Resolving the NE emission into three distinct streams
The masers in the NE emission are separated into three distinct streams identified with different symbols: dots, triangles and squares (see Fig. 5a). Each stream traces a sinusoid in the plane of the sky (see Fig. 5c), which is the signature for a spiral motion. The identification of the maser features belonging to each of the three sinusoids can be done with good confidence. First, we note that the masers at elevation z ≥ 90 au can be unambiguously divided into two streams, dots and triangles, on the basis of the very different V LSR of nearby emissions. Both dots and triangles at z ≥ 90 au have R decreasing with z, suggesting that they could trace the descending portion of a sinusoid, and one can argue that the corresponding ascending portion would be traced by masers with z < 90 au. Next, we note that the group of masers with 30 au ≤ z ≤ 60 au draw an arc, and, since they cannot be part of a sinusoid together with other masers at larger R, they have to trace a third sinusoid with maximum radius close to the radius, R apex = 17.2 au, of the apex of the arc. We have also verified that the blue-shifted cluster of masers located at elevation z ≈ 20 au (inside the dashed rectangle of Fig. 5a) cannot be adjusted within this third sinusoid. We have fitted the expression:  This method allows to appreciate the overall spatial structure traced by maser emission with very different brightness levels.
fixing the maximum velocity to ω R apex and fitting the zero level ω R apex − C, with the amplitude C as a free parameter. Fig. C.1 shows that, apart a minor fraction of scattered results, the radius estimated from the sinusoidal fit decreases regularly with ω, and attains a value consistent with the observations, 17.5 ± 2.5 au, for ω = 2.0 ± 0.5 km s −1 au −1 , thus effectively reducing the uncertainty to 0.5 km s −1 au −1 .

Figs. 4b and 5b
show that the maser V LSR are linearly correlated with R in the SW and NE-1 streams. For the masers belonging to the NE-2 stream, the measurement scatter from the linear fit of V LSR versus R is considerable (with large fit errors, see Table 1) which appears to be too large with respect to the observed maser V LSR . In the following, we investigate the physical reason for the observation of well defined sinusoidal patterns despite the co-occurrence of significant velocity scatters.
In a stationary, axisymmetric MHD flow, the two fundamental equations of motions 30,34,33 linking velocity and magnetic field along a field line are: where V p and V φ , and B p and B φ are the poloidal and toroidal components of the velocity and magnetic field, respectively, ρ is the gas mass volume density, ω K is the Keplerian angular velocity at the launch point, R is the rotation radius, and k is the "mass load" of the wind, expressing the fixed ratio of mass and magnetic fluxes along a given magnetic field line. Since V φ − ω K R is the toroidal velocity in the reference frame co-rotating with the launch point, Eqs. D.1 and D.2 lead to the well-known result that the velocity and magnetic field vectors are always paralell in the co-rotating reference frame.
where ω is the angular velocity of the trajectory at radius R, the two equations above can be combined into: where we have used the definition of ω e in Eq. 6. We can define the magnetic field helix angle α B = arctan(B φ /B p ), which is the angle with which a helical field line winds around the jet axis.
We have seen that the observation of well defined sinusoids in the plane of the sky requires that the rotation radius R, and the ratio of the effective angular velocity ω e on the streaming velocity V z keep constant (see Eq. 7). Since the poloidal velocity V p is equal to V z if R is constant, Eq. D.3 implies that α B does not vary along each of the observed maser streams. However, the reverse argument holds too. If the magnetic field is sufficiently stable and has a constant helix angle (that is, a helical configuration), the motion along the field line (in the co-rotating reference frame) requires a constant value of the ratio ω e /V z . That preserves the maser sinusoidal pattern despite the concomitant presence of a relatively large velocity scatter. From Eqs. D.3 and 7, with V p = V z , we have: Using the values reported in Table 1  to the latter across a few 10 au (see Fig. 2a). The interpretation of the water masers in IRAS 21078+5211 as external shocks at the cavity wall of a wide-angle wind faces also the difficulty to account for the huge maser V LSR gradients, ≥ 40 km s −1 over ≤ 10 au, transversal to the NE stream at elevation z ≥ 90 au (see Fig. 5a). If the masers, instead of delineating two streams at different velocities, traced shock fronts (seen close to edge-on) at the wind wall, it would be very difficult to explain such a large velocity difference for almost overlapping emissions.
We conclude this discussion on alternatives to the disk wind interpretation showing that it is improbable that we are observing multiple outflows from a binary system. Inspecting   Table 1), determined at larger elevations. That can be interpreted as evidence that the gas closer to the Alfvén point rotates faster than that at higher elevations, in agreement with the MHD disk-wind theory. As described in Appendix A, each maser feature is a collection of many compact spots for which we have accurately measured position and could indicate that the internal V LSR gradient is influenced by non-rotation terms, as, for instance, the increase of the streaming velocity with the elevation. In any case, assuming that also for these two features the internal V LSR gradient is dominated by rotation, we obtain a consistent average value for the in-feature angular velocity of ω = 2.4 km s −1 au −1 .
We stress that the internal motion of the maser features does not contribute to the features' V LSR used in our analysis of the maser velocity pattern on larger scales (see   Table 2, where higher streaming velocities correspond to smaller launch radii. The southern features are distributed in three separated clusters (see Fig. 2a ratio, which we take as 20 times the critical (collapse-preventing) value 27 and corresponds to a relatively weak initial magnetic field. A constant value of the opacity of 1 cm 2 g −1 was used to model the gas and dust, as well as an initial dust-to-gas mass ratio of 1%.
We used an axisymmetrical grid of 896×160 cells in spherical coordinates, with the radial coordinate increasing logarithmically with the distance to the center of the cloud.
An inner boundary of 3 au was set up, inside of which the protostar is formed through accretion. No flows are artificially injected from the inner boundary into the collapsing cloud.
The simulation starts with an initial gravitational collapse epoch. After ∼ 5 kyr, enough angular momentum is transported to the center of the cloud to start forming an accretion disk that grows in size over time. Roughly at the same time, we observe the launch of magnetically-driven outflows. Magnetic pressure arising from the dragging of magnetic field lines by the rotating flow eventually overcomes gravity and seeds the formation of the outflow cavity, thrusting a bow shock in the process (Fig. H.1d), which propagates outwards as the cavity grows in size. Previous observations of IRAS 21078+5211 have uncovered the presence of a bow shock located at distances of ≈ 36000 au from the forming massive YSO 26 . The initial launch of the magnetically-driven outflows provides a possible formation mechanism for the observed bow shock and in return, the propagation of the bow shock provides an estimation for the age of the system.
In the simulation, the protostar reaches a mass of 5.24 M (a value in the expected mass interval from observations) after 13.84 kyr of evolution. We estimate that the bow shock has propagated to a distance of ≈ 30000 au at that time, roughly in line with the observations. At the same time, the accretion disk has grown to about 180 au in radius, in agreement with the observational estimates 26 , as well. The data reveal a magnetocentrifugally launched jet, in a similar way as reported by the literature 4,13 , however, we see that the launching region of the jet is narrowed by the ram pressure of the infalling material from the envelope and the presence of a thick layer of the accretion disk which is vertically supported by magnetic pressure. Material transported from large scales through the accretion disk reaches the launching region, located at z 100 au (see Fig. H.1b), where the centrifugal force is stronger than gravity, and simultaneously, where the flow We provide this comparison as an example that a magneto-centrifugally launched jet around a forming massive star yields a consistent picture with the observations. However, the coincidences should be taken with caution, as different combinations of initial conditions may yield similar results, which means that a wider and deeper investigation is needed in order to determine the conditions of the onset of star formation from the observational data.