© The Authors 2023

1 Introduction

During the formation of massive stars, collimated bipolar outflows are launched, which become a marker of star formation and a source of mechanical feedback into their natal environment. They transport angular momentum into large scales, enhancing accretion through their accompanying disk. There is mounting observational and theoretical evidence, especially in recent years, that highly collimated protostellar outflows are driven by magnetic processes.

The presence of collimated jets around young massive stellar systems is well documented: Purser et al. (2021, 2016), for example, used radio continuum observations to build catalogs of candidates of jets around massive protostars, and McLeod et al. (2018) were able to even detect a jet around a massive young star from an extragalactic source. Guzmán et al. (2010) and Sanna et al. (2019) performed radio observations of jets around massive protostars for the star-forming regions IRAS 16562-3959 and G035.02+0.35, respectively. Synchrotron emission has been observed coming from the massive protostellar jets HH 80-81, Cep A HW2, and IRAS 21078+5211 (see Carrasco-González et al. 2010, 2021; Moscadelli et al. 2021) in regions not close to the protostar (at least a few hundred astronomical units away). The origin of this nonthermal emission has been attributed to a recollimation of the jet due to magnetic forces (Carrasco-González et al. 2021). In Carrasco-González et al. (2021), two components of the magnetically driven outflows are observed: a fast jet (with speeds of at least 500 km s⁻¹) and a wide-angle wind. Substructures with proper motions have been observed in the outflow ejecta (Obonyo et al. 2021), indicating nonsmoothness in the ejection process.

A new generation of observations that combine data from new or upgraded (sub)millimiter interferometers such as the Northern Extended Millimeter Array (NOEMA), the Atacama Large Millimeter/submillimeter Array (ALMA), the Very Long Baseline Interferometry network (VLBI), and the Jansky Very Large Array (JVLA) have enabled an unprecedented resolution of the innermost ~100 au around the forming massive star, such that the direct study of the dynamical processes that intervene in the launching (Moscadelli et al. 2022) and collimation (Carrasco-González et al. 2021) of the jet is now possible. However, previous numerical studies of massive protostellar outflows (for example, Rosen & Krumholz 2020; Machida & Hosokawa 2020 or Mignon-Risse et al. 2021; a more detailed literature review is offered in Sect. 5) have not been able to spatially resolve the launching mechanisms of the fast component of the magnetically driven collimated jets. The study by Kölligan & Kuiper (2018) has shown that a self-consistent launching of a fast magnetic outflow driven by the magneto-centrifugal mechanism was possible.

Thanks to observations of water masers done in Moscadelli et al. (2022), the innermost ~ 100 au region where the jet is launched can now be observed with a resolution of up to 0.05 au. In that article, we presented a simulation of a forming massive star surrounded by a disk-jet system that was able to reach those resolutions and reproduced the observed jet kinematics reasonably well. That simulation was part of a larger set of 31 simulations in total, and subjected to a deeper analysis of the physical processes and evolution of the disk and the magnetically driven outflows.

In a previous article (Oliva & Kuiper 2023, hereafter Paper I), we focused on the study of the effects of magnetic fields in the dynamics of the disk. We found that the disk can be distinguished into two vertical layers supported by thermal and magnetic pressure, respectively. Additionally, we extensively studied the size of the disk and the impact of magnetic braking under several star-forming environments. In this article, we focus on the dynamics of the magnetically driven outflows: their launching processes, acceleration, propagation, and termination. The use of the full parameter space we introduced in Paper I enables us to deepen our understanding of how those processes develop, and how the initial conditions for the gravitational collapse influence the resulting outflow, its structure, and propagation.

In Sect. 2, we summarize the setup and methods we used in Paper I. Section 3 deals with a dynamical analysis of the launching, acceleration, propagation, and termination of the magnetically driven outflows, focused on the fiducial case of our simulation catalog. In Sect. 4, we study how the conditions for the onset of gravitational collapse affect the protostellar outflows produced by the forming massive star. In Sect. 5, we offer a literature review and a comparison of our results with previous numerical studies. Section 6 contains a summary and the major conclusions of our work. Finally, we examine the numerical convergence of our results in Appendix A.

2 Methods and parameter space

Here, we summarize the methods and configuration of the simulation catalog already presented in Paper I. For a more in-depth discussion of the equations solved, and the solution methods used, we point the reader to that reference.

We simulate the formation of a massive star by following the gravitational collapse of a slowly rotating cloud core of 0.1 pc in radius, threaded by an initially uniform magnetic field, until the formation of a disk-jet system around the massive (proto)star. At the onset of gravitational collapse (initial configuration), the cloud core has a mass M_C, distributed along a density profile $\rho \propto r^{{\beta }_{\rho }}$ (r being the spherical radius), a rotational-to-gravitational energy ratio ζ, an angular velocity distribution according to $\text{Ω\hspace{0.17em}}\propto R^{{\beta }_{\text{Ω}}}$ (with R the cylindrical radius), and a uniform magnetic field parallel to the rotation axis whose strength is determined by the mass-to-flux ratio $\overline{\mu }$ in units of the critical (collapse-preventing) value. We used the Machida et al. (2007) resistivity model (both with the full temperature and density dependence and with a fixed temperature to save computational power and to allow for a direct comparison with the models by Kölligan & Kuiper 2018), and two other cases for comparison: neglecting the resistivity (ideal MHD) and neglecting magnetic fields.

For modeling the weakly ionized gas, we used the methods of magnetohydrodynamics with the code Pluto (Mignone et al. 2007), with additional modules for Ohmic resistivity (with the model by Machida et al. 2007), self-gravity (Kuiper et al. 2010), and the radiation transport of thermal emission from the dust and gas (Kuiper et al. 2020) with the gray flux-limited diffusion approximation. We assume axisymmetry and equatorial symmetry in the simulations. The inner boundary of the simulation domain, which we refer as the sink cell, models the (proto)star. Matter that flows into the sink cell is considered as accreted by the protostar and cannot escape. The magnetic field lines are prolonged across the sink cell with a zero-gradient boundary condition. No flows are artificially injected into the cloud from the protostar; the outflow launching occurs self-consistently within the computational domain. The accretion disk is massive and is expected to form spiral arms and fragment (Oliva & Kuiper 2020). In order to model the gravitational torque from the spiral arms in the azimuthal direction, we use shear α-viscosity, following the setup and results of Kuiper et al. (2011). Inside the outflow cavity, the Alfvén speed increases significantly due to the low densities and strong magnetic fields. In order to avoid extremely small numerical time steps, a limit to the Alfvén speed was used (a detailed discussion of its effects on the outflows is offered in Appendix A.2).

Convergence of the results was studied in the parameter space by using five different grids of increasing resolution, as well as several sizes for the inner boundary (which we refer as the sink cell, and which models the forming protostar). Grid × 16, the highest resolution one in our series, consists of 896 logarithmically spaced cells in the radial direction, and 160 linearly spaced cells along the polar angle. This translates to a minimum grid cell size is 0.03 au at 3 au from the forming massive star and 10 au at r = 1000 au. The rest of the grid sizes in the series are obtained by dividing by two the number of cells of each direction. Figure 1 presents a schematic summary of the parameter space considered in this study. The values of each parameter corresponding to the fiducial case are highlighted in blue. To simplify nomenclature, we use the following convention: modifications of the fiducial values are written followed by the name of the corresponding grid. For example, $\overline{\mu }=5\left[\times 8\right]$ means a cloud of a mass-to-flux ratio of 5, with the rest of the parameters as in the fiducial case run on the grid ×8.

3 Physical processes in the outflows

This section is focused on a detailed examination of the outflows produced in the fiducial simulation of our series. We compute several force terms from the magnetohydrodynamical equations and compare their relative strength with the aim of understanding the physical processes that intervene in the launching, acceleration, propagation and termination of the different components of the magnetically driven massive protostellar outflows.

3.1 Magneto-centrifugal launching

The low-density, bipolar jet cavity starts to appear close to the protostar shortly after the formation of the disk at around t ~ 5 kyr, in a process described in more detail in Sect. 3.7, which results in an initial bow shock being thrust out. Once the high-speed component of the outflow has established, (and similarly to Kölligan & Kuiper 2018) we find evidence for the magneto-centrifugal mechanism described in Blandford & Payne (1982) inside of the cavity, although with some key geometrical differences which we explain below.

