A Unified Model for Bipolar Outflows from Young Stars: Kinematic Signatures of Jets, Winds, and Their Magnetic Interplay with the Ambient Toroids

Kinematic signatures of the jet, winds, multicavities, and episodic shells arising in the unified model of bipolar outflows developed in Shang et al.\ (2020), in which an outflow forms by radially directed, wide-angle toroidally magnetized winds interacting with magnetized isothermal toroids, are extracted in the form of position--velocity diagrams. Elongated outflow lobes, driven by magnetized winds and their interplay with the environment, are dominated by extended bubble structures with mixing layers beyond the conventional thin-shell models. The axial cylindrically stratified density jet carries a broad profile near the base, across the projected velocity of the wide-angle wind, and narrows down along the axis with the collimated flow. The reverse shock encloses the magnetized free wind, forms an innermost cavity, and deflects the flow pattern. Shear, Kelvin--Helmholtz instabilities, and pseudopulses add fine and distinctive features between the jet--shell components, and the fluctuating jet velocities. The broad webbed velocity features connect the extremely high and the low velocities across the multicavities, mimicking nested outflowing slower-wind components. Rings and ovals in the perpendicular cuts trace multicavities at different heights, and the compressed ambient gap regions enrich the low-velocity features with protruding spikes. Our kinematic signatures capture the observed systematics of the high-, intermediate-, and low-velocity components from Class 0 to II jet--outflow systems in molecular and atomic lines. The nested shells observed in HH 212, HH 30, and DG Tau B are naturally explained. Outflows as bubbles are ubiquitous and form an inevitable integrative outcome of the interaction between wind and ambient media.


