A Unified Model for Bipolar Outflows from Young Stars: Apparent Magnetic Jet Acceleration

We explore a new, efficient mechanism that can power toroidally magnetized jets up to two to three times their original terminal velocity after they enter a self-similar phase of magnetic acceleration. Underneath the elongated outflow lobe formed by a magnetized bubble, a wide-angle free wind, through the interplay with its ambient toroid, is compressed and accelerated around its axial jet. The extremely magnetic bubble can inflate over its original size, depending on the initial Alfv\'en Mach number $M_A$ of the launched flow. The shape-independent slope $\partial{}v_r/\partial{}r=2/3t$ is a salient feature of the self-similarity in the acceleration phase. Peculiar kinematic signatures are observable in the position--velocity (PV) diagrams and can combine with other morphological signatures as probes for the density-collimated jets arising in toroidally dominated magnetized winds. The apparent second acceleration is powered by the decrease of the toroidal magnetic field but operates far beyond the scales of the primary magnetocentrifugal launch region and the free asymptotic terminal state. Rich implications may connect the jets arising from the youngest protostellar outflows such as HH 211 and HH 212 and similar systems with parsec-scale jets across the mass and evolutionary spectra.


INTRODUCTION
Astrophysical jets are ubiquitous in accretion systems, from young stars to supermassive black holes (e.g., Ray & Ferreira 2021;Romero 2021). They are energetic and highly collimated and remain so over large physical scales and dynamic ranges. While the detailed physics of their processes may change from one to another, jets display similarities in their general characteristic features. Although the dynamical mechanisms for their large-scale acceleration and collimation remain enigmatic and unsettled, the generality of the phenomena and fundamental physical processes from the birth to death of stars across the mass spectrum is recognized.
signatures Shang et al. 1998), from the innermost edge of the protostellar disk, by interacting with the stellar magnetosphere (Shu et al. 1994a;. Standard theories on the formation and propagation of astrophysical jets consider the initial launch and evolution of freely propagating flows in accretion-ejection systems (e.g., McKee & Ostriker 2007;Pudritz et al. 2012;Ray & Ferreira 2021). In the context of protostellar jets, the initial launch refers to the magnetocentrifugal BP mechanism from the original radial location in the disk, followed by passing through the critical points of sonic, Alfvénic, and fast-magnetosonic transitions before reaching the asymptotic state and collimation of streamlines (e.g., Königl & Pudritz 2000;Shu et al. 2000). The flow velocities reach the terminal values determined by the conversion of the rotation and magnetic energy obtained from accretion into the kinetic energy of the gas (e.g., Blandford & Payne 1982;Shu et al. 1994a;Spruit 1996). Beyond the fast-magnetosonic transition, in the asymptotic regime, Shu et al. (1995) demonstrate that the initial gradual acceleration continues to reach its ultimate terminal velocity and flow collimation logarithmically provided the physical scale is large enough.
Protostellar jets are intimately connected to largescale outflows, the presence of which is the outcome of the underlying driver interacting with the ambient environment (e.g., Bally 2016;Bachiller 1996;Lee 2020). Highly collimated outflows surrounding narrow molecular jets with extremely high velocity (EHV) are commonly associated with deeply embedded Class 0 protostars. The influence and feedback of the ambient environment on the overall system need to be included as an integral part of the theoretical treatment, especially for the earliest phase of star formation. Shang et al. (2020, 2023, hereafter Paper I andPaper II) advance the physics of such interplay for a jet-bearing hydromagnetic wind interacting with its ambient envelope. They conclude that outflows formed as wind-driven bubbles should be ubiquitous and are an inevitable integrative outcome of the interplay. Within the magnetized bubbles, the compressed and shocked regions can be thick and extended, and nested multicavities can form naturally as part of the process. Recent high-angularresolution and -sensitivity observations of nested-shell structures in several young outflows such as HH 212, DG Tau B, and HH 30 indeed support the formation mechanism of interconnected multicavities.
Within the integrative platform, a new secondary regime of magnetic acceleration of the flow can be identified, which appears to be a prominent feature arising from the magnetic interplay within confined shocked zones. The unique acceleration mechanism explored in this work differs from the primary one that launches and accelerates from the base region of the free magnetocentrifugal wind. We emphasize this specialized form of acceleration is most salient inside the context of very strongly magnetized self-similar elongated bubbles that take the functional form adopted.
We focus on this new phenomenon in this Letter because of its importance and generality in the context of our current framework. We highlight the ubiquity of energetic jets from astrophysical systems, their formation theory, and where our work will improve the theoretical picture. In Section 2, we review the theoretical framework advanced in Paper I and the formulation established for this series of works. We summarize the definitions of relevant variables and derive characteristic properties of the acceleration. In Section 3, we demonstrate the numerical examples and their consistency with the analytic expectations. We discuss and summarize the significance and implications of the findings in Sections 4 and 5.

THEORETICAL CONSIDERATIONS
We advance an integrated framework of outflows unifying the jet, the wide-angle wind, and a magnetized nonspherical (spindle-shaped) bubble elongated along the jet axis in Paper I. Contrary to the expectations of a thin-shell model with self-similar expansion (e.g., Shu et al. 1991;Koo & McKee 1992) and latitude-dependent velocity variations, the bubbles formed in this integrated framework possess extended internal fine structures that arise due to the magnetic interplay of the magnetized wind and the ambient medium.
We demonstrated in Paper I the physics of the presence of complex structures when a hydromagnetic wind interacts with an ambient toroid. These outflows are the interaction between an ambient medium and an advanced state of an asymptotic wind. The hydromagnetic wind can be simplified as a free constantvelocity radial flow with a strong density stratification of 1/ 2 . Through the interaction with the ambient medium, nested structures form surrounding the free wind, first enclosed by a reverse shock as an inner cavity, followed by an extended region filled with compressed wind, then covered by a layer of compressed ambient medium and compressed ambient poloidal field. Between the compressed wind and the ambient medium lies the tangential discontinuity, which is subject to substantial shear and unstable to Kelvin-Helmholtz instabilities (KHI). When the wind is toroidally magnetized, the compressed and shocked regions can become very extended and nonuniform due to the nonlinear growth of the KHI, further enhanced by vorticity generated by magnetic forces. The magnetic forces feed back into the flow velocity, resulting in magnetic pseudopulses in the extended compressed-wind region and giving an impression of bow shocks along the jet axis.
In the current framework, we find time-dependent structures arising in a compressed-wind region confined by shocks, which are self-similar with respect to the general features of outflow morphology and main kinematics. The elongated bubble establishes its self-similarity prior to the onset of the second acceleration through decay in the pressure of a toroidal magnetic field. This is different from a hypothetical magnetic tower mechanism, which would have to operate from the base of the launching process near the disk. In this work, we focus on the time-dependent self-similar conditions of the confined wind, which amplify the acceleration given by the decay of the magnetic energy term. We also study the increase of outflow size in interaction conditions out of a steady state for dynamical expansion.
This phenomenon of secondary acceleration occurs beyond the region of compressed wind with complex fine structures of KHI and magnetic pseudopulses, but before the terminal forward shock, within the zone bounded by the reverse and forward shocks. This zone of new physics, where the jet accelerates further beyond the terminal velocity, overlaps with the "tip" region in Paper I. In Figure 1, we illustrate our schematic configuration for the second acceleration. In this figure, we focus on the regime of physics discussed specifically in this work, different from those considered in Papers I and II, for the magnetically dominated bubbles (1 < M A 6) and their structures. We will discuss where this newly observed magnetic acceleration occurs in the configurations and their implications.
The foundation and the associated computational and analytical methodologies utilized are shown in Paper I, and their kinematic features in Paper II. For consistency, we follow the same methodology as previous works in the following subsections.

The Magnetized Wind
We consider a magnetocentrifugally driven wind, dominated by the wound-up toroidal component well beyond the transition region to reach its final asymptotic state, such as one considered in Shu et al. (1995) (for a comparison, see the end of this subsection). The toroidal velocity of the wind declines steadily at large distances due to the conservation of angular momentum. The free wind reaches approximately a poloidal flow with a purely toroidal field frozen inside, b p = 0 and v φ = 0, with mass density ρ = ρ wind and gas pressure p = a 2 wind ρ.
Within the spherical coordinates in axisymmetry, the velocity (v r , v θ , v φ ) and magnetic field (b r , b θ , b φ ), a simplified set of MHD equations was derived in Paper I. A summary of variables and symbols is given in Table  1.
In a sufficiently cold wind, steady-state equations can be written as Equation (24)-(27) in Section 2.2 of Paper I, which admit solutions that give these three constants v r , ρ 2 , and |b φ | for p = 0. These properties give the cylindrically stratified density profile ρ ∝ −2 , an exactly radial flow of velocity v = v rr , and magnetic field b = b φφ , for the freely propagating wind before it encounters the ambient medium.
For an axisymmetric b φ magnetic component, a current function C = − b φ , proportional to the total current carried in one hemisphere, can be defined and helps in computing the poloidal current density j p and magnetic force f c ( Table 1). The force term C 2 /2 behaves similarly to a pressure f c = −(ρ 2 ) −1 ∇ C 2 /2 . In a cold free wind, the value of C stays exactly constant, and relevant gradients vanish.
In this work, for mathematical convenience, we definẽ ρ ≡ ρ 2 ,b ≡ b φ = −C, and their ratio independent of for a free wind. This ratio µ symbolically defines a relative ratio of the magnetic field to mass density.
The steady-state free wind takes the form of constants of the problem: where D 0 , b 0 , and v 0 are constants of density, magnetic field, and velocity, respectively, defined at the inner radial boundary of the problem. Their numerical values are given in Table 1. The Alfvén speed v A = b φ /ρ 1/2 =b/ρ 1/2 is constant in the free-wind region and equal to its inner radial BC value v A0 ≡ b 0 /D 1/2 0 . In other regions of the flow, downstream from the free wind, the value of v A can and does differ from the BC and free-wind value v A0 . The free wind has an initial Alfvén Mach number M A = v 0 /v A0 , for the wind velocity v 0 and the initial Alfvén speed in the wind v A0 : Additionally, and are set exactly in the boundary condition. These values of M A , v 0 , wind density ρ w , and D 0 fix the scale and value of the toroidal magnetic field in the wind. Within the cold free wind, C = C 0 and µ = µ 0 . This overall character applies to strongly magnetized wind launched from a sufficiently close vicinity to the innermost portion of the disk, irrespective of the origin of its magnetic field, provided that the proper jet density and b φ can maintain their profiles for the features constructed and remain sufficiently radial v r v θ . The launch region also needs to be much smaller than the physical scales within which the very elongated bubble has reached its self-similar phase (Paper I; Paper II).
The initial free-wind flow in this work is assumed to have reached the superfast asymptotic regime in steady state before interacting with an ambient medium and becoming enclosed by the reverse shock. The free asymptotic steady-state X-wind studied by Shu et al. (1995) continues the magnetocentrifugal acceleration process beyond its fast point, considering carefully the inner θ boundary conditions and the expected magnetic behavior, including the slow logarithmic decay of magnetic energy toward more considerable distances. In its native format, Shu et al. (1995) provide an asymptotic description in which magnetocentrifugal effects proceed together with magnetic torque along each wind streamline (summarized in Equation (29) of Paper I through C X β X , where C X decays as 1/ log(r), despite C X sharing a similar role to the C-function adopted here). When such a wind is used as a BC to study outflows, it is equivalent to cutting off the magnetocentrifugal process at a specific scale corresponding to a particular value of C X . The Alfvén Mach number M A X is defined equivalently in Equation (30) of Paper I, showing the effect of C X β X , the magnetic Maxwell torque per unit of mass flux. The regime of 1 < M A 6 in this work is magnetically dominated, as explained there, and the C X β X term is relatively large, allowing a smaller M A . The wind magnetization labeled as M A in this series reflects a parameterized set of wind configurations in space and time. The resulting magnetically dominated self-similar bubble driven by a highly magnetized wind is out of steady state and dynamically expanding for our purposes.

Equations of Self-Similarity
We write down a simplified form of equations for the free-and compressed-wind regions, ignoring terms due to gas pressure, b p , v θ , and v φ . Under those approxima-tions, the equations of time evolution for v r , ρ, and b φ simplify to Equations for density and magnetic field can be written (usingρ,b, and D t Experimentally observed in our numerical results reported previously (Paper I; Paper II), self-similarity applies to the majority of the flow, either free or compressed, except for the exact KHI features at any instant of time.
Here we note the property that for any quantity A = A(x) behaving (to sufficient approximation) selfsimilarly as a function of x ≡ r/v 0 t, its total derivative in a region of self-similar acceleration can be estimated as while in steady state D t = v r ∂ r . Equation (11) also defines the velocity v E ≡ v r − r/t, pervasive in timedependent self-similar flows needing to consider total derivatives. Both steady-state and self-similar formulas agree in a steady-state region such that r v r t and thus v r ≈ v E . It is the case for the wind in the (usually small) innermost cavity enclosed by the reverse shock. The observation of the tip region shown in Paper I as discussed in Appendix B is now connected to this full set of selfsimilar equations.
In Appendix C we examine these self-similar equations and obtain an important property, Equation (C11), linking the magnetic quantity µb =b 2 /ρ = b 2 φ /ρ = v 2 A to the hydrodynamical results. That equation can be fulfilled in an uninteresting manner, ∂ r v r = 0, as in the free wind, or more interestingly as or v 2 A = v 2 E as most probably can happen in an acceleration region. Equation (12) gives an algebraic functional form for the behavior of b φ in the acceleration region. It is, however, not to be utilized in a free-wind region.

The Acceleration
The whole process of acceleration is under the control of the radial force, Equation (6). The acceleration originates in the magnetic force term which converts toroidal magnetic field to radial velocity and can be rewritten as The accelerating effects of this force are amplified by the ∂ t v r term of Equation (6). The magnetic force term in the acceleration region can be simplified to µ∂ rb ≈ ∂ r (µb), the first term on the right-hand side of Equation (13), by noticing that in that region, µ is nearly constant, This motivates the definition of a function behaving in the acceleration region as a specific enthalpy function for the magnetic force, which helps define a Bernoulli constant in Equation (18) below. The existence of a function such that −∇h = j × b/ρ is special. When present, the generation of vorticity by magnetic fields ∇ × (j × b/ρ) becomes zero (curl of a gradient), consistent with the acceleration region but not with the pseudopulses (Paper II). It may be illustrative to compare this magnetic specific enthalpy function with that of an ideal gas. Here we have a specific enthalpy µb = b 2 φ /ρ = v 2 A equal to twice the magnetic pressure per unit mass b 2 φ /2ρ, and equal to the square of the Alfvén speed. This relation between specific enthalpy, pressure, and wave speed could compare symbolically with that present in an ideal gas with pressure p, sound speed (γp/ρ) 1/2 , and specific enthalpy (p/ρ)[γ/(γ − 1)] with γ = 2. This analogy (short of a complete equivalence) works for the radial force component −∂ r h and its relation with gas pressure and wave (Alfvén) speed. In Appendix F we have extended it into a simplified model for the radial force, able to reproduce some of the results regarding radial velocity and position behavior of the acceleration region.
The momentum equation can be fulfilled if either v E = 0 (only at the end of the acceleration region after h has been spent), or the gradient of velocity takes the value independent of M A . In self-similar variables, the acceleration, is a constant. This functional form for the gradient of v r , ∂v r /∂r, could be "experimentally" checked against the strongly magnetized wind regimes for the phenomenological context below.

Bernoulli Constant and the Maximum Velocity
Inside the acceleration region, a Bernoulli function can be defined by combining three terms where the last term is obtained by using Equation (16) to integrate ∂ t v r ≈ −(r/t)∂ r v r . This quantity is therefore a Bernoulli constant (∂ r H = 0) within the acceleration region, beginning at r i and ending at r f . At the beginning of the acceleration region, and at the end because the magnetic strength has been spent to zero at the end of the acceleration region. The initial specific magnetic enthalpy h i , is calculated at r i , where h i = µ ibi , at the intersection extrapolated down through the 2/3 slope line down to the v r /v 0 = 1 line (shown in Figure 2 below in Section 3.2; for a detailed derivation, see Appendix D.).
As a simple convenient approximation, however, we Equation (D12) implies that v f ≤ 3v 0 ; in combination with Equation (21), this result requires M A ≥ 1. This is expected because sub-Alfvénic flows require a different treatment. The accelerated self-similar v f and r f found in Equation (21) can be substantially larger than the hydrodynamic momentum-conserving values, now analytically quantified. We demonstrate this analytic result in our numerical computations.
With regard to the validity of the formulation as presented, the trans-Alfvénic M A = 1 limit is simply symbolic but not realistic because the basic assumption of the wind having reached a superfast asymptotic state has been included in the boundary condition. A full study including wind launching will be implemented into this framework in future work.

NUMERICAL DEMONSTRATIONS
Simulations setups are described in detail in Papers I and II. We summarize key ingredients here for the computational data utilized in this work.

Numerical Results
The parameter space contains wind cases for seven values of M A = 1.2, 1.5, 2, 3, 6, 18, and 30 for v 0 = 100 km s −1 . For each of the free winds, there are four values of n = 1, 2, 4, and 6, and two levels of ambient magnetization α b = 1 and 0, at 6400 × 2016 resolution. Here we check the validity of our derivations of the formulation for the strongly magnetized winds of M A = 1.2 to M A = 30, most relevant to the regime considered in this work. The cases M A = 1.2, 1.5, 2, 3, and M A = 18 are newly performed runs using otherwise the same setup and parameter space as those of M A = 6 and M A = 30 reported in Papers I and II.
We note that self-similarity is established and maintained in a wind-blown bubble as long as the ambient medium [Equation (A2)] and the steady wind [Equation (A1)] follow the same 1/r 2 pattern, as explained, e.g., in Appendix A of Koo & McKee (1992) (see also Shu et al. 1991). This result is basically independent of the θ profile as the r and θ directions are independently separable. In this tip region, the background tapering n = 0 toroid (isothermal sphere) profile dominates over the n > 0 configurations as n becomes bigger.
The exploration of numerical results in this section shows that the acceleration effect is rather invariant and follows the theoretical results of Section 2.

The Velocity Gradient
We derive the velocity gradient in the self-similar acceleration zone in Equation (17). We first check if the functional forms for the gradient ∂v r /∂r match the numerical results. The top panel of Figure 2 shows the collection of cases from all n values at θ = 0. • 0446 near the jet axis for the maximum radial velocity v r . The constant velocity gradient is observed to be followed by all the cases at this same θ angle. Each of the v r curves starts from a horizontal portion full of oscillations from the magnetic pseudopulses, followed by the rise of v r /v 0 with a constant slope of 2/3 to a maximum value before ending just before the terminal forward shock. For cases with wider openings in n = 4 and 6, the smoother acceleration regions are less affected by magnetic pulses introduced by the interplay. For the stronger magnetization with smaller M A values, the acceleration often starts below the v r /v 0 = 1 level directly from a rebound from a noticeable, deep pulse.
Along the M A values shown, the second apparent acceleration has contributed to the increase of velocity from around the 10% level for the intermediate M A range, 18-30, about 1.75-2 for M A = 2-3, and can shoot up to an approximate factor of 2-3 close to M A = 1-1.2. All the cases are below the symbolic limit of v r /v 0 = 3, implied by Equation (21), when M A = 1. This approximate trend of v f /v 0 is shown in the top panel inset of Figure 2.
To inspect how the constant slope holds in the 2D acceleration zone, we sample the radial profiles of all the n and M A cases of interest for the θ angles. Self-similar acceleration is observed and reproduced in the 2D space, not limited to the jet axis, in the bottom panel of Figure  2. The peak velocities reached for different M A values, and θ values are naturally mingled on each of the panels labeled by n, and they all appear to follow the same slope and same trend. The fluctuations in v r /v 0 due to the magnetic pseudopulses are noisier as expected due to the smaller openings in n = 1 and, to some extent, in n = 2. They are virtually nonexistent in n = 4 and n = 6 in the acceleration zone. This indicates that the self-similar condition applied in the derivation is robust in the whole 2D domain of acceleration.

Very Magnetic Bubbles
For the primary expected behavior of b φ dropping to zero near the tip of the acceleration zone, we naturally look for the variableb and normalize to the value b 0 at the boundary. Figure 3 summarizes the situation using the local C normalized by b 0 . The variables −C/b 0 = b/b 0 should be 1 at the inner radial boundary and remain so throughout the uncompressed free-wind region. This quantity traces the strength of |b φ | to very small values. The locations where the transitions into the acceleration zone take place are marked by loci of v r /v 0 = 1.
We observe howb decreases toward the tip of the acceleration region with the value of M A . The acceleration region becomes significantly larger in the extremely magnetically dominated bubbles represented by small M A values, approaching up to z/(v 0 t) ≈ 3 in the axial direction. They demonstrate the characters as magnetic bubbles driven by very strong magnetic pressure. How the decay ofb powers the acceleration at the far end of the compressed-wind region is directly shown in Figure  4 through the variation ofb versus v r /v 0 .
The newly defined specific magnetic enthalpy in Equation (14) is shown in Figure 5, in the form of h/h 0 , where The profiles of µ are shown in Figure 7 and Appendix C, and they indeed vary very little in the smooth regions. Through the Bernoulli function H and ∂ r H = 0, the acceleration starts from v r /v 0 = 1, which often has a higher specific enthalpy than h/h 0 = 1 for most M A 6. The function h = µb is always an important parameter of this flow (equal to the square of the local Alfvén speed). However, in this work, it is a specific magnetic enthalpy function only in the acceleration region and with the condition that µ ≈ constant. Regions of magnetic vorticity generation cannot have a magnetic acceleration term derivable from an enthalpy function.
The color contours ofb trace the magnetic pseudopulses in the compressed wind from right across the reverse shock all the way up to the uppermost gray contour ofb/b 0 = 1. The magnetic pulses form bow-shaped compressedb and give the impression of aligned bow shocks along the jet axis, most evident in M A ∼ 18-30 especially for the smaller n = 1 and n = 2 range. Also worth noting is howb traces the KHI at the interface of the compressed wind and compressed ambient media on the sides.
The contrast of the KHI patterns at the interfaces of interplay suggests how the local conditions support the growth of different modes. However, these patterns all vanish similarly toward the very top portions of the bubbles. The KHI modes are suppressed in the upper acceleration zone despite their different local growth environments. This is consistent with the inability to generate vorticity because the magnetic force can be expressed as the gradient of a quantity [Equation (14)] in the acceleration region, similar to the condition arising in the free-wind zone.
The regions where µ is constant cannot generate magnetic vorticity. The regions where µ varies help to trace the presence of KHI regions with magnetic vorticity generation, enabling pseudopulses. After that region stops, it is possible to enter the acceleration region, with h = µb now behaving as a specific enthalpy function. The decay ofb along the acceleration region, as shown in Figure 3, powers the acceleration. A detailed comparison of the α b = 1 and α b = 0 panels of these figures shows more prominent pseudopulses for α b = 0, especially for larger or smaller n, due to the presence of KHI. This has the effect of delaying acceleration. Horizontal or slightly V-shaped contours of v r /v 0 = 1 for α b = 1 are often replaced for α b = 0 by sharply Vshaped contours for which a larger z value is needed to enter a smooth acceleration.
One salient feature from the extremely magnetized bubble regime (M A 3) is the significantly expanded   . The 2D spatial profiles of log 10 (h/h0), a specific magnetic enthalpy within the acceleration region, for winds of MA = 1.2, 1.5, 2, 3, 6, 18, and 30 with n = 1, 2, 4, and 6 toroids and ambient poloidal field strengths of α b = 1 (left) and α b = 0 (right), in which h0 ≡ v 2 0 /MA 2 . The black contours are loci of vr/v0 = 1 and the brown contours are loci of h/h0 = 1. The spatial axes are labeled in units of v0t, with horizontal axes exaggerated by factors as in Figure 3.
compressed-wind region overtaking almost the entire lobe volume. The whole bubble structure is powered by a wind that is launched from the disk by, e.g., an X-wind or an inner disk wind, a process taking place in the smaller scales. The large acceleration region (like any other part of the compressed-wind region) is powered by the launching region, but it is not part of it, distinct in this regard from a tower model. This is true even for those extremely magnetized (barely superfast) winds, for which most of the lobe volume is occupied by the apparent second acceleration region transferring magnetic to kinetic energy.

DISCUSSIONS
We discuss the theoretical and observational implications of our new findings presented in this Letter. We also discuss its generality, limitation, and applicability to other systems.

Signature of the Apparent Acceleration
The apparent second acceleration is a prominent signature of the bubble driven by the strongly toroidally magnetized wind under the self-similar elongated bubble environment. It makes a clear kinematic feature, especially on a position-velocity diagram or through proper motions.
In the kinematic signatures of the elongated magnetized bubbles explored in Paper II, there indeed appears a new velocity feature in the strongly magnetized systems of M A = 6 in which the jet velocity centroids can further shoot up in position-velocity diagrams made of column densities. This is the acceleration signature that motivates the study of this work. PV systems with a larger M A display smaller tips with fluctuating velocity centroids without shooting up.
This strong jet acceleration appears prominently from the top portion of the compressed-wind region as a linear increase of the velocity centroid if proper excitation of the local density is present. This is a strong effect in contrast to the mild fluctuations around the original ejection velocity caused by the magnetic pulses. In Figure 6 (a) of the PV diagrams of column density, the butterfly-shaped distribution of the column densities stretches out owing to acceleration in both jet portions compared to those of larger M A values near the bottom panels. The velocity peaks bend out toward the larger absolute velocity, and the wiggles of the magnetic pseudopulses fade moving up the z-axis. The bending of the velocity centroids outward starts from lower heights as M A goes lower. The column densities in the acceleration zone selected with mixed material are significantly reduced as shown in Figure 6 (b), indicating that this zone is powered with the primary jet/wind material. This acceleration has an exact slope of 2/3 in our presented figure panels, which are at 45 • . We note it is a welcome coincidence for our selected inclination angle. The viewing angle i changes the velocity along the line of sight by a factor of cos(i), and the position on the plane of the sky together with the proper motion by a factor of sin(i). Both of these factors equal each other and 2 1/2 = 0.7071 for our choice of 45 • , canceling the angle factors in this case (tan(45 • ) = 1) when the gradients are computed. This explains the special angle and slope exhibited by the out-bending jet velocity centroids in Figure 6. Hence, the apparent shifts of jet velocity centroids inferred from the PV diagrams give a direct indication of the magnetic acceleration, corrected for the individual inclination angle of the system (see Section 4.3).
This apparent acceleration at large scale is a feature that increases with the strength of the wind magnetization. This effect is to be observed in appropriate systems that can properly manifest the phenomenon. For extremely magnetized systems, this acceleration can increase v r to ∼ 2-3 times the terminal velocities and expand the outflow lobe to ∼ 2-3 times the originally expected sizes. For extremely strong small M A values, the acceleration zone will dominate most of the jet proper, leaving a very small cavity enclosing the free wind relative to the volume of the outflow bubble.
When the system exhibits such behavior at a moderate 20-60% level, the effect could otherwise be interpreted as a change of jet activities, especially when combined with intrinsic system episodicities caused by sporadic ejections. It is expected that for moderately to weakly magnetized systems, the effect is relatively weak and confined to a tiny tip or entirely nonexistent. Absence or no detection in these systems suggests magnetization in the underlying wind is likely not strong enough for the apparent acceleration to be visible.
In situations where actual pulses are present, as happens in real systems, similar signatures would apply for the formation of the acceleration zone and each of the pulses after the first pulse. The first acceleration zone should be located farthest downstream of the most extended compressed-wind region filled with pseudopulses. If the nonsteady wind follows the same 1/ 2 profiles, it will maintain this self-similarity discussed in Sections 2 and 3. The nonsteady wind will be self-shocked and self-similar because the factors cancel out. The nonsteady nature will not affect the first pulse that determines the outermost acceleration zone; however, each of the pulses can establish its own acceleration zone and its compressed-wind region bounded by each new pair of internal reverse-forward shocks. Despite their similarity Figure 6. (a) PV diagrams of column density (rescaled for log NH contrast) produced by winds with MA = 1.2, 1.5, 2, 3, 6, 18, and 30 (left to right) for the strongly magnetized α b = 1 (top) and nonmagnetized α b = 0 (bottom) toroids of n = 1, 2, 4, and 6 (top to bottom) at an inclination angle of 45 • . The column density is integrated for material with vp > a ambient . The spatial position is convolved with a Gaussian profile of 0.01v0t for smoother appearance as in Paper II. in physics and structures, each new zone will be interrupted by new pulses and new layers of shocked cavities. Within the interrupted acceleration zones, the velocities may not be able to develop to complete growth factors indicative of the underlying M A values. Detailed structures of the dynamically large pulsed ejections deeply connected to accretion or other system perturbations will be explored and reported in future publications.

Implications
We have demonstrated a potentially efficient acceleration mechanism for an extremely toroidally magnetized jet system with an extended intershock region between the reverse and forward shocks. The shocks are oblique locally and part of a very elongated self-similar bubble. This mechanism may be able to operate in similar extended inter-shock structures formed by configurations meeting the conditions outlined.
The mechanism outlined gives a new angle to understand the observed jet activities. Up to now, most variations in jet velocities inferred from proper motions or from the PV diagrams have been interpreted as the episodicity of the ejections, most likely tied to the underlying accretion. Systematic faster proper motions from knots located farther out have been simply interpreted as faster ejection velocities in the long past. Given the physics demonstrated in this work, such systems may be our candidates for manifesting the phenomena of apparent acceleration. However, due to a potential lack of shocking pulses in the top portion of the extended acceleration region, the identification of a proper detectable radiative mechanism may be needed. What systems to search for the signatures shall be our next challenge in future works.
We observe the jets formed and sustained by very strong toroidal magnetic fields remain stable for the large physical scales, despite the presence of KHI at the interface with the ambient medium and within the initial compressed and shocked regions. The top (tip) portion of the acceleration zone, however, appears to be KHIfree and stable. Vorticity seems to be suppressed, which has prohibited the KHI. This finding can help us understand how parsec-long jet systems can remain stable without being disrupted or fragmented by the instabilities.
M87 is known to have a helical field (Pasetto et al. 2021), which is consistent with Faraday depolarization at the projected axis and the edges of the conical jet. It also has an extended region with KHI filamentary patterns (Lobanov et al. 2003;Hardee & Eilek 2011), which coincides with the twisting of the toroidal fields (Blandford et al. 2019) within the conical jet. The M87 jet has an extended lobe-like far jet beyond the conical jet in the VLA image in Figure 1 of Pasetto et al. (2021), which strongly suggests the possibility of an underlying bubble structure. The knots in the conical jet and at least the beginning of the far jet (knot A) have been proposed to be produced by KHI, again sharing support for our present framework in the compressed-wind region in the regime dominated by magnetic pressure. Within this Newtonian acceleration picture, partial spending of (e.g.) 3 4 of the available specific magnetic enthalpy (leaving 1 4 of h i ≈ h 0 ) would correspond to an intermediate speed v m /v 0 = 1 + 1/M A < v f /v 0 = 1 + 2/M A , suggestive of a situation of transition between jet regions. The kind of acceleration mechanism and regions formed with the extremely magnetized bubble by toroidal magnetic configuration can be relevant to M87 and perhaps other black hole jet systems dominated by a toroidal magnetic field. The correlation with the measured Faraday depolarization and helical field structure may not be a coincidence.

Protostellar Candidate Systems
Similar systems of very magnetized bubbles relevant to the regime addressed in this work already exist in the literature.

HH 212
HH 212 is recognized as a candidate system exhibiting magnetized bubble features, as demonstrated in Section 8 of Paper I and Section 9.5 of Paper II. It demonstrates the signature multicavities of nested shell-like structures and velocity profiles.
Apparent long-range acceleration has also been observed and reported in HH 212 based on historical proper-motion measurements (Lee et al. 2022), which is otherwise deemed difficult for most traditional magnetocentrifugal-launching models of free winds. Figure 2 of Lee et al. (2022) shows the proper motions of SiO knots from ∼ 100 to 1000 au and much larger distances (up to ∼ 16000 au) along the jet axis (Claussen et al. 1998;Lee et al. 2015;Reipurth et al. 2019). There, an apparent acceleration is seen over the distance from 50 au to 3000 au, with the velocity increasing from ∼ 50-60 km s −1 at 40-50 au to ∼ 140 km s −1 at 3000 au, approximately a factor of ∼ 3 increase over a scale of more than two orders of magnitude (Figure 2 of Lee et al. (2022)). This increase by a factor of ∼ 3 from 50 km s −1 to 140 km s −1 at 3000 au is beyond the increase predicted by the asymptotics in Shu et al. (1995).
HH 212 is not only an excellent testbed of the multicavities as predicted by our framework (Paper I; Paper II) but also a lab for multiepoch large ejections of multiscale nested self-similar bubble structures. For a better understanding, the first segment of the apparent large shell up to the scale of ∼ 1200 au, is followed by an earlier ejected shell of ∼ 16000 au. The historical record of the proper motions should be grouped according to their coeval nested shell structures formed with different histories of bubble formation.
In the recently established shell as part of the selfsimilar bubble of 1200 au, we obtain the velocity along the line of sight (on the order of 10 km s −1 ) from the PV and as proper motion (on the order of 100 km s −1 ) for knots from both sides located at positions around 500-1000 au from Figures 2 and 5 in Lee et al. (2022). For the edge-on system with i ≈ 84 • in HH 212, we obtain the gradients based on proper motions on three accelerating knots (located respectively at y ≈ 590, 710, and 810 au). The gradients are shown to be constant and compatible with a dynamical time t ≈ 20 yr obtained from Equation (16). A simple visual inspection of the PV diagram shows a straight line connecting the velocity centroids, giving the same constant gradient as for the proper motion once we adjust for the viewing angle i with a factor ≈ 10 ≈ tan(84 • ), an intuitive result that is confirmed in the more careful analysis, leading to the same gradient and dynamical time of t ≈ 20 yr.
The segment of three accelerating knots with the constant gradient on the scale of 500-1200 au appears to follow a segment of knots oscillating around a constant centroid of ∼ 55 km s −1 . This region corresponds to the compressed region within the 1200 au bubble cavity, where the apparent acceleration occurs. Out to an even larger scale of ∼ 16000 au, Lee et al. (2022) also find that velocity further away remains roughly constant, which in our interpretation belongs to the compressed-wind region of the bubble cavity established by an earlier event. These large nested cavities appear self-similar to each other, and they belong to bubble cavities formed at different epochs of large ejections. This is consistent with the large-scale pseudopulses filling the previous generation of compressed-wind region outward to the previous terminal shock (see the discussion in Paper II), and that real large pulses by ejection events establish a new bubble nested within the previous one. We will follow up with further theoretical explorations in a later work.

HH 211
The Class 0 system HH 211 is one of the youngest sources driving a highly collimated molecular jet and an outflow (Gueth & Guilloteau 1999;Hirano et al. 2006;Palau et al. 2006;Lee et al. 2007). It is also the only known Class 0 jet surrounded by a toroidal configuration of a magnetic field established by SiO line polarization (Lee et al. 2018) at the time of writing. Hence, it is a desirable testbed candidate for our theoretical configurations with a confirmed magnetic field structure.
The signature of this apparent second acceleration can be observed and extracted from existing published data from the Submillimeter Array (SMA). In the PV diagrams of CO and SiO in Hirano et al. (2006), Palau et al. (2006), and Lee et al. (2007), the redshifted emission traces a shift of the velocity centroid increasing from 25 km s −1 to 35 km s −1 beyond ∼ 10 in SiO J = 5 -4, SiO J = 8 -7, CO J = 2 -1, and CO J = 3 -2. In Figure 3 of Hirano et al. (2006), the shift is continuous and follows a straight line in SiO J = 5 -4. In Figure 3 of Palau et al. (2006) and Figure 5 of Lee et al. (2007), SiO J = 8 -7 also traces the same shifted velocity centroid. CO J = 3 -2 in Figure 5 of Lee et al. (2007) traces a straight line following the trend in SiO J = 8 -7 and SiO J = 5 -4. These emission knots are identified between the SiO and H 2 knots RK5 to RK7 that are ∼ 3. 5 apart, with a rough increase of ∼ 4.5 km s −1 . Combined with the inclination angle i ≈ 80 • (Hirano et al. 2006) and a distance of 321 pc (Ortiz-León et al. 2018), an apparent acceleration is clearly shown to have a constant slope, compatible with the gradient of Equation (16) with a dynamical time of approximately 110 yr. We await further analysis of ALMA data in multiple SiO and CO transitions in future works.

Radio Continuum Jets
Radio jets from more massive protostellar systems show increasing velocities measured via proper motions of knots going away from the source. In recent observations, several massive protostars show indications of jets and wide-angle wind cavities. This suggests that some similar jet-driving mechanisms may operate deep down at the base of these sources. Rodríguez-Kamenetzky et al. (2022) reported a resolved deeply embedded collimation zone ( 100 au) at ∼ 15 au resolution of an intermediate-mass Class 0 protostar (∼ 3 M and 100 L ; Hull et al. 2016) in Serpens, known as the SMM1. Its radio jet consists of a central elongated thermal source (C) and two external lobes (NW and SE, 3000 au apart) (e.g., Rodríguez et al. 1989;Curiel et al. 1993;Rodríguez-Kamenetzky et al. 2016). NW and SE have been known to travel at tangential velocities 200-300 km s −1 in opposite directions (Rodríguez et al. 1989;Rodríguez-Kamenetzky et al. 2016), and the internal knots tend to have velocities ∼ 100 km s −1 . Periodic velocity variations due to tidal interactions, a binary companion with an eccentric orbit, or a past FU-Ori outburst event have been proposed as causes of such velocity trend (e.g., Rodríguez-Kamenetzky et al. 2022). However, a factor of 2 ve-locity increase can be observed on each side separately: SE N/(S1 & S2), and NW/(N3 & N4), using the information on the knots compiled in Table 3 of Rodríguez-Kamenetzky et al. (2022). This implies a highly magnetized flow with an equivalent M A ∼ 2 that can be inferred using the formulation derived in Section 2.4 with Equation (21). SMM1 is one of the few Class 0 sources with radio jets and has molecular cavities. It also has nested shell structures from CO J = 2 -1 (Hull et al. 2016;Tychoniec et al. 2019Tychoniec et al. , 2021, in widening cavity layers and decreasing velocity components (EHV, fast, slow) surrounding the radio jet observed by e-MERLIN (Rodríguez-Kamenetzky et al. 2022). This is consistent with the kinematic and morphological picture of a wind-driven elongated magnetized bubble developed in Paper I and Paper II. On the physical scale of 100 au, a collimation zone is resolved with an ionized stream at 60 au within a narrow cavity of ∼ 28 au near the source. This structure is likely the RS cavity enclosing the primary free Xwind (or an X-wind-like inner disk wind) launched from ∼ 0.4 au. Our current work adds additional information on velocity ratios based on the second acceleration and an inferred strength of wind magnetization of M A ∼ 2 from our framework. The higher luminosity of 100L for 3 M possibly helps with the ionization and illumination of the radio knots. Similar systems would help reveal more massive protostellar systems and their connections to magnetized jets through the proper-motion measurements of knots located far from the embedded sources through a different observational approach using molecular lines.

SUMMARY
We advance a new mechanism of magnetic acceleration within a self-similar highly elongated bubble driven by a strongly toroidally magnetized super-fast wind resulting from the magnetic interplay with the ambient medium. This bubble features a pinched tip along the jet axis generated by a cylindrically stratified density and magnetic field profiles.
This dynamical acceleration operates far beyond the initial launch and acceleration of the usual primary magnetocentrifugal mechanism and the free asymptotic state. The system enters a nonsteady dynamical phase, in which the toroidal wind magnetic field is converted into the increase of radial velocity toward the upper tip of the elongated bubble. The acceleration can maintain a shape-independent constant slope, reaching an M Adependent maximum velocity, determined at the initial ejection. An ideal system will demonstrate the correlation between the decrease of toroidal magnetic field strength and the increase of velocity toward the terminal shock far beyond the reverse shock.
Predicted observational signatures may be identifiable from systematically stronger and faster flow/jet velocity in the past, or otherwise unrecognizable in cases with a weaker jet magnetization. Three protostellar sources are discussed and show features compatible with the present mechanism. The constant velocity gradient of HH 211 and HH 212, a prediction of this mechanism, has been observed in PV for both sources and in proper motion for HH 212. The acceleration to larger velocities (by a factor of 2) in SMM1 and HH 212 can fit the terminal velocity of this mechanism. Additionally, the source M87 has basic physics potentially conducive to this mechanism (toroidally magnetized wind and jet), and we tentatively propose that it can contain a relativistic analog of the presented Newtonian mechanism.
The authors would like to thank the anonymous referee whose comments and suggestions have significantly improved the presentation of this work. The authors acknowledge grant support for the CHARMS group under Theory from the Institute of Astronomy and Astrophysics, Academia Sinica (ASIAA), and the National Science and Technology Council (NSTC) in Tai For the calculations carried out in this work, we adopt numerical values and notations consistent with the numerical simulations in Paper I and Paper II. The formulation, expressed in Lorentz-Heaviside (LH) cgs units for the magnetic field vector b, is consistent with the Zeus-TW code (Krasnopolsky et al. 2010). The constructed free wind ρ w and the tapered toroid ρ a are where R(θ) (a dimensionless density function; and also φ(θ), a dimensionless flux function) is found from the solutions of the toroids (Li & Shu 1996;Allen et al. 2003). The initial magnetic field in the toroids is purely poloidal, and it has the value b = α b b T (r, θ), where α b is a scaling constant, multiplying the toroid field b T computed using φ(θ). The density constants adopted for the wind D 0 , the tapered toroid D S , the density scale factorD, and the toroid magnetic scale α b have the numerical values listed in Table 1 below. A two-temperature equation of state (Wang et al. 2015) is built in for the pressure p and sound speed a, depending on the wind mass fraction f = ρ wind /ρ, and the ambient sound speed a ambient = 0.2 km s −1 and the wind sound speed a wind = 0.6 km s −1 : Most of this article focuses on the cold and essentially unmixed wind regions of the outflow, such that f ≈ 1, a ≈ a wind , and b p and p are minor or negligible terms in the force equations. In Paper I, we noted that a peculiar pointy-nose type of structure can arise in the (small) tip region, that the flow is nearly radial (v r v θ ), and that values of b p are nearly zero. These are similar to the assumptions adopted for the free wind. The transport equations for mass ρ and b φ are similar.
Hence, the radial momentum equation can take a simplified form, ignoring terms in v θ and v φ , as in a nearly radial wind: Ignoring the gas pressure term, the force equation is simplified to Inside the tip region, the magnetic force term is able to smoothly accelerate v r , while b φ smoothly decreases along the z direction. The observation that the tip region follows the expected behavior of the simplified force equation motivates the generalization in Section 2. The behavior of µ is explored in more detail in Appendix C.

C. PROPERTIES OF µ
In this appendix we study in some detail the quantity µ, defined in Equation (1), and utilize it to find some results about magnetism in the self-similar regime of this flow. We first investigate the variations of µ =b/ρ by computing the total derivative quantity D t µ. Using the usual formula for the derivative of a ratio, D t µ =ρ −2 (ρD tb −bD tρ ) . (C6) Using now Equations (9)-(10), the last two terms of which cancel to zero, resulting in D t µ ≈ 0 (C8) in a region for which the approximate Equations (9)-(10) can be applied. Figure 7 shows the variations of µ in space as found in numerical computation. We can find that consistent with Equation (C8) the quantity µ is nearly constant in most of the smooth regions we consider. That includes the self-similar laminar acceleration regions where D t µ ≈ v E ∂ r µ ≈ 0. However, in those parts of the compressed regions where a nonlaminar v θ is present due to interplay and vorticity generation, µ is a rapidly varying function tracing the pseudopulses driven by magnetic forces as shown in Figure 7.
We utilize now these properties of µ as a tool to investigate the variations of the magnetic field within the acceleration region. Substitutingρ andb into and rearranging terms in Equations (6)-(8), and making use of the property of Equation (11), we obtain (v r − r t )∂ r v r + µ∂ rb = 0 (C9) We multiply Equation (C9) by (v r − r/t), and multiply Equation (C10) by µ, reorder the terms, then take the difference of the two resulting equations: This result can be fulfilled in two manners: either by having its second factor ∂ r v r equal to zero as it can happen in the free wind, or by fulfilling Equation (12) in the main text, thus linking magnetic properties in its left-hand side to hydrodynamic properties in its right-hand side, as it can happen in the acceleration regions which are the focus of this work. Figure 7. The 2D spatial profiles of the quantity log 10 (µ/µ0) for winds of MA = 1.2, 1.5, 2, 3, 6, 18, and 30 with n = 1, 2, 4, and 6 toroids and ambient poloidal field strengths of α b = 1 (left) and α b = 0 (right), in which µ0 ≡ v0/MA √ D0. The black contours are loci of vr/v0 = 1, and the brown contours are loci of −C/b0 = 1. The spatial axes are labeled in units of v0t. The horizontal axes have been exaggerated by factors of 3.5, 2.3, 1.4, and 1.0 for n = 1, 2, 4, and 6, respectively. Figure 8. The 2D spatial profiles of number density (nH, the left halves) and column density (NH, the right halves) for winds of MA = 1.2, 1.5, 2, 3, 6, 18, and 30; toroids of n = 1, 2, 4, and 6; and ambient poloidal field strengths of α b = 1 (left) and α b = 0 (right). The spatial axes are labeled in units of v0t, with horizontal axes exaggerated by factors as in Figure 7.