A brief note on the visualizations that make use of directed contours throughout the article: an understanding of the launching mechanisms for outflows requires the simultaneous visualization of several intervening forces, morphological features, and other quantities that carry dynamical information. With this end, and taking advantage of the axial symmetry of our models, we present here a system of plots containing contour lines which aided the interpretation of our results and the drawing of our conclusions. The first instance of such visualizations is given in the right panel of Fig. 2. Each solid contour line is followed by an adjacent shadow of the same color, which indicates the direction on which the inequality listed in the legend holds true. Throughout the article, we refer to both the contour and its corresponding shadow as a directed contour¹.

First, we perform a dynamical examination of the launching mechanism of the fast jet. Figure 2 displays a typical scenario during the magneto-centrifugal phase of the simulation. The left panel provides a joint view of the gas mass density (background color map), poloidal velocity (as a vector field with arrows) and magnetic field (as field lines in solid color), and the right panel presents the directed contours that separate the fulfillment of conditions on two quantities of dynamical interest. The blue directed contour indicates the equilibrium between the (corotating) centrifugal acceleration $${\left(a_g\right)}_R=\frac{{\upsilon }_{\varphi }^2}{R},$$(1)

and the cylindrical radial component of the gravitational acceleration $${\left(a_c\right)}_R=\frac{GM\left(r\right)}{r^2}\sin \theta ,$$(2)

where M(r) is the mass enclosed in a spherical radius r (including the mass of the protostar), and which constitutes an approximation to the use of the full gravitational potential of stellar gravity and gas self-gravity. In the shadowed blue region, the gravitational force dominates over the centrifugal force along the cylindrical radial direction. A disadvantage of the use of contours to visualize a ratio of forces is the poor handling of regions where equilibrium holds, as opposed to regions with clear dominance of a force over another. This is the case in the thick and thin layers of the disk, for which the contour is no longer adequate. The equilibrium of forces in the disk was studied in detail in Paper I: in a nutshell, both the thick and thin disks are roughly Keplerian, with thermal and magnetic pressure causing small local regions of sub- and super-Keplerianity, but with the material always being transported inward.

The yellow directed contour delimits the Alfvén surface, that is the surface where the speed of the plasma is lower than the Alfvén speed $${\upsilon }_A=\frac{\left|B\right|}{\sqrt{4\pi \rho }},$$(3)

where the magnetic field strength $\left|B\right|$ is expressed in Gauss. For the purpose of our analysis, we use the corotating frame of reference, which means that we take only the magnitudes of the poloidal components of the velocity and magnetic field². Outside of the Alfvén surface, perturbations to the originally vertical magnetic field (in this case, caused by the gravitational collapse) occur faster than the restoring force of magnetic tension, resulting in the typical hour-glass morphology of gravitational collapse of a gaseous cloud in rotation (for an observational example, see Beuther et al. 2020). The contrary happens inside the Alfvén surface: the gas flow is forced to follow the magnetic field lines thanks to the restoring magnetic forces. The sub-Alfvénic regime is reached because of the strong magnetic fields close to the protostar (caused by the gravitational drag of the magnetic field) and the low density in the cavity.

The simulations presented include magnetic diffusion in the form of Ohmic resistivity, which acts mostly in the thin and thick disks (see Paper I) but its action is severely limited in the cavity due to its low density (and therefore, low resistivity). This implies that in the cavity, the plasma is essentially flux-frozen. This can be readily seen while comparing the directions of the velocity field and magnetic field lines in the cavity as presented in Fig. 2. We note, however, that while the plasma is flux-frozen in the simulations, this might not be fully the case in reality because other forms of magnetic diffusion (ambipolar diffusion and the Hall effect) are not considered in this study.

We can summarize the magneto-centrifugal mechanism as seen in our results in the following way. Consider a parcel of gas that has fallen in from large scales onto the thick disk. It is then transported through the accretion disk until it reaches the cavity wall. If the parcel enters the cavity, it experiences (a_c)_R >(a_g)_R, and the parcel would tend to move cylindrical-radially outward in the absence of magnetic fields. However, the parcel enters simultaneously the sub-Alfvénic regime, which means that the flow must follow the magnetic field lines resulting in the parcel being launched radially outward. The acceleration of the parcel inside the cavity is detailed in Sect. 3.2.

The model by Blandford & Payne (1982) was proposed in the context of accretion disks surrounding black holes, that is, an infalling envelope was not included, which is necessary in the study of the early stages of protostellar jets. Even though we find general agreement in the launching mechanism of the high-velocity jet, we observe that the jet launching area is not located on top of the full accretion disk. The ram pressure from the infalling envelope, and the effect of magnetic pressure on the disk (the existence of a thick layer in the disk) constrain the formation of the cavity to a conic region in the innermost scales, which means that the launching surface of the jet is located on the portion of the cavity wall that is in contact with the incoming material from the thick accretion disk. The densities found in the envelope and the thick disk impede the flow to become sub-Alfvénic except in the cavity, close to the launching area. The Alfvén surface does not occupy the full cavity, but it is restricted to a narrow region close to the rotation axis and reaches sizes of ~ 1000 au, after which the magnetic field becomes too weak. Beyond the Alfvén surface, the flow is collimated in a process described in Sect. 3.6.

Fig. 1 Schematic overview of the parameter space covered by the simulation series.

Fig. 2 Magneto-centrifugal mechanism for the fiducial case in ×8.

3.2 Acceleration of material in the jet

In this section, we continue the discussion of the trajectory of a parcel of plasma as it accelerates throughout the cavity and is expelled into larger scales. The right panel of Fig. 3 provides a view of the gravito-centrifugal equilibrium and the Alfvén surface as discussed above, but also presents a directed contour (in purple) that delimits the region where the magnetic field lines are mostly wound by rotation. In the launching area (marked as a), the magnetic field lines are mostly poloidal, and a (sub-Alfvénic) parcel of plasma is then forced to follow a mostly radial trajectory until it reaches the Alfvén surface.

To see this explicitly, we take the energy equation of the plasma that, ignoring radiation transport and considering a stationary state, reduces to $$\nabla \cdot \left[\left(E^g+E^K+E^{\text{th}}+P\right)\upsilon +cℰ\times B\right]=0,$$(4)

with the first terms inside the parenthesis being the gravitational, kinetic and thermal energy densities, and ℰ being the electric field. The presence of the gravitational energy density as stated is possible because it does not change significantly over time during the acceleration process. Given that in the cavity we can ignore Ohmic dissipation, cℰ = −υ × B. As discussed, the sub-Alfvénicity of the flow close to the launching area causes $B\Vert \upsilon$, which implies cℰ × B = −(υ × B) × B = 0 and allows us to conclude from Eq. (4) that $$E^g+E^K+E^{\text{th}}+P=\text{constant}$$(5)

along the flow in the cavity. This expression enables us to study the process of acceleration of material in the jet, illustrated in the left panel of Fig. 3. The background map represents the radial velocity (for which the color scale ignores the negative values), and the directed contour highlights the region where the azimuthal component of the velocity is higher in magnitude than its poloidal value. Outside of the disk and in the outer parts of the cavity, the velocity turns from mostly azimuthal to mostly poloidal. As the parcel is forced to follow the magnetic field lines from the position marked as a to b in Fig. 3, it travels to larger r, where the gravitational energy is lower. Then, according to Eq. (5), the kinetic energy of the parcel has to increase, given our result from Paper I that the contribution of the thermal energy and pressure are negligible. The parcel accelerates υ_r ~ 0 to ~400 km s⁻¹ along the spherical-radial direction. As a consequence of this process, the highest velocities of the jet at small scales are not found close to the rotation axis.

Once the parcel of plasma leaves the Alfvén surface (mark c in Fig. 3), it enters the region where the magnetic field lines are mostly wound. The trajectory of the flow becomes helical as the flow partially follows the magnetic field lines, but an additional upward acceleration is experienced by the parcel because of magnetic pressure (more details of this process are offered in Sect. 3.4). The jet now propagates away from the protostar and into larger scales of the cloud core.

A more thorough examination of the jet speeds across the whole simulation series reveals that speeds of at least 100 km s⁻¹, and over 1000 km s⁻¹ are reached in certain regions. However, the maximum speed reached by the jet is not a fully converged value for the reasons discussed in Appendix A.

The most direct observational evidence of the magneto-centrifugal mechanism and the acceleration of the jet as described here is provided by the recent observations of the most massive young stellar object (M_⋆ = 5.6 ± 2.0 M_⊙) in the star-forming region IRAS 21078+5211 (Moscadelli et al. 2022, 2021). In Moscadelli et al. (2022), 22 GHz water maser observations (which trace internal shocks of the gas), were able to uncover the helical trajectories of the plasma as it is magneto-centrifugally ejected. In that article, we used the simulation corresponding to the fiducial case with grid × 16 to produce a snapshot comparison when the protostellar mass is similar to the observed value, and found that the geometry and kinematics of the transition region where the flow abandons the Alfvén surface and becomes helical coincide very well with the observed masers. Furthermore, the high velocities of the water masers observed in IRAS 21078+5211 (~200 km s⁻¹) are fully compatible with the values found in the numerical simulations.