INTRODUCTION
The study of jets and outflows as signposts of star formation started with the serendipitous discoveries of the Herbig-Haro (HH) objects (Herbig 1951;Haro 1952 (Schwartz 1975). These properties led Schwartz (1978) to hypothesize a stellar source driving a wind (a stellar wind) of 100 km s −1 and a mass-loss rate of ∼ 1 × 10 −5 to 1 × 10 −6 M yr −1 , at least for a brief pe-Historically, a variety of models have been developed to describe the HH objects in terms of the interaction between a stellar supersonic wind and the ambient medium (Schwartz 1978;Norman & Silk 1979;Cantó & Rodríguez 1980;Cantó 1980;Smith 1986;Königl 1982). The framework of Shu et al. (1991, hereafter SRLL) connects the bipolar CO outflows with the momentum input from a (stellar) wind , and the density distribution of the ambient cloud core that gives birth to the star. Conservation laws of mass and momentum give a simple conceptual solution of a self-similar thin shell of the outflows whose shapes are determined by an angle-dependent function. These wind-driven momentum-conserving thinshell models descending from the early stellar wind hypotheses of off-axis HH objects establish the bases of the modern wide-angle wind paradigm of the molecular outflows (e.g., Shu et al. 1991;Masson & Chernin 1992;Li & Shu 1996a;Matzner & McKee 1999;Lee et al. 2001).
Nowadays, HH objects are identified with the jet knots on the axes, forming the basis of the jet-driven bowshock models. The connection between the formation of HH objects and the collimated jets from young stars was established by the linear alignment of HH 46 and 47 (Bok 1978;Schwartz 1977a,b). The linear alignment of HH 47A, HH 46, and HH 47C indicated high collimation in the bipolar flow (Dopita et al. 1982). Mundt & Fried (1983) found emission features of the HH objects resembling those of the HH 46/47 jet, leading to the identification of HH objects as "jets" from young stars. The subsequent discovery of HH 34 as a collimated jet pointing toward a series of bow shocks (Reipurth et al. 1986) indicates the connection of a jet and aligned bow shocks. Many similar collimated jets associated with HH objects from young stars (see the review by Reipurth & Bally 2001) have bow shocks appearing in the highly pointed jet-like momentum output carried in the gas (e.g., Raga & Cabrit 1993;Masson & Chernin 1993;Stone & Norman 1993a,b;Lee et al. 2000;Ostriker et al. 2001;Offner et al. 2011;Offner & Arce 2014;Offner & Chaban 2017).
For the shocked shells formed when a spherical wind runs into an ambient medium, Koo & McKee (1992a,b) introduced the structure and evolution of a spherical wind-blown radiative bubble. Two shocks arise in the interaction of a wind (radiative, adiabatic, or in between) with the ambient medium: a forward (ambient) shock, and a reverse (wind) shock. The shocked wind and the shocked ambient medium are separated by a contact discontinuity. The nature of bubbles adds a new dimension onto the structures formed during the interaction of the media.
In the framework depicted above, Shang et al. (2006) developed a unified model of bipolar outflows integrating the wind-driven and jet-driven paradigms. (Paper I, hereafter Paper II) considered an extension including bubble structures. In Shang et al. (2006), the wind is dominated by the toroidal component of the magnetic field, as expected from the asymptotic regime of an X-wind demonstrated in Shu et al. (1995) and Shang et al. (1998), and a magnetocentrifugal wind arises from the innermost regions of the disk (Li & Shu 1996a;Krasnopolsky et al. 2003;Shang et al. 2007). The outflows form by the interaction of these winds with the ambient medium.
The ambient environment is characterized by a series of magnetized, isothermal "toroids" representing molecular cloud cores before the onset of dynamical collapse (Li & Shu 1996b;Allen et al. 2003a,b;Galli et al. 2006). Analytic formulations (Shu et al. 1991Li & Shu 1996b) can serve good purpose for HD thin-shell approximation (Matzner & McKee 1999) using one toroid density profile (Lee et al. 2000). Full 2D MHD simulations with series of toroid configurations were carried out in Shang et al. (2006), followed by Wang et al. (2015Wang et al. ( , 2019, and Shang et al. (2020), using the Zeus family of codes (Stone & Norman 1992a,b;Krasnopolsky et al. 2010). The morphology correlates well with the systematic evolution of the outflow shapes (Arce & Sargent 2006;Velusamy et al. 2014) from the Class 0 to II stages of evolution.
The self-similar bubble nature of the outflow lobes in Paper I can only be revealed when the numerical resolution is sufficiently high. The fine structures are bounded by reverse and forward shocks, separated by a tangential discontinuity (TD). The wind-toroid interface is unstable to the Kelvin-Helmholtz instability (KHI) due to substantial shear, leading to the formation of extended magnetized mixing layers within the outflow lobes. Magnetic forces in the compressed toroidally magnetized wind can generate vorticity, magnetic pseudopulses, and nonlinear patterns, forming a feedback loop to the growth of structures. These become key morphological features belonging to an elongated magnetized wind-blown bubble, which maintains self-similarity in time.
Molecular outflows carrying both the jet and wind characters in the earliest phases are natural testbeds of the unified theoretical framework. The youngest Class 0 sources with elongated molecular shells surrounding the highly collimated molecular SiO or CO jets are identified in the submillimeter to millimeter wavelengths (nonexhaustively including HH 211, HH 212, L1448C(N), IRAS 04166+2706, CARMA-7, L1157, Serpens SMM1, and Ser-emb 8(N); see, e.g., Jhan & Lee 2021, Lee et al. 2021, Hirano et al. 2010, Wang et al. 2014, Plunkett et al. 2015b, Podio et al. 2016, Tychoniec et al. 2019, and the references therein). These sources usually possess distinct extremely high velocity (EHV) associated with their highly collimated molecular jets, and the expected low-velocity (LV) cavity. In Class I phase, similar phenomena can be identified with bright optical or nearinfrared jets giving the high-velocity component (HVC) surrounded by wide LV molecular cavity walls (e.g. HH 111, HH 46/47, L1551, DG Tau B, and HH 30;Lefloch et al. 2007, Zhang et al. 2019, Stojimirović et al. 2006, Zapata et al. 2015, Louvet et al. 2018. In T Tauri phase, high-velocity (HV) jets or microjets in bright atomic forbidden lines are mostly visible without obvious thick cavities, except for some flat-spectrum sources (Adams et al. 1988;Yokogawa et al. 2001) such as FS Tau B, DG Tau A, and HL Tau (Liu et al. 2012;Agra-Amboage et al. 2014;Takami et al. 2007). Additionally, in some active T Tauri stars, additional mysterious peaks in LV atomic oxygen [O I] have been known to exist (see, e.g., Hartigan et al. 1995;Simon et al. 2016;Banzatti et al. 2019).
High-resolution and high-sensitivity instruments from optical-infrared (such as the Very Large Telescope and Gemini) to submillimeter (such as the Atacama Large Millimeter/submillimeter Array, ALMA) wavelengths have revealed additional finer structures in those very young sources with collimated jet phenomena from Class 0 to II sources. Accumulating evidence of multiple molecular shells organized in nested layers is found in these youngest sources in ALMA observations, which have been described as "onion-like" or as "nested cones", coincidentally, in evolutionary phases spread from Class 0 to II (e.g., Lee et al. 2021;de Valon et al. 2020;Harsono et al. 2021;Güdel et al. 2018). These nested-shell structures also come with accompanying peculiar velocity features distinct from the high-velocity jet or the conventional wide low-velocity cavity. Together with the previously unidentified low-velocity [O I] peaks in T Tauri sources, some trends of systematics appear in the correlated phenomena. Yet a systematic explanation for these connected features has not existed in the literature until this work.
The purpose of the present paper is to document the rich kinematics displayed by the model developed in Paper I, which goes beyond the basic picture expected out of the conventional simple high-velocity jet and lowvelocity shell combinations (e.g. Bally 2016). The patterns of interplay generate their own kinematic signatures differentiable from the smooth bulk flows in the position-velocity diagrams. We document these kine-matic signatures as evidence of outflows being elongated bubbles by figure panels demonstrating the corresponding patterns step by step, for the parameter space explored (e.g. in Secs. 7 and 8).
Our work here offers a self-consistent systematic explanation to the unexplained nested-shell velocity structures as multicavities naturally arising in a magnetized bubble as a result of interplay as exemplified in our framework. Our framework also offers a natural foundation to the ubiquity of molecular outflows as magnetized bubbles and occurrence of their internal structures. These structures arise with their distinct kinematic features. We highlight such connections in Section 9.2. We note that the newly released image of the Class 0 protostar L1527 from the James Webb Space Telescope (JWST) lends the most spectacular support to the multicavity structures expected out of a magnetized wideangle wind interacting with its ambient medium.
The organization of the paper is as follows. Section 2 illustrates the theoretical framework to investigate the kinematics of the proposed unified model. Section 3 gives the basic data processing setups adopted and utilized. Section 4 constructs the internal structures of the 2D thin hydrodynamic elongated momentum-conserving shell with locally oblique shocks formed by the wind and toroid density profiles. Section 5 describes the internal structures of the nonspherical elongated hydromagnetic bubbles. Section 6 traces curves in position-velocity space representing the outer shape of hydrodynamic and hydromagnetic models of outflow lobes. Sections 7 and 8 give results for the velocity profiles and flow patterns of the cases explored in Paper I and their position-velocity diagrams of column densities (PVDCDs), respectively, parallel and perpendicular to the outflow axes. In Section 9, we discuss the physical and observational implications, and we summarize our findings in Section 10.

THEORETICAL FRAMEWORK
In this section we review the theoretical concepts underlying the unified jet and wide-angle wind framework, as well as the analytical and computational methodologies adopted to investigate the properties of the proposed model (Shang et al. 2006;Wang et al. 2015;Shang et al. 2020). A system of toroidally magnetized wideangle winds with a stratified profile of density ρ ∝ −2 and magnetic field b φ ∝ −1 can be formulated in the absence of rotation, for an exactly radial cold flow of constant velocity, with a set of simplified equations (which correspond to a wind in equilibrium in the cold limit). This system of winds can interact with the ambient system of singular isothermal toroids, which may be poloidally magnetized. These toroids represent the end states of quasi-static hydromagnetic evolution followed by gravitational collapse (Li & Shu 1996b;Allen et al. 2003a,b;Galli et al. 2006).
Through the magnetic interplay between the inner wind and the series of ambient toroids, large-scale outflows can form. We highlight the conceptual considerations of the classical momentum-conserving shell model of molecular outflow as presented in SRLL, the inclusion of magnetic pressure developed in Paper I, and the theoretical nature of the spherical wind-blown bubbles as presented in Koo & McKee (1992a,b). We develop an analytic bubble-within-shell model to help illustrate conceptually the internal structures of the thin hydrodynamic momentum-conserving shell. We generalize this conceptual approach to treat the shapes of the elongated hydromagnetic bubbles with thick interiors in the integrated unified framework.
The overall wind characters considered in this work apply in general to magnetocentrifugal winds (as Xwinds and disk winds) launched from the innermost portion of the disks, where the launch regions are much smaller than the size scales before which the bubble establishes its self-similarity. However, the Xwinds have natural angle-dependent asymptotic behaviors  and fanning streamlines underlying the strongest wide-angle wind at the base. They possess inner "hollow cones" occupied by coronal magnetic fields (Shu et al. 1997) from the inner boundary interacting with the stellar magnetosphere, which require very high-angular resolution beyond that sufficient to resolve the simple jet features.

Physics of the Magnetized Wind
The governing equations utilized in this work are those of time-dependent axisymmetric MHD. They can be simplified to study a steady-state magnetized wind through a force balance similar to the asymptotics of a cold X-wind  at the mid-to-large horizontal scale ( ) far away from the launching region, leading to equations fully derived in Paper I for a cold radially directed flow of constant superfast velocity, and briefly summarized in Appendix A.
For a magnetocentrifugal wind coming out from the innermost region of a disk, only a small region of the disk or simply a very small rotating disk is required for the launching. Such a wind is dominated by a woundup toroidal component well beyond a transition region. The toroidal velocity of the wind declines steadily at large distances due to angular momentum conservation so that the wind flow becomes approximately poloidal (v φ = 0) with a purely toroidal field frozen inside (b p = 0), mass density ρ = ρ w , and gas pressure p = a 2 w ρ. For a completely cold wind, p = 0, the simplified steady-state Equations (A8)-(A11) admit an exact solution with a radial wind of constant velocity v = v rr , cylindrically stratified density ρ w ∝ −2 and force-free magnetic field b = b φφ . This basic radial wind justifies boundary conditions that are lightweight options for use in numerical simulations of flow propagation and interaction with the ambient medium. Based on this solution, the wind is set up on the inner radial boundary at r = r in as where = r sin θ is the cylindrical radius from the axis. Alternatively, we define or,ρ w ≡ r 2 ρ w , as part of the wind boundary condition setup, and also in free-wind regions. The boundary condition value b 0 gives the scale of the wind magnetic field. That value is parameterized using the Alfvén Mach number at the boundary condition,

Magnetic Forces in the Wind
The forces due to the b φ component, dominant or exclusive within the wind, can be expressed under axisymmetry within a current function formulation, as follows. We define a current function C ≡ − b φ , proportional to the total current carried in one hemisphere of the wind. This function can be used to find the current density and force term due to an axisymmetric b φ . The current density j p and the magnetic force can be computed using C, leading to a force function L = C 2 /2 behaving in a way similar to a pressure in that its gradient helps to obtain the force per unit mass For a free wind, the value of C stays exactly or nearly exactly constant. In the free-wind regions, the relevant gradients vanish, and the wind is current-free and force-free.
Interactions of the wind with the ambient environment gives rise to nonzero gradients of C and L, and induces currents and magnetic forces.
The compressed wind region allows not only density but also magnetic field b φ to be compressed, generating gradients of L = (b φ ) 2 /2. The respective magnetic force components f r,c = − 1 ρ 2 ∂ r L and f θ,c = − 1 rρ 2 ∂ θ L , zero in the free wind, can be generated inside the bulk of the compressed wind, leading to corresponding variations in v r and v θ , and accelerations. The variations of C and L distinguish the free wind from the compressed wind across the reverse shock boundary and compressed ambient material across the wind-ambient interface when the wind dynamically interacts with the ambient toroid.
The compressed b φ can also generate vorticity ω φ according to the torque Magnetic pseudopulses facilitate and sustain mixing across the boundary through the strong KHI feedback (see Section 5.3).

Initial Ambient Medium: Isothermal Toroids
We base our initial ambient medium on the ALS toroids as in Allen et al. (2003a, ALS), with minor additions and scalings, explained in this subsection. In the initial ambient region, v φ = 0 and b T φ = 0, and the magnetic field of the toroids is purely poloidal, taking the same form as Equation (A6), with an angular dependence given by a magnetic flux function φ(θ) as below: (6) where R(θ) and φ(θ) are dimensionless angular functions defined through the solution of an ordinary differential equation (ODE) given in Li & Shu (1996b), and a a is the isothermal sound speed of the ambient gas. The linear sequence given by n = 0 labels a separate value of the mass-to-flux ratio in dimensionless form, and also measures the degree of magnetic support against selfgravity. The n = 0 solution of the ODE R n=0 (θ) = 1 is a constant value representing a spherically symmetric ambient medium. For n > 0, the ambient medium has an opening region featuring very low-density values, more and more open and less and less dense as n increases.
This function, which depends on the parameter n ≥ 0, is subject to the normalization condition π/2 0 R(θ) sin θ dθ = 1 + n/4, and a set of boundary conditions that includes the limiting property R(θ) ≈ a 0 (n) sin n (θ) for sufficiently small θ → 0, making the toroids more open at the axis as n increases. This functional form plays a role for the hydrodynamic shell to be shown in Section 4.2 and Appendix D.
The initial ambient medium is set up as minor additions and scalings given in Equations (7)-(8) of the ALS toroids defined in Equation (6). The initial ambient density and magnetic field are given as where ρ S ≡ D S r −2 , in which D S is a constant equal to either 2.25 × 10 11 g cm −1 , or zero (defining either of the "tapered" or "untapered" toroids). The term in Equation (8) defining the initial magnetic field is derived from Equation (6), and the components of b T are found using Equation (A7) multiplied by a scaling constant α b , set to be 1, 0.1, or 0. Equation (7) is used alternatively throughout this work as motivating the definitionρ a ≡ r 2 ρ a , particularly relevant to the initial conditions, and to unperturbed portions of the ambient medium. The value D S represents the addition of a spherical density (equal to a tenuous n = 0 toroid) where D S is much smaller than a 2 a /2πG, the density scale defined for a regular toroid [without the angle factor R(θ)]. Almost all of the runs set α ρ = 1, but a few runs lack an ALS toroid in their initial setup and have α ρ = 0 instead. We often refer to the initial and later ambient medium of our typical run as a toroid (in a general sense), because it is derived from the ALS toroid setup.
Equation of state -The gas pressure p = a 2 ρ is calculated from the two-temperature equation of state (Wang et al. 2015) by a sound speed a based on the wind mass fraction f = ρ w /ρ and a 2 = f a 2 w + (1 − f )a 2 a , where ρ w and ρ are wind and total densities. The sound speeds for the wind and the ambient medium, a w and a a are given in Table 1.

Shapes of Momentum-Conserving Thin Shells
SRLL builds upon the assumption that a thin momentum-conserving (MC) swept-up shell of radius r MC ≡ r s can form by the input momentum of the wind, moving out with shell velocity v MC ≡ v s = dr MC /dt exactly radially along each local piece of the shell where no distribution of matter happens along the θ-direction. This can be parameterized by an angle-dependent ambient density distribution function Q(θ) and a momentum input function by the wind P (θ), defined in units such that the rate of change in mass M and momentum where v w is a radial wind velocity, v MC is the radial velocity of the momentumconserving shell model, v s is any generic radial shell velocity (numerical, observational, or analytical). The equations as in SRLL represent a hydrodynamic angledependent MC thin shell.
In real systems, the shell can acquire mass from both sweeping up the ambient medium and from collecting wind material: The shell velocity v MC is not strongly dependent on time, so the rate of momentum per steradian can be approximated as Collecting the terms in ρ w on the left-hand side, and those in ρ a on the righthand side, and canceling the radial factors allows us to rewrite the formula as which can be used to find the shape of an MC thin shell in the hydrodynamic limit. The two sides of Equation (11) can be interpreted as pressure balance, respectively, representing the inner and outer ram pressures in a local frame comoving with the radial shell velocity v s = v MC (θ), with the hydrody- The consideration from Equations (10)-(11) links the SRLL MC formulation to the expression adopted in Koo & McKee (1992b) for a spherical bubble driven by a constant-velocity wind, at each individual angle θ. This forms the basis of the integration of an elongated angledependent bubble inside the MC shell as realized in Paper I. The thin shells of the hydrodynamic MC model presented here represent the shape of the outer contour of an outflow lobe. Outflow shapes of magnetized cases are examined in Section 2.3.2 below.

Shapes of Magnetized Outflow Bubbles
Magnetized outflow bubbles can be explored using the methods of the previous subsection, once more with the objective of representing the shape of the outer contour of the outflow lobe through one-dimensional curves, analog to Equation (12), but now for a magnetized case. We adopt here the notation v MS instead of the specific v MC for generality. For magnetized winds, the wind magnetic pressure b 2 φ /2 is added to the inner ram pressure, modifying Equation (11) in a way similar to Draine (1983), and resulting in equation although the ram pressures are not isotropic (nor do the magnetized winds always lead to thin shells). With the presence of a poloidal magnetic field, a poloidal magnetic pressure b 2 p /2 term can be added to the right-hand side of Equation (13), as in Königl (1982): On the other hand, the poloidal magnetic field is strongly amplified in a region of compressed ambient field b p (Figure 1, and also Königl 1982), and that needs to be taken into account through an efficiency factor e p (θ). An efficiency factor β(θ) is multiplied to the v w term, in order to treat velocity behaviors inside the magnetized cocoon deviating from the strict thin-shell approximation, geometry effects, and shell thickness. Combining these terms and factors, the pressure balance equation can be written as Solutions of this quadratic Equation (15) allow estimates and fits of the shape function r MS = v MS t.
Inside the magnetized cocoon enclosing the outflow bubble, the magnetic field b p is amplified with respect to its initial value in the toroids due to the exclusion of poloidal field from the wind region. This amplification can be estimated by conservation of magnetic flux and considering cross-sectional areas characterizing features of the magnetized compression, and it is the main effect producing e p . The balance of the total ram pressure with that of the poloidal field seems to be a good predictor of the outermost edge of the b p -compressed ambient material for n = 1 and n = 2 cases of α b = 1, and all n's for α b = 0.1. The updated shape curves are illustrated in Figure 10.
Most of our models show strong self-similarity as demonstrated phenomenologically in Paper I, and they follow the principle of momentum conservation and pressure balance. The updated formula (Equation (15)) using β(θ) and e p fits is used for tracing outflow shapes and production of the analytic PV-space curves for the magnetized outflow shape curves in Section 6.2. Details of the fits are presented in Appendix B.

Internal Structures of Nonspherical Bubbles
In the context of a wind-blown bubble, Koo & McKee (1992a) considered the dynamics of a bubble blown by a constant-velocity wind, and their evolution in a uniform medium. Koo & McKee (1992b) generalized the formulation to a power-law density distribution ρ a (r) ∝ r −kρ driven by a constant-velocity wind, and the evolution of hydrodynamical momentum-conserving bubbles as thin shells in which both the shocked wind and ambient medium are radiative. A special case of r −2 dependence of the ambient density profile, such as that in SRLL for a singular isothermal sphere, is discussed in their Appendix A, in which the bubble as a thin shell expands with a constant velocity, as in Section 2.3.1.
The outflow shapes in the hydrodynamic limit result from wind and ambient (toroid) density profiles as obtained in Paper I. They are elongated (nearly momentum-conserving) bubbles driven by a wind in interaction with a toroid, both of them angle-dependent. For a wind constant in time at each angle, self-similarity is achieved (see Appendix A and Equation [2.10] of Koo & McKee 1992b) when the system evolves beyond the influence of initial injection, as shown in Figures 10 and  14 in Paper I. Within a hydrodynamic bubble, a reverse (wind) shock (RS), a contact discontinuity (CD), and a forward (ambient) shock (FS) form from inside out when steady injection of mass and energy into the ambient medium occurs. The CD separates the shocked wind and shocked ambient medium in a 1D sphere, and the two shock fronts separate the physical space into four respective regions of distinct physical properties. This basic 1D spherical bubble structure is best illustrated in the grid of models in Appendix A of Paper I with the relevant equations of state (two-temperature, and γ = 1.0001). In the situation where the shocks are oblique, shown in Appendix C of Paper I, the contact discontinuity becomes a TD, and significant shear can develop and build up across the shocks.
Paper I reveals the 2D structures of highly elongated, nonspherical hydromagnetic bubbles through the exploration of a large parameter space. Instead of the simple 1D RS-CD-FS structures present in the spherical bubble, the hydrodynamic nonspherical bubble formed with substantial shear across the now very oblique RS-FS surfaces along the TD.
"Multicavities" form associated with the RS-FS in the nonspherical magnetized bubbles and are evidently delineated by the respective shock boundaries. The primary wind is confined by the innermost RS cavity near the base around the jet in the strongly magnetized wind. Complex structures are shown to develop inside the outflow lobes, contrary to the hydrodynamical thin-shell models.
The interplay allows fairly complex structures to originate from simple setups. These interfaces and their crossing are demonstrated in the variations of the Cfunction, wind fraction f , vorticity ω, and the θcomponent of the velocity in Paper I. The correlations with f reveal the extended mixing regions confined by the RS on the wind side and the FS on the ambient side. Figures 2 and 3 show a glimpse of the structures formed based on simulations from the 6400 × 2016 set for the 100 km s −1 wind. The morphology of these hydromagnetic bubbles appears very elongated and fairly collimated, and shows some "spindle"-like appearance. They appear to be "prolate"-shaped.
In the following sections, we demonstrate the conceptual development of forming hydrodynamic RS, CD, and FS layers from a nonspherical prolate-shaped bubble. We construct the structures of the hydrodynamic bubbles analytically by locally "plane-parallel" shocks in 1D and 2D, and implement laminar layers of shear. The concept is then generalized to physically thick and locally oblique structures for the magnetized mixing layers. We construct analytic properties in positionvelocity space at incremental stages to form the baselines of signatures. These analytic models and kinematic signatures will be studied and juxtaposed together with those extracted from the MHD simulations. The traces of the underlying processes with their respective signatures can be thus identified.

Simulations
The setups of numerical simulations are described in detail in Paper I. We summarize key ingredients here for the numerical data utilized in this work, and their associated numerical values can be found in Table 1.
The parameter space originally explored in Paper I is huge. The wind is parameterized with seven M A = 6, 30, 60, 90, 180, 600, and ∞ values, and two velocities of v 0 = 50 and 100 km s −1 (a total of 14 wind combinations explored). Initial ambient medium is parameterized with α ρ = 1, four toroid configurations n = 1, 2, 4, and 6, and three levels of ambient toroid magnetization α b . Adding the four no-toroid configurations with α ρ = 0 but without n and α b , a total of 4 × 3 + 1 = 13 ambient configurations were explored. For each of these 14 × 13 parameter combinations, two numerical resolutions of 3200 × 1008 and 6400 × 2016 are available in the archive.
We focus on the v 0 = 100 km s −1 velocity and the higher-resolution 6400 × 2016 cases in this work. An additional set of hydrodynamic (M A = ∞, α b = 0) cases is performed for v 0 = 100 km s −1 at 6400 × 2016 resolution, with v θ = 0 to imitate an equivalent 1D bubble in r-coordinate extended to 2D at every θ by independent variation of θ. This set serves as our reference for smooth v θ = 0 hydrodynamic models without the onset of the KHI.

Synline Data Post-processing
The numerical data are post-processed through the tool, Synline, developed in our data pipeline locally by the team. Synline is a package generating synthetic observables from two-dimensional (2D) axisymmetric simulation data. The 2D data blocks in desired properties such as the number density and velocity are first rotated about the symmetry axis, interpolated and mapped onto a new three-dimensional (3D) Cartesian uniform grid. Inclination of the system to the plane of the sky is incorporated during the construction of the new 3D data cubes. The synthetic observables are produced by performing line-of-sight integration of the emissivities of local cells at the optically thin limit.
In the Synline package, one has the ability to calculate the emissivity by a selection of line transition and its associated level populations with constants out of a database included in the package, under the non-LTE assumption with temperature as an input field. This final module of Synline is not used in our work, as we generate only column density maps without calling the emissivity routines. The reason for this is that our simulations do not provide a point-to-point map of the local temperature, which would be necessary to estimate the emissivity.
Column density maps are generated without calling the emissivity routines by directly summing up the number densities along the line of sight with the proper inclination angles. A position-position-velocity (PPV) data cube in column density is generated by binning the local line-of-sight velocities at an (adjustable) interval of 0.06 km s −1 . A thermal line profile is applied according to the local sound speed a 2 = f a 2 w + (1 − f )a 2 a for the ∼ 100 K wind and the ∼ 10 K toroid in our setup (see Paper I). Subsequent production of the PV diagrams of column density (PVDCD) follows by manipulating the data cubes parallel and perpendicular to the outflow axis in the projected plane of sky.
We stress once again that the resulting diagrams that illustrate the column densities mapped in the PV space are not directly analogous to the PV diagrams as resulting from observations, which are based on surface brightness of emission lines. The line excitation depends strongly on the temperature, and the resulting surface brightness depends on the emission line selected, and on its appropriate radiative transport: for forbidden lines, it results just as a sum of the emissivities of the volume elements along the line of sight, while for permitted lines, the full radiative transport has to be taken into account, as the emission is mediated by absorption and re-emission at Doppler shifts corresponding to the local radial velocity.
The self-similar nature of the evolution as demonstrated in Paper I allows a self-similar presentation of the kinematic features. The PVDCDs are shown in the self-similar units normalized to v 0 t and v 0 , respectively, at a representative epoch of 100 yr. The spatial position is convolved with a Gaussian profile of 0.01v 0 t for a smoother appearance. Different criteria can be applied prior to the density integration to extract different physical conditions from the data cubes. The range of the (a) For α b = 1 cases, the absolute values of C function, normalized to b0, are shown on the left halves and number density nH, rescaled and spanning 4.5 dex, on the right halves in individual panels. Both quantities are shown in logarithmic scales and filtered with a threshold of wind mass fraction f > 10 −4 . The -axis has been exaggerated to show detailed structures. The exaggeration factors are 4.7, 2.7, 1.5, and 1.0 for n = 1, 2, 4, and 6, respectively.   (a) For α b = 1 cases, the azimuthal vorticity |ω φ | in logarithmic scales is shown on the left halves and wind mass fraction f on the right halves of individual panels. The -axis has been exaggerated to show detailed structures. The exaggeration factors are 4.2, 2.6, 1.5, and 1.0 for n = 1, 2, 4, and 6, respectively.   Figure 3 (a), but for α b = 0 cases. The -axis has been exaggerated to show detailed structures. The exaggeration factors are 4.2, 2.6, 1.5, and 1.0 for n = 1, 2, 4, and 6, respectively. wind mass fraction 0 < f < 1 is used to identify the contribution from the compressed wind and compressed ambient material confined by the RS and FS. The criterion, v p /a a > 1, is implemented as the outer boundary across the FS.

STRUCTURES OF 2D ELONGATED HYDRODYNAMIC BUBBLES WITH LOCALLY OBLIQUE SHOCKS
We establish in this section the basic structures of the 2D momentum-conserving hydrodynamic bubbles using semianalytic procedures. The procedures are novel in the construction of analytical expressions for the shell thickness across the RS, CD, and FS, the extension to 2D with locally oblique shocks, and development of a shear profile within the post-shocked regions. This whole set of procedures gives the reference thin-shell models of 2D hydrodynamic nonspherical bubbles with desired shear profiles without the turbulence induced by the KHI nor the thick post-shocked regions. The resulting kinematic signatures will serve as a baseline for detecting those generated by the KHI and the thick extended magnetized shocked-wind region.
The obliqueness angle -We define the obliqueness χ locally in terms of the angle that the radial wind makes to the normal to the shape curve r s (θ): χ = arcsin(n ·θ) = arcsin(−r s /(r 2 s + r 2 s ) 1/2 ), (16) defined so that minimum obliqueness (such as in a spherical shape curve) corresponds to χ = 0 • , and maximum obliqueness χ = 90 • corresponds to a locally radial tangent to the shape curve. In our most usual case, both in simulations and in analytical models, the obliqueness χ is in the first quadrant, allowing us to write it as arccos(n ·r). Occasionally in simulations there are a few cases where the normal is closer to the axis than the radial line, moving that χ < 0 into the fourth quadrant at that place. The fact that typically χ > 0 implies that the equatorial radius of these outflows is narrower than the axial radius, giving them an overall elongated shape.

Construction of Smooth Hydrodynamic Models
We lay out the conceptual development in four progressive steps to demonstrate how the structures are constructed. For completeness and for reading as a unit, the entire package of the procedure is described in more details in Appendix C. Here we outline the four steps in sequence and their respective methodology.
We begin with the time-evolving position of the MC shape, which traces a curve in space for each of the n = 1, 2, 4, and 6 nonmagnetized toroids, with an input wind of M A = ∞. Procedures to generate these curves are given in Section 6 and Appendix C, and curves are shown in the upper-left panel in Figure 22 of Paper I. This is referred to as Step 1 in Appendix C, and in Equation 12 of this work. In Step 2, the locations of the reverse and forward shocks are constructed into the 1D flow solutions with the proper obliqueness in the 2D flow without the complexity encountered in the real numerical simulations. We use these solutions as the clean and smooth background references. The local obliqueness χ(θ) at the position r s (θ) of the momentum-conserving curves found in Step 1 can be incorporated into the construction process.
A hydrodynamic 1D Riemann problem can be defined, having its left-side state with the density, projected velocity and sound speed of the wind, and its right-side state the with ambient density and sound speeds, with a velocity of zero. When solving this Riemann problem, an EOS with γ = 1.0001 is adopted for simplicity and for its close resemblance to our two-temperature equation of state, as shown in Appendix A of Paper I. The solution to the 1D Riemann problem gives the three discontinuities (or signal speeds) and is adopted to construct the locations of the RS, CD, and FS. The respective shock surfaces are approximated locally as two plane-parallel waves to allow for such choice of methodology. Depending on how the propagation of the shocks is calculated, three models of thickness can follow with various levels of fitting and projection of vectors (the normal, radial, and hybrid models; see Step 2 of Appendix C). The shaded areas in Figure 4 show the compressed (shocked) regions between the positions of RS and FS obtained by the Riemann problem, on top of the "no-v θ " (v θ = 0) numerical smooth reference calculation.
In Steps 3 and 4, we utilize the analytic formulation of the jump conditions and behaviors of the parallel oblique shocks, obtained previously from Paper I, and construct analytically a new family of shear profiles in Step 4.
Step 3 requires detailed construction of shock velocities in different frames of reference, and in their respective normal and tangential directions. Therefore, we build upon the understanding of parallel oblique shocks and their respective jump conditions at the RS and FS, as shown in Appendix C of Paper I. The detailed transformations among the frames of reference in RS and FS, pre-shock, post-shock and fixed frames, velocities normal or tangential to the shock surfaces, are constructed and implemented in Step 3 of Appendix C. In Step 4, an analytical expression further completes the local velocity due to shear for an arbitrary point residing within the post-shock region bounded by the RS and FS. This step extends the knowledge on the total shear across the plane-parallel oblique shocks derived in Appendix C of Paper I. We construct the shear profiles for any point between a pair of locally parallel oblique shocks, once the end points on the RS and the FS can be identified. Figure 23 shows the shear profiles for the shocked region in self-similar coordinates. The patterns across the RS of the deflected flow directions from the original incoming radial wind are evidently demonstrated.

Obliqueness of Thin Hydrodynamic Shells
The obliqueness angle χ has positive values (in the first quadrant) in our usual case in whichn is further away from the axis thanr, and negative values in the opposite case, which can happen only if r s > 0 (not in our set of momentum-conserving examples). Hence, it is safe to define the shape curve S as that of the momentum-conserving model, with positions r s = r MC and normal vectorsn =n MC as in Steps 1 and 2 of Appendix C. The limiting behavior of χ near the axis (small values of the coordinate angle θ), near the tip of the outflow, is a peculiar property of the detailed hydrodynamic momentum-conserving models.
For the obliqueness at the tip, we show in Appendix D the derivation of the small θ behavior with the density ratio δ = ρ w /ρ a . Using Equation (16), one can obtain the obliqueness for small angles in Equations (D33)-(D36), and apply them to both the untapered and tapered toroids with the background scaled density D S (see Equation (9)).
For the tapered toroids, as we have adopted primarily throughout Paper I and this work, in which D S > 0, one reaches a sin χ that is independent of n. For our usual parameter values as shown in Table 1 of Paper I, this leads to a modest value ∼ 47 • with our usual winds and setup. Such behaviors can be observed in Figure 9, and in the limiting values of χ in Figure 22 of Paper I. The numerical values of χ indeed approach the same ≈ 47 • for the constructed MC curves, and numerical calculations without v θ (v θ = 0) and with regular, unrestricted v θ simulations. A scan of the n-independent χ values for a broad range of D S can be seen in Figure 24. Figure 25 shows a collection of detailed resultant MC curves of hydrodynamic bubble shapes associated with the varying D S and their obliqueness in the zoomed-in tips. As D S decreases, the medium is less dense and resists the wind less. For n > 0 and D S = 0, the opening region of the toroid is so light that it allows the shell velocity v MC to approach its maximum value v 0 t = r MC /t. The zoomedin tips show this locally flat (χ = 0) limiting shape for small θ values, within the toroid opening for each n case.
As D S increases, the shell velocity decreases, the shapes become narrower, and acquire a pointy tip (χ > 0) at the axis θ → 0 as shown in Figure 24. The n = 0 toroid is different because it has no opening, even at θ = 0, leading to a shape that is always narrow with a pointy tip. We note that the match of behaviors of χ at small (including θ → 0) angles with the opening angles at the bottom indeed decides the overall shapes of elongated bubbles.

INTERNAL STRUCTURES OF THE
ELONGATED NONSPHERICAL HYDROMAGNETIC BUBBLES

Multicavities Revealed
We reviewed the conceptual development behind the "multicavities" in Section 2.3.3 of the physical framework. The primary wind and the ambient nonperturbed toroids are vorticity-free. The vorticity ω is generated behind the shocks and further enhanced by magnetic forces. The C function varies abruptly across the shocks and boundaries into the nonmagnetized regions, and v θ becomes nonzero and varies in the post-shocked regions. Figures 2 and 3 illustrate the structures and the sense of "multiple" cavities out of the thick and extended mixing region and the magnetized ambient medium. An innermost cavity forms inside the loci of the reverse shock, in which the primary wind is unshocked and unperturbed. The outermost cavity forms with the magnetic cocoon by the surrounding ambient poloidal field and the compressed ambient medium exterior to the compressed field. Such layered structures could be interpreted as the layered shells from episodic ejections. The nonlinear growth and coalescence of the KHI modes further complicate the apparent presentation of the structures as a mixture of large and small density concentrations resembling multiple large and small shells. These structures are not expected of the thin-shell models. A magnetized outflow is an elongated bubble filled with internal structures, not an empty cavity.
We illustrate in Figures 5 and 6 density distributions of internal structures pertaining to the hydromagnetic bubbles within the outflow lobes surrounded by a moderate to strong ambient magnetic field. These figures reveal more details at a resolution ∼ 2×2 higher than that of Figure 18 of Paper I, and scaled to the self-similar coordinates. The column and number densities juxtaposed naturally cast the visual impressions of nested "multicavities" from the inner to outer cavities and the apparent extended shells mimicking the episodicity. The visual impression of the outermost cavities, however, only extends to partial lengths of the outflow lobes. The contributions to column density pass through the density- depleted gap region filled with the ambient poloidal field lines, which is wider in horizontal scales for the larger n values. This feature is enhanced at higher resolution, and forms the basis of the interpretation of an extended low-velocity outflowing gas from the base of the outer outflow lobes. The newly detected structures around HH 212 (reported in Lee et al. 2021; see Section 9.5 below) may be an example of the manifestation.

Velocity Structures of the Multicavities within the Magnetized Bubbles
Here we reveal the velocity profiles and patterns appearing in the internal structures of the compressed layers in the elongated hydromagnetic bubbles. As demonstrated using plane-parallel hydrodynamic oblique shocks in Section 4.1 and Appendix C, substantial shear is present across the shocked region bounded by the RS and FS. The hydrodynamic shear Figure 5. Distribution of the number density (left) and the column density integrated over regions with vp > aa (right) for magnetized ambient media with α b = 1. The rescaled logarithmic number density has a range of 4.5 orders of magnitude, while the also rescaled logarithmic column density has a range of 4 orders of magnitude. The ambient bp field lines have been drawn on top, as in Figure 4 of Paper I. The -direction of these panels is stretched to reveal more enlarged views of the detailed structures. The exaggeration factors are 4.2, 2.6, 1.5, and 1.0 for n = 1, 2, 4, and 6, respectively. Figure 6. Distribution of the number density (left) and the column density integrated over regions with vp > aa (right) for magnetized ambient media with α b = 0.1, both rescaled and shown as in Figure 5. The -direction of these panels is stretched to reveal more detailed structures. The exaggeration factors are 4.2, 2.6, 1.5, and 1.0 for n = 1, 2, 4, and 6, respectively. profiles for smooth analytical models are illustrated in Figure 23. The directions and magnitudes of velocity vectors vary as the flows cross the reverse shock surfaces by following the jump condition at the RS. This property has significant observational implications for wind diagnostics.
In Figures 7 and 8, we track flow properties in the extended compressed regions between the RS and the FS. Our setup of a constant-velocity primary wind allows an easier tracking of the change of directions across the oblique shocks. The generation of vorticity through magnetic forces, oblique shocks, and nonlinearly grown modes of the KHI in the toroidally magnetized winds introduces local perturbations to the bulk flow. However, the overall patterns of bent and deflected flow vectors relative to the original radial free-wind are self-evident. The resultant lobe-scale flow patterns appear as emerging from near the base of the outflows, moving around the widest waist, then converging to the top. This apparent bulk flow motion can give an impression of a secondary wind rising on top of the "disk" or "disk atmosphere" as if it were an extended disk wind (EDW), which is launched from a few to tens of astronomical units from the inner launch loci of the primary wind. The inference of a separate launch region of tens of astronomical units by the reduced post-RS speed and deflected direction is a risky application of the connection between the terminal velocities and the launch radii, and a misuse of the formula derived in Anderson et al. (2003). Further discussion of the conceptual misunderstanding continues in Sections 9.2 and 9.4.
Theoretically, these complex nested velocity shells simply occur naturally in elongated toroidally magnetized wind-blown bubbles as part of the physical mechanisms. A separate EDW launched slightly outside of and surrounding the primary wind is not required to generate the range of the extended intermediate velocities. The spatial correlation of observed PV diagrams with the illustrated column density maps provides a hint for the nature of the occurrence of the different velocity components between the jet and the very low shell velocities near the outflow base (see Section 9.2).

Magnetic Pseudopulses
Here we revisit the generation of "pseudopulses" in the very extended compressed wind region. These pseudopulses were called magnetic pulses in Paper I.
In a toroidally magnetized compressed wind medium, vorticity can be generated by magnetic forces and fed back into the already grown KHI modes. Likewise, the enhanced vorticity drives the already pulsed b φ with stronger amplitude, allowing the magnetic field to com-press and relax further in the z-direction. This feedback loop produces self-generating pseudopulses cascading the magnetic energy in this region. Compounded by density stratification near the jet axis, it creates the impression of pulsed jets in the compressed wind region.
In a magnetized medium, the vorticity, ω = ∇ × v, can be generated by magnetic forces: where the terms on the right are the baroclinic term and the effects of magnetic forces. Inside the unmixed wind, where b p is zero or tiny, Equation (5) shows the dominant term (due to b φ ) of this magnetic torque, The magnetic accelerations and the induced variations of v θ = 0 in turn feed into the induction equation allowing acceleration and deceleration in the rdirection, which then causes the variation of v r , leading to the presence of pseudopulses in the radial direction. They leave even stronger marks at the higher-resolution (6400 × 2016) set. The connections of f r to v r and f θ to v θ are shown in Figure 26 in Appendix E.

LOBE SHAPE CURVES IN PV SPACE
We start with producing curves, one-dimensional loci of the MC shapes, in axisymmetric position space. Then we project those MC curves onto PV space, without line-of-sight integration. Subsection 6.2 continues this for magnetized MS models.

Momentum-Conserving Curves in PV Space
The self-similar MC curves in r-θ axisymmetric position space are computed in Section 2.3.1, Equation (12). Projecting them onto axisymmetric -z position coordinates, we build MC curves in -z space or plane: Usually the coordinates θ and are restricted to positive values. However, mirror axisymmetry of the MC curve in the 2D position space ( -z or r-θ) allows an extension of and θ into negative values. Following the curve defined by Equation (18) for a range such as −180 • ≤ θ ≤ +180 • parameterizes a continuous curve representing the two lobes, with lobe tips located at θ = 0 • and 180 • . Inside our n-range, the equatorial wind/ambient density ratio is very small, giving the curves a narrow waist or "neck" at the equator, located at θ = ±90 • . The −90 • ≤ θ ≤ +90 • range suffices to parameterize a single lobe. These vectors are normalized to the nominal wind speed v0. They are plotted on top of the number density on left, along with the column density integrated over regions with vp > aa shown on the right, both rescaled and shown as in Figure 5. Both the -direction of these panels and -component of the vectors are stretched, with the same n-dependent exaggeration factors as in Figure 5, to reveal more detailed structures. The ambient bp field lines are shown as in Figure 5.  Figure 6. Both the -direction of these panels and -component of the vectors are stretched, with the same n-dependent exaggeration factors as in Figure 6, to reveal more detailed structures.
We now construct the MC curves in PV space. By defining a line of sight with inclination angle i, the MC curves of Equation (18) can be projected onto the PV plane (velocity v, position p) as showing perfect self-similarity. The inclination angle i = 0 • corresponds to a bipolar outflow axis aligned along the line of sight, with |θ| < 90 • corresponding to the blueshifted side (negative v values), and 90 • < |θ| < 180 • corresponding to the redshifted side. The redshifted side can be alternatively parameterized with i = 180 • and |θ| < 90 • . A bipolar outflow axis lying on the plane of the sky corresponds to i = 90 • . Figure  9 illustrates the relationship of the inclination angles, base opening angles, and the toroid n values, showing a consistent trend of shapes and orientations with respect to the evolution in the n values. Figure 9 and Equation (19) show that when using self-similar units (v 0 for velocity and v 0 t for position) a change of inclination angle i corresponds to a simple rotation of the MC curves in PV space. We also note that the mathematical shapes of the curves traced by Equations (18) and (19) are identical. Both of these remarkable properties are related to self-similarity. The mathematical identity and rotation symmetry allow to use Figure 9 not only to represent curves in PV space as labeled, but also the curves in -z space. These MC space curves (both in PV and -z spaces) are elongated oval-shaped contours, symmetric across an equatorial line, and also around an axis connecting the two apices (each apex touching the tips of the jets). The series of oval curves follows the shapes of the individual lobes labeled by the toroid index n, and is further connected to the curve family with varying M A discussed in Section 6.2. For typical (not too small) values of n, the base of the curves near |θ| ∼ 90 • has a "neck" that gets very close to the origin of coordinates, but the curves are smooth and do not cross each other before moving into other quadrants of the plane. This situation can be seen in the insets of Figure 9.
The base region of the outflow can be described in terms of the MC curves near |θ| = 90 • as shown in the insets. A useful parameter to describe that base region is an opening angle o between the two tangent (or secant) lines to opposite sides of the lobes. The inflection point of the MC curves provides a maximum opening value, which can be used to parameterize the base of the outflow at very low heights (top-right inset of Figure 9). Alternatively, a tangent or secant connected to the larger scale of the curves (bottom-left inset) may tell more about the overall outflow structure and observed images. The opening angle based on the idealized MC curves in self-similar units has the same value in either -z coordinates, or on PV space for any inclination angle i. This is because of self-similarity making the curves mathematically identical between these two spaces and responding to a change in i as a rotation in PV space.

Magnetized Outflow Shape Curves in PV space
We consider now the regimes of hydromagnetic winds interacting with magnetized or nonmagnetized ambient toroids. When the winds are toroidally magnetized, new phenomena arise with the thick and extended compressed regions, the KHI and pseudopulses feedback loop with generation of vorticity. This motivates us to compute magnetized outflow shape curves for different M A values using Equation (15).
Equation (15) gives the MS (magnetized outflow shape) of a hydromagnetic wind carrying a b φ interacting with a possibly magnetized toroid carrying a b p . This is a magnetized analog of Equation (11). Equation (15) is exactly as self-similar and mirror-symmetric as the hydrodynamic model of Section 6.1, and therefore equations similar to Equations (18) and (19) can be written, sharing thus in symmetry properties such as the identity of mathematical shapes traced in PV and -z spaces, and changes of inclination angle i being equivalent to rotations of the curves in PV space. Figure 10 summarizes such MS shape curves. The variation of the contours of MS curves with different n values is evident for the entire M A range. Prominent dependence of shape curves on the magnetic parameters is salient for the more magnetized cases, the M A 30 range for the wind, and α b = 1 for the medium. Detailed properties derived from Equation (15) are given in Appendix B.

THE POSITION-VELOCITY DIAGRAMS
In this section, we build signatures of each featured component of the integrated hydromagnetic outflow bubble in position-velocity diagrams of column density (PVDCD). We present the parallel PVDCDs, in which the position y represents the distance along the projected outflow axis.
We note that by producing the column density maps in PV space, it is difficult to recover identical information from real observations without including specific tracers and their line excitation. Different tracers probe different physical conditions, and tracers are subject to thermochemistry and chemical evolution, whose detailed conditions are beyond the scope of current work.
We have assumed an inclination angle of 45 • for PVDCD maps in this work. At this angle, the signa- tures of each of the components are separated the farthest. However, in real systems, the signatures can be embedded in different inclination angles from 0 • to 90 • , and they will create different impressions. Modeling of real systems requires the knowledge of their inclination angles and of their physical conditions, and comparison should include a calculation of the emissivity maps. We leave such an effort to a future work.

The Hydrodynamic No-v θ Models
In this subsection, we build four smooth reference PVDCDs of the "no-v θ " (v θ = 0) calculation runs for the hydrodynamic M A = ∞ and α b = 0 cases, one for each of our four n values. The four runs represent at every θ angle a 1D radial bubble. This series with a restricted v θ value will be contrasted against the cases with regular v θ simulation treatment in subsections that follow.
The no-v θ models are expected to preserve the basic features of the dual jet plus thin-shell structures without the complication from the presence of the KHI, despite the presence of shear. The central jet is the densest part of the underlying wind with a −2 profile, which straddles a velocity centroid with broad wings extended to the positive and the negative ranges. The divergence of the radially pointed flow naturally produces a velocity spread projected onto the line of sight, as demonstrated in Figures 1 and 2 of Paper I. The flow near the base has close to equal projections toward and away from the observer at a full self-similar velocity |v/v 0 | = 1, giv-  . Position-velocity diagrams of column density (PVDCD) of the no-v θ (v θ = 0) models for n = 1, 2, 4, and 6 toroids (left to right) at inclinations of 90 • , 60 • , 45 • , 30 • , and 1 • (top to bottom). Column density is rescaled logarithmically for the 4 orders of magnitude that are covered. The unshocked ambient material with vp < aa has been filtered out while producing the PVDCDs. Momentum-conserving curves are overlaid as white dashed lines.
ing rise to a total velocity width of ∼ 2 near the base. Moving up the length of the jet, the projected velocity along the line of sight decreases, and the radial velocity vectors narrow down to the vicinity of the projected velocity centroid, causing the widths to decrease.
The PVDCDs of the smooth reference no-v θ models at different inclination angles i with the full coverage of both lobes in all four quadrants are depicted in Figure 11. These panels illustrate the variation of the PVDCD patterns with n and the inclination angle i, along with the thin-shell MC curves in PV space. The four quadrants are needed to show the positive and negative projected velocities across the hemispheres due to inclinations. Clearly the two jet-shell components are distinctly visible. The densest part of the wide-angle wind maintains well the jet appearance at the projected velocity v 0 cos i while the radial part of the flow is projected around the jet velocity centroid, giving a width of 2v 0 sin i at the base. The width varies with the thinshell opening confined by the ambient toroid, leaving fast narrowing of the width for the small n but gradual narrowing for the larger n.
The divergent radial flow with the density profile of −2 allows the integrated column density to drop off approximately as −1 . The MC curves in PV space, as in Section 6.1, tracing the thin-shell cavities naturally lay on top of the projected jet-wind features in the PVDCD of these no-v θ models. The bottom portion of the MC curves coincides with the base of the thin-shell cavities very well. These cases remain well-traced by the MC curves. It is expected that MC thin shells would become more conical for larger n and more collimated for smaller n, because the toroid density functions for larger n are also further away from being spherical. The face-on situation at 1 • , equivalent to the situation of a cone with small to large opening angles moving at v 0 , is shown at the bottom row of Figure 11. The small to large opening angles of the axisymmetric cones give the vertical width in y/v 0 t located at v/v 0 = ±1.

Hydrodynamic Simulations with Regular v θ
To enter the more realistic regime in which shear is allowed to induce KHI and other complexities, we first make a comparison between hydrodynamical simulations with regular v θ and those calculations with the nov θ (v θ = 0) constraint. The comparison of the two sets of models is made for the hydrodynamic cases, having control parameter values M A = ∞ and α b = 0. Figure  12 illustrates the differences the presence of KHI makes for one inclination angle (45 • ). Significant changes can be seen clearly in the bottom panels. The changes in the distribution of column density with the position and ve-locity are considered to come from the physical effects of active shear within the thin layer between the reverse and forward shocks.
Within models with active shear (a contrast with the no-v θ models in Section 7.1, which do have intense shear, but "inactivated" to induce a KHI), patches of featherlike structures extend across the ellipsoidal region of the cavity walls between the two wide-angle jet features. This action smears the thin-shell features traced by the MC curves. These patches arise from the large spread of velocity across the post-shock regions, produced by the large "spikes" or the large vortices aligned along the spikes present in the HD models illustrated in Appendix D of Paper I. The fuzzy threads protruding the MC curves trace the locally parallel shock fronts making different obliqueness angles across the shocked shells from the interiors of the physical cavities, through the outer surfaces, then into the ambient medium.
The velocity vectors moving along the shell surfaces are compared in Figure 13 between the no-v θ (left) and regular v θ (right) cases. The velocity vectors from the simulations are shown along the positions on the theoretical MC curves. Variations of the velocities can be seen as a result of the development of KHI. In the nov θ calculations, the resulting shells follow more or less the MC curves, and the velocity variations appear to be smooth. In the simulations with regular v θ , the MC curves pass through the instability patterns. The vectors follow zigzag paths, resulting in fluctuating directions and magnitudes.

Kinematic Signatures of Outflows Formed by Magnetized Winds
For the synthetic PVDCDs of outflows made with column densities, we integrate densities along different individual lines of sight and collect the binned velocities. These PVDCD panels are shown in Figure 14 for the nonmagnetized ambient media and in Figures 15 and 16 for the respective weakly and strongly magnetized media. In this subsection, we demonstrate signatures produced by magnetized winds, free or compressed, from their interactions with the ambient media.
The PVDCDs derived from the hydromagnetic winds constructed by our models contain two intrinsic components: the jet and the wide-angle radial velocity vectors, and the structures induced by their interactions with the ambient medium. For the free-wind models, as shown in Figures 1 and 2 of Paper I, the jet collimated with the density profile ρ ∝ −2 appears at the velocity peaked as v 0 cos i, as the peak of emission near the base, surrounded by the broad wings produced by the radial velocity vectors of the wide-angle wind. The outlining con-  tours of the column density imitated emission (CDIE) originating from the highest and lowest projected velocities move up with the length of the jet and converge gradually near the tip of the jet. The basic feature of the outflow cavity produced by the wide-angle portion appears near the origin of the PV coordinates. They vary in physical extent from n = 1 to n = 6, in their respective spread from the tip of their respective MS cone that forms from the revolution of the MS curves.
These jet plus wide-angle signatures are clearly distinguishable from the MS cone on the PV plane, similar to the MC analogs in the PVDCD of Figure 12 The interactions of the magnetized wind with ambient media add the thick and extended compressed wind regions to the characteristics of the PVDCDs, as threads connecting the PV origin and the jet features. The jets with wide-angle wind features are confined under smaller triangular-shaped areas following the boundaries of the RS. (We refer to this triangular-shaped area as the "RS cavity".) On top of the jet tips, the jet portion of the compressed winds near the jet axes is associated with clear but wiggled velocity centroids around the original values. The wiggled jet features trace the influence of the magnetic pseudopulses in the jet portion. These kinematic features are evident in Figures 14, 15, and 16. The CDIE distribution appears extended in rectangular or butterfly-like areas between the systemic and the projected jet velocities as extended intermediate to low velocities, which are linked to the nested velocity features discussed in Section 9.2.
Small to large patches and threads of velocity patterns across the jet-shell trace the underlying shear and instabilities within the shocked regions. This feature is known to follow the occurrence of the KHI from Sections 7.1 and 7.2, and Figure 12. It increases in PV coverage as the wind magnetization grows stronger until M A ≈ 30. The patches in M A 600 spread wider toward the jet features. The patches in the smaller n cases are more horizontally intertwined than the larger n cases, while the n = 6 case has less of this effect due to their wide-opened base cavities.
The presence of these features is consistent with the origins of KHI in the extended compressed wind regions. The boundaries between the MS cone and those of the jet plus wide-angle wind are blurred; however, the cavity formed by the RS is distinctly traceable and confined by these intertwined patches on the PVDCD. The patchy appearance of these features suggests the likely association with mild KHI shocks. In real PVs, they may appear as fainter, scattered emissions with different local excitation conditions.
Another noticeable feature associated with the jet component arises on the PVDCD for a strongly magnetized wind at large positions. This apparent acceleration exists beyond a certain y/v 0 t from the top portion of the compressed wind region, for the strongly magnetized M A = 6 and the moderately magnetized M A = 30 cases. Some tiny tips remain all the way to M A = 60. We note that this increase of velocity only occurs at the outer boundary of the compressed wind region, which is far beyond the RS enclosing the free-wind propagation region. The nature of this feature is to be explored in depth in Paper III.

Signatures of the Magnetized Mixing Regions
Here we highlight the kinematic features that arise from the magnetic interplay. These include structures originating in the interacting compressed wind and ambient regions, which cover a substantial volume of mixing and entrainment, the jet fluctuations caused by the pseudopulses, the patterns of the nonlinear KHI, the multicavities, and the RS cavity enclosing the primary wind with broad velocity widths near the base of the jet. Naturally, the resultant PVDCDs would contain the imprints of these ingredients in their respective regions that alter the simple smooth thin-shell systems.
The signatures of the compressed magnetized winds appear on the PVDCD maps as extended threads of CDIEs connecting the outskirts of the "jet" with wiggling velocity centroids and the "sporadic" patches around the origin and the y-axis, confining the broad velocity feature at the jet base across the RS cavity. These are nonexistent in Figure 11. These features arising from the compressed regions with mixing can be viewed on the PVDCD maps filtered with the f values between 0 (pure ambient material) and 1 (pure wind material). With this procedure, the gas of the primary wind and the ambient toroid is subtracted from the contributions to the PVDCDs maps. One can probe directly the PV patterns contributed by mixing. Figure 17 is one such map for the strongly magnetized ambient medium α b = 1, with both the pure wind f very close to 1 and pure ambient material f very close to 0 subtracted. The   panels capture the spatial and spectral distributions of mixing resulting from the interplay.
These patterns suggest that identifying the patchy or feather-like structures by their additional occurrence in the PVDCDs compared with the reference ones is appropriate and consistent with their proposed origins in the internal structures of magnetized elongated bubbles. The intermediate velocity components (IVCs) on the PVDCDs may appear as part of these extended threads of CDIEs when local excitation conditions meet the requirements for the tracing molecular or atomic lines for the jets and winds in Section 9.2. These can further connect to the even lower-velocity components (below).
An extended patch of low-velocity features is naturally present in Figures 14, 15, 16 and 17 throughout the parameter space. Because of the conventional knowledge of the MC thin-shell models in the HD regime, the occurrence of low-velocity components (LVCs) near the origin of the PV coordinates is always expected. The LVCs near the origin of PV coordinates are thus considered to be part of the traces for the MC or MS contours. Distributions of the low-velocity profiles outside of the elliptical contours of the MC/MS curves along the y-axis increase with the magnetization of the winds, which are more evident for the n = 1 and 2 systems and more in the α b > 0 than the α b = 0 ones. The magnetized gap regions appear spur-like along the y-axis for the strongly magnetized ambient α b = 1, indicative of the vortices aligned to the lobe surface confined by the ambient poloidal cocoons.
We note that the disappearance of the observed increase of velocity in the jet component in Figure 17 suggests the phenomenon is connected to the intrinsic wind, free of ambient materials. That is, this acceleration phenomenon is not caused by mixed or entrained material. We note this apparent acceleration is connected to the needle-like appearance at the tip in the 2D density distribution in our Figures 2, 3, and 5-8, at the far end of the post-RS regions before hitting the FS. The needles were attributed to magnetic forces that accelerate v r while b φ smoothly decreases inside the tip region. The last paragraph of Section 5.2 in Paper I gives a first account of such possibility, which we will explore further in Paper III.
In the magnetic interplay, pseudopulses are generated in the post-RS compressed wind regions. This is evident in PVDCDs as oscillations in jet velocities downstream behind the RS cavities but upstream of the FS. They appear and operate cooperatively with the generation of vorticity through magnetic forces as described in Section 5.3 and in Figure 26 shown in Appendix E. When the fingers grow sufficiently large around the RS cavity as they advect, they generate the impression of multiple converging "zigzag" shells mimicking episodicity along the cavity as shown in Figures 5 and 6. In the extended compressed wind region, the KHI fingers advect and coalesce, and larger interconnected fingers form, as part of the feedback loop in Equation (17). When the large fingers merge near the jet proper, the feedback loop operates across the whole width of the lobe, leading to perturbations of density and velocity on the jet. Variations in density and the radial velocity component v r advected, projected onto the PVDCD, will appear as small quasiperiodic wiggles around the jet velocity centroids propagating toward the tip. Such variation in v r along the axis is mild, and occurs with shorter periodicity, extended along the upper portion of the whole outflow lobe. As the actual patterns depend on local properties as discussed in Section 8 of Paper I, asymmetry could occur across the lobes.
Real pulses are varying incoming wind velocities or mass-loss rates from the base wind launching region, likely caused by episodic accretion and ejection in the disk or magnetic star-disk interaction, which will form their respective projected velocity centroids on PV diagrams. The expected amplitudes and frequencies of the episodic ejections may be different from those of the pseudopulses. The new ejections will make multiple larger direct shifts in the projected jet intensity peaks near the base as distinct knots or blobs, rather than occurring post-RS at some large distances. The shapes and velocity patterns are best manifested in Figures 7  and 8 of Wang et al. (2019), where the ejected knots can jump around distinct velocity centroids, interact and coalesce to form large-scale patterns, as in the case of IRAS 04166+2706 (Tafalla et al. 2010). We note that in a real series of pulses, a new reverse shock cavity forms near the base with each single ejection, and due to the density profile, the perturbations appear strongest on the axis rather than on the wide-angle portion. Details of the physics will be explored in depth for the coverage of parameter space at high numerical resolutions in follow-up publications.

THE TRANSVERSE POSITION-VELOCITY DIAGRAMS OF COLUMN DENSITIES
In this subsection, we explore the PVDCDs from cuts made perpendicular to the jet axes. Such PVDCDs are particularly interesting in probing the kinematic features caused by asymmetry in flow kinematics, especially those due to rotation. Despite the essential absence of rotation, the internal structures of the axisymmetric elongated magnetized bubbles can carry peculiar kinematic signatures revealing the nested layers of multicavities in the transverse PVDCDs generated.
The bubble signatures in these perpendicular PVD-CDs serve as nonrotating baselines to avoid misinterpretation of the observational patterns. We note that the magnitude of poloidal velocity dominates the major features explored here, and the presence of rotation will not alter these results. On these scales, effects of rotation from launch and those of the collapsing rotating toroids will be small. These may likely complicate the signatures with small systematic skews or twists in the final patterns in real maps. The patterns displayed by the PVDCDs vary according to the internal structures traversed and sampled along the line of sight with each of the cuts. The maps are sensitive to the loci of the cuts on the jet axes in spite of the systematic self-similarity.

Perpendicular PV Maps Made with Multicavities
We illustrate the situation using the cases of strongly magnetized α b = 1 and nonmagnetized α b = 0 ambient media with cuts made at three different positions in the self-similar coordinate z/v 0 t with widths of 0.05v 0 t. The line-of-sight inclination is chosen perpendicular to the jet axis. Each cut passes through different parts of the outflow structure. The jet feature is concentrated to = 0 and v = 0 as a red dot surrounded by concentric shadows around it.
A cut passing through the lower portion of the outflow structure at z/v 0 t = 0.05 encounters the densest part of the multicavities, which include the compressed ambient material that ends with the FS of largest -extent as the elongated narrow rings along the position , and the compressed wind regions that form ovals gradually extending horizontally along the v/v 0 direction. The cut also passes through the lines of sight that will penetrate the lowest part of the primary RS cavity that contains the fast-moving wide-angle free wind with largest velocity dispersion. This is reflected in the largest velocity width (in v/v 0 ), as shown in Figure 18 (a) (for α b = 1) and (b) (for α b = 0). The very faint CDIE at this largest width at ±v/v 0 gives the sense of the free wind.
The cut made at z/v 0 t = 0.25 (Figure 18 (c) and (d)) has less dispersion in v/v 0 as the lines of sight usually traverse the upper portion of the free wind, which has mostly upward-pointing velocity vectors, and always traverse the compressed wind regions. The velocity vectors are deflected after they pass through the reverse shocks. The traversed regions cover a significant fraction of the multicavities, which are the sites of the instabilities, the compressed wind, the compressed ambient medium, and outer FS. The ring contributed by the compressed ambient medium is elliptical in , and ap-parently wider in v than at the lower z-cut. Another ellipse, corresponding to the compressed wind, is even more extended in v, and may be visible for larger n cases.
Moving further up in z/v 0 t, the lines of sight traverse the middle of the outflow lobe at z/v 0 t = 0.5 (Figure 18 (e) and (f)). At this height, the lines of sight pass through the pseudopulses, wind-bubble edge, multiarcs, and instabilities, integrating them into the generated PVDCDs. Ovals corresponding, respectively, to the compressed ambient medium and compressed wind are present in our images. For larger n, both ovals drop in CDIE with respect to lower heights, even more so for the ambient oval, which is also stretched in thedirection. For smaller n, the two ovals become virtually indistinguishable.
The outer edge of the compressed wind is visible in the PVDCDs as an oval region. At lower heights, some of the free wind is present in our cuts, and its feature in our PVDCDs is often narrow in and very extended in the v-direction, and typically corresponds to low (sometimes very low) CDIE values. In many practical cases, only the higher CDIE portions or features will be observationally detectable.
One signature that emerges when one compares the α b = 1 and α b = 0 cases is the feature of the "gap" region threaded by the compressed strong ambient poloidal field. The prominent CDIE peaks protrude out on the -axis for α b = 1, but are absent for the α b = 0 cases. This feature appears around the zero velocity. It reflects the compressed ambient regions that are cushioned by the ambient poloidal field near the base, projected toward and away along the line of sight. This is most evident in the configurations n = 4 and n = 6, for all M A values; although for smaller n, it may be difficult to distinguish it from the wind ovals.
The restriction to brighter CDIE parts still makes the compressed wind and ambient features valuable for PV observations, applicable both to systems with and without rotation. For example, Figure A.2 of Louvet et al. (2018) shows the expected characteristics in a series of transverse PV diagrams of the edge-on system HH 30. These panels, if rotated 90 • to show velocity axis horizontally and position axis vertically, shown from y = 0. 35, then to y > 1 , and up to y ≈ 2. 5, appear to follow the general features of the height variations of compressed materials shown in the different cuts of transverse PVDCDs for the relevant n 6. They vary from filled oval shapes to relatively unfilled ringlike shapes as one examines from the bottom to about the mid-height of the outflow lobe. Figure 18. Position-velocity diagrams of column density (PVDCD) perpendicular to the jet axis at a location of z/v0t = 0.05. The PVDCDs are shown at an inclination angle of i = 90 • for (a) α b = 1 (top panels) and (b) α b = 0 (bottom panels) cases with n = 1, 2, 4, and 6 toroids (from left to right) and MA = 6, 30, 60, and 90 winds (from top to bottom). The unshocked ambient material with vp < aa has been filtered out while producing the PVDCDs. The rescaled logarithmic column density covers 4 orders of magnitude. For each n-column, an associated exaggeration factor n magnifies to show detail structures ( n = 7 for n = 1, 2, and n = 5 for n = 4, 6). PV coordinate axes are shown in black dotted lines. Figure 18. PVDCDs perpendicular to the jet axis at a location of z/v0t = 0.25, shown at an inclination angle of i = 90 • for (c) α b = 1 (top panels) and (d) α b = 0 (bottom panels) cases with n = 1, 2, 4, and 6 toroids (from left to right) and MA = 6, 30, 60, and 90 winds (from top to bottom). The unshocked ambient material with vp < aa has been filtered out while producing the PVDCDs. The rescaled logarithmic column density covers 4 orders of magnitude. For each n-column, an associated exaggeration factor n magnifies both and v to show detail structures ( n = 1.6 for n = 1, 2, and n = 1 for n = 4, 6). PV coordinate axes are shown in black dotted lines. Figure 18. PVDCDs perpendicular to the jet axis at a location of z/v0t = 0.5, shown at an inclination angle of i = 90 • for (e) α b = 1 (top panels) and (f) α b = 0 (bottom panels) cases with n = 1, 2, 4, and 6 toroids (from left to right) and MA = 6, 30, 60, and 90 winds (from top to bottom). The unshocked ambient material with vp < aa has been filtered out while producing the PVDCDs. The rescaled logarithmic column density covers 4 orders of magnitude. For each n-column, an associated exaggeration factor n magnifies both and v to show detail structures ( n = 1.6 for n = 1, 2, and n = 1 for n = 4, 6). PV coordinate axes are shown in black dotted lines.

A Schematic Two-Shell System
We further clarify the situations conveyed in the perpendicular cuts made at different heights by a schematic illustration in Figure 19 for the expected appearance of PV density diagrams for a multiple shell. We start from two shells by drawing two nested cones with poloidal velocity higher for the inner one, and sketch the expected PVDCDs for an imaginary slit perpendicular to the flow axis. Panels (a) and (b) in Figure 19 illustrate in 3D the two shells simplified as cones, cut by a slit at constant height z and viewed along the x-axis, where the resulting cut is shown in 2D position space in panel (b). The two curves produced by the cut are represented by two circles, inner and outer, labeled as 1 (blue) and 2 (red). Here, the position coordinates on these circles are simply called x = cos φ and y = sin φ. The velocity component of this nonrotating flow in the plane of the panel is v = v r sin θ + v θ cos θ, often v ≈ v r sin θ for the flows considered here, v r v θ . Projecting the circles in panel (b) along the line of sight into the position-velocity space leads to the situations in panels (c) and (d) of Figure 19. Here the position coordinate (labeled as y) is y = sin φ, and the velocity coordinate is v x = v cos φ, which give the mathematical expression of an exact ellipse with the corresponding semiaxes and v . The ellipses would cross each other as in panel (c) or not, as in panel (d), depending on whether the velocity component v 2 of the outer red cone is smaller or larger than that of the inner blue cone, v 1 . The function v ≈ v r sin θ can grow with θ because of the sine factor (leading to case (d)), but it can also decrease with θ when the poloidal velocity decreases fast enough from the inner to outer parts inside the compressed flow region (leading to case (c)), a typical effect in a compressed wind bubble. Both cases may be observed in perpendicular PVDCDs, for different slits, cases, or features. The inner (blue) and outer (red) cones may be separate features, in which case the final PVDCD is either as two crossed ellipses • • in panel (c), or as two concentric ellipses • • or • • in panel (d).
We also consider the alternative case in which the two cones are the inner and outer end of a continuous structure of cones, where v is either monotonically decreasing with , or monotonically increasing. Panels (e)-(h) in Figure 19 explore different cases, with a green dashed ellipse exemplifying the cones intermediate to the blue and red boundaries. The overall shape summing the set of intermediate ellipses is indicated with light-gray dashed lines of various shapes. If v decreases with , the intersecting blue and red ellipses will be located as in panel (c). The intermediate ellipses may be rather small in PV space, leading to an overall shape enveloping the crossed blue, green, and red ovals in panel (e).
The intermediate ellipses may be rather large, leading, in panel (g), to the overall shape of a bulging oval enveloping the set of large intercrossing ellipses. Panel (f) shows the intermediate case for which the envelope of the ellipses takes the shape of a rhombus. For the situation in which v increases with , the blue and red ellipses will be located as in panel (d) in Figure 19. The intermediate ellipses would also be noncrossing and concentric, leading to the overall shape of a set of concentric ellipses shown in panel (h). Within the overall shapes described in (e)-(h), a hollow region may exist, large or small according to the behavior of the set of ellipses at small values of the PV coordinates.
Examples of these overall shapes are present in different features and panels of Figure 18. For example the rhombus shapes are present in Figure 18(c) and (d) for most panels with n = 2 and nearly all panels for (d). The intersecting ellipses or ovals are present in Figure  18 In real observations, some of these ellipses, crosses, and rhombuses may be seen only in parts due to natural limitations of the slit length or width, and of the excitation conditions of different parts of the source. For instance, the shapes like parentheses '( )' in Louvet et al.
(2018) may be explained as ellipses clipped by limits of position or velocity coverage. Inclinations different from 90 • may also introduce complexity of shapes due to projection effects. The presence of v φ components in either wind or ambient material also leads to a more complex picture, such as in the ellipses tilted with respect to the PV axes shown in the diagrams of Appendix B in Tabone et al. (2020), or other factors such as asymmetry in the flow structures.

OBSERVATIONAL IMPLICATIONS AND DISCUSSIONS
In what follows, we discuss the theoretical interpretations and the implications of our results, in the context of existing models and of recent observational results.

The Thin-Shell Models
We start with the thin-shell models. In our Paper I and this work, we have demonstrated important internal structures as bubbles exist inside the outflow lobes. This holds true for the thin shells, which were the stateof-the-art for the hydrodynamic regimes of the smooth model assumptions.
In generating the hydrodynamic MC curves, we utilize a hydrodynamic deduction based on mass and momentum conservation inside a thin shell, but for a general ratio of wind-ambient mass densities. Despite the different deduction, the first appearance of Equation (12) Figure 19. Schematic illustrations of transverse PVDCDs. Panel (a) is a 3D sketch of a transverse slit cut through an outflow lobe having two nested cones (representing outflow features). Blue and red circles mark the intersections of 3D cones with the plane cut. In this figure, the outflow axis is along z, the line of sight is along x, and y is the position axis of the transverse PVDCDs. In panel (b) the two cones are seen on the plane of the cut. The inner and outer cones are, respectively, cut at the cylindrical radii 1 and 2, with velocity components v 1 and v 2. The transverse PVDCDs of (b) are shown in two cases: (c) v 1 > v 2 (two intersecting ellipses) and (d) v 1 < v 2 (two nested nonintersecting ellipses). Panels (e)-(h) show the results of adding a hypothetical cone (green dashes) representing a wide feature spreading from blue to red (details in Section 8). The addition of the green series leads from (c) either to (e) crossed ovals, (f) a rhombus, or (g) a single oval envelope of intercrossing ellipses, and then from (d) to (h) concentric ellipses. The total sum of the blue-green-red series is indicated in light-gray dashes.
is notably shown with Equation (12) in Schwartz (1978), starting from the discussion below Equation (7) in Mc-Kee & Cowie (1975), and was nonetheless attributed to McKee & Cowie (1975). Likewise, Equations (13)-(17) in Lee et al. (2001), in which wind mass is included into SRLL, lead to the same expression Equation (12) as discussed above. These thin-shell MC curves are often good predictors of the simulated lobe shape curves (including some of the cases in which the wind-bubble shell is not so thin). In this work, we utilize them to construct and to compare with HD models in Sections 2.3.1, 4, 6.1, and Appendix C. We also use them as a baseline to compare with simulations of PVDCDs for the HD cases in Sections 7.1 and 7.2.
These thin shells with the locally oblique structures we implemented in Section 4 serve as the HD limit for bubbles with smooth and laminar flow across the shocked regions bounded by the RS and FS without the onset of KHI. These treatments are fundamentally different from driving an outflow by entrainment, and the same applies to all of our MHD models in this work.

The High-, Intermediate-, Low-Velocity Components
In Section 1 and the Introduction of Paper I, we review the observational evidences of the jets and outflows from the deeply embedded protostars to the optical jets of the revealed T-Tauri systems, based on which our theoretical framework has been built and developed (references therein). Notably there are those very collimated molecular jets and outflows with prominent SiO and CO jet emissions from Class 0 protostars (e.g., HH 211, HH 212, L1448C(N), and IRAS 04166+2706), while losing the SiO emissions before entering the Class I phase. The wide cavities, the highly collimated shapes, and large terminal and internal bow shocks establish the foundations of the outflow paradigms.
Kinematically, the features of the wind-driven and jetdriven paradigms of the protostellar outflows are best summarized in the velocity components. The jet-shell signatures are identified as the extremely-high-velocity (EHV) and low-velocity (LVC) components with the jet signatures associated to the EHV, and the thin shells to the LVC, respectively. Between the EHV (∼ 100 km s −1 , used interchangeably with the HVC) and the LVC, there exists the standard-high-velocity (SHV) component (a few to ∼ 20-50 km s −1 , also as the IVC here), and extended LVC (a few kilometers per second down to systemic velocity) that are often present in Class 0/I sources (e.g., Bally 2016;Bachiller 1996). These are what may link to the signatures of the magnetic interplay identified in this work (Section 7). All of these components appear in our predicted parallel and transverse PVDCDs as signatures of jet-bearing elongated magnetized bubbles.
Recently, evidence has been reported of another phenomenon, the so-called "slow molecular winds" arising in Class 0/I sources at scales of 2000 au, and apparently "rotating" in the same sense as their gaseous disks. They are often inferred to root at radii ∼ 10-100 au from the protostars, e.g. the review of Pascucci et al. (2022). They are often associated with cavities or shells nested like onions. These systematic characters match well with the morphological and kinematic features of the elongated magnetized "bubbles", especially the multicavities in Paper  Some T-Tauri sources also have visible multicavity structures surrounding the central optical/IR jets, by CO/H 2 nested cones, with "nested" velocity structures. We note the presence of the nested multicavities in the source DG Tau B (de Valon et al. 2020) in Paper I, and connect them with the image contrast in the column density maps (Figure 18, Paper I). It displays onionlike cavities inside-out, with a gap potentially filled with magnetic field lines, as a more evolved version from that shown in HH 212 (see Section 9.5). Similar structures exist in DG Tau A (Agra-Amboage et al. 2014; Güdel et al. 2018) andHH 30 (Louvet et al. 2018). The slower molecular cones surrounding the jet may reveal the compressed wind as their flow directions deflected (Section 5.2) and velocity magnitude reduced after passing the RS. Combined with the multicavity morphology, the nested jet-slower cone-cavity layers do support the structures of magnetized bubbles produced by a jetbearing (unseen) wide-angle (X-)wind from the innermost regions.
It may not be surprising that the magnetized outflow bubble phenomena can last into the Class II phase. Remnant cavities in H 2 are known to exist, around T Tauri stars, such as FS Tau B (Liu et al. 2012) and T Tau (Saucedo et al. 2003;Beck et al. 2008). Multiple velocity components do exist in these more evolved systems, also identified as the high-velocity component (HVC, 100 km s −1 ), and the low-velocity component (LVC, 30 km s −1 ) in optical/IR forbidden lines (Bachiller 1996;Bally 2016). The LVC as a whole is usually interpreted as the presence of a lower-velocity "(extended) disk wind" launched farther away from the central source (see further discussions in Sections 9.3 and 9.4; Hartigan et al. 1995;Kwan & Tademaru 1995;Shang et al. 2007). The HVCs are connected to the microjets of the T Tauri stars. However, the origins of the LVCs and their connections to the HVCs have remained evasive and controversial. Our results here can provide insights.
In recent years, some LVCs, seen in atomic [O I] λ6300, can be further decomposed into two distinct components (e.g., Rigliaco et al. 2013;Simon et al. 2016;Fang et al. 2018;Banzatti et al. 2019), one broad (BC), and one narrow (NC). The two-peak profiles of the LVCs are more frequently associated with higher accretion rates (Banzatti et al. 2019). The BCs and NCs are seen as two distinct kinematic features in the systems, often referred to as "slow winds" (Pascucci et al. 2022) in a similar fashion as those from the Class 0/I sources mentioned earlier. The BCs can have highest velocities of 30-40 km s −1 , and have broader widths than the NCs. Conventional application of Anderson et al. (2003) has led to interpret the BCs as arising from inner (< 0.5 au), and outer as 0.5-5 au, in disk radii, with the former being faster (see discussion in Section 9.4).
Like the Class 0/I counterparts, these broad LVCs may be the kinematic manifestations of the more evolved (larger n) magnetized bubbles, but in different atomic tracers. The appearance of LVC-BC correlates strongly with the HVC. As disks disperse, in [O I], the HVC fades away eventually as LVC-BC becomes nondetectable (Simon et al. 2016;Ercolano & Pascucci 2017;McGinnis et al. 2018;Banzatti et al. 2019). The LVC-NC persists. The connection to the nested velocity structure appears to persist in this phase. If the HVC is identified with the jet, the identified broad LVC(-BC) may arise from the compressed wind across the RS from the bubble interior. The LVC-NC persists as the compressed and entrained ambient material all the way out to the FS.
The DG Tau A system is extensively studied in the optical and near-infrared wavelengths for its multiple velocity components and line profiles.  (Hartigan et al. 1995;Bacciotti et al. 2000;Agra-Amboage et al. 2011;Liu et al. 2016), and the [O I] LVC can be further decomposed into a BC and an NC (Simon et al. 2016;Banzatti et al. 2019). Morphologically, these ve-locity components are encompassed by an even lowervelocity wind (of 10 km s −1 ) traced by H 2 2.12 µm (Takami et al. 2004;Beck et al. 2008;Agra-Amboage et al. 2014), with an apparent trend of decreasing vertical extent, followed by the CO J = 2 − 1 emission with a speed of ∼ 3.1 km s −1 (Güdel et al. 2018). These evidences suggest they manifest the multicavities.
In our framework, we can reproduce the kinematic features seen in systems from Class 0, I, and II phases as different velocity components in tracers ranging from molecular lines to forbidden lines. The known observational features are consistent with the natural signatures produced by jet-bearing wind-driven magnetized bubbles demonstrated in this work.
Our framework offers a different but self-consistent interpretation of the inferred presence of outflowing gas with slower velocities confined within multicavities, the so-called "slower (molecular) disk winds," actually as evidence of the compressed bubble structures revealed.
In Sections 9.3 and 9.4, we discuss conventional assumptions and approaches adopted to explain the observational data, and their limitations and uncertainties. The magnetic interplay of the primary jet-bearing X-wind or an X-wind resembling inner disk wind (IDW) from the innermost regions, with the ambient media and, or with an outer Extended Disk Wind (EDW) would lead to extensive shocked and compressed regions, making the hypothesis of existence of an outer EDW less solid.
Note -The new JWST image of L1527 is one that captures the large KHI fingers generated by magnetic interplay when a magnetized wide-angle wind interacts with its surrounding ambient toroid. The structures resemble those produced in the evolved cavity formed with a moderately magnetized wind (M A 60) interacting with an n 6 toroid in a mildly magnetized (α b ∼ 0.1) ambient medium. These structures originate deeply from the outflow base as part of the magnetized bubble, grow, and advect to large scales with the underlying flow.

The Disk Wind Scenarios
Historically, the general physics of magnetocentrifugal launch of winds from disks is well recognized as proposed in Blandford & Payne (1982), the BP mechanism. This is the original definition of the term "disk winds" adopted for the context in star formation of the earliest phases (Königl & Ruden 1993;Königl & Pudritz 2000), and in reference to the X-wind models (Shu et al. 1994a(Shu et al. ,b, 2000.
Nowadays, the term "disk winds" as presented has become analogous to phenomena of outflowing gas from the general protoplanetary disks or disk-like structures of all scales and of all evolutionary stages. They facilitate transport of angular momentum from the disk to the wind that allows or helps disk accretion. This term has become too loosely defined to convey specifics of basic physical mechanisms.
The BP mechanism applies in principle to all radii on the disk from the inner to the outer parts, because the effective gravitational potential (gravitational plus centrifugal) is able to accelerate gas along inclined poloidal field lines at any given disk radius (Figure 1 in Blandford & Payne 1982). The effective potential for gas attached to a poloidal field line has a gradient that accelerates gas in two regions between the equator and an angle of 60 • from it, allowing the gas to launch, while for angles closer to the vertical than 30 • , the effective potential has the opposite gradient sign and does not allow magnetocentrifugal launching. Once launched, the mechanism allows gas to accelerate due to the projection of the centrifugal force onto the poloidal field line attached to the gas (typically with some additional push by projected gradients of gas pressure and of b φ magnetic pressure). These processes can happen in principle at any radius, either inner or outer, including in particular the inner disk radius. Gas in open field lines of the acceleration region leaves the disk into larger values and eventually forms part of disk winds of different categories discussed below.
Field lines that are loaded with matter (as allowed by gas-magnetic attachment and field line angle) can happen at any location. Disk-wind scenarios are classified according to the launching location along the disk. The range of launching locations can be extended and the mass loading be distributed along that range, leading in this case to an extended disk wind (EDW) scenario. Alternatively, mass loading can be either limited to a narrow inner range of locations or have a mass loading concentrated toward the inner edge of the disk, leading in both cases to an inner disk wind (IDW) scenario (including the wind portion of the X-wind model as a particular case). Both scenarios are equally able to participate in the BP mechanism. The general EDWs possess a radius-dependent profile of their terminal veloci- (Blandford & Payne 1982;Shu et al. 1995) where the ratio ( A / 0 ) depends on the specific model. This ratio, the terminal velocity and their respective density profiles are related to the local mass loading onto each of the individual streamlines which characterize a given disk-wind model. Different from the EDWs, the IDWs have finite inner and outer disk radii, launched from small inner disk regions that can reach high escape velocities to account for those observed in jets, and usually identifiable with the HVC.
The X-wind models launch a disk wind from a region located close to the inner edge of the disk, connected also to the stellar magnetosphere (Shu et al. 1994a;. Therefore, X-wind models launch the disk wind from the deepest part of the gravitational potential in the system and can drive jets reaching the highest terminal velocities. Note that both the X-winds and the IDWs in their original forms have their respective velocity profiles as shown in Shang et al. (1998), Shu et al. (2000), Shang et al. (2000Shang et al. ( , 2010, Liu & Shang (2012), and, e.g., Krasnopolsky et al. (2003) and Anderson et al. (2005). The details of the velocity profiles pertaining to the launch will be reflected in the associated HVC. The formulation shown in Paper I is a simplified solution for lightweight computation for otherwise complex physics and expensive cost for the focus of the larger-scale bubble structures they form.
In IDWs, asymptotic behavior for large z with vertical isodensity contours is regularly observed, as seen in Krasnopolsky et al. (2003), Figure 6, and in Anderson et al. (2005), Figure 12 (dash-dotted line). Simulations having mass loading concentrated toward the inner radius approach nearly vertical isodensity contours for all z, and in that regard, they resemble X-wind disklaunching models.
In the context of magnetized gravitational collapse, magnetic braking can slow down rotation of material in an infalling core through the emission of torsional Alfvén waves to outer parts of the environment (Basu & Mouschovias 1994;Krasnopolsky & Königl 2002), in some circumstances leading to values of v φ problematically small for large disk formation (the "magnetic braking catastrophe"; e.g. Mellon & Li 2008, 2009Krasnopolsky et al. 2010). The mechanism operates through poloidal field lines connecting an outer environment having slower angular rotation to a magnetically braked inner part. In principle, this can include inner parts of collapsing cores, infalling pseudodisks, and also parts of the disk proper.
At a later evolutionary stage, when the disk can evolve self-consistently, a new generation of application of the EDW in the context of "wind-driven accretion" is gaining popularity for the consideration of angular momentum transport in the disk. The magnetorotational instability (MRI) can generate the anomalous viscosity required to transport angular momentum within the disk by the α-disk model for accretion. The alternative mechanism utilizes the capabilities of the magnetic torque in EDW to carry angular momentum out of this disk in place of the MRI. These models are considered with dif-ferent ingredients of the nonideal MHD terms for their respective physical mechanisms (e.g. Bai 2016). Alternatively, these EDWs are connected to the spontaneous formation of structures (Suriano et al. 2017(Suriano et al. , 2018(Suriano et al. , 2019, leading to varieties of rings and reconnection islands in the planet-making disks. These lighter EDWs share some common features that differentiate them from the original cold EDWs, as they require help from thermal effects and large magnetic pressure for the successful launch at larger radii. Other flows within the new "disk-wind" category can also be launched thermally by photoevaporation, "photoevaporative winds" (Hollenbach et al. 1993(Hollenbach et al. , 1994Gorti et al. 2009Gorti et al. , 2015, from the extended radial regions farther away from the central source out of a transitional or debris disk at a much later evolutionary stage. These photoevaporative flows (mostly features of flow close to the upper disk layers) are not large-scale winds in the same sense as IDW or EDW. Additionally, because they are not magnetized winds, they cannot extract angular momentum from the disk.
Regardless of the underlying wind-driving mechanism, the basic structure of free wind, wind shock, shocked wind, contact, shocked ambient medium, and free ambient medium, should apply.

Inferences of Slower Disk Winds Launched from an Extended Disk Region
In this subsection, we discuss the theoretical basis, validity, and limitation of the conventional approach to infer the presence of extended disk winds (EDWs). We will make cautionary remarks for the inferred launch radii of an EDW at large, and the presence of the socalled "slow molecular winds".
The standard formulation of the magnetocentrifugal winds follows a number of conservation laws, which, in axisymmetric steady state, lead to a number of conserved quantities, in the theory of the BP mechanism (Blandford & Payne 1982). Applying the conserved quantities at relatively large distances within the wind would allow, in principle, to infer knowledge about conditions at the launching point, as shown in Anderson et al. (2003), provided strict requirements are fulfilled. These requirements, especially regarding accurate data on rotation velocity, have often been ignored when the methodology of Anderson et al. (2003) has been utilized. The methodology is based on the "fieldline angular velocity" Ω ≡ (v φ − b φ v p /b p )/ , and the Bernoulli constant H ≡ v 2 /2+h+Φ g −Ω v φ (e.g., Shu et al. 1994bor Lovelace et al. 1986), both conserved along field lines of an ideal axisymmetric steady-state pristine wind. Under those conditions, the values of the Bernoulli constant H are equal at the observing and launching points P ∞ and P 0 . A number of simplifying assumptions is then applied (cold wind allowing setting enthalpy h = 0, Keplerian value v K = (GM * / ) 1/2 = (GM * Ω 0 ) 1/3 for both Ω and v φ ( v p ) at the launching point, and negligible Φ g at the observation point), together writing an easily solvable cubic equation giving Ω 1/3 0 and the location as 0 ≈ (GM * /Ω 2 0 ) 1/3 . Another assumption of validity is that the input data must be purely from the free wind, without involving the ambient medium directly. This is because the ambient data are not connected to the conserved quantities. That is, the region where the data are used and extrapolated, must be far enough from the region, which has strong interaction with the surrounding medium. Avoidance of the known shocked region should be exercised. The interaction may break axisymmetry of rotation confusing the measurements of rotation needed to apply Anderson et al. (2003) if wrong assumptions are taken.
Application of Anderson et al. (2003) requires resolved pristine-wind data so as to trust values of ∞ and v φ,∞ . Therefore, the use of data from under-resolved sources (Pascucci et al. 2022) can lead to uncertain values of 0 , depending on the details and purpose of data management. Since typically v φ,∞ v p,∞ , distinguishing these two components in order to measure v φ,∞ can be sensitive to the estimate of source inclination i, and to minor asymmetries in the source. Interaction and interplay between the magnetized wind and the ambient medium can reach much further inside, up to the RS. All models, without exception, are subject to the RS-limitations in comparing a pure-wind model against observational data unavoidably in interaction with the environment. This not only applies to the interactions with the ambient media but also to the interactions between the inner and outer winds if they coexist.
Observational data obtained between the RS and the FS can be too severely modified by the complex interplay to show the pristine conserved quantities along streamlines that are necessary to apply Anderson et al. (2003). This can be observationally very demanding, as toroidally magnetized winds often have very thick shells moving the RS close to the origin, and presenting vortex structures caused by magnetization and by oblique shocks in parts of the wind region. Time variabilities can also break the assumptions for the valid application of Anderson et al. (2003). New shocked structures will alter the pristine environment required for the validity of the approach.
For the reasons presented above, identifying the observed shells with an EDW appears to be doubtful. On the contrary, the model presented in this work could justify the observations without the need of assuming a separate wind launched from the outer regions of the disk. We illustrate these points further with the ALMA observations obtained for the HH 212 in Section 9.5.

The HH 212 System
In Paper I, the extensively studied HH 212 is inferred to be a candidate for a wide-angle wind interacting with a somewhat opened environment of 2 n 4, given the shape of the opening near the base, and the presence of an SiO jet and CO cavities. The molecules SO and SO 2 appear to also originate from the same location where the SiO jet appears. Lee et al. (2018) resolved the SO emission into a collimated jet and a wide-angle rotating shell. They appear to be spatially surrounding the highly collimated, cylindrically narrow SiO jet. The slower SO emission appears to be tracing the wide-angle wind, with SO 2 , connected to its inner densest core surrounding the highly excited SiO J = 8 − 7 transition (Lee et al. 2018).
In Lee et al. (2021), additional shell structures are revealed in SiO and SO within 1400 au of the central protostar at 13 au resolution. In their Figure 1, the relative relationship between SiO and SO is shown, in which the SO is indeed aligned in jet with SiO but distributed in a wide-angle way more broadly around the core near the base, up to ±3 km s −1 from the systemic velocity. The spatial distribution confirms the suggestion in Paper I. Furthermore, the inner jet core is clearly enclosed by an extended shell, further surrounded by a wider more opened thick butterfly-like SO shell at 100 au scale (shown in their Figure 2). The cross-connected spatial morphologies in SiO and SO illustrate the nested-shell or multicavities identified in Paper I, and as summarized in our Figures 2(a), 3(a), and 5. The thick SO shell protruding from the base corresponds to the "halfwall" seen in the lower half of the column density passing through the "gap" region threaded by the compressed ambient magnetic field lines in Figure 5.
These features match well with the predicted structures from our framework, and serve as a manifestation of the physical mechanism elucidated. The nested SiO and SO emission structure of HH 212 is seen as an inner elongated shell surrounded by an X-shaped cavity with outward pointing velocity vectors. This situation is manifested in Figure 7 of this work. However, in Lee et al. (2021), such structure, as shown in their Figure 2, is interpreted as evidence of a jet interacting with a surrounding disk wind launched from 4 au all the way to 40 au, following a previously reported slow wideangle wind in SO and SO 2 by Tabone et al. (2017). Lee et al. (2015) showed prominent "multicavities" of CO J = 3 − 2 "nested shells" in channel maps (their Figure 1) and integrated intensity maps (their Figure 2). Lee et al. (2017) revealed a flattened envelope traced by HCO + J = 4 − 3 in the absolute velocity range ∼ 0.5 to ∼ 3.0 km s −1 away from the systemic velocity (their Figure 1).
Our work can accommodate all of these observed features without the need for the presence of a separately launched disk wind from 4-40 au, in one unified framework. The structure can be interpreted as part of the multicavities interior to the compressed magnetized ambient envelope (the cocoon). Their "inferred" SO/SO 2 disk wind shown in Figure 3 of Lee et al. (2021) in fact can be interpreted as the compressed wind outside the RS cavity, whose flow velocities in the post-shock region have been reduced and deflected by the crossing of the RS at different local oblique angles, showing an optical illusion of direct launch from the disk surface. The outwardly pointing but decreasing velocities with values much reduced than those of the primary wind, extended to the compressed gap region threaded by strong ambient field lines, will result in a range of apparent launch by using the formulation of Anderson et al. (2003), as discussed in Section 9.4. The outer HCO + labels the infalling material outside the influence of the expansion, which represents the background infalling envelope.
We note that the overall SiO parallel PV diagram of HH 212 as reported in Lee et al. (2017) Figure 5 of Lee et al. (2022), the velocity at the base of the SiO emission is broad, with an MS-type shell outlining the compressed wind across the RS-enclosed cavity, similar to our predictions. In Figure 1 of Lee et al. (2021), the relative positions of the SO HVCs and LVCs are shown with respect to SiO. The SO HVC traces the primary wind interior to the reverse shock, while the SO LVC traces the post-shock emission. Indeed, the radial wide-angle wind is seen through the detection of the wide velocity dispersion at the very base of the jet in Figure 1 of Lee et al. (2022), as predicted in this work. 9.6. Extended Mixing Regions in CARMA-7 The source CARMA-7 (C7) reported in Plunkett et al. (2015b) is a Class 0 protostar that has 22 noted knotlike episodic features in 12 CO. Some of these knots ap-pear in the high-velocity portion of the PV diagram, and the rest appear in the low-velocity portion. According to Figure 2 of Plunkett et al. (2015b), the low-velocity points are located along the "collimated" axis near the systemic velocity, and the high-velocity ones are those located near the velocity maxima. The high-velocity points are near the base of the "embedded" unresolved blue-shifted and red-shifted jets. As we demonstrated in this work, the high-velocity component at a ∼ 90 • inclination angle traces the contribution from the wideangle wind flowing along the line of sight, and the near systemic low velocity on the contrary arises from the density-collimated wind streamlines moving along the jet axis in C7 due to the effect of the inclination angle. This effect can be properly identified by comparison with its image in Figure 1 of Plunkett et al. (2015b). The overall outflow itself is very compact and jet-like near the base, indicating its youth in evolution (Paper I).
The reconstructed emission in Figure 2 of Plunkett et al. (2015b) is consistent with the PVDCD features within this work in that a cavity feature resembling the image appears in both the PV and PP spaces, enclosing the triangular-shaped emission on the PV at the baseline of the offset. The high-velocity portion of the PV appears as zigzag patches tracing the contour of a cavity reminiscent of the MS curves (Section 6) projected at the same inclination angle on the upper half (beyond ∼ 6 as shown in Figure 2 of Plunkett et al. 2015b). The jet portion of the data was missing because of the combined effects of 12 CO self-absorption and interferometric filtering (Plunkett et al. 2015a). Nonetheless, the velocity dispersion broadens toward the source (within ∼ 6 ) clearly, especially on the blue side, suggesting that the emission originated in a wide-angle wind, which gives rise to a velocity with a relatively large magnitude and with a projected direction along the line of sight at 86 • . The magnitude of 15-20 km s −1 is already significantly reduced from that of a fast free wide-angle wind that could have a velocity of ∼ 100 km s −1 as in a fast jet. Therefore, this emission mostly originates from the compressed wind to be become visible in 12 CO on top of a faint potentially high-velocity background. A lobe-like contour encompassing the systemic velocity is also traceable, in both of the physical lobes and in PV, revealing the nature of multicavities.
The apparent episodicity of C7 supports the magnetic interplay, as pointed out in Paper I. The presence of an extended region of entrainment was suggested in Figure  3 of the Extended Data in Plunkett et al. (2015b). The image shown in Figure 1 of Plunkett et al. (2015b) can be compared with the moderately magnetized wind cases in Figures 5 and 6. Within these complex structures, den-sity enhancement and induced variation in velocity components in the multicavities can produce zigzag boundaries as part of the manifestation of the KHI within the extended mixing layer. The images separated in HV and LV panels appear nested with the HVC enclosed by the LVC. The PV patterns of C7 without the central jet can also be interpreted as originating from the multicavity of the mixing and gap structures surrounding it.

SUMMARY
In this work, we illustrate the kinematic features of the unified model of jet-and wide-angle-wind-driven outflow described in Paper I and further developed here. The description makes use of synthetic Position-Velocity Diagrams of the column density (PVDCDs) formed in directions parallel and transverse to the elongated outflow, and explored in a wide parameter space. We demonstrate different contributions and expected patterns in controlled steps. We start from basic analytic formulations with simple locally parallel oblique shock structures, followed by hydrodynamical numerical thin-shell calculations, producing reference smooth PVDCD signatures. We proceed to the subsequent production of full PVDCD patterns from HD and MHD simulations that allow the development of KHI. We demonstrate the corresponding changes in kinematic signatures generated in full simulations from HD to MHD, covering the parameter space adopted for different strengths of magnetization in winds and ambient toroids of various mass-to-flux configurations. The produced gallery can serve as a guide for observational comparisons for diagnosing the presence of outflows as magnetized elongated bubbles.
Our results suggest that bubble structures dominate the outflows driven by magnetized winds, and their interplay with the environment. These structures are thick and extended, contrary to the conventional expectations arising in thin-shell models. One particular objective of the figure gallery is providing insights into the unavoidable bubble structures that result, with the jetbearing inner disk winds possessing a strong wide-angle component. Real observations, if not very well-resolved, will typically show information connected to the largescale nested-shell structure, and if better resolved, will reveal information in their respective multicavity layers and velocity components. Thus, observations require careful analyses to separate the contributions from individual signatures. These will have profound implications for outflow phenomena based on similar physical mechanisms of magnetized bubbles driven by magnetocentrifugal winds, as their basic principles hold in general.
Outflow bubbles are the ubiquitous, natural, and inevitable result of interaction between a wind and an ambient medium. They will be present, no matter the launching mechanisms, launching locations, or evolutionary stages. Interaction between the wind and the ambient medium provides the shape by solving the shock jump equations. The jet core is the densest part of the wind launched sufficiently close to the innermost region. Wide-angle wind is seen through the interacting shock layers, usually in the form of a compressed rather than pristine wind; however, it is often mistaken as a separate slower extended disk wind component. The more revealed and more resolved bubbles that formed in an environment with a much wider toroid of larger n 6 can give similar identification when signatures of the compressed interacting wide-angle wind are revealed in optical lines.
Here we summarize our highlights: 3. The basic jet-shell structures as seen in PVDCDs are shown in the hydrodynamic limit using different n values of toroids. The kinematic features are shown in PVDCDs with a 45 • inclination angle. This angle only serves as a demonstration where the contribution from the jet and the wide-angle wind can be separated the farthest to be identified in the PV plane.
4. The reverse shock has its own kinematic presence as the boundary across the change of behavior in the velocity widths of the wind projected on the PVDCD.
5. The crossing of the reverse shock from free to compressed wind can produce deflection or rotation of flow velocity vectors that give an impression of the existence of a different layer of flow. This effect of deflection can imitate the behavior of a disk wind that may be interpreted and extrapolated to originate from somewhere inside the shell but outside of the innermost launch point of the primary fast wind. Because of the reduced poloidal flow velocity and a modified direction that is less radial compared to the original, the extrapolated location on the Keplerian disk may be farther out. The deflected flow direction inside the compressed wind could be misinterpreted as a separate slower "disk wind" originating from an outer region.
6. Signatures of KHI, compressed wind and compressed ambient media, magnetized or not, strongly or weakly, are directly revealed through the PVDCDs. They are juxtaposed for comparisons and identification in figure panels in this work for source-to-source comparisons. The notable presence in fact appears patchy and episodic as extended threads of emission features. These could be alternatively interpreted as episodicity in ejections, out of a smooth and steady flow.
7. Magnetic pseudopulses are generated in the compressed wind regions. Pseudopulses appear and cooperate with the KHI fingers in an interlinked feedback loop of magnetic generation of vorticity. Jet pseudopulses are particularly evident in PVD-CDs as fluctuating jet velocities on top of the RS cavities, that appear projected onto the PV diagrams as mild quasiperiodic wiggles around the velocity centroids propagating along the axis toward the tip. Pseudopulses can be distinguished from real wind pulsations, which create their own respective velocity centroids and have reverse shocks near the base, and by the association of the pseudopulses with magnetized KHI and magnetized vorticity generation within compressed wind.
8. The high-velocity component is easily identified as the EHV jet component, while the lower-velocity components can arise due to the interaction of wind with its surroundings, as the compressed ambient media, vortices aligned along the cavities. Multiple velocity peaks can arise that have broad and extended features in between the extremely high jet velocity, and the traditional low-velocity shell or the extremely low velocity near the systemic velocity reflecting the complexity of the flow patterns.
9. The velocity components arising from the multicavities follow a nested onion-like structure from the inside-out, displaying a profile of the various components as the extremely high (EHV), (standard) high or high velocity (SHV/HVC), intermediate velocity (IVC), the broad (BC) and narrow (NC) parts of the low velocity (LVC), and the very low velocity.
10. The presence of the broad low-velocity profiles found with atomic [O I] forbidden line (Pascucci et al. 2022) in T-Tauri stars can also be attributed to originate in the compressed wind resulting from the interplay with a less-dense ambient environment.
11. Transverse PVDCDs afford identification of compressed wind/ambient features and establish a nonrotating axisymmetric baseline in the comparisons between simulations and real data.
12. Our predicted model images and PVDCDs can provide an explanation to the observed nested bubbles and their associated velocity components in sources such as HH 212, HH 30, and DG Tau A and B. This demonstrates the power of the methods and models utilized.
13. Apparent acceleration can arise in the tip region when the toroidal magnetic field in wind decreases toward the tip. This can be observed in the PVD-CDs for small enough M A values. The underlying physics will be followed in Paper III.
14. We review assumptions and caveats needed for the application of Anderson et al. (2003) to find the launch radii from the measurements of rotation in the flows. In this work we indeed demonstrate that the bubble interiors include compressed regions arising from the interaction with extended shocks. The applicability of the Anderson et al. (2003) formulation appears thus to be questionable in such winds, as here the condition of laminar ideal flow is not verified.
15. We classify our results as the interaction between a toroid and a wide-angle X-wind or X-wind-like IDW, giving rise to the specific form of density profiles that must be cylindrically stratified, at all vertical heights along the jet axis, and a wide-angle wind. Such winds belong to the X-wind Class. The X-wind naturally displays the strongest feature and has the widest angle in streamline distribution in the base when the launch region narrows to the neighborhood of the inner disk edge.
16. Regardless of the underlying wind-driving mechanism, the basic bubble structure of free wind, wind shock, shocked wind, contact, and shocked ambient medium should apply.
The successful extended comparison with outflow sources from Class 0 to Class II phases, such as HH 212, HH 30, and DG Tau A and B demonstrated in our image and position-velocity maps has suggested a clear trace of bubble signatures existing in the known observational data. Our framework provides a self-consistent and integrated understanding of the enigmatic presence of the lower-velocity components of the embedded and revealed outflows, which had not been properly explained in the literature since their discovery nearly 30 yr ago until our work presenting them as arising in the multicavities formed in the magnetic interplay. Our gallery of results can serve as a guide for observational probes for the physical and dynamical situations inside the evolution of jet-outflow systems from Class 0, I, to II phases in molecular to atomic forbidden lines. Krasnopolsky, R., & Königl, A. 2002, ApJ, 580, 987,  opening of the initial toroid R(θ) ∝ sin n (θ) for θ → 0 1, 2, 4, 6 α b ambient poloidal field strength factor Equation (8) 0, 0.1, 1 αρ initial toroid density factor Equation (7) 0 (without toroid), 1 (with toroid) D S tapering density factor for the toroid Equations (7,9) 2.25 × 10 11 g cm −1 (or zero) other setup parameters and chosen units rout position of the outer radial boundary 10 5 au G Newton's constant 6.67428 × 10 −8 cm 3 g −1 s −2 1 au approximate astronomical unit 1.5 × 10 13 cm dependent parameter of wind b0 wind toroidal magnetic field constant Equation (1) v0 √ D0/MA angular functions for the initial conditions of toroids and boundary conditions of winds R(θ) or Q(θ) ALS toroid density distribution function This system of steady-state equations (Equations (A8)-(A11)) admits (for the p = 0 case) solutions such that the three quantities (v r , ρ 2 , |b φ | ) are constant. These extremely simplified equations can be useful to describe the steady-state portions of wind, and have been used to design the boundary conditions for the wind at r = r in , a region in which steady-state can be achieved, despite the fact that the model is not completely cold in that region, it does feature zero rotation and b p = 0 for our wind model.

B. SEMIANALYTIC SHAPE FITTING OF THE MAGNETIZED BUBBLES
We adopt the following procedure to obtain the analytic expression for outflow shapes, leading to the final production of Figure 10. The algebraic basis of the procedure is the pressure balance equation (Equation (15)), a quadratic equation in the velocity v MS (θ) that gives the shape r MS (θ) = tv MS of our self-similar outflow fitting.
The outflow is surrounded by a magnetic cocoon (Figure 1) if b p is substantial, and can possess a tip if b φ is very strong, for M A 30 (see PVDCDs for M A = 6 cases, and stronger ones in Paper III). Pressure balance equations can be written in all cases, and for small magnetization values, the resulting shapes do not deviate much from those obtained by the unmagnetized momentum-conserving Equation (12) (see also Shu et al. 1991;Koo & McKee 1992b;Lee et al. 2001).
In Paper I, the ratio η = (b 2 p /2) / (ρv 2 p ) of local poloidal magnetic pressure to local total ram pressure is defined to measure the relative strengths of these quantities. Contours of constant η are estimators of the location of the edge of the outflow (e.g. η = 0.5 as experimentally adopted in Figure 14 of Paper I).
The α b = 1 cases contain substantial magnetic poloidal pressure, represented in Equation (15) as a term e p b 2 p /2. Here e p represents the amplification of the initial poloidal field due to compression. Figure 1 shows how the wind regions drive out their initial b p , compressing it into a fairly narrow magnetic cocoon surrounding the wind. Magnetic flux conservation allows us to estimate e p utilizing the cross-sectional area ratio between the cocoon+wind divided by the much narrower cocoon alone. This effect depends both on the parameter n and θ along the outflow shape, but for the purposes of computing Figure 10 we consider it as only a function of n (e p = 40, 20, 9, and 8 for n = 1, 2, 4, and 6, respectively).
The factor β(θ) adjusts the wind velocity term in Equation (15) for effects that either increase (> 1) or decrease (< 1) its strength. The tips for M A < 30 are fit by locally increasing β at small θ. Deceleration inside the magnetized extended structures formed by compressed wind and ambient medium decreases β. Thickening of these same structures due to their magnetization and KHI is represented as a modest increase in r MS and β. The adopted β values (Table  2) combine these reductions and increases, often staying very close to unity.
Equation (15) is a quadratic equation in v MS . For our experimental values of β(θ) and e p , the equation discriminant can become negative if ρ a < ρ w (small θ values) and b 2 φ /2 is large (small M A values). This occasionally happens in small tip portions of Figure 10, and for such cases, the v MS (θ) values are carefully extrapolated from adjacent regions without complex roots, obtaining excellent fitting results.
The selected numerical values of β(θ) are summarized in Table 2. Monotone cubic Hermite interpolation of β(θ) is followed for all other intermediate θ values.

C. CONSTRUCTION OF HYDRODYNAMIC SHELL STRUCTURE
In this appendix, we devise four steps to demonstrate structures formed by an oblique bubble approximated by a pair of locally parallel oblique shocks. The detailed procedures are laid out for the conceptual development. Table 3 compiles symbols and quantities used here. C.1.
Step 1: The Momentum-Conserving Shell The position of the momentum-conserving shell found analytically has been previously demonstrated to well capture the behaviors of the thin shells. The first step of the development is obtaining the solution of the momentum-conserving shell.

C.2.
Step 2: Locations of the Reverse and Forward Shocks In this step, we add the structure of an oblique bubble to the position r S of the momentum-conserving shell. We construct the locations of the reverse and forward shocks to the thin-shell curve.
We start by constructing a system of unit vectors. The unit vectors tangent and normal to the MC shell arê The unit vectors in the r-and θ-directions can be projected onto the unit vectors in the and z-directionsẑ and , resulting inr = sin θˆ + cos θẑ, andθ = cos θˆ − sin θẑ. The normal vector can also be similarly written aŝ n MC = n MC ˆ + n MC zẑ . The derivative dr MC /dθ ≡ṙ MC can be found numerically by finite differences at high θ resolution.
We illustrate the schematic configuration in Figure 20. The construction needs to choose a unit vectord at each point P MC = r MCr located along the MC shell. We will use either of the choicesd =n MC ord =r. The radial incoming wind velocity v w is projected alongd obtaining a set of values v d = v w (r ·d). A hydrodynamic 1D Riemann problem is then set, having, as its left-side state, this velocity v d with the density and sound speed of the wind, and, as its right-side state, the ambient density and sound speeds, with a velocity of zero, as follows: left side : ρ L =ρ w /r 2 s , a L = a w , p L = a 2 w ρ L /γ, v L = v w (r ·d); right side : ρ R =ρ a /r 2 s , a R = a a , p R = a 2 a ρ R /γ, v R = 0 .

(C14)
A value of γ = 1.0001 is set for the equation of state of the 1D Riemann problem, which is solved by the usual analytical methods, such as, e.g., in Chapter 4 of Toro (2009). This value of γ = 1.0001 has been shown to closely reproduce the results of two-temperature numerical simulations in Paper I. polynomial of order n used to fit the reverse shock position p n FS (zFS) = FS polynomial of order n used to fit the forward shock position Model choices normal modeld =p =nMC, f pd = 1 radial modeld =p =r, f pd = 1 hybrid normal modeld =nMC,p =r,f pd = (d ·p) −1 Construction of fields (Pg = P, zg, g ) generic point, with its cylindrical coordinates (Pproj, zproj, proj ) projection of Pg onto the RS, with cylindrical coordinates (P fproj , z fproj , fproj ) "projection" of Pg  The solution of this 1D Riemann problem has three discontinuities, left, middle, and right, with speeds s L , s M and s R . A second unit vectorp is now chosen, equal to either ofd orr, and a factor f pd is defined to be either equal to 1 or to (d ·p) −1 . The modeled positions of the reverse and forward shocks of the axisymmetric bubble (RS, FS) are then constructed as The construction completed by Equation (C15) allows an oblique 2D axisymmetric flow to be represented by a set of 1D flows having a direction of motion determined by the choice ofd (either normal or radial).
The left and right initial states (Equation (C14)) of those 1D flows are built from the wind (on the left side) and the toroids (on the right side), with the wind velocity projected in the direction ofd and applied to the left side. Based on these initial states, the discontinuity structure of the 1D flows at an arbitrary time t sim is found analytically, and it is expressed as the three speeds s L , s M , and s R , corresponding to the reverse shock, CD, and forward shock of the 1D flows. These three speeds are now used to construct locations for the RS and FS of the 2D bubbles. The middle speed s M is made to correspond to the location P MC of the MC shell, and the left and right discontinuity speeds are made to correspond to the RS and FS of the bubble. These correspondences are applied in the direction of the vector p, and with a scaling factor equal to the product of f pd and time: these notations are chosen for generality, but more frequently are applied simply asp equal tod, and the factor f pd equal to 1.
The respective shock velocities for the RS and FS are estimated as which is based on a time derivative of Equation (C15). For relatively thin shells, the approximation v RS = v FS = v MCr can also be applied. The positions P RS and P FS are often not used directly, but via a polynomial fit procedure, which gives as a result polynomials p n (z) (n = 8 as of now) such that the point FS = p n FS (z FS ) is on the fitted forward shock and RS = p n RS (z RS ) is on the fitted reverse shock.
We have chosen the following configurations for the constructions of the positions of the reverse and forward shocks: (a) The normal modeld =p =n MC , f pd = 1: represents the 2D bubble flow as a set of 1D flows moving in the direction normal to the MC shell, and constructs the RS and FS positions based on that, motivating the particular choice ofd =p, and f pd = 1 values; (b) The radial modeld =p =r, f pd = 1: represents the 2D bubble flow also as a set of 1D flows, but now moving in the radial direction. It is observed to lead to less separation between RS and FS than the normal model. The radial model is expected to be connected to the results of the v θ = 0 numerical models, which are also a set of radial flow; (c) The hybrid normal modeld =n MC ,p =r, f pd = (d ·p) −1 : is intended as a mathematical simplification of the normal model, with 1D flows modeled to be in the normal direction, but with their results applied in the radial direction with the help of a projection factor f pd approximately restoring the intended normal direction of the 1D flows. Such a hybrid approximation is closer to representing the normal model in regions of the bubble having small curvature and small thickness.
The shock surfaces constructed by these models will be utilized in Step 3. Figure 21 summarizes the constructed profiles for the positions and the thickness.

C.3. Construction of the Field of Shear Velocities
In the following two steps (3 and 4), we will construct a field of velocity vectors for the whole flow, both inside and outside the shell of the outflow bubble. With this construction, the flow structures across the RS and FS can be analyzed analytically in the post-shock regions between the RS and FS.
Basically, starting from an arbitrary generic point P = P g given in cylindrical coordinates as P g = z gẑ + gˆ , its gas velocity v gas (P g ) ≡ v(P g ) is estimated in a sequence of calculations as follows, using the fitted curves of RS and FS found in Step 2. Two trivial cases are checked first: if g > p n FS (z g ) then P g is well inside the toroids with v(P g ) = 0, and if g < p n RS (z g ) then P g is inside the pure-wind region with v(P g ) = v 0r . Otherwise, P g is inside the shell of the bubble, located between RS and FS.
Step 3 treats this case by finding two representative projection points for each P g , one at the RS and another one at the FS, as follows.

C.4.
Step 3: Projections onto the shock surfaces First, the point P g is projected onto the RS at the place closest to it. If d is the distance between P g and a point on fitted RS curve then the minimum of d 2 (as a function of z) would correspond to the projection point. This minimum is an internal extremum of a smooth function, and so it fulfills the usual zero-derivative criterion d(d 2 ) dz = 0 Therefore we need the real root of the polynomial which finds the point P proj = z projẑ + projˆ . We need also a representative point on the FS. We choose to find it as the intersection of the FS with a line connecting P g to the projection point P proj (located on the RS). The coordinates of the intersection point P fproj can be called z fproj and fproj ≡ p n FS (z fproj ). We find this intersection point analytically by first computing the slope m of the line m = ( proj − g )/(z proj − z g ) , (C20) Figure 21. Top: Density ratios of the wind and ambient material, used to construct the momentum-conserving shell solutions for n = 1, 2, 4, and 6 toroids. Bottom: Properties of the 2D hydrodynamic bubbles constructed from (a) normal model and (b) radial model by Equation (C15) in Step 2, and (c) the v θ = 0 numerical model. The upper half of the panels shows the radial shell positions as a function of θ, with the shades indicating the regions between the reverse shocks and forward shocks (whose positions have been stretched by 2% for better visualization). The r-axis is nonlinearly stretched as a square root. The lower half shows the respective shell thickness measured radially between the RS and FS, the RS and the CD, and the CD and the FS, with a logarithmic r-axis. and then finding the place where that line cuts the FS, using the polynomial fitting curve. That gives a polynomial equation whose relevant real root gives the value of z fproj . Once the two representative points P proj and P fproj have been identified, we find their pre-shock and then their post-shock velocities. The pre-shock (gas) velocities in the fixed frame (ff) are, at P proj and P fproj , v ff pre (RS) = v 0r and v ff pre (FS) = 0 .
The usual formulas to obtain post-shock velocities are given in frames comoving with the shocks and a simple transformation provides that, as follows. The shock (pattern) velocities have been estimated above (Equation (C16)) as The numerical values of these vectors are stored in arrays, which are then interpolated in z(θ). Velocity subtractions such as provide the transformations of the pre-shock gas velocities to the frames of the reverse and forward shocks (rsf and fsf). In order to obtain the post-shock gas velocities, it is useful to separate the pre-and post-shock velocity components, respectively, tangential and normal to the shocks. That requires knowing the respective unit vectors tangential and normal to the shock. Simple differential geometry using the fitting polynomials gives the tangent and normal unit vectors at P proj and P fproj , which we may call here (T RS ,n RS ) and (T FS ,n FS ).
These unit vectors are now used to obtain the tangential and normal pre-shock velocities in the shock frame. The tangential components are unchanged before and after the shocks in the shock frames, whereas the normal component is altered based on shock jump conditions. The effective (normal) shock Mach numbers are now known, and applying formulas for oblique shocks such as Equations (94) and (133)  which give the dependence of the ratios of normal components on the effective Mach numbers. Combining now the tangential and normal post-shock gas velocities in the shock frames gives the total velocity vectors in the shock frames, which can be simply transformed to the fixed frame of reference by adding back the respective shock velocities (v RS and v FS ), yielding as a result the post-shock gas velocities in the fixed frame at the two representative points. C.5.
Step 4: Construction of the Shear Profile The final step is to construct and identify the local velocity field for any point P g located between the RS and FS. The two points P proj and P fproj with fixed-frame post-shock velocities as illustrated in Figure 22 are determined in Step 3. Those post-shock velocities are known as v ff post (RS) and v ff post (FS). The distances d RS and d FS from P g to the projection points P proj and P fproj are also known. Simple interpolation gives the constructed estimate for the shear velocity at P g : thus completing the identification of velocity profile for any location in the post-shocked region. The resulting shear velocity field is illustrated in the panels of Figure 23. In this step, we build upon the understanding of parallel oblique shock structures presented in Appendix C of Paper I for the constructions of the analytic procedures here. There we had derived an analytical expression for the total shear across the post-shock region (without information of the structures inside), bounded by the RS and FS. Here we utilize the analytic formulations of the jump conditions and behavior of the oblique shocks, and construct analytically a whole new shear profile. We deploy methods of interpolation between two newly computed estimates of the locations of the RS and FS. Our new procedures are novel in giving analytical constructions for the shear velocity profile and for the thickness of the shells that possess 1D radial RS, CD, and FS structures, while extended to 2D in θ. Because of the shift of focus, a new set of notations different from those in Appendix C of Paper I for the basic schematics of the parallel oblique shocks has been adopted in this work. Figure 23. Shear velocity panels, illustrating the results of Appendix C. The panels are in deeply zoomed-in scales, able to show details of the velocity structure inside the shell as constructed with the models presented in Appendix C. Columns, left to right: n = 1, 2, 4, and 6. Rows, top to bottom: top, middle, and base of the outflow. RS (red), CD (green), and FS (blue) drawn as solid lines. Velocity vector scales of each panel given by the free wind velocity v0 in the left portions before the RS.

D. SMALL-ANGLE OBLIQUENESS FOR THE MOMENTUM-CONSERVING THIN-SHELL CURVES
Under the momentum-conserving model, the shape of the outflow shell is controlled by two functions of θ: the wind velocity v w (θ) and the ratio δ(θ) of wind to ambient densities. For this work, the wind velocity is adopted to be v w (θ) = v 0 . The density ratio results from dividing a wind function ρ w =ρ w (θ)/r 2 over an ambient function ρ a =ρ a (θ)/r 2 . We consider the wind function to have the form ρ w = D 0 (sin θ) −k /r 2 , k ≥ 0. The ambient function of this work is the sum of two terms (Equation (7)): a main term ρ T (r, θ) = (a 2 a /2πGr 2 )R(θ) (Equation (6)) representing the toroid solutions of Li & Shu (1996b), parameterized by n ≥ 0, and a minor term ρ S (r) = D S /r 2 with D S ≥ 0.
The momentum-conserving shape function corresponding to this δ is r(θ) = v 0 t(1 + δ −1/2 ) −1 . Obliqueness can be obtained from this shape function and the equation sin χ = −r / r 2 + r 2 = −(r /r)/ 1 + (r /r) 2 , where the ratio r /r ≡ (1/r)(dr/dθ) needs to be obtained for a given θ value. The equations up to this point in the appendix are valid for all θ. However, henceforth, in this appendix, we will focus on the small-angle limit θ → 0. For small θ, the toroid solutions behave as ρ T ∝ (sin θ) n . The D S > 0 cases are referred to as the tapered toroids, and the D S = 0 cases are the untapered toroids. In the limit of small θ angles, the density ratio is δ −1 = ρ a /ρ w = a 2 a 2πG a 0 (n) D 0 (sin θ) n+k + D S D 0 (sin θ) k , where the ambient sound speed a a sets the scale for the toroids, and a 0 (n) is a constant arising from their solutions. Unless D S is set to zero, the dominant term (for small θ) in the above expression has a power law of the form (sin θ) k . Using this value of δ, the ratio r /r for small values of θ is obtained as r /r = − δ 1/2 2 (1 + δ −1/2 ) −1 cos θ (n + k) a 2 a 2πG a 0 (n) D 0 (sin θ) n+k−1 + k D S D 0 (sin θ) k−1 , which can be used in combination with Equation (D28) to obtain values of χ(θ) valid for small θ. We need now to establish the leading form of these equations as θ → 0, starting from Equation (D29). For the untapered toroids with D S = 0, only the first term remains, and therefore it leads. For n = 0, even when tapered, both terms are of the same order in sin θ, and our choices of parameters are such that if n = 0 we have that a 0 = 1 and (a 2 a )/(2πG) D S , making the tapered and untapered toroids effectively coincide in that the leading term is again the first one. In the remaining scenario, for tapered toroids with n > 0, for a sufficiently small value of θ it is the second term that dominates: independent of n for a sufficiently small θ. Therefore, the leading term always has the form where for D S = 0 we have p = n + k, and for D S > 0 we have p = k, and the coefficient c is equal to either a 2 a 2πG (a 0 /D 0 ) (for the untapered n > 0 toroids, and for n = 0) or (D S /D 0 ) (for the tapered toroids with n > 0).
We now apply this result to our cases. For the n > 0 untapered toroids with our usual winds, we have that p = n + k with k = 2, and therefore the obliqueness is given by sin χ = [1 + (4/c)(n + 2) −2 (sin θ) −n ] −1/2 , an expression in which for n > 0 the sine terms dominates the RHS of Equation (D34), leading to sin χ ∝ θ n/2 , which has a limiting value of zero at θ = 0 for the untapered toroids with n > 0. For the tapered toroids, and for the n = 0 case, both terms are of the same order, leading to a nonzero limiting value of sin χ for θ → 0, equal to sin χ = [1 + (1/c)] −1/2 , depending on c, which is a function of problem parameters. This equation is illustrated in detail in Figure 24. Figure 25 shows the outflow shapes at all angles, with its top row focusing on very small angles to help in examining the θ → 0 limit of the inclination angle χ. A variety of D S values are shown, including tapered cases (D S > 0 such as in our usual simulations) and untapered cases (D S = 0). Naturally due to their less-dense ambient medium, outflows of untapered cases can proceed faster, and they are the larger outflow curves in each panel of Figure 25. For any untapered n > 0 case, it has been shown (from Equation (D34)) that lim θ→0 χ = 0, and therefore the tangent to the outflow shape at that point θ = 0 • is horizontal, parallel to the disk and perpendicular to the jet axis. The approach of the curves to their horizontal tangents is seen in the n > 0 top panels of Figure 25, most clearly for n > 1, but nonetheless true for any n > 0. However, approaching the horizontal tangent limit requires smaller values of θ → 0 for n = 1 than for, e.g., n = 2, and it would require even smaller θ values if we want to consider n = 0.5. For the tapered cases D S > 0, and also for n = 0, the tangent is oblique even at θ → 0, following Equation (D35). Consistent with that, all of the curves with D S > 0 and in the n = 0 column of Figure 25 have nonhorizontal tangents at all θ, including θ = 0 • , for which the inclination is shown in Figure 24.
We have adopted a wind with a density profile (Equation 1) ρ ∝ (sin θ) −k with k = 2, with a formal divergence at θ → 0 • , compatible with the simulations because of the implicit tapering provided by angular resolution, which keepṡ M finite and well controlled (Shang et al. 2006;Wang et al. 2015). In addition to that, here we briefly discuss the k = 0 Figure 24. The blue curve shows the small-angle (θ → 0) limit of obliqueness χ (ranging between 0 • to 90 • , horizontal black dotted lines) as a function of the density factor D S . This limit is computed using Equation (D35) with the tapered n > 0 cases c = D S /D0 and D0 = 2.0025 × 10 11 g cm −2 . The n > 0 tapered HD obliqueness approaches ∼ 47 • (horizontal green dashed line) for D S = 2.25 × 10 11 g cm −2 (vertical orange dashed line). The horizontal red dashed line is the θ → 0 limiting value of the untapered n = 0 case (∼ 88 • ) with c = (a 2 a /2πG)/D0. The tapered χ can be seen to approach the value of the untapered n = 0 case as D S ≈ 10 15 g cm −2 .
This expression has discontinuities at n = 0 and n = 2. The k = 0, n = 2 case has similar algebraic behavior to the k = 2, n = 0 case, consistent with the definition p = n + k.

E. MAGNETIC FORCES AND PSEUDOPULSES
The panels of Figure 26 are shown here to illustrate the effects and connections between the velocity and force components.

F. SHAPES OF PERPENDICULAR PVDCD DIAGRAMS
In this appendix, the diagram shapes sketched in Figure 19 are listed and then compared to the PVDCD shapes of Figure 18 computed using Synline. This continues (in more detail) the comparison of perpendicular PVDCDs, which was summarized with a few examples in Section 8.2.
• Intersecting ellipses or ovals • • as in Figure 19(c) or (e) are present in many cases, including Figure 18  . Shapes of the HD MC curves for n = 0, 1, 2, 4, and 6 toroids, for variation of the density factor D S of 1000, 100, 10, 1, 0.1, 0.01, 0.001, and 0 relative to the current setup of D S = 2.25 × 10 11 g cm −2 . The tips change from pointy to flat as D S decreases, except for n = 0. The lower row shows the whole outflow shape, and the upper row shows the zoomed-in tips.
• In this set of simulations, the concentric ellipses • • as in Figure 19(d) and (h) are present in many of the densest features (see below). The overall shape of Figure 18(e) (n = 1-2,M A = 6) is also present in Figure 19(h), probably due to the CDIE cutoff level isolating the densest features in this case.
In Figure 18(e) (n = 1,M A = 30), the overall shape is as that of Figure 19(e) due to a small high-velocity portion; otherwise, it would be as that of Figure 19(g).
Here we mention a few of the different shapes of the high CDIE features, which are the possible result in realistic observations, depending on excitation conditions and instrument sensitivity. Figure 18(c) (n = 1-2,M A = 6) contains high CDIE concentric ellipses as in Figure 19(h): the thicknesses ∆y and ∆v x of these sets of ellipses are connected to the slit thickness ∆z. Dense CDIE features shaped as in Figure 19(h) are also present in Figure 18