Fig. 3 Acceleration of material along the magnetic field lines.

Fig. 4 Cavity wall ejection.

3.3 Cavity wall ejections

The cavity wall contributes to both inflows and outflows episodically. When a section from the cavity wall enters the region where centrifugal forces dominate over gravity, the section starts to move outward, forming a toroidal prominence, as illustrated in Fig. 4. This prominence grows quickly and develops a neck, thanks to the material infalling from larger scales that provides an inward flow that narrows the neck. The densities present on the cavity wall are sufficient for magnetic diffusion to take place, allowing magnetic field lines to reconnect. This in turn allows the protuberance to decouple from the rest of the cavity wall and it is ejected through the cavity thanks to the magneto-centrifugal mechanism.

The solid pink directed contour in the right panel of Fig. 4 shows the regions where magnetic pressure dominates over thermal pressure, that is, it delimits the region where $${\beta }_p\equiv \frac{P}{P_B}=\frac{E^{th}}{E^B}<1.$$(6)

The contour confirms that the thin disk is thermally supported while the thick disk is magnetically supported in the vertical direction, as it was concluded from Paper I. The densest parts of the cavity wall become thermally dominated, as magnetic diffusion becomes high enough in those regions and then magnetic energy is dissipated into thermal energy. We also display the contour for β_p = 0.01, which temporarily includes the pro-tuberance and shows that some magnetic energy is dissipated into thermal energy as the magnetic field lines are reconnected. After the ejection event takes place, the magnetic field lines are quickly reconfigured in the cavity wall in a similar way to Fig. 2. We remark, however, that the ejection mechanism for the cavity wall ejections is the magneto-centrifugal mechanism and not the energy released by magnetic diffusion during reconnection. The reconnection of magnetic field lines provides a way for the mass to be released into the cavity.

The cavity wall ejections provide not only a mechanism for episodic ejection, but also a mechanism for episodic accretion, as the infall from the envelope onto the protostar is temporarily halted while the ejection event takes place. The observed timescale for the formation and ejection of the protuberance is ~ 10 yr, with the reconnection event happening in ≲ 1 yr, although those values increase with the elapsed time as the disk grows in size.

We remind the reader that the simulations are performed under the assumption of axisymmetry, and therefore, a three-dimensional study that resolves similar scales would be required to know the role that a turbulent flow would have on ejected material from the cavity wall. For example, turbulent reconnection might make cavity wall ejections more frequent, asymmetric and smaller in size, instead of episodic and large in size, as we find. Additionally, the inclusion of ambipolar diffusion and the Hall effect may affect the reconnection process or not allow the formation of the prominence at all.

3.4 Magnetic tower flow

We observe a different kind of magnetically driven outflow which we refer to as the magnetic tower flow, and corresponds to the outflow mechanism described in Lynden-Bell (2003). As time progresses, the dragging of magnetic field lines by rotation increase magnetic pressure until its gradient is eventually larger than the local gravitational force. At large scales, outside of the Alfvén surface, this force drives the magnetic tower flow.

An example of the tower flow is given in Fig. 5, taken from late stages (t = 28 kyr = 0.59t_ff) and large scales (the cavity reaches the radius of the cloud). Again, the purple directed contour delimits the region where the magnetic field lines are wound by the rotational drag (i.e., where $\left|B_{\varphi }\right|>\left|B_p\right|$). We observe that the magnetic field lines are mostly toroidal in most parts of the cavity at large scales, with the exception of the regions close to the rotation axis and also inside of the Alfvén surface, where the magnetic field is mostly poloidal.

The z-component of the momentum equation, considering a steady-state flow, reads $$\nabla \cdot \left(\rho {\upsilon }_z\upsilon -\frac{1}{4\pi }{B}_{z}B\right)+{\partial }_z{P}_t=\rho a_{gz},$$(7)

where the divergence ∇·ƒ of a each vector field f(R, z) is ∇·ƒ = (1/R)∂_R(Rf_R) + ∂_zf_z. Given that the magnetic field lines are strongly toroidal in the tower flow, we can neglect the contribution from magnetic tension, which only involves the poloidal components. The low thermal pressures at large scales and in the cavity allow us to consider only the contribution of the magnetic pressure to the total pressure gradient. With these approximations, Eq. (7) becomes $$\nabla \cdot \left(\rho {\upsilon }_z\upsilon \right)=\rho a_{gz}-\frac{1}{8\pi }{\partial }_z{B}^2\cdot$$(8)

The z-component of the gravitational force is negative in the upper hemisphere. In general, the magnetic field strength increases toward the center of the cloud because it is being dragged there by gravity, and so, the z-component of the pressure gradient is generally negative in the same region. Therefore an outward momentum flow is obtained if $$\left|\frac{1}{8\pi }{\partial }_z{B}^2\right|>\left|\rho a_{gz}\right|.$$(9)

The condition of Eq. (9) is calculated and shown in Fig. 5 as the nonshadowed part of the orange directed contour. That directed contour was computed with the spherical-radial components of both quantities for simplicity, given that z ~ r for θ → 0. Indeed, in the tower flow, the magnetic pressure gradient is larger than gravity, and since B is toroidal in the same region, we conclude that the tower flow is driven by the winding of magnetic field lines by rotation.

The speeds obtained in the tower flow are slower than the magneto-centrifugally driven jet. Typical values are on the order of 10 km s⁻¹. The morphology and velocity structure of the tower flow can vary significantly with the initial conditions chosen for the cloud core. This is explored in more detail in Sect. 4.

Fig. 5 Magnetic tower flow.

Fig. 6 Magnetic broadening of the cavity due to the tower flow.

3.5 Magnetic broadening of the outflow cavity

We observe (Fig. 6) that the tower flow increases its thickness over time. As the accretion disk grows in size with time (cf. Oliva & Kuiper 2023), so does the region of the cloud where the magnetic field lines are wound by rotation. Magnetic pressure keeps increasing in the cavity, strengthening the tower flow and broadening it as a consequence. Additionally, the slow depletion of the envelope means that the ram pressure from large scales in the cloud core decreases over time, allowing further growth of the cavity. The broadening of the tower flow offers an additional component for the earliest outflow broadening (before the driving of HII regions) around early O and B-type stars (see Beuther & Shepherd 2005). According to our results, very collimated outflows are indicative of the earliest stages of star formation (fiducial case t ~ 10 kyr or 0.2t_ff). A first broadening comes as a result of the broadening of the tower flow (t ≳ 20 kyr), and we expect a second broadening due to the growth of the HII region (t ≳ 100 kyr, see Kuiper & Hosokawa 2018) before other forms of radiative feedback (e.g., stellar winds and radiation pressure) take over.

3.6 Jet recollimation

At the launching region, the jet is accelerated along the spherical-radial direction, which means that the highest speeds are not reached along the rotation axis. However, at scales of a few thousand astronomical units along the propagation of the jet, and fully outside of the Alfvén surface, we find evidence that magnetic hoop stress and the ram pressure from the infalling envelope exert forces that recollimate the flow and yield high radial velocities at the rotation axis. In Fig. 7, we observe those recollimation points marked as a, b, and c, where the cavity appears narrower. The left panel of Fig. 7 presents the positive radial velocity, which peaks particularly strongly around the recollimation points a and b. The right-hand panel of the figure provides a visualization of the forces involved in the recollimation. Similarly to the previous analysis, the region where the magnetic field lines are wound is shown in the shadowed area of the purple directed contour. Additionally, the green directed contour displays the regions where the cylindrical-radial component of the Lorentz force density $$F_{\text{BR}}={\left(J\times B\right)}_R/c=\frac{1}{4\pi }{\left[\left(\nabla \times B\right)\times B\right]}_R$$(10)

is negative (i.e., directed inward; J is the current density) and in magnitude larger than the cylindrical-radial component of the centrifugal force density F_cR, that is, the regions where the Lorentz force tends to bring the flow closer to the rotation axis. The overlapping of the green and purple contours shows that the inwardly directed Lorentz force constitutes an example of magnetic hoop stress.

Recollimation in massive protostellar outflows has been observed in the form of jet knots in, for example, IRAS 21078+5211 (Moscadelli et al. 2021) and the radio jet HH 80–81 (Rodríguez-Kamenetzky et al. 2017). In those systems, nonthermal synchrotron emission originated in the knots is observed, where the particles are accelerated until they reach relativistic speeds. We find velocities on the order of ~1000 km s⁻¹ in the jet material, however, the highest speed of the ejecta cannot be determined from our results because of the Alfvén limiter we set up in order to be able to avoid extremely small time steps. Even though we cannot fully recover the kinematics that produce the nonthermal emission observed in jet knots, the recollimation mechanism we observe provides an explanation and a model for the production of the knots themselves. For a concrete example, the recollimation point a in Fig. 7, located at around 1000 au away from the massive protostar, and where such high speeds are seen in the simulation data, coincides reasonably well in position with the nonthermal emission lobe (JVLA observations at 5 cm) observed in IRAS 21078+5211 (Moscadelli et al. 2021, Fig. 5). The recollimation of the outflow, which restricts the width of the cavity, is also visible in the shape of the radio lobes at scales of ~10 000 au, taken from SO molecular transitions which trace external shocks of the flow (e.g., coming from the cavity wall; see Moscadelli et al. 2021).

Three-dimensional effects could potentially intervene in the production of jet knots and the propagation of the outflow. These effects should be studied with three-dimensional simulations of a disk-jet system, where the assumption of axisymmetry is relaxed. Material from fragments formed in the accretion disk could potentially be ejected along the jet cavity, producing jet knots or bow shocks (however, this has not been reported in three-dimensional simulations in the literature; see Sect. 5). Additionally, episodic accretion might cause small precession in the jet axis, as suggested by the different regions of the red and blueshift in the ejected material at large and small scales in IRAS 21078+5211 (Moscadelli et al. 2021, 2022).

Fig. 7 Recollimation of the jet. The color scale for the radial velocity is saturated at υr = 500 km s−1. For a fuller picture of the speeds reached over time, we refer the reader to Fig. 10.

Fig. 8 Production of the bow shock that gives rise to the jet cavity.

3.7 Bow shock and the formation of the cavity

Given that the magneto-centrifugal mechanism requires the existence of the cavity, the natural question to ask is how is this cavity formed in the first place. The formation of the cavity is captured in Fig. 8. The disk starts to appear shortly before the instant of time depicted (where it is only a few astronomical units in size). In the thick disk, the magnetic field is dragged by rotation and gravity, and magnetic pressure increases rapidly, until it eventually overcomes gravity and the outflow is launched. This is evidenced with the fact that the magnetic field in the cavity (excluding the region surrounding the rotation axis) is toroidal, and that the magnetic pressure gradient is larger in magnitude than gravity in the spherical radial direction. The velocities achieved in this stage are on the order of 10 km s⁻¹, so they are in line with the values observed in the tower flow stage. The formation of the cavity produces a bow shock, marked in the figure with the label a, which propagates outward as the cavity expands over time. At t ~ 5.4 kyr in the fiducial simulation for the × 16 grid, the cavity is broad enough such that the magneto-centrifugal mechanism is able to start from the cavity wall, as described in Sect. 3.1. Then, the velocities of the jet increase to more than 100 km s⁻¹.

Bow shocks are often observed in protostellar outflows (e.g., Moscadelli et al. 2021; Lefloch et al. 2015). Even though we expect bow shocks to be originated by multiple causes (for example, episodic accretion), the formation of the cavity provides a production mechanism for bow shocks that are highly symmetrical and appear distant in the ejected material. Moreover, the position and propagation velocity of the bow shock could serve as a crude estimate for the age of a given protostellar disk-jet system, as we suggested in Moscadelli et al. (2022). In that particular case, we estimated from the fiducial simulation on grid ×16 that the propagation of the bow shock is of about 3300 au kyr⁻¹, and with this, that when the protostar is ≈5.24 M_⊙, the bow shock has propagated at around ~30 000 au. This coincides reasonably well with the observed position of the bow shock at ≈36 000 au away from the forming massive star.

Fig. 9 Effects of magnetic braking on the magnetically driven outflows.

3.8 Magnetic braking: Effects on the outflows

The rotational dragging of magnetic field lines creates magnetic tension, which tends to locally remove angular momentum very efficiently from the small scales in a process called magnetic braking (Galli et al. 2006). We discussed magnetic braking in the disk region in Paper I, and showed that its effects are stronger as time passes and the magnetic field lines are progressively more wound. When enough angular momentum is removed from the innermost parts of the disk and cavity wall (which in our fiducial case happens at t ~ 15 kyr), those regions lose centrifugal support and the flow becomes mostly infalling. Figure 9 shows an extensive morphological and dynamical analysis of the situation. The blue directed contour (gravito-centrifugal equilibrium) reveals the lack of centrifugal support for the inner parts of the disk (mark a) and the cavity wall (mark b), which corresponds to the magnetically braked flow region.

The magnetic field is mostly toroidal in the thin and thick layers of the disk, however, in the magnetically braked region of the disk it becomes mostly poloidal, following the flow. As time progresses, magnetic tension continues to transport angular momentum outward from the braking radius of the disk until the full depletion of angular momentum in the braked region. When this happens, sections in the low-density cavity are observed to counter-rotate (with the corresponding sign reversal in the toroidal magnetic field; see black directed contour in Fig. 9). The counter-rotation is intermittent in the early stages. Counter-rotation in jets has previously been studied as a theoretical possibility (Sauty et al. 2012; Cayatte et al. 2014; Staff et al. 2015), and some observational evidence has been gathered (see, e.g., Coffey et al. 2011; Louvet et al. 2016) although it is not conclusive (Tabone et al. 2020).

The removal of angular momentum from the innermost scales interferes with the magneto-centrifugal mechanism as discussed in Sect. 3.1. The center and right panels of Fig. 2 show that the lack of angular momentum in the inner parts of the thick disk and cavity wall (marker b) impedes the magneto-centrifugal launching of material despite being in the sub-Alfvénic regime. As a consequence, the cavity becomes narrower close to the protostar. An outflow driven by magnetic pressure (marker c) can still be launched from the innermost parts of the cavity (cf. orange contour, that shows dominance of magnetic pressure over gravity in the vertical direction). Despite the loss of angular momentum in the innermost scales, the magneto-centrifugal mechanism and cavity wall ejections are able to continue from higher points in the disk (marker d), where the conditions of centrifugal support and sub-Alfvenicity are still satisfied. The tower flow, on the other hand, remains unaffected as the rotation of the disk continues to wind the magnetic field lines outside of the Alfvén surface, and therefore, magnetic pressure continues to be generated over time.

In this work, we did not include the effects of ambipolar diffusion and the Hall effect (see also the discussion in Sect. 4.2). In particular, we expect ambipolar diffusion to have the following effects. Because ambipolar diffusivities increase with decreasing density (Marchand et al. 2022) gravitational collapse would be promoted specially at its earliest stages. However, at the number densities found in the accretion disk (10⁹ ≲ n ≲ 10¹² cm⁻³), the diffusivity curves from Marchand et al. (2022) show that the ambipolar and Hall effect diffusivities behave similarly to Ohmic resistivity. Because of this, we would not expect a strong impact in the disk physics, although additional diffusivities might cause the thick layer of the disk to be thinner because of the decreased magnetic pressure. From Paper I, we expect higher diffusivities to delay the effects of magnetic braking (see also Masson et al. 2016 for an example study in the low-mass case where the interplay between ambipolar diffusion and magnetic braking is present).

Inside of the cavity, however, the situation is different. Even though ambipolar diffusivities increase with low densities and strong magnetic fields (such as inside the cavity), all magnetic diffusivities (Ohmic resistivity, ambipolar diffusivity and the Hall diffusivity) strongly drop when the material becomes ionized (Marchand et al. 2022; Machida et al. 2007). The ejected material from collimated outflows from massive (proto)stars has been found to be shock-ionized (Moscadelli et al. 2021; Rodríguez-Kamenetzky et al. 2017), so we expect our treatment of the jet physics in the cavity to be correct, especially at early times when magnetic braking effects are not yet present.

As the main focus of this article is the study of the mechanisms that mediate the production, propagation and termination of outflows, we point out that magnetic braking is nevertheless a mechanism capable of terminating a jet. Magnetic braking could then be the mechanism responsible for setting the timescale of collimated protostellar jets, and regulating the expansion of HII regions, although this conclusion must be confirmed by the addition of more physical effects to the model.

Fig. 10 Morphology and kinematics of the outflows for different initial conditions, using the grid ×8. The velocities shown as a function of time correspond to the maximum radial velocity within an angle of 22.5° from the rotation axis, and the logarithmic scale has a cutoff of ±1 (velocities in the interval −1 ≤ υr ≤ 1 km s−1 are shown as zero).

4 Variation of the outflow properties with initial conditions

In this section, we utilize the analysis from Sect. 3 and the results from other simulations in our series in order to investigate how the physical processes of the outflows are impacted by the conditions of the natal environment of the protostar. For comparisons between simulations, we use three kinds of graphics, which we present next.

Each panel of Fig. 10 displays the morphology of the outflows at early and late times (with a corresponding adjustment to the spatial scale) accompanied by a velocity profile as a function of time. The velocity profile is taken along a cone of 22.5° around the rotation axis, using the maximum value of the radial velocity at a given radial distance, and it permits us to study the kinematical signature of the jet and its propagation through the cloud core. In order to produce an effective morphological comparison, we used only the results for the high resolution grid (×8) because coarser grids in our series are not able to fully resolve all the scales of interest (see Appendix A), although the basic physical processes are still present.

Figure 11 exhibits a comprehensive overview of the temporal evolution of the dynamical processes discussed in Sect. 3, namely, the formation of the cavity (leading to a bow shock), the magneto-centrifugal mechanism, magnetic braking and the magnetic tower flow. The time of launching of the bow shock is simply computed as when υ_r > 0 for the first time at a reference altitude of z = 100 au. The interval of presence of the magneto-centrifugal jet requires additionally that (a_c)_R > (a_g)_R and ${\upsilon }_p<{\upsilon }_{Ap}$ at the same altitude. Magnetic braking, in the other hand, is graphed by using the magnetic braking radius computed in Paper I. Finally, the presence of the tower flow is fixed by simultaneously enforcing the conditions $\left|B_{\varphi }\right|>\left|B_p\right|$, ∇_rP_B > (ρa_g)_r, and υ_r > 0 at an altitude of z = 5000 au.

Finally, Fig. 12 aims to aid in the understanding of the effects of the jet in their environment with time, using only the data from the high-resolution series ×8. We calculate the momentum flux (panel c) of the jet through a hemisphere of radius 15 000 au, which gives an idea of the mechanical feedback that the magnetically driven outflows impart to their natal environment at larger scales. In terms of mass, panel d presents the accretion-to-ejection rate for the fiducial case, computed using the mass flux through spheres at different scales. The time average of the accretion-to-ejection rate at large scales in the cloud core remains within 10%. Both quantities are characterized by a strong variability over time, due to the episodic ejection events (specially at late times) caused by the cavity wall ejections; for this reason, a time average is necessary to study their overall behavior. A smoother quantity is the total radially-outward momentum contained in a hemisphere (i.e., the total momentum expelled by one branch of the jet), and it is presented in panels a and c.

Fig. 11 Comparison of the processes present in the jet for different initial magnetic field strengths, rotation profiles, rotational energies and resistivity models. For all the simulations except MC = 150 MC [×8], the free-fall timescale is 52.4kyr. The gray boxes denote that the simulation was not run until 30 kyr.

Fig. 12 Mechanical feedback from the outflow. Panel a: momentum ejected by the magnetically driven outflows within the cloud core, for the high-resolution series on grid ×8. Panel b: time-averaged momentum flux through a hemisphere of radius r = 15 000 au for the same simulation series. Panel c: total momentum ejected by the magnetic outflows, normalized to the mass of the cloud core, as a function of the fraction of the free-fall timescale elapsed, for two different masses of the cloud core. Panel d: ejection-to-accretion rate measured using spheres of three different radii.

4.1 Mass of the cloud core

Morphologically, the fiducial case (M_C = 100M_⊙) and the higher mass cloud core (M_C = 150 M_⊙) are very similar, except that the cavity is “pinched” in the former case compared to the latter, which indicates that the jet recollimation does not occur in the same place for both simulations. The velocity profile along the rotation axis shows different stages of evolution of the system in the fiducial case: the initial gravitational infall before t ≈ 5 kyr, the formation of the cavity and the launching of the magneto-centrifugal jet (velocities on the order of 100 km s⁻¹ after t ≈ 5 kyr). The velocity profile does not show the acceleration of the jet on the smallest scales because it occurs radially outward (Sect. 3.2) and therefore it is not enclosed in the conical region used to compute the profile. However, the profile does capture the recollimation of the jet at scales of 10²–10³ au (Sect. 3.6), where the velocity reaches values on the order of 1000 km s⁻¹.

All the dynamical mechanisms related to the outflows (Fig. 11) are present in both simulations, although they start earlier in time for the more massive case. The same can be seen when examining the total linear momentum carried out by the outflows (Fig. 12a), with the additional remark that the mechanical feedback from the magnetic outflows is stronger for more massive cloud cores. In Paper I, we concluded that the disk dynamics scale with the mass of the cloud core, when examining the wider range of cloud core masses of 50–200 M_⊙ using the lower resolution grid ×4. Motivated by that result, we scaled the curves of the higher resolution results for the linear momentum by dividing by the mass of the cloud and expressed the time as a fraction of the free-fall timescale (Fig. 12c). We find that the momentum scales well with the mass of the cloud, and so, despite the nonscalable nature of the thermodynamics and magnetic diffusivities of the system and its self-gravity, the results reported in this article can be reasonably used for cloud cores in a wider interval of masses, possibly including low-mass star formation.

Taking into account the respective free-fall timescales for the two clouds (52.4 kyr for M_C = 100 M_⊙ and 42.8 kyr for M_C = 150 M_⊙), we corroborate the scalability of the outflow dynamics as follows. The bow shock and subsequent magneto-centrifugal jet start at t = 5.2 kyr = 0.1 t_ff for M_C = 100 M_⊙, which corresponds to t = 4.2 kyr for M_C = 150 M_⊙ and is in agreement with the respective simulation result. Similarly, the tower flow reaches the height of z = 5000 au at t = 0.18t_ff in the fiducial case and at t = 0.16t_ff in the more massive cloud core. We had checked the scalability of magnetic braking with the fraction of free-fall timescale in Paper I already.

The curve for the outflow total linear momentum is explained as follows. At t ~ 0.1t_ff, the magneto-centrifugal jet is launched, and the total linear momentum increases as the ejecta propagates throughout the cloud core. The outflow reaches a distance of 0.1 pc at around 0.2t_ff. Once this happens, the total momentum reaches a maximum value, with variations caused by changes in the recollimation zone and episodic ejection events (cavity wall ejections). We see from Fig. 10 that the high velocities associated with recollimation and jet launching become less frequent at late stages, because of the effects of magnetic braking. This, together with the slower tower flow (which broadens the cavity at late stages), cause the total linear momentum to decline slightly over time for t > 0.35t_ff.

4.2 Magnetic field

Our parameter space for the magnetic field responds to two different assumptions on the environmental conditions for massive star formation. In low-mass star formation, supercriticality has been observed in prestellar cores not yet under gravitational collapse, with measured magnetic field strengths on the order of a few tens of microgauss (see, e.g., Ching et al. 2022; Crutcher 2012), which mean mass-to-flux ratios on the order of a few. In the case of massive star formation, the precollapse conditions of the magnetic field are not well-known, and because of this, we assume cases of lower and higher magnetization. The lower magnetization case assumes that the same order of magnitude of the magnetic field strength observed in low-mass star formation is also present in the massive case and it constitutes our fiducial case. The higher magnetization case assumes that a similar mass-to-flux ratio as observed in the low-mass star formation would be present in the high-mass scenario.

In what follows, we compare the outcome of the fiducial run (normalized mass-to-flux ratio $\overline{\mu }=20$, which corresponds to an initial magnetic field of 68 µG) and a simulation with a stronger initial magnetic field ($\overline{\mu }=5$, corresponding to 0.27 mG). Stronger magnetic fields cause stronger and earlier magnetic braking as already discussed in Paper I, and the additional support from magnetic pressure allows for the formation of a larger accretion disk. The outflow cavity is narrower at early times (Fig. 10) with strong magnetic fields, which can be attributed to the additional magnetic pressure which makes the observed thick layer of the disk thicker, constraining the cavity close to the protostar.

In Sect. 3.7, we presented our finding that the outflow cavity is initially formed by magnetic pressure, launching a bow shock in the process. A comparison of the points and bars that indicate the bow shock and the magneto-centrifugal phase in Fig. 11 confirms that, with more magnetic pressure, the cavity is able to form earlier in time, followed by the jet. However, the magneto-centrifugal jet is shorter-lived in the case with $\overline{\mu }=5$, (cf. its velocity profile, Fig. 10). Magnetic braking for a cloud core with $\overline{\mu }=5$ acts earlier in time compared to when $\overline{\mu }=20$. The extraction of angular momentum in the inner disk becomes very strong after t ≳ 15 kyr, which impedes the magneto-centrifugal mechanism to continue. Additionally, the magnetic braking starts to act in the thick disk as well (but not in the thin disk, see Paper I) turning the flow into sub-Keplerian and the tower flow cannot form because the magnetic pressure derived from the rotational dragging of magnetic field lines is not enough to overcome gravity. The magnetic pressure originated from the dragging of magnetic field lines by gravity, on the other hand, makes the thick disk thicker. We conclude from those results that very strong magnetic fields have the power to terminate magnetically driven outflows. Because of the lower jet speeds, the momentum output (Fig. 12a) of the outflows is considerably weaker in the case of stronger initial magnetic fields.

Because ambipolar diffusivities increase with magnetic field strength, we expect the higher magnetization case to be the most affected by the noninclusion of ambipolar diffusion in the present study. However, we remark that we still see the formation of a disk-jet system, and that magnetic braking, dependent on the initial magnetic field strength, would set a timescale for the formation and duration of the magneto-centrifugal phase of the outflows. For example, for low-mass star formation (where the $\overline{\mu }$ is on the order of a few) the jet duration (as a fraction of the freefall time) could be shorter. The interplay between jet termination, magnetic braking and angular momentum transport should be examined in future studies.

4.3 Angular momentum content

We performed simulations of various ratios of rotational to gravitational energy (ζ = 0.01, 0.04) and initial distributions of angular velocity $\text{Ω}\propto R^{{\beta }_{\text{Ω}}}$, where β_Ω = 0 (solid body rotation) and β_Ω = β_ρ/2 = −0.75 (steep initial rotation profile, as we considered in our simulations of a fragmenting accretion disk in Oliva & Kuiper 2020). The fiducial case is β_Ω = 0, ζ = 0.04. In the high-resolution series, we only investigated the case where the rotational-to-gravitational energy ratio was kept the same but the angular momentum is initially more concentrated in the center of the cloud (steep rotation profile instead of the solid body assumption). In Paper I, we found that the accretion disk becomes larger thanks to the greater initial availability of angular momentum close to the protostar, and that the thick disk is flatter. As a result of a larger accretion disk for β_Ω = −0.75, we find an initially wider tower flow (Fig. 10), which translates into a more cylindrical outflow in both small and large scales compared to solid body rotation.

Both the velocity profile in Fig. 10 and the total linear momentum (Fig. 12a) indicate that the propagation speed of the bow shock and the jet material is slower for the steep angular velocity profile compared to solid body rotation, although the recollimation speeds are similar in both cases. From Paper I, we know that the mass of the protostar does not increase the same way in both cases. For a steep initial angular momentum profile, the early formation of the accretion disk slows down the accretion rate onto the protostar. As an example, at t = 10 kyr, M_⋆ = 3.1 M_⊙ for the fiducial case, but only 1.77 M_⊙ for β_Ω = −0.75. Even though the magnetic field lines are dragged more by rotation with more availability of angular momentum, and therefore magnetic pressure should increase, the lower radial mass flux onto the protostar translates into a reduced dragging of magnetic field lines by gravity, which means that magnetic pressure actually decreases and the propagation of the outflow becomes slower. This difference can also be seen in the formation of the accretion disk: while in the fiducial case the magnetic outflows and the disk are formed almost simultaneously, in the case with β_Ω = −0.75, the jet is launched with a delay with respect to the formation of the accretion disk (cf. Fig. 11 and Fig. 8 of Paper I). The differences in protostellar mass are also responsible for the delay in magnetic braking (see Paper I).

Overall, the total outward momentum reaches similar values for both cases. After examining the full parameter space including the results for the lower resolution series ×4, we find that even though the initial angular momentum can cause differences in the launching and propagation of the jet, all the dynamical processes seen in the fiducial case are also seen in cloud cores with a lower ratio of rotational to gravitational energy.

4.4 Ohmic resistivity

As a final exploration, we look at the effect of resistivity in the dynamical processes present in the outflows (Fig. 11). We found in Paper I that the presence of resistivity delays magnetic braking compared to the ideal MHD case (the delay is more noticeable when comparing against the protostellar mass as opposed to time). Apart from our standard case, that considers resistivity only as a function of density, we included the full resistivity formula by Machida et al. (2007), which also depends on temperature and can cause an increase of a factor of 10–100 compared to our standard case. Additionally, we included a simulation where Ohmic resistivity is artificially high by a factor of 10. These four simulations provide insight to what happens when the resistivity model is varied, that is to say when more magnetic diffusivity is added. From Fig. 11, we see that the jet and the tower flow are launched slightly later for the case with higher resistivity. This is expected, because the jet launching mechanism depends on flux freezing inside of the Alfvén surface and the dragging of magnetic field lines by rotation. However, Fig. 11 shows the presence of all processes in all cases.

5 Comparison to previous numerical studies

In Paper I, we presented a comparison of our results for the accretion disk dynamics with several studies in the literature. Here we continue the comparison of our results for the dynamics of the magnetic outflows against the results found in the same studies.

Kölligan & Kuiper (2018) performed simulations of a forming massive star starting from the collapse of a cloud core with magnetic fields, but considering an isothermal equation of state (in contrast to our treatment of radiation transport with the flux-limited diffusion approximation). We base our fiducial case on the setup proposed in that study, except that they use a rotational to gravitational energy ratio of 2% (our fiducial case is of 4%). The study by Kölligan & Kuiper (2018) was the first to be able to clearly distinguish a fast, magneto-centrifugally launched jet from the slower magnetic tower flow in the context of massive star formation. They found evidence of the magneto-centrifugal launching of the jet, however, they applied a series of analytical criteria from the classic bead-on-a-wire jet theory from Blandford & Payne (1982; for which the effects of the infalling envelope are missing) to identify the jet launching area both at early and late times (when it moves upward). Our explicit identification of the Alfvén surface and the gravito-centrifugal contour (cf. Fig. 2) enabled us to identify the launching region more accurately for all times and determine the effects of the ram pressure of the infalling envelope on the geometry of jet launching, which cause departures from the classic theory. Kölligan & Kuiper (2018) performed a convergence study using grids with increasing resolution, but their highest resolution grid is equivalent to our ×2 grid, which means that in this study we are able to resolve more structures in the plasma. The insufficient resolution of the outer regions in the cloud core in their convergence study meant that the authors were not able to find convergence in several quantities, including the jet duration and the momentum transported by the magnetic outflows. In contrast, we find outflow momenta that tend to converge for the two highest resolution grids (×8 and ×16), and the minimum cell size required for the adequate resolution of the tower flow at large scales (see Appendix A.1). In Kölligan & Kuiper (2018), the artificial mass generated by the Alfvén limiter grew to unrealistic limits for some simulations. The higher resolution grids used in our study and fine tuning of the Alfvén limiter allow us to keep the artificial mass generated by the limiter to negligible levels (cf. Appendix A.2).

5.1 Studies with ideal MHD

The study conducted by Banerjee & Pudritz (2007) consisted in 3D ideal MHD simulations of a slowly rotating Bonnor-Ebert sphere of 168 M_⊙ on 1.62 pc a combination of cooling and radiation diffusion for the gas thermodynamics and using the code FLASH with an adaptive mesh refinement (AMR) grid. Despite only obtaining a sub-Keplerian disk (probably due to strong magnetic braking due to the lack of magnetic diffusivity), they observe early pseudo-disk-driven magnetic winds of typical velocities of 4 km s⁻¹, whose launching is credited to the tower flow mechanism.

A setup more comparable to ours is used in Seifried et al. (2012; 100 M_⊙ within a radius of 0.125 pc and initial solid body rotation), where the effects of different mass-to-flux ratios was investigated. They find outflow velocities of ~5 km s⁻¹, which are both very collimated (weakly magnetized runs) and poorly collimated (strong magnetization). The authors point to the existence of a “fast” outflow component launched by the magneto-centrifugal mechanism, which is detected by a launching criterion based on the energy equation (cf. Eq. (4)). However, the velocities of the fast outflow are still on the order of 10 km s⁻¹, which we credit to the coarse grid used (3D AMR grid with a minimum cell size of 4.7 au) and the unconstrained effects that magnetic braking causes when no magnetic diffusivity is considered. They also observe an outflow driven by magnetic pressure, although the authors do not call it a tower flow, even though the whole cavity is dominated by the toroidal component of the magnetic field. The ejection-to-accretion rates at large scales are somewhat larger than our findings (they find Ṁ_ej/Ṁ_acc ≲ 0.4).

The study of Myers et al. (2013) considers a denser cloud core (300 M_⊙ within a radius of 0.1 pc) and $\overline{\mu }=2,$, but supersonic turbulence instead of rotation. They include a gray flux-limited diffusion approximation to radiation transport. Even though the authors report on the existence of magnetically driven outflows of up to 40 km s⁻¹, they do not discuss their origin. A similar approach was taken by Rosen & Krumholz (2020), but with masses of the cloud core of 150 M_⊙. However, the coarse spacial grid used (minimum cell size of 20 au) impedes them to resolve the launching mechanisms of the magnetic outflows, adopting a subgrid outflow injection approach instead. This means that we cannot directly compare our results with theirs.

5.2 Studies including Ohmic resistivity

The three-dimensional studies by Matsushita et al. (2017) and Machida & Hosokawa (2020) use the same Ohmic resistivity model by Machida et al. (2007) as we do. However, while we use treatment of radiation transport based on the flux-limited diffusion approximation, both studies model the thermodynamics of the dust and gas with a barotropic equation of state. Matsushita et al. (2017) start from a Bonnor-Ebert sphere of several masses that range from 32 to 1542 M_⊙ and radius 0.28 pc which rotate as a solid body and that are threaded with an initially uniform magnetic field determined by normalized mass-to-flux ratios of $\overline{\mu }=2,5$ and 10. The authors report on outflows of magnetic origin associated with the outer disk, but do not discuss their origin in detail. The speeds obtained in the outflows are on the order of 10 km s⁻¹. Machida & Hosokawa (2020) start from a similar setup, but expand the parameter space of the cloud core masses (from 11 to 525 M_⊙) and normalized mass-to-flux ratios (from 2 to 20). Additionally, they introduce an enhancement factor to the density profile to promote gravitational collapse and their cloud core radius is 0.2 pc. The authors classify their outflows into successful, delayed and failed. They observe protostellar outflows failing to evolve and collapse by the strong ram pressure when a massive initial cloud is weakly magnetized; a similar preliminary finding was reported in Matsushita et al. (2017) with the conclusion that only strongly magnetized cloud cores can drive protostellar outflows. This trend is opposite to our findings: we see that a relatively weak initial magnetic field $\left(\overline{\mu }=20\right)$ is sufficient to produce a magneto-centrifugally launched jet, while strong magnetic fields launch shorted-lived jets because of the increased action of magnetic braking over time. In Matsushita et al. (2017) and Machida & Hosokawa (2020), the origin of the outflows is not thoroughly discussed, but given the magnitudes of the velocities obtained (a few tens of kilometers per second), it is reasonable to think that they are launched by magnetic pressure and not the magneto-centrifugal mechanism. The propagation of the failed outflows reported in Machida & Hosokawa (2020) are reported to be stopped by ram pressure from the infall, a direct consequence of the interplay between gravity and magnetic forces. We have seen both in this article and in Paper I that the initial density and angular momentum distribution are crucial to determine the features of the accretion disk and the magnetically driven outflows. For this reason, we speculate that the differences in the initial cloud mass distribution (they consider higher cloud masses which means stronger ram pressure) and the launching mechanisms for the outflows could be responsible for the contradictory conclusions obtained in our study and Matsushita et al. (2017) and Machida & Hosokawa (2020).

5.3 Studies including ambipolar diffusion

Finally, we turn our attention to studies that include the effects of ambipolar diffusion but no Ohmic dissipation. Mignon-Risse et al. (2021) found magnetically driven outflows of velocities on the order of 30 km s⁻¹ in simulations that start from a cloud core of 100 M_⊙ within a radius of 0.2 pc (roughly equivalent to M_C = 50 in our setup), a density profile increasing toward the center of the cloud but with a plateau, solid body rotation, a rotational-to-gravitational energy ratio of 1% and a mass-to-flux ratio of 5. Cloud cores with subsonic and supersonic turbulence were also considered in that study. The authors use the code RAMSES with an AMR grid of maximum resolution of 5 au, including radiation transport using the flux-limited diffusion approximation and gray irradiation from the protostar(s). Commerçon et al. (2022) start from a similar setup, but consider additional configurations with a rotational-to-gravitational energy ratio of 5% and a normalized mass-to-flux ratio of 2. The characteristics of the outflows obtained in both studies correspond mostly to our wide tower flow rather than the fast jet. In Mignon-Risse et al. (2021), the authors determine that the origin of the outflows is indeed magnetic (as opposed to produced by irradiation from the protostar), and that they are driven by the magnetic pressure gradient. The authors argue that the magneto-centrifugal mechanism might be responsible for some of the high velocity flow present in their simulations, but do not perform an in-depth confirmation of this possibility. We show that very high resolution is needed close to the massive protostar in order to obtain a fast jet. On the other hand, their consideration of ambipolar diffusion might mean that the early stages of gravitational collapse and the late stages of the outflow dynamics are better represented in their simulations for low values of the mass-to-flux ratio.

5.4 Sink particles and the launching region

When comparing our results against studies that make use of subgrid sink particle algorithms, such as the ones mentioned in this section with the exception of Kölligan & Kuiper (2018), it is important to take into account that those algorithms usually require the definition of an accretion radius. For example, the accretion radius in Commerçon et al. (2022) is four times the minimum cell size, that is to say ~20 au. This means that the dynamics of the gas are only expected to be fully consistent with the system of equations beyond that radius. The use of a sink cell of 3 au in radius and a logarithmically spaced grid permitted us to resolve the launching region of the magneto-centrifugal jet in a more consistent way in the regions very close to the protostar. Additionally, the definition of an accretion radius when using sink particles affects any comparisons of the effects of magnetic braking, because we observe magnetic braking only in the innermost tens of astronomical units in the disk, within the size of the accretion radius of several of the studies above.

The minimum cell sizes of the simulations run in grids ×8 and ×16, of 0.06 and 0.03 au respectively, have enabled us to directly compare our results to the VLBI observations of IRAS 21078+5211, which have a resolution of 0.05 au (Moscadelli et al. 2022). This level of detail, in both observations and simulations, had not been obtained so far. The recent advances in observational techniques for studying protostellar jets evidences the need for future numerical studies to accurately describe the processes in the innermost ~ 100 au from the forming massive star.

6 Summary and conclusions

We have paved the path toward a more complete theory of magnetically driven massive protostellar outflows by performing simulations that self-consistently generate them, by analyzing their dynamics and evolution over time, and by relating the conditions of the onset of gravitational collapse with their production and propagation into larger scales. We performed a set of 31 axisymmetric magnetohydrodynamical simulations, including Ohmic dissipation as a nonideal effect, and radiation transport of the thermal emission and absorption of the dust and gas. We summarize our findings as follows:

• The outflow cavity is originally formed by the increase in magnetic pressure as a product of the dragging of the magnetic field by the gravitational collapse and rotation. In the process, a bow shock is thrust outward. Once the cavity is formed and the Alfvén velocity increases inside, the magneto-centrifugal mechanism starts.

• The fast (≳100 km s⁻¹) and highly collimated component of the outflow is launched by the magneto-centrifugal mechanism. Ram pressure from the infalling envelope and the thick layer of the accretion disk constrain the geometry of the launching region into a narrower configuration compared to the predictions of classical jet theory.

• The jet material is accelerated along the magnetic field lines while its flow is sub-Alfvénic. Once it leaves the Alfvén surface, it enters a helical trajectory along the cavity.

• We observe episodic ejection events originating in the innermost parts of the cavity wall that are in contact with the thick layer of the accretion disk.

• Magnetic pressure originating from the rotational dragging of magnetic field lines by the accretion disk drives a slower (~10 km s⁻¹) magnetic tower flow that broadens over time.

• Magnetic hoop stress and the ram pressure from the infalling envelope recollimate the jet in regions that are, in principle, compatible with jet knots observed in protostellar jets, where the material reaches relativistic velocities and nonthermal radiation is produced.

• Magnetic braking modifies the innermost regions of the disk and the jet, narrowing the cavity, and it is capable of terminating the magnetically driven outflows.

• Our results indicate that the ejected momentum scales with the mass of the initial cloud core and its free-fall timescale.

• The morphology and dynamics of the magnetically driven outflows are affected by the initial conditions for the gravitational collapse. More centrally concentrated distributions of angular momentum, which produce a larger disk, drive a more cylindrical jet. Stronger initial magnetic fields launch jets and produce narrower jets that last shorter periods of time.

By assuming axisymmetry and taking maximum advantage of a grid in spherical coordinates that is able to resolve the phenomena of interest on both large and small scales, we were able to obtain unprecedented detail in the inner ~100 au around the protostar. This region is uncovered by a new generation of observational evidence (Carrasco-González et al. 2021; Moscadelli et al. 2021, 2022), which brings new challenges for theoretical predictions on massive protostellar outflows.

Acknowledgements

We thank Richard Nies for his contributions to the analysis of part of the dataset at the early stages of the project. A.O. acknowledges financial support from the Deutscher Akademischer Austauschdienst (DAAD), under the program Research Grants – Doctoral Projects in Germany, and complementary financial support for the completion of the Doctoral degree by the University of Costa Rica, as part of their scholarship program for postgraduate studies in foreign institutions. R.K. acknowledges financial support via the Emmy Noether and Heisenberg Research Grants funded by the German Research Foundation (DFG) under grant no. KU 2849/3 and 2849/9.

Appendix A Convergence of the results

A part of the aims of this study is to determine numerical parameter setups that produce consistent results for the disk and outflow dynamics. We present here a selection of the results in search for numerical convergence.

Appendix A.1 Resolution

From a comparison between the results of Kölligan & Kuiper (2018) and studies that use AMR Cartesian grids, it was clear that resolution of the jet launching region is extremely important for capturing the magneto-centrifugal mechanism and obtaining velocities on the order of ~ 1000 km s⁻¹. In the resolution study we performed, however, we observed several other critical phenomena to be properly resolved by spatial grid.

A proper study of the propagation of the outflow is only possible if the large scales of the cavity are properly resolved. Figure A.1 shows a comparison of the outflow cavity at large scales for all the grids in our investigation. Grid ×16 produces an outflow cavity that is well resolved at large scales. The cavity contains the tower flow at large scales, and its thickness changes smoothly with distance. The results for grid ×8 show most of the features of grid × 16, although the thickness of the cavity does not change as smoothly as in simulation ×16 for early times, even though it does for late times when the tower flow broadens. The outer parts of the cavity are resolved by ~ 6 grid cells. For the grid ×4, geometry of the cavity changes qualitatively. In smaller scales, we observe the flow to be restricted to lobes, while in outer scales the cavity is narrower. A closer look to the outflow cavity for late stages (where the flow is dominated by the tower flow) is offered in Fig. A.2. Beyond z ~ 6 000 au, the cavity is not properly resolved (only 2-3 grid cells) and the toroidal nature of the magnetic field at large scales cannot be seen. The same is true for grids ×2 (where the lobes are even smaller) and ×1 (where the propagation of the outflow itself is not properly represented any more.

The total outflow momentum (Fig. A.3) is strongly hit by the resolution of the outflow cavity, but also by the minimum density resolved in the simulation. A higher resolution grid is able to produce lower densities as fine structures in the cavity are distinguished instead of being represented by their average density. Lower densities produce a weaker linear momentum. The curves for the total momentum in grids ×4, ×8 and ×16 seem to agree in order of magnitude, which is attributable to a better agreement in the densities and velocities obtained. However, the curves also show a peak at early times only present in simulations × 1, ×2 and ×4, which may be due to the improper resolution of the outflow cavity at large scales.

The narrowing of the cavity wall with the effects of magnetic braking was discussed in Sect. 3.8. We note that, when this width is not properly resolved, the launching point of the magneto-centrifugal mechanism quickly moves outward over time until the jet is completely terminated. This is the reason why in Fig. 11, we only show the results for grids ×4 and ×8, where the cavity is partially and properly resolved at all times, respectively.

Fig. A.1 Convergence of the jet cavity with resolution.

Appendix A.2 Alfvén limiter

As discussed in Paper I, we use a varying density floor such that the Alfvén speed of a given grid cell does not increase beyond a fixed limit. This Alfvén limiter creates artificial mass as a result. We performed a parameter scan with different values of the maximum allowed Alfvén velocity with the following values: ${\upsilon }_{\max }^A=200$, 500, 1000, 2000 and 5000 km s⁻¹. We monitored the total artificial mass generated in each case, and established a maximum acceptable value of 1 M_⊙, that is, 1% of the mass of the cloud core, for the low-resolution exploratory runs. For grid ×2, we observed that this limit was surpassed only for ${\upsilon }_{\max }^A=200$ and 500 km s⁻¹. For this reason, we chose ${\upsilon }_{\max }^A=2000\text{\hspace{0.17em}km\hspace{0.17em}}{\text{s}}^{-1}$ as the fiducial value for grids × 1, ×2 and ×4. In grids ×8 and × 16, however, another parameter scan revealed that ${\upsilon }_{\max }^A=1000\text{\hspace{0.17em}km\hspace{0.17em}}{\text{s}}^{-1}$ yielded sufficiently low artificial masses (well below 0.005 M_⊙), while saving computational costs, and therefore this value was chosen for the two highest-resolution runs. As a consequence, the artificial mass created by the Alfvén limiter has a negligible impact in the case of the high resolution runs we used for the main parameter scan in this study.

The Alfvén limiter is mostly active in the launching region of the jet (close to the massive protostar) and around the recolli-mation regions. Therefore, we cannot conclude on the maximum velocities observed in a protostellar jet. This limitation of our setup, which is based on classical physics, leaves theoretical room for the relativistic velocities required for the nonthermal emission observed in recollimation regions in protostellar jets (see, e.g., Moscadelli et al. 2022 and Rodríguez-Kamenetzky et al. 2017).

Fig. A.2 Convergence of the tower flow with resolution at large scales.

Appendix A.3 Sink cell size

Finally, we wanted to check that our choice of a sink cell size does not fundamentally interfere with the physical processes described in this article. To this end, we performed a parameter scan with three values of the radius of the sink cell: 1 au, 3 au (our fiducial case) and 10 au. Due to computational costs, we only display the results for the first 15 kyr of evolution in Fig. A.4. We find that the magneto-centrifugal launching of the jet and the formation of the bow shock occur earlier with a smaller sink cell. However, all the processes described in Sect. 3 are present in all the simulations, and the launching of the jet seems to happen still around t ~ 5 kyr for the two smallest sink cell sizes.

Fig. A.3 Convergence of the total jet momentum with resolution.

Fig. A.4 Convergence of the observed dynamical processes of the magnetic outflows as a function of the size of the sink cell.

References

All Figures

	Fig. 1 Schematic overview of the parameter space covered by the simulation series.
In the text

	Fig. 2 Magneto-centrifugal mechanism for the fiducial case in ×8.
In the text

	Fig. 3 Acceleration of material along the magnetic field lines.
In the text

	Fig. 4 Cavity wall ejection.
In the text

	Fig. 5 Magnetic tower flow.
In the text

	Fig. 6 Magnetic broadening of the cavity due to the tower flow.
In the text

	Fig. 7 Recollimation of the jet. The color scale for the radial velocity is saturated at υr = 500 km s−1. For a fuller picture of the speeds reached over time, we refer the reader to Fig. 10.
In the text

	Fig. 8 Production of the bow shock that gives rise to the jet cavity.
In the text

	Fig. 9 Effects of magnetic braking on the magnetically driven outflows.
In the text

	Fig. 10 Morphology and kinematics of the outflows for different initial conditions, using the grid ×8. The velocities shown as a function of time correspond to the maximum radial velocity within an angle of 22.5° from the rotation axis, and the logarithmic scale has a cutoff of ±1 (velocities in the interval −1 ≤ υr ≤ 1 km s−1 are shown as zero).
In the text

	Fig. 11 Comparison of the processes present in the jet for different initial magnetic field strengths, rotation profiles, rotational energies and resistivity models. For all the simulations except MC = 150 MC [×8], the free-fall timescale is 52.4kyr. The gray boxes denote that the simulation was not run until 30 kyr.
In the text

Fig. 12 Mechanical feedback from the outflow. Panel a: momentum ejected by the magnetically driven outflows within the cloud core, for the high-resolution series on grid ×8. Panel b: time-averaged momentum flux through a hemisphere of radius r = 15 000 au for the same simulation series. Panel c: total momentum ejected by the magnetic outflows, normalized to the mass of the cloud core, as a function of the fraction of the free-fall timescale elapsed, for two different masses of the cloud core. Panel d: ejection-to-accretion rate measured using spheres of three different radii.

In the text

	Fig. A.1 Convergence of the jet cavity with resolution.
In the text

	Fig. A.2 Convergence of the tower flow with resolution at large scales.
In the text

	Fig. A.3 Convergence of the total jet momentum with resolution.
In the text

	Fig. A.4 Convergence of the observed dynamical processes of the magnetic outflows as a function of the size of the sink cell.
In the text