Issue 
A&A
Volume 668, December 2022



Article Number  A78  
Number of page(s)  29  
Section  Interstellar and circumstellar matter  
DOI  https://doi.org/10.1051/00046361/202141316  
Published online  09 December 2022 
Modeling the CO outflow in DG Tauri B: Sweptup shells versus perturbed MHD disk wind
^{1}
Univ. Grenoble Alpes, CNRS, IPAG,
38000
Grenoble, France
email: catherine.dougados@univgrenoblealpes.fr
^{2}
Observatoire de Paris, PSL University, Sorbonne Université, CNRS UMR 8112, LERMA,
61 avenue de l’Observatoire,
75014
Paris, France
^{3}
Departamento de Astronomía, Universidad de Chile,
Casilla 36D,
Santiago, Chile
^{4}
Instituto de Radioastronomía y Astrofísica, Universidad Nacional Autónoma de México,
Antigua Carretera a Pátzcuaro 8701, ExHda. San José de la Huerta,
C.P. 58089,
Morelia, Michoacán, México
Received:
14
May
2021
Accepted:
20
September
2022
Context. The origin of outflows and their exact impact on disk evolution and planet formation remain crucial open questions. DG Tau B is a Class I protostar associated with a rotating conical CO outflow and a structured disk. Hence it is an ideal target to study these questions.
Aims. We aim to characterize the morphology and kinematics of the DG Tau B outflow in order to elucidate its origin and potential impact on the disk.
Methods. Our analysis is based on Atacama Large Millimeter Array (ALMA) ^{12}CO(2–1) observations of DG Tau B at 0.15″ (20 au) angular resolution. We developed a tomographic method to recover 2D (R,Z) maps of vertical velocity V_{Z} and specific angular momentum j = R × V_{ϕ}. We created synthetic data cubes for parametric models of winddriven shells and disk winds, which we fit to the observed channel maps.
Results. Tomographic analysis of the bright inner conical outflow shows that both V_{Z} and j remain roughly constant along conical surfaces, defining a shearlike structure. We characterize three different types of substructures in this outflow (arches, fingers, and cusps) with apparent acceleration. Winddriven shell models with a Hubble law fail to explain these substructures. In contrast, both the morphology and kinematics of the conical flow can be explained by a steady conical magnetohydrodynamic (MHD) disk wind with footpoint radii r_{0} ≃ 0.7–3.4 au, a small magnetic level arm parameter (λ ≤ 1.6), and quasi periodic brightness enhancements. These might be caused by the impact of jet bow shocks, source orbital motion caused by a 25 M_{J} companion at 50 au, or disk density perturbations accreting through the wind launching region. The large CO wind mass flux (four times the accretion rate onto the central star) can also be explained if the MHD disk wind removes most of the angular momentum required for steady disk accretion.
Conclusions. Our results provide the strongest evidence so far for the presence of massive MHD disk winds in Class I sources with residual infall, and they suggest that the initial stages of planet formation take place in a highly dynamic environment.
Key words: stars: individual: DG Tau B / stars: formation / ISM: jets and outflows / stars: protostars
© A. de Valon et al. 2022
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the SubscribetoOpen model. Subscribe to A&A to support open access publication.
1 Introduction
Understanding the origin of protostellar flows is a key element to our full comprehension of the star formation process. Protostellar flows come in two components: high speed collimated jets and slower, often less collimated winds and outflows. We focus here on the slow molecular outflows which are traditionally associated with the earlier stages of star formation. However, they have also recently been detected around more evolved Class II systems (e.g., Pety et al. 2006; Louvet et al. 2018; FernándezLópez et al. 2020). Despite their ubiquity, the exact origin of molecular outflows, their link to the highvelocity jets, and their impact on the young forming star and disk are still crucial open questions.
Two main paradigms are currently considered. The first traditional model describes these slow outflows as sweptup material, tracing the interaction between an inner jet or a wideangle wind with the infalling envelope or parent core. These models have been mainly used for interpreting outflows from Class 0 and I stars, which are still surrounded by massive envelopes (Zhang et al. 2019; Shang et al. 2020; Lee et al. 2000). However, on small scales (less than a few 1000 au), recent observations have revealed rotating molecular outflows to originate from well within the disk at all evolutionary stages from Class 0 to Class II (e.g., Launhardt et al. 2009; Zapata et al. 2015; Bjerkeli et al. 2016; Tabone et al. 2017; Hirota et al. 2017; Louvet et al. 2018; Zhang et al. 2018; Lee et al. 2018, Lee et al. 2021; de Valon et al. 2020). These observations suggest an alternative paradigm by which these slow molecular outflows, at least at their base, would trace matter directly ejected from the disk, by thermal or magnetic processes. In support of this interpretation, the flow rotation signatures are consistent with an origin from disk radii r_{0} ≃ 1–50 au (see references above), where Panoglou et al. (2012) have shown that magnetic disk winds could remain molecular.
These two paradigms imply different evolutions for the disk. Jet and winddriven shell models predict that an important mass is swept up from the envelope, impacting the reservoir of matter infalling onto the disk. Disk wind models predict an extraction of mass from the disk and, in the case of magnetohydrodynamic (MHD) models, an extraction of angular momentum which could drive disk accretion (Bai et al. 2016). Evaluating the contributions of each of these mechanisms to the slow molecular outflow emission requires highangular resolution studies of the molecular outflow base. The flux of angular momentum extracted by the rotating molecular wind is estimated in only two sources so far, HH30 and HH212, and it is found sufficient to drive disk accretion at the current observed rate in both cases (Louvet et al. 2018; Tabone et al. 2020).
Recent highangular resolution observations also reveal striking signatures of multiple CO shells in a few sources (Zhang et al. 2019; FernándezLópez et al. 2020). Under the classical paradigm where CO outflows trace sweptup shells, they would require short episodic wind and jet outbursts every few 100 yr. Characterizing and understanding the origin of these variabilities could bring critical insights into the star formation dynamics.
We present here an analysis of the DG Tau B CO outflow, based on recent Atacama Large Millimeter Array (ALMA) observations at 0.15″ resolution by de Valon et al. (2020) (hereafter DV20). DG Tau B is a Class I 1.1 M_{⊙} protostar located in the Taurus cloud (≈ 140 pc) and is associated with a bipolar atomic jet (Mundt & Fried 1983) and a strongly asymmetric CO outflow first mapped by Mitchell et al. (1997). The bright redshifted CO outflow lobe displays a striking bright and narrow conical shape at its base. Zapata et al. (2015) detect rotation signatures in the same sense as the disk. The ALMA observations by DV20 clearly confirm rotation in the bright inner conical redshifted lobe and show that it is surrounded by a wider and slower outflow. Residual infall signatures are detected at opening angles ≥70°, almost tangent to the disk surface. In addition, DV20 report striking substructures in the CO channel maps at different lineofsight velocities, reminiscent of the nested layers recently identified in HH46/47 by Zhang et al. (2019) and suggesting variability or interaction processes. The exquisite levels of detail provided by these new ALMA observations provide a prime opportunity to distinguish between sweptup and disk wind origins.
In Sect. 2, we recall the main properties of the DG Tau B outflow and characterize the three types of substructures visible in the channel maps. In Sect. 3, we present a modelindependent analysis of the inner conical outflow component, which allowed us to retrieve 2D maps of the expansion velocity V_{Z} and specific angular momentum j = RV_{ϕ}. We compare these overall properties with parametric models of winddriven shells (in Sect. 4) and disk winds (in Sect. 5). We discuss our results and their implications for the origin of the CO outflows in DG Tau B in Sect. 6. Section 7 summarizes our conclusions.
Fig. 1 Substructures in the DG Tau B redshifted CO outflow. Left panels: ^{12}CO channel maps at selected lineofsight velocities. The white dashed line traces the θ = 17° outer limiting cone of the inner outflow as defined by DV20. Right panels: transverse PV diagrams across the flow averaged over a slice of ∆Z = 0.2″. Black dotted lines indicate the outer limits of the conical outflow at the specified height. The triangle (a) and the two circle symbols (c and d) highlight the onaxis height of one arch and two different cusps, respectively. The square symbol (b) is located in the extended outer flow. The symbols are represented both on the channel maps and PV diagrams. 
2 Summary of outflow structure
Figure 1 summarizes the main properties of the DG Tau B redshifted CO outflow as identified in DV20. A narrow, limb brightened conical outflow is visible in the channel maps at (V − V_{sys}) > 2.0 km s^{−1}. Its opening angle decreases from 17° at (V − V_{sys}) = 2.8 km s^{−1} until 12° at (V − V_{sys}) ≥ 5.0 km s^{−1}. The sheerlike velocity gradient across the conical layer is best seen in transverse positionvelocity (hereafter PV) cuts (Fig. 1), where the flow width clearly narrows down at higher velocity, up (V − V_{sys}) > 5.0 km s^{−1}. The conical outflow is surrounded by a slower and wider outflowing component visible at (V − V_{sys}) < 2 km s^{−1}. This outer flow is visible in PV diagrams as an extended pedestal with a shallower velocity gradient (see Fig. 1).
Another striking property of the DG Tau B redshifted outflow, revealed by the ALMA observations in DV20, are brightness enhancements visible in channel maps. These various substructures, illustrated in Fig. 1, can be classified in three types:
Bowshaped intensity enhancements are visible at low velocities (V − V_{sys}) = 0.88–1.51 km s^{−1} (see Fig. 1, left panels). We refer to these substructures as arches. The radial extent of the biggest arch is larger than the inner conical outflow, implying that this arch is at least partially formed outside of the conical outflow. The arch seems to increase in height with increasing velocity although this phenomenon is difficult to quantify due to the limited spectral sampling.
At intermediate velocities, from (V − V_{sys}) = 2.46–3.42 km s^{−1}, thin quasi vertical lines (see white arrow in Fig. 1) are visible inside the conical outflow, close to the edge. They are almost vertical at (V − V_{sys}) = 2.46 km s^{−1} and more open at higher projected velocity until they become almost tangent to the edge of the conical outflow at (V − V_{sys}) = 3.42 km s^{−1}. We refer to these substructures as fingers.
At high velocities, from (V − V_{sys}) = 3.1 to almost 7 km s^{−1}, multiple Ushaped structures are visible inside the conical outflow. We refer to these substructures as cusps. The contrast of these cusps is maximal at (V − V_{sys}) = 4.37 km s^{−1} and decreases with increasing velocity. The cusps show a signature of apparent acceleration: their projected distance from the source increases with increasing projected velocity.
Table 1 lists the characteristics of the main arches and cusps. We identify four arches (A0 to A3) and six cusps (U0 to U5). On the channel map at (V − V_{sys}) = 1.19 km s^{−1} we derived the maximal height of each arch on axis (at δx = 0) and the maximal radial extension. We divided these two values to derive the arch aspect ratio. The cusps are also characterized from the channel maps at 3.73 ≤ (V − V_{sys}) ≤ 5.32 km s^{−1} (see Fig. B.1). At higher velocities, the cusps could not be characterized because the outflow signaltonoise ratio (S/N) decreases drastically. Moreover, the region at δz < 2.2″ was not studied because the cusps locations are complex to identify due to overlapping structures. We derived the cusp reference height onaxis on the channel map at (V − V_{sys}) = 4.37 km s^{−1}. The apparent acceleration of each cusp, in (″)/km s^{−1} listed in Table 1, was obtained by measuring the average spatial shift of the cusp between two consecutive channel maps (taking as error bar the rms dispersion between measurements in different channels).
Internal discrete structures are also visible in transverse PVdiagrams as pseudoellipses (see Fig. 1). The top and bottom of the ellipses seem to match with respectively the top of some arches and bottom of some cusps (see Fig. 1). This potentially indicates that cusps and arches are linked to the same phenomenon. We present a modeldependent study of these ellipses in Sect. 5.
Characteristics of the observed arches and cusps.
Fig. 2 Principle of tomographic reconstruction method. Top panel: 3D representation of the two reference systems used: the outflow (X, Y, Z) in black and the observer (δx, δy, δz) in yellow. The colored circles illustrate the emitting rings defined by five parameters (see text). The colored dots trace the locations at ϕ = 0 and ϕ = π along the emitting rings. Bottom panels: transverse PV diagrams at δz = 2.5″ across the flow axis and averaged over a slice of ∆z = 0.2″. In the left panel is shown the schematic projection of the colored rings and corresponding colored dots in the PV diagram. Because of the flow inclination, the rings are in fact projected at slightly different heights. Their real projection is studied in Appendix D.1. The white dots in the right panel illustrate the outer limits of the PV diagram. The red curve shows the polynomial fitting of the two edges. 
3 Tomography of the inner conical outflow
In this section, we develop a modelindependent method that allows us to recover the dynamics and the morphology of the inner conical outflow component. This method assumes that the outflow is axisymmetric. We later discuss possible departures from axisymmetry and their implications on the analysis conducted here.
3.1 Method
We followed Louvet et al. (2018) who modeled the outflow of HH30 at a given vertical offset by an emitting ring with radius R and extended their method to take into account the inclination of the outflow. For this purpose, we defined the outflow and the observer reference systems (see Fig. 2). On the outflow reference system, Z is defined by the outflow axis and X is tangent to the plane of sky. The observer reference system is defined by δz the projection of the outflow axis onto the plane of the sky, δy, the lineofsight direction, and δx = X, in the plane of the sky. The inclination of the outflow i is then defined by the angle between δy and Z. In the case of edgeon disks such as HH30, the two reference systems are identical. We modeled one layer of the outflow at a specified height Z by an emitting ring of radius R and azimuthal angle ϕ (see Fig. 2) with X = R cos ϕ and Y = R sin ϕ. For each ring the velocity components are defined in cylindrical coordinates with: V_{R}(R, Z), V_{Z}(R, Z) and V_{ϕ}(R, Z) (see Fig. 2). Hence, an emitting ring is defined by 5 parameters: Z, R, V_{r}(R, Z), V_{z}(R, Z) and V_{ϕ}(R, Z).
The observational coordinates on a positionpositionvelocity (PPV) data cube are defined by the projection of the outflow on the plane of sky (δx, δz) and the projected velocities on the line of sight V_{los} = −V · e_{y} with redshifted velocities considered as positive. This depends on R, Z and ϕ as: (1) (2) (3)
A transverse PV diagram corresponds to a pseudoslit of the data cube perpendicular to the flow axis. This corresponds to a solution of Eqs. (1)(3) with δz = cst. In the case of edgeon flows, a ring traces a perfect ellipse in the PV diagram. A fit of these ellipses give complete information about the morphology and dynamics of the outflow as shown by Louvet et al. (2018).
In the case of an inclined outflow such as DG Tau B, different rings overlap on the transversal PV diagrams (see Fig. 2 left). Hence it was not possible to fit them individually. However constraints on some of the ring parameters could be derived from characterizing the outer limits of the PV diagram. The radius of the ring corresponds on the first order to δx_{max} = δx(ϕ ≈ 0, π). In addition, the projected velocities at the edge of the ellipses V_{los}(δx_{max}) allow one to recover both V_{Z}(R, Z) and V_{ϕ}(R, Z) from the following equations: (4) (5) (6) (7)
By consequence, characterizing the outer limits of the PV diagram along multiple heights allowed us to recover a 2D map of the expansion velocity V_{Z} and specific angular momentum j = RV_{ϕ}. In the following, we used the inclination derived from DV20 at i = 117° ± 2. To characterize the outer shape of the transverse PV diagram at a given δz, we derived the maximal projected velocity for each value of ±δx. Numerically, we computed the gradient of the emission profile at a fixed δx and localized its maximum. We also derived an uncertainty on V_{los} which is found to vary in the range 0.05 to 0.2 km s^{−1}. We used in the following a mean value of 0.1 km s^{−1}. Figure 2 illustrates our method to determine V_{los} on the edges of the PV diagram. This procedure failed at low radii (or high velocities) because of the low S/N and the almost vertical profile that generated high uncertainties on the velocity estimate. The determination of the velocity was also limited by our spectral resolution of 0.3 km s^{−1}. We fit the two edges for each PV diagram with a polynomial curve (see Fig. 2). We also applied a Gaussian filter with δz using a standard deviation of 0.16″.
Fig. 3 Tomographic maps of V_{Z} (left panel) and specific angular momentum j (right panel) in the outflow referential. The black dashed line traces the conical fit of the region V_{Z} = 10–11 km s^{−1}. The white (resp. black) contours in the left (resp. right) panel show V_{Z} contours. The red dashed lines indicate the heights of the two extrema in specific angular momentum. 
3.2 Results
Using Eqs. (5), (4), (6), and (7), a tomographic map of V_{Z}(R, Z) and V_{ϕ}(R, Z) could be recovered. We show the specific angular momentum RV_{ϕ}(R, Z) instead of the rotation alone as this is more meaningful in the understanding of the dynamics. Figure 3 shows the resulting tomographic map of V_{Z} and RV_{ϕ} in the outflow referential. The tomography efficiently traces the conical shape visible on the channel maps. Curves of constant V_{Z} trace conical surfaces with semiopening angles varying from 12° for the highest velocities to 17° for the lowest velocities. V_{Z} radially decreases from ≈ 14 km s^{−1} to 5 km s^{−1}. This range of velocities is conserved until at least Z = 1200 au.
The specific angular momentum derived from the tomographic study varies from 0 to 140 au km s^{−1}and is consistently in the same sense as the disk rotation. At Z < 500 au the specific angular momentum increases with radius from ≈30 au km s^{−1} in the inner radius to ≈ 100 au km s^{−1} on the outer radius. The specific angular momentum is also roughly constant on conical lines of constant V_{Z} until Z ≈ 500 au (see Fig. 4). Our average value around 70 au km s^{−1} is consistent with the previous estimate of DV20 of ≈65 au km s^{−1}.
Two extrema in the specific angular momentum map can be observed at Z = 550–800 au and 900–1100 (see Fig. 4). At the lowest altitude, the specific angular momentum reaches zero while at the highest altitude the specific angular momentum increases up to j > 140 au km s^{−1}. These irregularities are also visible on channel maps. They correspond to regions where bumps in the cones are observed: toward δx < 0 at ≈4.5″, δx > 0 at ≈6.5″ (see the channel map (V − V_{sys}) = 4.37 km s^{−1} on Fig. 1). These bumps may be due to local radial displacements of the outflow axis, due for example to wiggling. We study the impact of small amplitude wiggling in Sect. 5.2.
Fig. 4 Specific angular momentum j along curves of constant V_{Z}. The corresponding range of VZ is shown in the box. The uncertainty of specific angular momentum was obtained by propagating the V_{proj} uncertainty. The red dashed line corresponds to the median value of each curve. 
3.3 Limitations and biases
In this section, we discuss the different biases and limitations of this tomographic study. Firstly, this study could not recover the radial velocity component V_{R} as it impacts mostly the size of the ellipse at x = 0 where all the ellipses are stacked. A modeldependent study to characterize this radial velocity will be achieved in Sect. 5. Furthermore, in order to apply this method, it is critical that the centers of the rings are not significantly displaced from δx = 0. Such displacements can be induced by a poorly estimated outflow position axis (PA) or by outflow axis wiggling. We derived the PA of the redshifted outflow in Appendix A at PA = 295° ± 1°. This value is in very good agreement with the disk rotation axis PA = 115.7° ± 0.3° determined by DV20. We also determined an upper limit of 0.5° for the wiggling of the CO outflow axis. We discuss in Sect. 5.2 the impact of possible lowamplitude wiggling on our results.
The different biases were also computed. We show in Appendix D.1 that assuming that the maximal radial extent corresponded to ϕ = 0, π was partially inaccurate, and could introduce a bias in the estimate of R, Z, V_{Z}(R, Z) and V_{ϕ}(R, Z). The effect of ellipse stacking and its effect on the estimate of the velocities was studied in Appendix D.2. We evaluate at <20% the potential bias in our estimate of the conical outflow dynamics. Our estimates of V_{Z} and specific angular momentum are overestimated and underestimated respectively (see Appendix D.2). We estimate the bias on R and Z to be respectively <1.5% and <3%, resulting in an error <3% on our estimate of the opening angle θ.
The highly asymmetric pedestal emission visible on the two sides of the PV diagram at (V − V_{sys}) < 2 km s^{−1} traces the outer region (see Fig. 1). DV20 show that this region is outflowing and surrounds the conical outflow. We did not attempt to apply the tomographic method to this region. From its morphology in the channel maps, we derived for this component an opening angle >30^{°}. Such a large opening angle produces a bias of ≈60% in the estimation of V_{Z} using our reconstruction method (see Appendix D). Moreover, this pedestal may potentially be explained by the top of one large ellipse, with the extremal region located at velocities < 1 km s^{−1}, absorbed by the medium. In that case, our reconstruction method is not applicable. By consequence, we did not apply our method for lineofsight velocities (V − V_{sys}) < 2 km s^{−1}. In the following, we investigate to which extent winddriven shells and disk winds can account for both the conical velocity stratification determined here and the striking substructures (arches, fingers, cusps) identified in Sect. 2.
4 Winddriven shell modeling
The traditional interpretation proposed for CO molecular outflow cavities around young stars is that they trace shells of ambient material swept up by a wideangle wind or by jet bow shocks (see Cabrit et al. 1997; Lee et al. 2000; Arce et al. 2007, for reviews). In this section, we investigate the simplest and most widely used model to interpret CO outflow observations, namely the winddriven shell (hereafter WDS) solution of Lee et al. (2000) where the shell is a parabola that expands radially in all directions with a velocity proportional to the local distance from the source (hereafter referred to as the Hubble law). Such a shell structure is predicted under a set of specific conditions in the wind and ambient medium^{1} (see Sect. 6 for details). This simple WDS model is recently shown by Zhang et al. (2019) to reproduce several features of the multiple CO shell structures at the base of the HH46/47 molecular outflow. Therefore it is natural to investigate whether the same WDS model can also reproduce the morphology and kinematics of the DG Tau B outflow, on smaller spatial scales.
Following Lee et al. (2000), the parabolic morphology and the radial Hubblelaw kinematics of the shell can be empirically described by two parameters, C and τ, through: (8)
where τ is the age of the shell, C is the inverse size of the parabola at θ = 45° (where Z = R = 1/C), and the product τC defines the shell expansion speed at each polar angle θ = arctan(R/Z) through: (9) (10)
The above equations always produce ellipses in both channel maps and transverse PV diagrams (Lee et al. 2000). This is a direct result of the assumed Hubble law, where the shell velocity vector is proportional to the position vector. A channel map at a given lineofsight velocity is then equivalent to making a cut through the shell at a given depth along the line of sight and this cut is shaped as an ellipse. Similarly, a transverse PV diagram has the same (elliptical) shape as a cut through the shell at the corresponding projected height. On channel maps, the ellipse is projected at increasing distances from the source with increasing velocity, due to the Hubble law. Similarly, on transverse PV diagrams, the mean velocity of the ellipse increases with the distance of the PV cut from the source (see Figs. 24 and 26 in Lee et al. 2000). In Appendix E, we derived analytical formulae for the center of the ellipse in channel maps as well as for its aspect ratio. Interestingly, we find that the ellipse aspect ratio only depends on the inclination i and is equal to 1/cos i for the classical model of Eq. (8) (see Appendix E). In the following, we attempt to fit with this model first the lowvelocity outer flow component, and then the bright conical outflow and its discrete structures.
Fig. 5 Comparison of the lowvelocity outer CO outflow with a classical WDS model. Left panel: ^{12}CO individual channel maps at different lineofsight velocities tracing the lowvelocity outer CO outflow. The white contours trace the model of a WDS of parabolic shape defined by C = 10^{−3} au^{−1}, dynamical age τ = 6000 yr and specific angular momentum j = 250 ± 50 au km s^{−1} with an inclination of i = 117°. Right panel: transverse PV cut at δz = 7.5″ averaged over a slice of ∆Z = 0.2″. The red ellipse traces the WDS model. (V − V_{sys}) units are km s^{−1}. 
4.1 Lowvelocity outer flow
The wide and lowvelocity outflow at (V − V_{sys}) < 2 km s^{−1} shows several properties suggestive of a “classical” parabolic WDS with a Hubble law. Its outer border in channel maps has a parabolic shape, and it exhibits a larger offset from the origin at higher lineofsight velocities (see Fig. 1 left panels). Although such WDS models do not usually consider rotation, we included rotation to properly fit the large leftright asymmetry observed in the channel maps. To reduce the number of free parameters, we considered that the specific angular momentum j is the same at all positions of the sweptup shell. C was then fixed by the global parabolic shape of the cavity, τ by its spatial shift between the channel maps, and j by its global leftright asymmetry.
Figure 5 (left panels) shows that the outer contour of the lowvelocity outflow and its increased spatial offset with velocity are well fit by a WDS obeying Eq. (8) with parameters C = 10^{−3} au^{−1}, τ = 6000 yr and j = 250 ± 50 au km s^{−1}. The rightmost panel in Fig. 5 shows that this shell model reproduces well the most extended, lowest velocity emission of the broad pedestal in transverse PV cuts; the predicted blueshifted emission from the front side of the shell falls very close to systemic velocity, consistent with its nondetection in our data. On the other hand, our assumption of a thin parabolic shell does not match the observed outflow thickness at high altitudes. This discrepancy is visible at δz ≈ 15″ on Fig. 5 where the observed width of the emissive outer layer is ≈3″, significantly larger than predicted by our model. The shell should actually have a thickness ≃3″ ≃ 500 au.
4.2 Conical outflow and discrete structures
In this section, we attempt to model the conical outflow and its discrete structures (arches, cusps, fingers, described in Sect. 2) by a stacking of several parabolic winddriven shells with a Hubble velocity law. Several qualitative features are suggestive of such a model: the loop shapes of the arches at low velocity are reminiscent of the ellipses predicted in channel maps (see Appendix E), the apparent acceleration of the cusps (increased altitude with increasing velocity) is reminiscent of the predicted Hubble law dynamics, and finally, in the conical flow studied by tomography, contours of constant τ = Z/V_{Z} follow quasiparabolic curves above Z ≃ 400 au (see Fig. 6). Hence we investigate below whether the conical outflow could be made of successive nested parabolic winddriven shells, where the apparent continuous aspect of the tomography would be an artifact of our limited spatial and spectral sampling, and the discrete structures (arches, cusps, fingers) would trace a few individual shells brighter than average.
Contrary to the slow outer flow modeled in Sect. 4.1, the leftright asymmetry in these faster flow regions is small. Therefore, we neglected rotation when fitting WDS models to the channel maps. We set the WDS axis inclination equal to the largescale disk inclination derived from ALMA studies (i = 63° ± 2°, DV20), leading to i = 180°−63° = 117° in the redshifted lobe.
Here, we find that the “classical” WDS model of Lee et al. (2000) encounters a major problem, as shown in the top row of Fig. 7: the predicted aspect ratio of ellipses in channel maps, A = 1/ cos i (see Appendix E) is too large (≃2.2). In order to match the observed aspect ratio of the arches (≃1.4), the inclination of the shell axis should be i ≈ 135° instead of i = 117°. In the WDS model, however, the direction of shell elongation is not arbitrary but must follow the direction of both highest wind density (traced by the axial jet) and lowest ambient density (traced by core flattening). Proper motions of jet knots in DG Tau B imply a jet axis inclination of i ≥ 65° for the blueshifted lobe (Eislöffel & Mundt 1998), hence i ≤ 115° for the redshifted lobe. This limit agrees within 2^{°} with the largescale disk inclination determined by ALMA (i = 63° ± 2°, DV20), which should follow the core flattening. Therefore, we can exclude a shell axis at i ≈ 135° as a solution to the ellipse aspectratio problem of the WDS model of Lee et al. (2000).
Since the ellipse aspect ratio in channel maps does not depend on τ nor C (see Appendix E) the only possibility to reduce it without changing the shell parabolic shape is to modify the shell dynamics. The maximum height of the ellipse is reached onaxis (δx = 0) where the projected velocity is greatly affected by the radial velocity component V_{R} (see Eq. (3) with ). To keep a small number of model parameters, we thus chose to add an ad hoc free parameter η that modified the radial velocity as: (11)
In Appendix E, we show that the ellipse aspect ratio in this modified model is set at A = (η tan^{2} i + 1) cos i. Since we wanted to reduce the aspect ratio, we needed η < 1. In other words, we needed a velocity vector that is more collimated (forwarddirected) than the radial shell expansion in the original WDS model of Lee et al. (2000). We refer to this ad hoc model as “modified collimated WDS”.
As shown in the second row of Fig. 7, the shape of the largest arch A0 at (V − V_{sys}) = 1.19 km s^{−1} is well fit by a modified collimated WDS with η = 0.6. The smallest two arches A2, A3 and the smallest cusp U5 are also well fit by a collimated WDS with η = 0.5 (Fig. 7, bottom row). The parameters of these bestfit solutions are listed in Table 2. The inferred shell dynamical times have typical intervals of ∆τ ≃ 300–750 yr, similar to those inferred by Zhang et al. (2019) in HH4647. The fit velocity values V_{Z} at θ = 14° (in the region of the conical flow) are also listed^{2}. Not surprisingly, their range of ≃7.5–16 km s^{−1} is similar to our “modelindependent” tomographic results for V_{Z} in the conical flow region.
The Hubble law in the WDS model predicts an everincreasing shell speed at higher altitudes, until it reaches the polar wind speed, which is ≃125 km s^{−1} according to the redshifted jet speed in DG Tau B (Eislöffel & Mundt 1998). In contrast, the maximum lineofsight velocity with detectable emission in our transverse PV cuts (averaged between the left and right sides) is found to stay roughly constant with altitude at V_{max} ≃ 8 ± 1 km s^{−1} (see Fig. 8b). Therefore, all successive winddriven shells should be truncated, or have their CO emission strongly suppressed, above the point where they reach V_{Z} = V_{max}/cosi ≃ 18 km s^{−1}. This velocity limit is close to the molecule dissociation limit ≃20 km s^{−1} in dense hydrodynamical shocks (see e.g., Wilgenbus et al. 2000). Therefore, the disappearance of CO emission above a certain speed might be explained by shockdissociation of ambient CO. Figure 8c shows that the corresponding truncation region for the bestfit WDS models in Table 2 has a rough conical shape with θ ≃ 9°. However, our ad hoc modified collimated WDS model meets two serious issues, detailed below.
The model predicts a full ellipse in each channel map (white contours in Fig. 7), which is not observed. In contrast, discrete structures highlight only a portion of ellipse, depending on the velocity range (see Fig. 7 and Sect. 2): the ellipse top at lowvelocities (arches), ellipse flanks at midvelocities (fingers), and ellipse bottom at highvelocities (cusps). We find that a transition from arches at low velocity to cusps at high velocity can only be obtained if emission is restricted to a range of heights from z_{min} to z_{max}, as illustrated by the colored contours in Fig. 7. Serious discrepancies still remain with observations, however: In the broadest shell, WDSA0, the extents of Arch A0 and Cusp U1 require inconsistent ranges of emitting heights, and the predicted “fingers” at intermediate velocity are much wider than observed (see blue and green contours in middle row of Fig. 7). In the smaller inner WDS, full ellipses are still predicted in intermediate velocity channels, which are not observed (see yellow contours in bottom row of Fig. 7). The same problems remain even if we adopt conical shapes for the shells instead of parabolae, hence the above discrepancies appear intrinsically linked to the assumed Hubblelaw dynamics.
Another serious issue is that the bestfitting value of η in our modified collimated WDS models is always close to 0.5 (see Table 2). A ratio V_{R}/V_{Z} = 0.5 R/Z corresponds to a velocity vector locally tangent to the parabola. Hence the shell is not expanding but stationary. The physical justification for the Hubble law in the WDS model, namely a shell expanding at constant speed over time (Shu et al. 1991), is then no longer applicable. If the emitting material is moving parallel to the shell, a velocity increasing in proportion to distance would require, instead, a constant accelerating force of unknown nature operating out to z = 3000 au, which is totally unphysical.
In summary, we find that only the outer faint, lowvelocity flow in DG Tau B can be reproduced with the parabolic WDS model with radial Hubble law proposed by Lee et al. (2000). In contrast, the bright conical outflow at mid to high velocity, although reminiscent of WDS models because of the apparent acceleration of its discrete structures, cannot be explained by such models, even when ad hoc modifications to the kinematics, emissivity range, and shape are introduced. The model faces important issues which seem intrinsically linked to the Hubblelaw dynamics. We discuss in Sect. 6.1 the implications of these results in the context of more general wind or jetdriven shell scenarios.
Fig. 6 Dynamical times (τ = Z/V_{Z}) derived from the tomographic map of V_{Z} are shown in color. The white lines correspond to a parabolic curve with 1/C varying from 30 to 135 au, in steps of 15 au. 
Fig. 7 Comparison of ^{12}CO channel maps at three different lineofsight velocities (color maps) with predicted ellipses for parabolic WDS (white contours). The model used on each row is sketched in the leftmost panel. Top row: classical model with radial Hubble expansion (see Eq. (8) and Lee et al. 2000). Middle and bottom rows: modified model with “collimated” expansion (see Eq. (11)). Green, blue, and yellow contours highlight specific height ranges, indicated in the first column. White dotted lines outline the radial boundary defined by the shell ellipses with increasing velocity. Model parameters are listed in Table 2. δx and δz units are arcseconds. 
5 Diskwind modeling
In this section, and alternatively to the WDS models considered in Sect. 4, we investigate a simple kinematical disk wind model for the DG Tau B redshifted outflow where the conical morphology visible in Fig. 3 would trace the trajectory of CO molecules ejected from the disk. Although we cannot derive V_{R} in a model independent way, from the external contours of transverse PV diagrams (see Sect. 3), we showed in Sect. 2 the existence of brighter elliptical structures visible on transverse PV diagrams. Assuming that they trace each a specific layer of the flow along which V_{Z} and V_{R} stay roughly constant with height, we fit these ellipses to derive both the V_{R} and V_{Z} components of the velocity (see Appendix C). Figure C.1 shows that the derived velocity directions are parallel to the conical contours of constant V_{Z}. This comforts our hypothesis that the trajectory of the outflow follow lines of constant V_{Z}. From this hypothesis, we derived the collimation and kinematics of the streamlines using the tomography and created a synthetic data cube of the conical disk wind outflow.
Fig. 8 Kinematic evolution along the flow axis. Panels a: positionvelocity diagrams averaged over a slice of width ∆Z = 0.2″ at three different δz positions along the flow. The red contours trace solutions of collimated WDS used to fit A2, A3, and U5 (third line on Fig. 7). Panel b: in black is represented the maximal velocity of emission of the PV diagram. In red is shown the center velocity of the ellipses on the PV diagrams. Panel c: the gray region highlights the conical outflow domain. The black lines represent the parabolic morphology of the three solutions. The red dots show for each solution the region where the projected velocity reach V_{max}. 
5.1 Steady diskwind model
We made the assumption that the flow is axisymmetric and that the matter has reached its terminal velocity and has a constant poloidal velocity along its trajectory. We fit this trajectory by a conical surface defined by an angle θ from the Zaxis and an anchoring radius r_{0,geo}. We extracted from the tomography the specific angular momentum, j = R × V_{ϕ}, along curves of constant V_{Z} (see Fig. 4) and derived a median value for each streamline. We defined the uncertainty of this value as the standard deviation of specific angular momentum. We also computed the poloidal velocity V_{P} = V_{Z}/cos θ.
Figure 9 represents the derived values of r_{0,geo}, θ and V_{P} × j for each streamline of constant V_{P}. We fit the variation of anchoring radius r_{0,con} with the poloidal velocity by a polynomial law. The variation of ejection angle θ was fit by a power law. V_{P} × j was taken constant at 570 ± 50 au km^{2} s^{−2} for all the streamlines. The fits were achieved using nonlinear least squares, the equations are as follows: (12) (13) (14)
Here r_{0} is in au, V_{P} in km s^{−1}, θ in degrees, and j in au km s^{−1}. We modeled the disk wind with axisymmetric conical streamlines with the dynamics and morphology laws derived in Fig. 9 and created a synthetic data cube of the conical outflow. We set the external and slower layer at V_{p} = 6 km s^{−1}, corresponding to the smallest value that could be mapped with our tomography (see Sect. 3. We set the internal, faster velocity at 20 km s^{−1}. The parameters in the velocity range V_{p} = 14–20 km s^{−1}, not covered by the tomography, were determined from an extrapolation of our fits (blue dotted line in Fig. 9). This extrapolation was done in order to describe the almostvertical highvelocity component not described by the tomography due to insufficient S/N. For each layer, we set the initial value of V_{ϕ}(R) assuming V_{p} × j = 570 au km^{2} s^{−2} (see Fig. 9). We assumed optically thin emission throughout the outflow, which is justified by the observed ratios of ^{13}CO/^{12}CO (see DV20). We did not consider a variation of emissivity with radius of ejection nor with height (see Appendix F). Proper modeling would require CO chemistry and temperature profiles, which is well beyond the scope of this paper. Projection and beam convolution effects were also taken into account.
Figure 10 shows synthetic channel maps and PV diagrams for our model compared with observations. The global morphology of the outflow at (V − V_{sys}) > 2 km s^{−1} as well as its variation with lineofsight velocity are well recovered as expected since we use the tomographic results to constrain the wind collimation and kinematics. This model does not attempt to describe the extended outflow surrounding the cone at low velocities (V − V_{sys}) = 1.15 km s^{−1} (see Sect. 2, Fig. 10). To describe completely this extended lowvelocity component, we would need to extrapolate the diskwind model at larger ejection radii. However, as this component falls at absorbed cloud velocities, we are not able to derive modelindependent constraints on the dynamics of this component. Proper modeling would require timeconsuming and uncertain parameter space exploration. We choose therefore to focus in the following on the conical outflow; nonetheless, the disk wind could be more extended radially than we describe with our current modeling.
Although effective to describe the global morphology of the outflow, our simple axisymmetric and steady diskwind model does not reproduce the different substructures identified in our observations: cusps and arches and the local deviations of specific angular momentum at Z ≈ 600 and 1000 au. In the following subsections, we discuss two small perturbations of our diskwind models which could explain the various substructures observed.
Fig. 9 Diskwind properties derived from the tomography along lines of constant V_{Z}. Red dashed curves show fits as a function of the poloidal velocity of the streamline V_{p}: the anchoring radius of the streamline r_{0,con} was fit by a thirdorder polynomial, the angle of the streamline with the flow axis θ by a power law and the product V_{P} × j by a constant value of 570 au km^{2} s^{−2}. Blue dotted lines represent the extrapolation used to model the highvelocity component that could not be mapped by tomography due to its low S/N. 
Fig. 10 Comparison of observations with steady diskwind model. Left panels: individual ^{12}CO channel maps computed from the global diskwind model at selected lineofsight velocities (top row) are compared to observations (bottom row). The color scale is the same for all the channel maps. Right panels: transverse PV diagrams at two positions δz along the flow and averaged over a slice of width ∆Z = 0.2″. The background grayscale image shows the observations and the red contours trace the predictions from the diskwind model. (V − V_{sys}) units are km s^{−1}. 
5.2 Wiggling of the flow axis
Although we do not detect a clear signature of wiggling in our data, we cannot exclude a small amplitude wiggling of the CO axis <0.5° (see Appendix A). A wiggling of the outflow axis could explain the variations observed on the specific angular momentum tomographic map. Indeed, in order to determine the specific angular momentum using Eq. (7), we assumed that the center of the layer is located at δx = 0 at all heights. If the center is shifted toward δx > 0 or δx < 0, the specific angular momentum computed with our method will be respectively higher or lower than the true value. This effect is more critical if the PV diagram shows a strong velocity gradient, which is the case for DG Tau B. In this section, we investigate this effect, and show that small amplitude wiggling may also create the substructures observed (cusps, fingers and arches).
We modifed the diskwind model presented in Sect. 5.1 to add a precession of the outflow axis. Each conical layer of the outflow precesses with an angle α and a precession period τ_{p}. This is an extension of the model developed by Masciadri & Raga (2002) for jets, modified to take into account the conical morphology of the outflow and inclination to the plane of the sky. We modeled both a prograde and retrograde precession. However, due to the small value of α in our models, the two models give very similar results. We present here results for the prograde model only. We first assumed that all diskwind layers precess with the same α and the same precession period τ_{p}. Due to the velocity shear across the outflow, the spatial period Λ_{p} = τ_{p} × V_{p} then varies between layers according to the poloidal velocity. We also investigated a precession model where the spatial period Λ_{p} is constant across all streamlines. A constant Λ_{p} corresponds to a variation of precession period as . We visually fit τ_{p} and Λ_{p} to best reproduce the location of the two extrema variations in the specific angular momentum map separated by ≃400 au. For each model, we computed synthetic data cubes and derived the specific angular momentum map using the same method used in Sect. 3 for the observations.
Figure 11 shows the resulting channel maps and the specific angular momentum maps for the two precession models (constant τ_{p} and constant Λ_{p}) with two different precession angles (0.2 and 0.5°) compatible with the upper limit derived for the CO outflow axis wiggling in Annex A. Precession models with constant τ_{p} = 400 yr successfully reproduce the channel maps morphology, in particular the cusps at highvelocity and arches at low velocities (see Fig. G.1). A best match to the intensity contrast is obtained for α = 0.5°. However, the resulting map of specific angular momentum j is not fully consistent with our observations. Indeed, as Λ_{p} is different for each layer, the perturbations of specific angular momentum are not localized at one specific height, as in the observations.
The modified model with constant Λ_{p} = 800 au for all layers better reproduce the positions of the two extrema at Z = 600 and Z = 1000 au in the specific angular momentum map. However, the cusps have a lower intensity contrast than observed, even with the maximum allowed wiggling angle of 0.5°. In addition, this model predicts clear detectable wiggling on the edges of the cone in channel maps, which is not seen in the observations. A model in between these two extremes, that is to say with α ≃ 0.5° and a precession period τ_{p} increasing more slowly than may better account for all observational properties.
A remaining discrepancy with observations is that none of the wiggling models reproduce the short spacing of ≃1.5″ = 200 au between the inner cusps, as well as the apparent increase of cusp separation with distance from the source (see Table 1, Fig. G.1), although this latter effect is mostly seen in the farthest cusp A0 and may result from a lack of sensitivity. The constant Λ_{p} model predicts a projected separation between the cusps of Λ_{p} sin(i), corresponding to ≃5″, while the constant τ_{p} model a twice smaller separation typically.
However, we stress that our wiggling models are probably too simplistic as they do not take into account the (magneto)hydrodynamical interactions between the layers. Masciadri & Raga (2002) have shown that simulations depart rapidly from analytical solutions in the case of jet wiggling due to precession. This difference is potentially even greater with a shearing outflow. Dedicated numerical simulations are required to fully test this scenario. Nonetheless, this model is a promising candidate to explain the variation of specific angular momentum along the DG Tau B outflow. We discuss in Sect. 6 possible wiggling mechanisms and their implications.
Fig. 11 Synthetic data cubes computed from the generic diskwind model with flow axis precession: with a precession angle of 0.2° (left panels) and 0.5° (right panels) and assuming a constant precession period τ_{p} = 400 yr (top row) or a constant precession spatial wavelength Λ_{p} = 800 au (bottom row). For each model are shown both a channel map at (V − V_{sys}) = 4.37 km s^{−1} and the corresponding tomography of j derived with the method described in Sect. 3. The color scale of j ranging from 0 to 140 au km s^{−1} is identical to the color scale of Fig. 3. 
5.3 Emissivity enhancements in the disk wind
In this section, we investigate an alternative model where cusps and arches are created by localized emissivity enhancements in the conical disk wind. We first derived their location from the tomography, assuming that they are axisymmetric, and then created a synthetic data cube to compare with our observations.
Figure 12 shows the projected velocity onaxis at δx = 0 for the front side and the back side () of each conical wind layer in our model, predicted from Eq. (3). We also represent the domain of lineofsight velocities where arches and cusps positions could be measured in channel maps as described in Sect. 2. The cusps observed at (V − V_{sys}) > 5 km s^{−1} could not be characterized in Sect. 2 because of low S/N.
Figure 12 shows that if the observed substructures are due to axisymmetric emissivity enhancements in the conical outflow, the cusps come from the back side () and the arches from the front side () of the enhanced ring. Figure 12 also shows that wind layers with poloidal velocities V_{p} > 9 km s^{−1} are not located on the arch domain, and are located on the low S/N cusp domain, which would require weaker emissivity enhancements in the fastest internal layers.
Each cusp was characterized by its projected height at δx = 0 at a specific projected lineofsight velocity V_{los}. The projected height of the cusp corresponds to an equation Z(R) derived from Eq. (2), with . Similarly, the projected velocity of the cusps corresponds to a poloidal velocity as shown in Fig. 12 and assuming V_{R}, a conical line of constant V_{Z} in the tomography. As a result, we can associate the cusp observed in a given channel map with a (R, Z) location on the V_{Z} tomography map. This location would be the intersection between the Z(R) equation derived from the projected height, and the conical line derived from the projected velocity.
Figure 13 shows the tomography of the conical outflow using equations from Fig. 9 for extrapolation at Z > 1200 au. The white dots correspond to the solutions (named SU0 to SU5) for the cusps heights on axis identified in the channel map at (V − V_{sys}) = 4.37 km s^{−1}. We were able to follow six different cusps in up to six different channel maps (see Fig. B.1), which allowed us to reconstruct the shape for some of these enhancements, shown as hatched areas in Fig. 13. As mentioned in Appendix B, the cusps are also visible at higher lineofsight velocities, but could not be characterized reliably due to their lower S/N. Consequently, the hatched area shown in Fig. 13 should be extended toward the inner streamlines. Moreover, density enhancements could also be present at δz < 2.2″ but were not identified by our procedure. The apparent acceleration of the cusps seen in channel maps can be readily reproduced, in this diskwind model, by axisymmetric emissivity enhancements that cross obliquely the flow streamlines. The upward shift of each cusp at higher velocity (apparent acceleration) is simply a result of the velocity shear across flow streamlines. It does not require a Hubblelaw dynamics in the underlying flow.
Interestingly enough, the density enhancements SU34 are located close to the two extrema of specific angular momentum, at Z = 600 and 1000 au, suggesting a possible link. We could not derive the shape of the emissivity enhancements for V_{p} > 9 km s^{−1}, due to the low S/N of the cusps at the corresponding projected velocities (V − V_{sys}) > 5.32 km s^{−1} (see Fig. 12).
We computed the dynamical age of these emissivity enhancements with . This value would give the true dynamical age of the enhancement if it is created by variability of ejection from the disk surface. The derived values of τ_{w} for cusps visible on multiple channel maps are shown in Fig. 14. The dynamical age is almost constant within each cusp with only a slight decrease with increasing poloidal velocities, especially for U3 and U4. The difference in dynamical ages between two successive cusps varies from 190 to 490 yr and roughly increases linearly with the age.
These solutions were determined using only the properties of the cusps onaxis (at δx = 0), corresponding to . In order to check if the reconstructed enhancements were consistent with the full cusp morphology in all channel maps, we developed a 3D model where we modified the emissivity profile of the synthetic diskwind model presented in Sect. 5.1. We multiplied the underlying DW emissivity profile by Gaussian components representing the emissivity enhancement. We used the mean τ_{i} values derived from Fig. 14 to determine the location of the enhancements along each wind streamline. In order to reproduce the observed width and intensities in the channel maps, we set the full width at half maximum of the Gaussian component at 23 yr and the maximal emissivity enhancement at G_{i} = 3. We introduced these enhancements only in the external layers of the outflow (V_{P} < 9 km s^{−1}), since we did not have constraints on their location at higher velocities. Figure 15 shows computed and observed channel maps at different lineofsight velocities. The contour of the modeled enhancement SU0 is also represented on top of the observations. This model successfully reproduces the morphology of the cusps as well as their apparent offset from source with increasing velocity. Interestingly enough, the model used to reproduce the cusps also matches the arches at low velocity, with the locations of the apexes consistent with our observations (see Fig. G.1). The intensity contrast is also roughly recovered with our G_{i} value, meaning that the outflow brightness is locally multiplied by three.
The channel map at (V − V_{sys}) = 1.51 km s^{−1} shows that only the top of the arches is reproduced. Indeed, the flanks of the largest arches are wider than the conical flow region modeled by our disk wind in Fig. 9. In order to fit completely the arches, we would need to extend the diskwind model to larger radii and opening angles and lower poloidal velocities as suggested in Sect. 5.1, and extend the density enhancements to these regions. However, due to cloud absorption at low velocities <2 km s^{−1}, we do not have modelindependent tomographic constraints on the streamline shape and kinematics in this slow external flow. Therefore we cannot determine a reliable solution for the emissivity enhancements producing these arch flanks.
Fig. 12 Lineofsight velocity onaxis (at δx = 0) of the model diskwind flow surfaces as a function of their poloidal velocity. Orange circles show the front side (), and green circles show the back side () of the diskwind flow surfaces. In gray, we show the lineofsight velocity range where the arches and cusps are observed (see Appendix B). In light gray, is represented the low S/N domain where the cusps could not be characterized. 
Fig. 13 Tomographic map of V_{Z} from Fig. 3 extrapolated from Z = 1200 au to Z = 3000 au using the conical wind model from Fig. 9. Gray hatched regions named SU05 represent the locus of emissivity enhancements required to explain the cusps onaxis positions at V − V_{sys} = 3.73–5.32 km s^{−1} (see Fig. B.1). White dots represent the solutions for the cusps at V − V_{sys} = 4.37 km s^{−1}. Red dotted lines indicate the height of the two extrema in specific angular momentum (from Fig. 3). The gray points located on the R = 0 axis represent the positions of jet knots observed by Podio et al. (2011). The red dots indicate the positions with associated errors of the ellipses identified in transverse PV diagrams (see Appendix C) and indicate that emissivity enhancements extend to the inner streamlines. 
Fig. 14 Dynamical evolution of the enhancements. Left panel: dynamical age of the enhancements SU05 as a function of V_{P}. The red dotted line represents the mean value. Right panel: dynamical age versus timescale between adjacent cusps (τ_{i} − τ_{i+1}). The uncertainties are determined by propagating the errors on the cusp positions. 
Fig. 15 ^{12}CO channel maps at selected lineofsight velocities for the steady conical diskwind model (bottom row) and observations (top row). White contours highlight the predicted emission from the modeled density enhancement SU0. The intensity color scale is the same for all channel maps. 
6 Discussion
In this section, we use the results of our parametric modeling of the flow (Sects. 4, 5) and tomographic study to critically discuss two possible origins for the smallscale redshifted CO outflow in DG Tau B: 1) a stacking of multiple shells sweptup by an inner wind (or jet), without any contribution from an extended disk wind (Sect. 6.1), 2) an extended disk wind (Sect. 6.2), with internal perturbations causing the observed substructures (Sect. 6.3). We find that our measurements of rotation put stringent constraints on each of these scenarios, and we discuss physical implications for the ejection process and its relation to the disk accretion process.
6.1 Stacking of winddriven shells
We discuss here a scenario where the redshifted CO outflow in DG Tau B can be accounted for by a stacking of multiple sweptup shells resulting from the interaction of an episodic wideangle wind (or jet) with the ambient medium. Such a scenario is recently proposed by Zhang et al. (2019) to reproduce the multiple shell structures at the base of the HH46/47 outflow. They find good agreement with the WDS model of Lee et al. (2000, 2001), namely a parabolic layer undergoing radial expansion following a Hubble law V ∝ r. Below we summarize the key successes and failures encountered in Sect. 4 by the same WDS model when applied to the various flow components in DG Tau B, and we show that our rotation measurements raise additional issues for this model in terms of cavity refilling.
6.1.1 Outer flow component
In Sect. 4.1, we find that the morphology and kinematics of the lowvelocity outflow can be well reproduced by the simple WDS model of Lee et al. (2000). We now examine whether a WDS origin is physically consistent with the large specific angular momentum inferred from our modeling, j_{outer} ≃ 250 ± 50 au km s^{−1}, in the same sense as the disk (see Sect. 4.1).
Studies of rotation signatures in prestellar and protostellar cores show that specific angular momentum stops decreasing with radius below scales ≤3000 au and becomes roughly constant (Ohashi et al. 1997; Gaudel et al. 2020). This “plateau” is interpreted as the region where infall motions start to dominate, and specific angular momentum is roughly conserved along streamlines. Depending on the object, the specific angular momentum in the “plateau” is ≈40–400 au km s^{−1} (0.2–2 × 10^{−3} km s^{−1} pc, Gaudel et al. 2020). Our estimate for the outer flow in DG Tau B, j_{outer} ≃ 250 ± 50 au km s^{−1} falls well within this range. In addition, infall signatures are identified around DG Tau B at large polar angles θ ≃ 70^{°} (DV20). No such signatures are seen at smaller polar angles, but rotational flattening predicts lower envelope densities there (Ulrich 1976), hence they might be too faint for detection. It thus appears promising to consider that infalling material might dominate the rotation in the outer CO layer.
Strictly speaking, the Hubble law assumed in the WDS model of Lee et al. (2000) is only valid for a static ambient medium with a 1/r^{2} radial density decrease^{3}. A rotating infalling envelope, in contrast, has a nonradial motion and a flatter radial density law ∝ 1/r^{1.5} (Ulrich 1976). However, the calculations of LópezVázquez et al. (2019) for such an ambient medium show that the WDS expansion remains quasi radial, except close to the midplane, and with almost constant speed after 200 yr. Our simple model in Sect. 4.1 thus remains roughly valid if the wind expands into an infalling envelope.
Using the shell rotation speeds computed by LópezVázquez et al. (2019), we expect a shell specific angular momentum close to that of the infalling material immediately ahead of it. Therefore, we consider in Fig. 16 the spatial distribution of specific angular momentum in a rotating, freefalling envelope from Ulrich (1976), with a centrifugal radius R_{d} = 700 au and central mass M_{*} = 1.1 M_{⊙} appropriate to DG Tau B (DV20). The predicted j values along the fit parabolic outer flow boundary are very similar to the observed one, j_{outer}. This detailed comparison confirms that an infalling envelope in DG Tau B, if present up to polar angles θ ≃ 30°, could provide enough specific angular momentum to explain the rotation in the outer CO layer.
However, this analysis is highly simplistic. First, we consider here that the specific angular momentum of the envelope is fully transferred to the entrained layer. This assumption gives an upper limit for the shell rotation velocity, as turbulent mixing would decrease its specific angular momentum. Secondly, the spherical and ballistic infall model used here (Ulrich 1976) does not take into account the effects of pressure gradients or the magnetic field. Pressure gradients could potentially increase the specific angular momentum of the infalling envelope at large polar angles, through “pushing” outer infalling streamlines toward the axis, while the magnetic field would decrease it due to magnetic breaking. Dedicated numerical simulations of the interaction of an infalling material with an inner wind component and taking into account all these effects are needed to fully test the entrainment scenario for the outer CO layer.
A serious issue with this interpretation, however, is the young inferred shell age, τ = V_{Z}/Z = 6000 yr (Sect. 4.1), much younger than the true age of DG Tau B. A first way out of this “short age problem” would be that the interface between wind and envelope in DG Tau B is not expanding, but static. A static shell may form when mixing between the wind and the ambient material is not instantaneous, as assumed in most WDS models, but very gradual. Shocked ambient material is then slowly entrained along the shell surface by the shocked wind in a thin turbulent mixinglayer. The static shell shape and mixinglayer properties were recently computed in the case of a freefalling rotating envelope by Liang et al. (2020). Using again R_{d} = 700 au and M_{*} = 1.1M_{⊙} for DG Tau B (DV20), the specific angular momentum in the mixinglayer is predicted to be au km s^{−1}, twice lower than estimated in the outer flow. Therefore, a static wind/envelope interface does not seem able to explain the rotating outer flow.
A second way out of the short age problem would be that the outer flow component does not trace the first wind encounter with the infalling envelope, but a more recent wind outburst from 6000 yr ago. To provide its high angular momentum material to the young shell, however, the infalling envelope should somehow manage to penetrate and “refill” all the older shells created by previous (unobserved) wind outbursts. Whether such an efficient cavity refilling by the envelope is physically possible on ≤6000 yr timescales is a difficult open question, well outside the scope of the present paper. As shown below, the issue of cavity “refilling” becomes even more acute when the WDS scenario is applied to the inner conical flow.
Fig. 16 Specific angular momentum map (in color) for the model of ballistic rotating infall (Ulrich 1976) with R_{d} = 700 au and M_{*} = 1.1 M_{⊙}. The angular momentum varies from 25 au km s^{−1} to 900 au km s^{−1}. The green hatched and red filled regions represent respectively the limits of the conical outflow derived from the tomography in Sect. 3 and the parabolic shape of the classical Hubblelaw WDS model fit to the external CO outflow. The white and red contours outline the infalling streamlines with a specific angular momentum similar to the conical outflow j = 40–100 au km s^{−1} (white), and to the WDS model j = 200–300 au km s^{−1} (red). The black hatched region outline the streamlines with initial θ_{0} = 70 ± 5° reproducing the infall signatures seen in DG Tau B (DV20). 
6.1.2 Inner conical outflow
In contrast to the outer flow component, we find in Sect. 4.2 that the morphology and kinematics of the inner conical outflow and its bright substructures (arches, fingers and cusps) cannot be reproduced by the model of parabolic WDS with radial Hubble law used by Lee et al. (2000) and Zhang et al. (2019), even after several ad hoc modifications. We identify two serious issues: (1) The observed aspect ratio of arches in channel maps at low velocity is significantly shorter than predicted by the original model (ellipse aspect ratio =1/cos i); it can be reproduced by a more collimated WDS model where the flow is parallel to the parabola; but that is no longer physically consistent with a Hubblelaw velocity field, which requires an expanding shell (Shu et al. 1991). (2) The shell models fitting the arches at lowvelocity and cusps at high velocity do not agree with the observed emission morphology in channel maps at midvelocities, predicting fingers that are broader than the cone, or full ellipses that are not seen.
We find that assuming a conical shell instead of a parabolic one, but keeping a Hubble law, still creates the same problems. They appear intrinsically caused by the Hubblelaw velocity field, regardless of the shell detailed morphology. Therefore, any model where the shell is expanding quasiradially and at nearly constant speed over time will fail to reproduce our observations. This, in particular, discards all models where the wideangle wind and the ambient medium share a similar powerlaw in r, and where they instantly mix in the shell (e.g., the models of LópezVázquez et al. 2019; Shang et al. 2006).
Simulations including magnetic field (Wang et al. 2015; Shang et al. 2020) and stationary solutions (Liang et al. 2020) show the formation of a shear layer along the shell more in line with DG Tau B. However the first model predicts a shell anchoring radius increasing with time, while the second has a shell anchored near the centrifugal radius R_{d} ≃ 700 au (DV20). This is inconsistent with the small observed anchoring radius of the DG Tau B conical outflow (≤50 au, DV20). Alternative models of sweptup shells exist involving an infalling sheet (Cunningham et al. 2005) or a jet instead of a wide angle wind (e.g., in Downes & Cabrit 2007), but they have no analytical solutions. Therefore, dedicated numerical simulations would be required to test them in DG Tau B.
In the following, we show that the specific angular momentum measured in the conical outflow by tomography, j_{cone} ≃ 40–100 au km s^{−1} (see Fig. 3), raises additional issues for the sweptup shell scenario. We first note that a wideangle “Xwind” cannot explain the observed rotation in the conical outflow; with a launching radius ≃0.05–0.1 au from the central protostar and a wind magnetic lever arm parameter λ = j_{X}/j_{Kep} ≃ 3 (Shang et al. 1998), its specific angular momentum is predicted at j_{X} ≈ 20–30 au km s^{−1}, a factor two to three times less than observed. In addition, the wind cannot dominate the sweptup shell mass unless it is slower than twice the shell speed^{4}, which in the present case would require V_{w} ≤ 15 km s^{−1} (see Table 2). This is inconsistent with a wind originating from close to the protostar. The “Xwind” model, for example, has V_{w} ≃ 150 km s^{−1} (Shang et al. 1998).
The low observed expansion speeds ≃6–14 km s^{−1} in the conical flow (see Fig. 3) is more consistent with jet bow shocks dominating the shell mass. A jet magnetic lever arm parameter λ ≃ 9 would then provide enough angular momentum. Such a scenario, however, cannot explain why regions of lower speed in the conical flow have inversely higher specific angular momentum (see Fig. 3). In a jet bow shock, lower speeds arise where more ambient mass has been sweptup. Assuming the ambient medium provides no angular momentum, the jet angular momentum would get more diluted, and the shell specific angular momentum would drop, instead of increasing. We conclude that the observed rotation in the conical flow cannot come from an inner wind or jet.
To reproduce the observed conical flow rotation in the sweptup shell scenario, we thus need an external medium with an important angular momentum. The infalling rotating envelope is an obvious candidate. However, the observed specific angular momentum in the conical flow is twice larger than predicted, in the same region, by the Ulrich infalling envelope model (see Fig. 16). In addition, it is unclear how infalling matter could penetrate and “refill” the space between the closely spaced shells producing the cone substructures, especially when the region immediately outside the conical flow is instead in outflow motion (see Fig. 5).
An alternative would be that the sweptup material originates from the rotating disk atmosphere, at radii R_{0} ≃ j^{2}/(GM_{*}) ≃ 2–9 au. The problem is then to refill the cleared cavities between shells with a “new” static disk atmosphere. We note that our three shell models fitting the substructures in the conical flow have remarkably identical expansion speeds within 1 km s^{−1} (see last column of Table 2 for WDSA2, WDSA3 and WDSU5). Assuming that the corresponding wind and jet outbursts were of similar strength, it implies that they met an identical ambient density ahead of them, hence the atmosphere refilling process should be extremely efficient. This appears difficult to achieve unless a largescale disk wind is present.
Realistic simulations of the interaction between an episodic inner wind or jet with an infalling rotating envelope and the disk atmosphere will be required to definitely exclude that sweptup shells with enough angular momentum and appropriate kinematics could be generated. However, at this stage and taking into account all the abovementioned difficulties, we do not favor this scenario as the origin of the smallscale rotating CO outflow in DG Tau B. In the next section, we discuss an alternative scenario where this rotating outflow traces a (perturbed) extended disk wind.
6.2 The diskwind scenario: constraints on the driving mechanism
We therefore favor the scenario in which the inner CO conical outflow traces matter directly ejected from the disk. We show in Sect. 5 that the stratified kinematical structure derived for the conical outflow is suggestive of a quasisteady diskwind. In this section, we discuss constraints on the driving mechanism. Disk winds come in different flavors, depending on the main physical mechanism responsible for driving the flow: pure thermal effects in photoevaporated disk winds (PDW; Alexander et al. 2014), cold magnetocentrifugal ejection (Pudritz et al. 2007) or a combination of the two processes in the socalled warm or magnetothermal disk winds (Casse & Ferreira 2000; Bai et al. 2016). In the following we refer to the two last classes of magnetized disk winds as MHD disk winds (hereafter MHD DW).
6.2.1 Photoevaporated disk winds
In photoevaporated disk winds, the high energy radiation (UV and Xrays) of the central accreting protostar heats the surface layers of the disk to high temperatures (10^{3}–10^{4} K) up to significant radial distances. Beyond the gravitational radius , where c_{s} is the sound speed in the upper disk surface layers, thermal energy exceeds the gravitational binding energy. Numerical simulations show that significant mass loss starts before r_{g}, at the critical radius r_{c} ≃ 0.1–0.2r_{g}. Consequently, matter is ejected and reaches terminal velocities of typically two to three times c_{s} (Alexander et al. 2014, and references therein). The exact properties of the wind depends on the dominant source of highenergy irradiation. ExtremeUV (EUV) heating creates an isothermal ionized layer on the disk surface with temperature T ≃ 10^{4} K (c_{s} = 10 km s^{−1}) which drives a fast wind (~35 km s^{−1}) with massloss rates 10^{−9}–10^{−10} M_{⊙} yr^{−1} (Font et al. 2004; Wang & Goodman 2017). Xray irradiation results in cooler and slower flows (c_{s} ≃ 3–5 km s^{−1}, v ≃ 15–20 km s^{−1}) but penetrates at higher densities and therefore can drive massloss rates up to a few 10^{−8} M_{⊙} yr^{−1} (Picogna et al. 2019). Non ionizing farUV (FUV) heating mostly drives slow massloss (v ≃ 18 km s^{−1}) from the outer disk regions.
The conical morphology of the DG Tau B CO outflow matches expectations from both selfsimilar PDW models by Clarke & Alexander (2016) and hydrodynamical simulations by Owen et al. (2011) and Wang & Goodman (2017). In such diskwind solutions, the specific angular momentum is conserved along the streamlines and is equal to the Keplerian value at the footpoint radius. The observed values of j in the CO conical outflow imply footpoint radii in the range r_{0} = 1.6 to 8.2 au. Figure 17 shows that these footpoint radii are lower (by a factor ≃three) than the ones derived from a straight extrapolation of the conical morphology on large scales. This requires a larger opening angle of the streamlines at the base, consistent with the selfsimilar models of Clarke & Alexander (2016).
The values of r_{0} are in line with expectations from EUV dominated PDW models. Indeed, for a 1M_{⊙} star and c_{s} ≃ 10 km s^{−1}, expected in EUV heated PDW, the gravitational radius r_{g} ≃ 10 au and significant mass loss starts at r_{c} ≃ 1–2 au. On the other hand, observed terminal velocities of V_{p} = 4–16 km s^{−1} require sound speeds ≤2–8 km s^{−1} at radii r_{0} ≤ 10 au, excluding EUVdriven models.
The large derived mass flux of 1.7–2.9 × 10^{−7} M_{⊙} yr^{−1} for the CO conical flow (DV20) within r_{0} ≤ 10 au however excludes FUV driven winds, which fail by at least one order of magnitude (Wang & Goodman 2017). We show below that it is also inconsistent with the latest Xray driven models. Indeed, recent Xray driven photoevaporation models of Picogna et al. (2019) predict mass loss rates up to 10^{−7} M_{⊙} yr^{−1} for stellar Xray luminosities L_{X} ≥ 10^{31} erg s^{−1}. However, these high L_{X} models also predict a stronger contribution of the outer disk regions to the total mass flux due to the increased penetration of Xray photons. Figure 9 in Picogna et al. (2019) shows that only 10% of the total mass flux originates from disk footpoint radii below 10 au. In summary, current PDW models fail to account for the combination of large mass flux and small footpoint radii of r_{0} = 1.6–8.2 au, derived for the DG Tau B CO conical flow.
Last but not least, the survival of CO molecules in such a wind is problematic. The full thermochemical computation of Wang & Goodman (2017) shows that in their fiducial models, CO survives only at the very base of the wind in an intermediate layer on scales z/r ≤ 0.6. However, the models of Wang & Goodman (2017) are EUV dominated and hence result in warm and fully ionized winds. Similar problems are expected in thermally driven winds launched from the inner disk, which require base temperature greater than 2000 K. Therefore, we conclude that pure thermal processes appear highly unlikely as the main driving mechanism of the DG Tau B CO conical wind.
Fig. 17 Launching radius of the streamline as a function of its poloidal velocity assuming: conical extrapolation for the geometrical radius (green circles), steady thermal ejection conserving angular momentum (orange circles) and steady cold magnetocentrifugal ejection (blue circles). See text for more details. We do not derive the launching radius for the largest poloidal velocities due to the large uncertainties. 
6.2.2 Magnetic disk winds
Disk winds driven by magnetic forces require a large scale poloidal magnetic field anchored in the disk. This largescale field exerts a torque on the rotating disk that both ejects matter and removes angular momentum from the disk (Blandford & Payne 1982; Pudritz et al. 2007). The strength of this torque is characterized by the magnetic lever arm parameter λ ≃ (r_{A}/r_{0})^{2}, where r_{A} is the poloidal Alfven radius and r_{0} the disk footpoint radius of the streamline. In principle, such disk winds can produce at the same time fast and collimated jets originating from the inner streamlines and much slower and less collimated winds originating from outer disk radii. The full kinematics and morphology of these solutions also depend on whether thermal effects are important in the launching regions. Numerical simulations show that the mass loss can be significantly increased when thermal effects are taken into account (Casse & Ferreira 2000; Bai et al. 2016). Such magnetothermal winds have low to moderate λ values but can extract significant mass and angular momentum from the disk.
Terminal velocities depend on both the footpoint radii r_{0}, λ and thermal effects. Under the assumption of steady magnetically driven ejection, the asymptotic values of V_{p} and V_{ϕ} are given by (Eqs. (4) and (5) in Ferreira et al. 2006): (15) (16)
where β encompasses all pressure effects, including thermal and turbulent Alfvén waves (see Ferreira et al. 2006). The streamline footpoint radius can be estimated from these equations assuming cold MHD ejection, ie. negligible thermal effects (β ≈ 0) following Anderson et al. (2003). Figure 18 traces the relationship between the mean values of V_{P} and j = RV_{ϕ} for the various conical layers of constant V_{Z} derived from the tomography. In this figure is also represented the parameter space (r_{0},λ) predicted by cold magnetocentrifugal diskwind models (Ferreira et al. 2006). The mean poloidal velocities and specific angular momentum coincide with a line of constant λ_{cold} ≃ 1.58 with footpoint radii r_{0,cold} = 0.7–3.4 au. If thermal effects play a dynamical role at the base of the wind, Eqs. (15) and (16) show that the values of λ and r_{0} derived under the cold assumption are respectively upper and lower limits. This effect is illustrated in Fig. 18 where we plot V_{P} and j from the warmest solution of Casse & Ferreira (2000) with λ = 1.9, at r_{0} = 1 au. We see that using the cold MHD curves would lead to overestimate λ ≈ 2. 5 and underestimate r_{0} ≈ 0.7 au. The cold assumption is only valid if β ≪ 2λ − 3 ≃ 0.2, which would require a cold disk atmosphere and no substantial wind heating. On the other hand, the low derived upper limit on λ ≤ 1.6 is consistent with warm MHD DW models (Casse & Ferreira 2000; Bai et al. 2016; Wang et al. 2019) or cold MHD DW from weakly magnetized disks (JacqueminIde et al. 2019).
The derived minimum footpoint radius of r_{in} ≥ 0.7 au for the CO wind streamlines is in good agreement with the thermochemical predictions of Panoglou et al. (2012), who demonstrate that CO molecules magnetically launched from footpoint radii ≥1 au survive in the case of accretion rates in the disk ≥10^{−7} M_{⊙} yr^{−1}. Similar results are obtained for warm magnetothermal wind solutions (Wang et al. 2019).
The streamline footpoint radii derived from the kinematics are significantly smaller than the radii obtained from direct conical extrapolation (Fig. 17), suggesting wider opening angle of the streamlines at their base. This is indeed expected in MHD DW solutions where streamlines originally follow a conical trajectory and recollimate on larger scales due to the hoop stress provided by the azimuthal Bfield. The constant opening angle of the streamlines observed out to Z = 3000 au suggests that the magnetic hoop stress drops rapidly above Z ≃ 50 au.
Contrary to pure thermal disk winds, MHD DW also account for the observed large mass flux in the DG Tau B conical CO flow. Wind mass loss rates in the range 10^{−8}–10^{−7} M_{⊙} yr^{−1} are predicted by the magnetothermal wind solutions of Wang et al. (2019), on the same order as the accretion rate in the underlying disk and increasing with disk magnetization.
We estimated the local ejection efficiency ξ defined as (Ferreira et al. 2006). We estimated ξ from estimates of the disk accretion rate at the inner launching radius of the CO outflow. Indeed, from the mass conservation across the disk region launching the conical CO outflow () and the definition of ξ, we derived the following expression: (17)
We estimated the accretion rate onto the central star by taking 10 % of the jet mass flux (Ellerbroek et al. 2013). Podio et al. (2011) estimate the red jet mass flux at 6.4 × 10^{−9} M_{⊙} yr^{−1}, giving a mass accretion rate onto the star of ≈6 × 10^{−8} M_{⊙} yr^{−1}. We took this value as a lower limit to . From the measured mass flux in the conical CO flow Ṁ_{DW} = 2.3 × 10^{−7} M_{⊙} yr^{−1} and the diskwind launching zone r_{out}/r_{in} ≈ 5, we then derived an upper limit on ξ ≤ 1. The mass flux in the conical wind is ~4 times the estimated accretion rate onto the star, implying that 80% of the mass accreting at r_{out} is being ejected before reaching the star.
If the transport of angular momentum in the disk is entirely provided by the torque exerted by the MHD DW, one expects in steady state the following relationship: ξ = 1/(2(λ − 1)). The upper limit derived above on ξ ≤ 1 translates into λ ≥ 1.5. This condition is compatible with our upper limit on λ ≤ 1.58 derived from the kinematics. Thus the CO mass flux combined with the constraints on launching radii and magnetic lever arm appear compatible with an MHD DW extracting all angular momentum required for the disk to accrete from r_{out} to r_{in}.
If the inner conical outflow is tracing a disk wind, then the outer parabolic outflow cannot be explained by the interaction between the envelope and an inner jet or Xwind. However, the outer outflow could be tracing the interaction between the envelope and outer diskwind streamlines located outside of the conical outflow (r_{0} > 4 au). The arches located at least partially outside of the conical outflow as well as the continuous aspect in the PV diagrams suggests an “intermediate” outflow located between the conical outflow and the outer parabolic surface. Alternatively, the outer flow could also be tracing directly the outer diskwind streamlines. The derived j and maximal velocity in the outer layer would indicate launching radii ≥30 au and a similar low λ value as for the inner cone (see Fig. 18). Moreover, Bai et al. (2016) show that some MHD DW solutions can accelerate until R ~ 100 r_{0}, possibly explaining the apparent acceleration of this component seen in our channel maps. The global outflow would then be a continuous MHD DW originating from 0.7 au to ≥30 au. Unfortunately, the morphology and dynamics of the potential disk wind originating from r_{0} > 4 au could not be studied in detail due to our limited spectral sampling and the absorption by the surrounding cloud or envelope. However, in that scenario, the origin for the difference of emissivity between the bright conical outflow and the faint outer flow is not clear but could reflect the radial distribution of magnetic field strength or surface density in the underlying disk.
Fig. 18 Specific angular momentum ( j = R × V_{ϕ}) versus poloidal velocity V_{P} for steady and axisymmetric MHD disk winds. The black symbols represent values derived from the tomography and averaged along lines of constant V_{Z}, with their uncertainty. Red curves show the expected relation from selfsimilar cold magnetocentrifugal disk winds with r_{0} varying from 0.5 to 100 au and λ between 1.52 and 2.3 from Ferreira et al. (2006). The green dot corresponds to a warm MHD disk wind solution from Casse & Ferreira (2000) with r_{0} = 1 au and λ = 1.9. The green arrow shows the path followed by the solution with an increase of r_{0}. The gray box represents the estimated specific angular momentum for the outer low velocity layer on Sect. 4.1, taking a maximal poloidal velocity corresponding to a height of Z = 3300 au in the outflow referential. 
6.3 Diskwind scenario: origin of perturbations
We discuss here the merits of different models for the origin of the substructures (arches, fingers, and cusps) in the conical flow, observed in the channel maps.
6.3.1 Perturbation by jet bow shocks
A first potential explanation for the bright substructures seen in the conical outflow is that the steady disk wind is perturbed by nested bow shock wings created by the propagation of the variable axial jet (see Fig. 19, scenario B). Perturbation of the inner streamlines of a rotating disk wind by a large jet bow shock is recently reported in the much younger system of HH 212 by Lee et al. (2021). In DG Tau B, this interpretation is supported by the similar spatial spacings between the inferred locations of the perturbations producing the substructures in the CO conical wind, and the axial jet knots identified in optical images, as shown in Fig. 13 (for z = 500–3000 au along the flow axis spacings range between 200 and 1300 au for the jet knots, 300–700 au for the overdensities). Indeed, in the jetwind interaction scenario, the overdensities trace the point of contact between each bow shock and the outer disk wind so they propagate along the interface at the bow shock propagation speed, which is similar to the inner jet knot propagation speed.
Although the optical knot observations are not synchronous with our ALMA CO observations, and jet knots move away from the source at the jet speed ≃150 km s^{−1}, the general pattern of knot spacing as a function of distance is set by the underlying jet variability properties, and thus will tend to remain similar at different epochs in a given jet. If, in addition, the jet undergoes lowamplitude wiggling (as frequently seen in young stars), perturbations to the disk wind caused by jet bow shocks would be slightly nonaxisymmetric, possibly explaining the apparent distortions in specific angular momentum along diskwind streamlines using tomography (see Fig. 4). Finally, this scenario might also explain the lack of recollimation of the conical diskwind streamlines, due to the additional internal pressure created by the jet driven bow shocks wings.
Hydrodynamical simulations of the interaction of a variable inner jet with a slower outer disk wind have been recently performed by Tabone et al. (2018). These simulations show the formation of a dense stationary conical layer closing down at the source, created by the stacking of jet bow shock wings in the disk wind. This shell exhibits local overdensities at the positions of individual jet bow shocks, illustrated by the red regions in Fig. 10 (right panels) in Tabone et al. (2018). These overdensities globally reproduce the observed shapes of the perturbations in the DG Tau B conical layer: for the perturbations closerin, the density map is dominated by the regions close to the bow shock apex which bend inward, while for the perturbations farther out the density map is dominated by the bow shock wings which bend outward. This simulation also shows that the conical dense shell mostly retains the velocity of the surrounding disk wind, because the shock is weak and oblique. Therefore, the simulations in Tabone et al. (2018) does not show the characteristic stratification in V_{Z} observed in the DG Tau B inner conical flow. The numerical simulations of Tabone et al. (2018) also reproduce the observed trend of similar spacings between the axial jet knots and the overdensities in the contact region (shown in red in their Fig. 10, right panels). However, the simulations are made in a simplistic configuration where the outer disk wind is assumed to have a constant vertical velocity of 40% of the jet speed, uniform density, and no rotation motion. More realistic simulations taking into account the velocity and density gradients across the outer disk wind, and including rotation and magnetic fields, are strongly needed to reliably test this scenario.
6.3.2 Possible origin of CO outflow axis wiggling
As discussed in Sect. 5.2, wiggling of the wind ejection axis is the only model explored here that can reproduce simply the variations of angular momentum observed along the DG Tau B outflow. From this analysis, we derived estimates of the wiggling period τ ≃ 400 yr and semiamplitude wiggling angle α ≃ 0.5°.
Wind axis wiggling can originate from the precession of the underlying disk angular momentum vector. Disk axis precession can be induced by a misaligned companion. In that scenario, orbital periods are expected to be significantly shorter than the disk axis precession period (Terquem et al. 1998), hence the companion would be located well within the disk. A companion in a misaligned orbit can open a gap in the disk, separating the dynamical evolution of the inner and outer disks (Zhu 2019). The inner disk then starts to precess with a period τ_{p} related to the orbital period of the companion by , where µ = m_{2}/(m_{1} + m_{2}) is the ratio of the companion mass to the total mass of the binary system, and assuming a small misalignment ip as suggested by the maximal wind precession angle α = 0.5° (see Eq. (27) in Zhu 2019). An additional constraint can be obtained by requiring that the semiamplitude wiggling angle α is dominated by the precession motion, which translates into the condition V_{0} ≤ V_{Z} tan(α), that is V_{0} ≤ 0.1 km s^{−1} where V_{0} is the orbital velocity of the flow source. Combining these two constraints, we derive a companion mass ratio µ ≤ 2.5 × 10^{−3} and binary separation: a ≤ 0.5 au. This separation is smaller than the launching radius of the CO disk wind (r_{0} ≥ 0.7 au, see Fig. 17). So the precessing disk launching the CO flow would be outside the orbit of the planetary mass companion, which is inconsistent with the scenario investigated here. Indeed, in such a scenario the outer disk is precessing on much longer timescales than the ones given by the formula above.
Nixon et al. (2013) and Facchini et al. (2018) investigate circumbinary disk precession around an inner misaligned binary system. In such circumstances the inner rim of the circumbinary disk can break from the outer disk and precess. However, large misalignments or massive enough companions are necessary for this situation to occur (Facchini et al. 2018). We show that this is not likely in the DG Tau B case. The inner binary system would truncate the disk at 1.5–1.7 times the separation a of the binary (Facchini et al. 2018). If we take an upper limit of r_{in} = 0.7 au for the inner circumbinary disk rim, to allow the launching of the CO flow, we get an upper limit of a ≤ 0.5 au for the binary separation, corresponding to an orbital period τ_{0} ≤ 0.35 yr using Third Kepler’s law and a total mass of 1 M_{⊙} for the binary system (DV20). With a small misalignment (less than a few degrees) between the outer disk and the binary orbital planes, suggested by the small wiggling angle of the CO outflow, a mass ratio of the companion µ ≤ 10^{−2} would be required to get a precession period of 400 yrs (see Eq. (4) in Facchini et al. 2018). Equation (2) in Facchini et al. (2018) then shows that such a low mass companion combined with a small misalignment will not break the inner circumbinary disk. Therefore this second precession scenario can be ruled out to explain the CO outflow wiggling.
So far we have investigated only precession as the origin of the wiggling of the CO flow axis. However, wiggling of the diskwind ejection axis can be also induced by the orbital motion of the CO outflow source in a binary system. The equations of motion of the ejected gas will be equivalent to the precession solution investigated so far. We follow the formulation developed by Masciadri & Raga (2002) and Anglada et al. (2007) assuming an orbital plane perpendicular to the outflow axis. From the third Kepler’s law of motion, , with τ_{o} = 400 yr and M_{tot} = 1 M_{⊙} we can derive the mean separation of the companion at a ≃ 50 au. On the other hand, from the semiamplitude of the wiggling (α ≃ 0.5°) the orbital velocity of the CO outflow source is constrained at: V_{0} = V_{CO} × tan(α) = 0.1 km s^{−1}, using an average velocity of ≃10 km s^{−1} for the CO outflow. This in turn gives the orbital radius of the outflow source around the center of mass of the binary r_{0} = 1.3 au and the ratio µ between the mass of the companion and the total mass of the system: µ = r_{0}/a = 0.025. Thus a brown dwarf or massive planetary mass companion located at ≃50 au (=0.35″) separation would be required to account for the observed wiggling of the CO outflow in the orbital scenario. Such a low mass companion could have escaped direct detection so far (Rodríguez et al. 2012). Strikingly the predicted companion separation is very close to an emission bump at r = 62 au detected in the continuum emission profile of the disk at millimetric wavelengths (de Valon et al. 2020; Garufi et al. 2020). However no clear gap is detected in the disk emission at this position which would be expected for such a massive companion.
Therefore, the precession scenario is excluded to account for the observed wiggling of the CO flow while the orbital scenario requires a brown dwarf or massive planetary mass companion at 50 au separation, which signature we do not clearly see in the disk yet. We also recall that our wiggling models cannot account for the observed variable separations between cusps in channel maps. So we conclude that although attractive to explain some of the substructures observed in the DG Tau B CO outflow, the interpretation of the wiggling scenario faces some difficulties. We discuss below the alternative model where substructures arise from axisymmetric brightness enhancements in the disk wind.
6.3.3 A variable disk wind
We show in Sect. 5.3 that the cusps, fingers, and a section of the arches can be explained by brightness or density enhancements in the conical outflow. The timescales between the density enhancements are typically a few hundred years (see Fig. 14). Unfortunately, these timescales cannot be directly compared to the ones observed in the DG Tau B jet due to the larger jet velocity and its fading brightness at large distances. However, such timescales are observed on younger molecular outflows. The cluster W43MM1 in Nony et al. (2020), CARMA7 in Plunkett et al. (2015) as well as HH46/47 in Zhang et al. (2019) show signatures of variability in molecular outflows with timescales between episodic events typically of a few hundred years. We may be witnessing similar variability in the DG Tau B CO outflow.
In the following, we discuss the possibility that these density enhancements are created by variability at the source in the diskwind launching regions. For the model presented in Fig. 15, we considered for the sake of simplicity that the density bursts take place simultaneously at all radii in the disk (that is, τ_{i} = cst for all layers). A more physical assumption would be to assume that the density burst propagates radially across the disk with a velocity V_{prop}. If we consider that the burst takes place close to the midplane (Z ≈ 0), the expression of the travel time for the density launched from a fixed radius r_{0} is: (18) (19) (20)
where t is the current time and t_{eject}(r_{0}) is the epoch of density ejection. Here, V_{prop} is positive when the density burst moves from the inner to the outer disk regions.
Figure 14 indicates that τ decreases with increasing poloidal velocity, which is consistent with a density enhancement moving from the outer to the inner regions of the disk. Figure 20 is a modification of Fig. 14 in which we transformed the poloidal velocity into the corresponding radius of ejection under the MHD diskwind hypothesis, using Eq. (8) from Ferreira et al. (2006) with λ_{ϕ} = 1.58 and β = 0. We did not take into account the uncertainty on the estimation of r_{0} (≃20%, see Fig. 17) as these uncertainties would make impossible the estimation of V_{prop}. Nonetheless, this study gave an estimate of the propagation velocity. For each cusp, we traced the profile τ(r_{0}) derived from the observed positions of the cusp apex in the different channel maps, as described in Sect. 5.3. Linear fits to these profiles with associated τ uncertainties are also shown. The slopes of these profiles is directly linked to 1/V_{prop}, and hence give the radial propagation velocities of the density enhancement at the origin. The average velocity over all cusps is V_{prop} ≈ −0.2 km s^{−1}. The same study could be achieved with PDW models. In that case, the radius of ejection were multiplied by λ^{2} (=2.5) and therefore V_{prop} values were also multiplied by 2.5, giving an average radial velocity V_{prop} ≈ −0.5 km s^{−1}. From uncertainties on the slopes, we derived uncertainties on V_{prop} values in the case of U3, U4, and U5. In the case of U0,1,2, a horizontal solution (V_{prop} = ∞) could not be excluded. For these last three fits, we derived the minimal radial velocities for the two opposite directions of propagation.
Figure 20 seems to indicate that the density bursts propagate from the outer to the inner regions of the disk at ≈V_{prop} = 0.2–0.5 km s^{−1}. The accreting density must also propagate toward the inner regions of the conical flow (r_{0} < 2 au) where the density enhancements were visible but not characterized. This propagation velocity corresponds to ≈0.01–0.04 V_{Kep}(r = 2–5 au). Such accretion velocities match expectations for radial surface velocities in MHD winddriven accretion (Riols et al. 2020). From this velocity and our estimate of the burst duration of ≃23 yr, we can constraint the radial extent of the burst propagating along the disk at ∆R ≈ 1 au. Figure 19, left panel, illustrates the proposed scenario.
Episodic density bursts propagating inward are obviously reminiscent of Fu Ori and Ex Ori type variable accretion events. Some Fu Ori have burst durations <30 yr (e.g., V1515 Cyg, or V1714 Cyg in Hartmann & Kenyon 1996). But periods are usually assumed to be ≈10^{4}–10^{5} yr. On the other extreme, Ex Ori have typically bursts of a few months and periods of few years. DG Tau B variability seems to be located between these two extrema. However, the mass accretion rate increase in Fu Ori type events is typically three to four orders of magnitude higher than suggested in the DG Tau B outflow by the moderate factor three emissivity enhancement during the bursts. Moreover, models of Fu Ori events predict a global accretion affecting the whole vertical structure of the disk during the high state. In magnetically accreting disk models, most of the mass is concentrated in the midplane, where the radial accretion velocity is ≃V_{Kep}/10^{3} (Riols et al. 2020). The high propagation velocity combined with the moderate emissivity enhancements derived in DG Tau B suggest that the accretion burst takes place locally on the disk surface and is less extreme than in typical Fu Ori phenomena. Such moderate accretion bursts could be due for example to residual infalling envelope material creating a shock wave when infalling into the disk (Hennebelle et al. 2017). It is important to note, however, that this axisymmetric model does not reproduce the local deviations observed in the specific angular momentum map (Fig. 3). Nonaxisymmetric perturbations would be required.
Fig. 19 Schematic scenarios describing the different components of the DG Tau B redshifted outflow in the diskwind paradigm (see text). 
Fig. 20 Dynamical timescales of the density enhancements responsible for the highvelocity cusps in the channel maps, as a function of the launching radius r_{0} of the streamline under the cold MHD diskwind hypothesis (assuming λ = 1.58). For each cusp, the red dashed line represents the best linear fit of τ(r_{0}) and its corresponding velocity. The one σ domain is shown in gray and the uncertainties of the fits are noted on the right. 
7 Conclusions
We present a detailed analysis and modeling of the ALMA ^{12}CO(2–1) observations of the DG Tau B redshifted outflow published in de Valon et al. (2020), with the aim to constrain its origin. Our main conclusions are as follows:
We identify three classes of discrete structures visible on the ^{12}CO channel maps: arches at low velocities, fingers at medium velocities, and cusps at high velocities. Both cusps and arches show apparent acceleration in channel maps.
We derived the 2D kinematics of the inner conical outflow using a tomographic method, assuming only axisymmetry of the outflow. We reconstructed 2D maps for both the expansion velocity V_{Z} and specific angular momentum j = R × V_{ϕ}. The inner outflow shows a striking V_{Z} shear with faster material closer to the flow axis. Lines of constant V_{Z} are conical from Z ≃ 50 au out to ≈1200 au. Specific angular momentum is roughly constant along those lines (except in two localized regions), and it increases outward inversely with V_{Z} from ≃40 to 100 au km s^{−1}.
The lower velocity external CO outflow shows a parabolic morphology, apparent acceleration, and a large specific angular momentum of j = 250 ± 50 au km s^{−1}. This suggests that it is tracing either a sweptup infalling and rotating envelope, or an extended disk wind launched from ~30 au.
The conical outflow and the discrete structures could not be described by winddriven shells with radial Hubble velocity laws. Such models fail to reproduce at the same time the morphologies of the observed structures (arches, fingers, and cusps) in the channel maps. Numerical simulations of the interaction between an episodic jet and wideangle wind with an infalling envelope are, however, required to confirm these conclusions.
Instead, the conical outflow global morphology and kinematics appear consistent with matter directly ejected from the disk. Constraints on the diskwind footpoint radii were derived at r_{0} = 1.6–8.2 au (resp. 0.7–3.4 au) in the limiting cases where thermal (resp. magnetocentrifugal) processes dominate. However, none of the current photoevaporated wind models can reproduce the large observed mass flux (2.3 × 10^{−7} M_{⊙} yr^{−1}) ejected from r_{0} < 10 au. In contrast, an MHD diskwind model with a constant magnetic level arm parameter λ ≤ 1.58 can reproduce – at the same time – the flow velocity and angular momentum, as well as the large mass flux if it extracts most of the angular momentum for disk accretion across the wind launching region. The low lambda value is consistent with recent models of warm or weakly magnetized MHD DW (Bai et al. 2016; JacqueminIde et al. 2019).
The wiggling of the flow axis may explain both the localized deviations of specific angular momentum and the morphology of the substructures in the conical flow. Orbital motion of the flow source in a binary of separation ≃50 au with companion mass 2.5 × 10^{−3} M_{⊙} can explain the inferred wiggling period and amplitude. Such a low mass companion could have escaped direct detection so far, but it should produce a gap signature in the continuum dust disk emission, which is not currently detected. In addition, the wiggling scenario fails to account for the variable separation between cusps in channel maps.
Alternatively, the substructures observed in the CO channel maps can be explained by a series of mild (a factor three) density perturbations in the wind launching region, propagating inward at a radial velocity of ≃0.2–0.5 km s^{−1}, consistent with the surface accretion flow predicted in MHD winddriven accretion. We derived a typical perturbation width of ~ 1 au and intervals of 200–500 yr between perturbations. Alternatively the conical morphology and local density enhancements might be explained by the interaction of inner jet bow shocks with the disk wind (Tabone et al. 2018). although, further numerical simulations are required to fully test this hypothesis.
The discrete structures that are increasingly observed in Class 0 and Class I outflows on larger scales have been usually interpreted in terms of nested shells swept up by an episodic inner wind. In contrast, we have shown in this paper that the substructures in the DG Tau B outflow appear best explained by density enhancements at the disk surface and that they propagate in a shearlike MHD disk wind. If confirmed, these results would directly demonstrate the link between accretion and ejection processes in embedded sources.
DG Tau B is a Class I protostar with a structured disk, infalling flows, and possible variable disk wind. Structures in Class I disks are assumed to trace early stages of planetary formation processes. Therefore, our results suggest that planetary formation is taking place in a very dynamic environment. The impact of such outflowing and infalling flows on the disk and its evolution remains an open question. Additional models and simulations of variable disk wind are needed in order to fully comprehend its impact on disk evolution and planet formation.
Acknowledgements
The authors would like to thank the referee, whose comments helped improve the quality of the paper. This paper makes use of the following ALMA data: ADS/JAO.ALMA#2015.1.01108.S, ADS/JAO.ALMA#2017.1.01605.S. ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ. This work was supported by the Programme National de Physique Stellaire (PNPS), the Programme National de Physique et Chimie du Milieu Interstellaire (PCMI) of CNRS/INSU cofunded by CEA and CNES and the Conseil Scientifique of Observatoire de Paris. FL acknowledges the support of the Marie Curie Action of the European Union (project MagiKStar, Grant agreement number 841276). L.A.Z. acknowledges financial support from CONACyT280775 and UNAMPAPIIT IN110618 and IN103921 grants, México. DM acknowledges support from ANID Basal project FB210003.
Appendix A Determination of the redshifted outflow PA
We show in Sect. 3 that the velocity difference between two radially symmetric positions, at a given altitude z in the flow, is directly linked to the rotation velocity of the outflow. However, this derivation is highly sensitive to the determination of the outflow radial center position. Because of the steep gradient in projected velocity versus radius observed in transverse PV diagrams, a small variation in radial position would cause an important variation in the velocity difference. Therefore, the determination of the position axis is a key concern in the characterization of the specific angular momentum of the outflow.
To accurately derive the flow center position as a function of the projected altitude, we adopted the following method, illustrated in Fig. A.1. At each projected altitude above the disk δz, we integrated the emissivity from (V − V_{sys}) = 6.64 km s^{−1} to (V − V_{sys}) = 9.77 km s^{−1}. This gave the radial emissivity profile of the highvelocity emission tracing the limiting inner cone. Indeed as shown in DV20, at large projected velocities, the emission traces an inner cone with almost constant radius. We derived the radial positions of the two peaks tracing the edges of the highvelocity component (See Fig. A.1) using Gaussian fitting. From the median value of these two positions, we derived the radial center position and its associated uncertainty. We applied this method after rotating the data cube with three different values of outflow PA (sampled around the disk PA), and the transverse PV diagram was obtained with a slice perpendicular to this axis. The effect of the different PV cut is completely negligible, as the considered variation of outflow P.A is small (less than 2 degrees).
Figure A.1 shows the derived radial offset as a function of the projected height for three values of the flow position angle (PA). With a PA of 296° and 294°, the radial offset increases consistently with a missestimation of the outflow PA of ± 1°. The global offset is minimized with a PA of 295 ± 1°. This is consistent with the atomic redshifted jet PA of 296° derived by Mundt et al. (1987) and the PA of the projected disk in the plane of the sky, determined at 25.7 ± 0.3° by DV20 and 24 ± 1° by Guilloteau et al. (2011).
The outflow rotation could potentially introduce a bias in this method. Indeed, at a given projected velocity, the two sides of a rotating ring fall at different projected radial offsets from the axis, inducing an artificial shift of the position centroid of the flow axis. However, since a rotating ring produces a tilted ellipse in the transverse PV diagram, centered at V_{Z} × cos(i), the radial shift has an opposite sense at velocities above and below V_{Z} × cos(i). By integrating emission over a broad range of velocity, we thus averaged out this effect. Moreover, if the variation of radial offset was due to flow rotation, the measured radial offsets should decrease with distance from the source, as the rotation velocity decreases (due to the conservation of angular momentum along conical streamlines). This is not consistent with our results in Fig. A.1 where the radial offsets stay constant (at our nominal PA) or increase linearly (at nonnominal PAs) with height. Fig. A.1 also shows no clear signature of wiggling in the highvelocity component of DG Tau B. This is consistent with the absence of wiggling in the atomic jet observed by Mundt et al. (1991). From the maximum error bars observed, we derived an upper limit for the wiggling angle of θ ≤ 0.5°. We consider in the following that the redshifted outflow PA is constant for all layers (velocities) at PA = 295°, corresponding to the PA derived here for the highvelocity emission.
Fig. A.1 Determination of the outflow axis PA. Top panels: Determination of the highvelocity component radial center as a function of projected heights above the disk for different outflow axis PA orientations. The black dotted lines show an angle of ±1° with respect to the central axis. Middle panel: ^{12}CO PV diagram perpendicular to the outflow axis. (V − V_{sys}) unit is km s^{−1}. Bottom panel: The black dashed profile shows the emissivity integrated from the PV diagram between (V − V_{sys}) = 6.64 km s^{−1} to (V − V_{sys}) = 9.77 km s^{−1}, this domain is indicated by the red dashed lines in the middle panel. The red curves show Gaussian fits used to derive the radial positions of the two edge peaks. The vertical black dotted lines in the bottom and middle panels indicate the positions of these peaks. The average of these two radial positions give the radial center at this height. 
Appendix B Characterization of the cusps
Fig. B.1 Characterization of the cusps. Left panels: Individual ^{12}CO channel maps at different lineofsight velocities illustrating the spatial evolution of the cusps. Six cusps, labeled U0U5, are clearly identified (U4 is only visible on three channel maps). We identify the projected heights of the bottom of the cusps (red circles) and their evolution in the channel maps (white dashed lines). Right panel: Longitudinal PV diagram at δx = 0″ averaged over a slice of ∆x = 0.2″. The velocity and height locations of the cusps derived in the channel maps are shown as red dots. The four orange crosses indicate the maximum onaxis projected velocity (back side) of the top of the fit four bright elliptical structures identified in transverse PV diagrams at the corresponding heights (see Fig. C.1). 
Appendix C Ellipse fits in the transverse PV diagrams
In this section, we discuss an alternative method to derive the radial component of the velocity, not constrained by our tomographic method. At a few positions along the flow, transverse PV diagrams clearly show elliptical structures nested inside the main emission (see Fig. C.1). We assumed that each of these ellipses traces one emitting layer of the outflow. If the outflow velocity does not vary drastically with height, one layer of the outflow of fixed radius will be projected as an inclined elliptical structure in the PV diagram (Louvet et al. 2018). The velocity width of the ellipse at x = 0 is directly linked to the radial velocity component of the layer. We fit ellipses by the naked eye to the structures observed. From these fits, we recovered the radial velocity V_{R} as well as V_{Z} at a few specific positions along the outflow.
Figure C.1 shows the V_{Z} profile of the outflow reconstructed with the tomographic technique and extrapolated until Z ≈ 2700 au. The red arrows trace the poloidal vectors derived from individual elliptical fits. The angle of the arrow is determined by the ratio between V_{Z} and V_{R}. The length and colors of the arrows correspond to V_{Z}. The V_{Z} values derived from the elliptical fits appear consistent with the estimates from the tomographic study. The poloidal velocity direction is also consistent with the conical lines of constant V_{Z} in the tomography. This suggests that the flow is indeed aligned with these conical lines.
We studied the relation between the ellipses visible in the transverse PV diagrams and the cusps visible in the channel maps. The onaxis maximal velocity of the fit ellipses was compared to the cusp location in the onaxis longitudinal PV diagram (see Fig. B.1). The (R,Z) positions of the fit ellipses were also compared to the cusp locations in the tomographic V_{Z} map (see Fig. 13). The ellipse located at δz = 1.1″ comes from a lower altitude region, not included in our cusp analysis because of crowding. The ellipses at δz = 2.1″ and 2.9″ seem to extend the U5 cusp, at higher and lower lineofsight velocities respectively. The ellipse at δz = 7.5″ is not so clearly associated with a cusp extension in Fig. B.1. It could be the high velocity extension of the U3 cusp, or associated with a fainter cusp located between U2 and U3, not included in our analysis. Moreover, Fig. C.2 indicates that the tops of the ellipses located at δz = 2.9″ and δz = 7.5″ are possibly consistent with a cusp structure which was not included in our study, due to their low S/N. Similarly, Fig. C.2 shows that the cusps U1 and U2 also coincide with internal structures in the transverse PV diagram, but the stacking of ellipses in the PV diagram makes the identification more confusing than in the channel maps.
The cusps could be precisely located only at moderate projected velocities <5.3 km s^{−1}, where they are sufficiently bright (see Fig. B.1), hence they probe outer streamlines; in contrast elliptical structures in transverse PV diagrams are best distinguished at high projected velocities > 5.3 km s^{−1} where they have less overlap with each other (see Fig. C.1), hence most of them probe inner faster flow streamlines. However, we see clear correspondences between faint cusps in channel maps and the onaxis highvelocity portion of some ellipses in transverse PV cuts, and viceversa, which demonstrates that they trace different portions of the same underlying substructures, extending across the whole conical outflow. Therefore, assuming a radial flow to deproject the cusp apparent positions seems fully justified.
Fig. C.1 Constraints on the radial velocity component V_{R}. Left panels: Transverse PV diagrams averaged over a slice of width ∆Z = 0.2″ at selected positions δz along the flow. In red are shown the ellipse fits. We show side by side the PV diagram with and without the fit for more visibility. Right panel: V_{Z} tomography of the outflow in the outflow referential, extrapolated beyond Z > 1200 au. The red arrows represent the velocity directions determined from the ellipse fits. The arrow colors trace the V_{Z} values determined from the ellipse fit. R and Z units are in au. 
Fig. C.2 Illustration of the correspondences between ellipses and cusps. Left panels: transverse PV diagrams at δz = 2.9″ and 7.5″ (from Fig C.1), with red dots indicating the onaxis projected velocities of cusps U1 and U2 at δz = 7.5″. Each of these two cusps corresponds in position and velocity to the maximum velocity (back side) of an ellipse in the PV cut. Right panels: channel maps at the same lineofsight velocity as the orange crosses in the left panels, showing that elliptical structures seem to correspond to fainter cusps in positionposition space, but falling outside of the velocity domain where cusp characterization was possible. 
Appendix D Biases in the tomographic reconstruction
We studied possible biases in the method presented in Sect. 3 to reconstruct poloidal maps of V_{Z} and j = R × V_{ϕ}. We first studied the bias introduced by projection effects considering one single layer of the outflow using the analytical solutions from Eqs. 1, 2, 3. To study the impact of beam convolution and multiple layers of the outflow, we also applied our tomographic reconstruction method to the synthetic data cube presented in Sect. 10 and estimated the difference between the reconstructed V_{Z} and j and the initial values of the model.
D.1 Single shell
In our tomographic study, we assumed that one shell will be projected as an ellipse in transverse PV diagrams. We also assumed that the projected velocities at the extrema radii allowed us to recover V_{Z} and V_{ϕ}. We determine in this section biases introduced by these two assumptions.
We assumed a conical shell of radius R(Z) in the outflow referential and with local opening angle θ, such as R(Z) = R_{0} + Z tan θ. For the DG Tau B inner conical outflow, the opening angle θ of the layers vary between 17° and 12°, increasing with decreasing velocities. The estimated opening angle of the lower velocity emission contributing to the pedestal is θ ≃ 30°. We defined the parameter θ_{v} corresponding to the angle of the velocity vector with the Z axis in the poloidal plane (V_{R} = V_{Z} tan θ_{V}). For a flow parallel to the conical surface, θ = θ_{V}. We also assumed a constant specific angular momentum (V_{ϕ} = j/R(ϕ)) as well as constant V_{Z} and V_{R} over the transverse slit width. Solving Eq. 2 with constant δz gave: (D.1)
In the case i = θ and i < θ, the cut of the conical outflow will be respectively a parabola and an hyperbola. Consequently, R(ϕ) will tend toward infinite. This is not the case for the DG Tau B conical outflow. We then implemented this solution into Eqs. 3,1 for multiple values of the inclination. The resulting transverse PV diagrams are represented in Fig. D.1 for θ_{V} = θ = 17°. The PV diagrams were generated using Eqs. 3, 2 with the solution from Eq. D.1. The difference with an ellipse increases with the inclination but is expected to be small for the DG Tau B case (i = 117°).
We assume in Sect. 3 that the radial edges of the ellipse correspond to ϕ = 0 − π. However, this is an approximation. The radial extrema of the ellipse correspond to solutions of the equation: . For a conical layer, the solution is obtained for: (D.2)
Fig. D.1 Biases in the tomographic reconstruction for one single conical shell. Top panels: Black curves show computed transverse PV diagrams for one conical layer with constant V_{Z}, V_{R}, and j seen at different inclinations i and with θ_{V} = θ = 17°. The red dashed curve traces an elliptical fit. Bottom Panel: Computed relative bias in V_{Z} (green curves) and j (orange curves) for one conical layer as a function of the inclination and for different combinations of θ and θ_{v}: solid curves: θ = θ_{v} = 17°, dotted curves: θ = 17° & θ_{v} = 30°, and dashed curves: θ = θ_{v} = 30°. 
This difference is small in our situation. However, this leads to bias in the estimation of R, Z, V_{Z} and V_{ϕ}. We determined the estimated value of V_{est,Z} and V_{est,ϕ} using Eqs. 6,7. We then computed the relative differences with the real values R_{real}, Z_{real}, V_{real,Z}, and V_{real,ϕ}: (D.3) (D.4) (D.5) (D.6)
Hence, R, Z, and j are systematically underestimated, while V_{Z} is overestimated. Figure D.1 shows the relative biases due to inclination and projection effects in the estimation of V_{Z} and V_{ϕ} for different values of θ and θ_{v}. In the conical outflow, where θ_{V} = θ ≤ 17°, the inclination bias is expected to be ≤10% in V_{Z} and Z, ≤ 1 % in j and R.
D.2 Multiple layers
In the previous section, we estimated the bias in the estimation of V_{Z}, V_{ϕ} due to projection effects for one single conical layer. However, the DG Tau B transverse PV diagrams shows a clear shearlike velocity structure suggesting a stacking of layers. We modeled this effect directly in the transverse PV diagrams by stacking the elliptical projections for conical layers of increasing radii at origin R_{0} with the same opening angle θ. The velocity shear in V_{Z} of the conical layers was defined by a shearing parameter . The specific angular momentum of each layer, j, was assumed constant with Z and vary between 20 and 90 au km s^{−1}, increasing with increasing radius, such as j × V_{P} was kept constant to mimic the DG Tau B observations.
Due to the stacking, the maximal velocity at ±R does not perfectly describe the velocity of the ellipse of radius R. This effect is larger when f_{sh} is small. Using a similar procedure as before, we derived the relative difference between the estimation of V_{Z} and j with the tomographic method and the input theoretical values. This was achieved for a range of f_{sh} values between 0.2 to 20 km s^{−1} au^{−1} and with the three (θ, θ_{V}) configurations studied in the previous subsection.
In order to study efficiently this bias, we defined the adimensional parameter , where we defined [V_{los}(±R)] = 1/2(V_{los}(R) + V_{los}(−R)). This adimensional parameter can be derived from the observations. We show in Fig. D.2 the predicted relative biases in V_{Z} and j as a function of this parameter x. Biases increase with increasing x values illustrating the effect of the velocity shear. We also show the distribution of observed x values computed at each (R,Z) position in the conical outflow. The x values are concentrated around ≃0.5, suggesting moderate biases in both V_{Z} and j for the conical outflow where θ_{V} = θ ≤ 17°.
However, this modeling did not include the effect of beam smearing nor the impact of our polynomial fitting method to describe the shape of the PV diagram. In order to study these effects, we directly applied our tomographic method to the synthetic data cube presented in Sect. 10 and determined the relative differences between the computed and the input V_{Z} and j values at each (R,Z) position along the conical flow. We show in Fig. D.2 the derived relative biases in V_{Z} and j as a function of the x parameter. Their distributions are broader and flatter than predicted, especially in j, likely due to the combination of velocity shear and beam smearing effects.
Fig. D.2 Biases in the tomographic reconstruction for a stacking of conical shells. Left panels: The gray area shows the limits of the transverse PV diagram predicted for a stacking of conical shells with a shearing factor f_{sh} (see text). In black is shown the ellipse corresponding to a single layer of maximal radius R. The blue and red dots correspond respectively to the velocity of the ellipse at ±R and the maximal velocity of the PV diagram at ± R. The red dots are the ones used in our method, while the blue dots correspond to the true measurements. The differences in velocities will lead to a bias in our j and V_{Z} estimation. The red line shows the slope of the PV diagram , see text for more details. Right panel: Plot of the relative bias in V_{Z} (green curves) and j (orange curves) due to the shearing aspect of the PV diagram as a function of the x parameter (see text). The solid, dashed, and dotted lines correspond to different configurations of θ and θ_{v}. The black histogram, labeled C.O., shows the distribution of the x parameter derived at different heights Z and radii R in the conical outflow. The relative biases corresponding to the a and b ellipses shown in the left panels are represented. The orange (resp. green) colored contours trace the relative biases derived from applying the tomographic method to the synthetic data cube. Contours outline 30 to 90 % of the distributions in step of 15 %. 
From this study, we estimate that in the conical part of the outflow (at (V − V_{sys}) ≥ 2 km s^{−1}), the tomographic method suffers from a relative bias ≤15% in the estimation of V_{Z} and ≤20% in the estimation of j. The X values extend up to ≈1.3 in the pedestal region. With an opening angle θ, θ_{v} > 30^{°}, a tomographic study of the low velocity component would suffer from a relative bias larger than 60% and 30% in the estimation of V_{Z} and j respectively.
Appendix E Winddriven shell analytical solutions
In this section, using Eqs. 1, 2,3 and 11, we derive an analytical solution for predicted channel maps in the case of the generalized WDS model introduced in Sect. 4. The WDS model is defined by three parameters: C, τ and n, see Eqs. 11. The projection on the plane of the sky (δx, δz) for the emissivity map at (V − V_{sys}) = V_{cm} could be recovered by solving Eq. (3) with V_{los} = V_{CM} and using Eqs. 11 for V_{Z} and V_{R}: (E.1) (E.2) (E.2) (E.4)
Equation E.2 has no solution in the case ζ > 1 (). This corresponds to the case where no emission is predicted at V_{los} = V_{CM}. Similarly, in the situation 0 ≥ ζ ≥ 1, only a fraction of ϕ will be projected such as sin^{2}ϕ > ζ. Developing Eq. E.4, we obtained the following: (E.5) (E.6)
We reformulated δz as z_{0} + ∆Z sin β with: (E.7) (E.8) (E.9)
Similarly, manipulating Eq. E.3, δx could be reformulated as δx = ∆X cosβ with the following: (E.10)
Hence (δx, δz) trace an ellipse of center (0, z_{0}) and aspect ratio: (E.11)
Therefore, in the classical WDS models with radial velocity vectors (η = 1), the aspect ratio of the ellipse on the channel maps only depends on the inclination. We represent on the channel maps in Figs. 5, 7, and 8 the two limiting ellipses computed at V − δV/2 and V + δV/2 to take into account the width of the channel map.
Appendix F Global model
In order to comprehend the impact of projection or convolution effects on the observations, we developed a locally axisymmetric code that allows us to simulate optically thin observations of simple outflow models. We took advantage of the axisymmetric hypothesis that permits to reduce the complexity of the model by defining at each height the radius R(Z) and the cylindrical velocities V_{z}, V_{r} and V_{ϕ} (see the schematic view of Fig. 2). The dependency between the height and the radius or the velocities vary with the model used. We then created at each height an emitting ring with azimuth parametrized with ϕ. As the morphology is axisymmetric, neither the radius nor velocities depend on ϕ. Under the assumption of optically thin emission, the emissivity at position (Z,ϕ) is proportional to the elementary volume dV = R(Z)dϕdRdZ. We added an additional variation of the emissivity with the height and radius as a powerlaw with parameters α and β respectively. Proper modeling of the emissivity would require the temperature and chemistry to be solved, which is well beyond the scope of this model. The positions of the outflow emission on the data cube (δx,δy, V_{proj}) were then defined by Eqs. 1,2,3.
Fig. F.1 Schematic representation of the method used to create synthetic data cube. 
We then created a data cube with the same spectral and spatial resolution than our observations, and placed on each point (δx,δy, V_{proj}) the emissivity I. Under the assumption of optically thin emission, we summed each emissivity corresponding to the same positions on the data cube. We set a step size of 1° for ϕ and a fraction of the spatial pixel for Z. We then convolved the data cube by a 2D Gaussian matching the spatial beam characteristics in order to fully simulate the ALMA observations. The code, written in Python 3, is publicly online^{5}.
Appendix G Channel maps of diskwind models
Fig. G.1 ^{12}CO observed (first column) and synthetic diskwind channel maps at selected lineofsight velocities. The second and third columns correspond to the two diskwind models with axis precession presented in Sect. 5.2. The last column corresponds to a diskwind model with axisymmetric density enhancements, as presented in Sect. 5.3. δx and δz units are arcseconds. 
References
 Alexander, R., Pascucci, I., Andrews, S., Armitage, P., & Cieza, L. 2014, in Protostars and Planets VI, eds. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning, 475 [Google Scholar]
 Anderson, J. M., Li, Z.Y., Krasnopolsky, R., & Blandford, R. D. 2003, ApJ, 590, L107 [Google Scholar]
 Anglada, G., López, R., Estalella, R., et al. 2007, AJ, 133, 2799 [Google Scholar]
 Arce, H. G., Shepherd, D., Gueth, F., et al. 2007, in Protostars and Planets V, eds. B. Reipurth, D. Jewitt, & K. Keil, 245 [Google Scholar]
 Bai, X.N., Ye, J., Goodman, J., & Yuan, F. 2016, ApJ, 818, 152 [Google Scholar]
 Bjerkeli, P., van der Wiel, M. H. D., Harsono, D., Ramsey, J. P., & Jørgensen, J. K. 2016, Nature, 540, 406 [Google Scholar]
 Blandford, R. D., & Payne, D. G. 1982, MNRAS, 199, 883 [CrossRef] [Google Scholar]
 Cabrit, S., Raga, A., & Gueth, F. 1997, in HerbigHaro Flows and the Birth of Stars, 182, eds. B. Reipurth, & C. Bertout, 163 [NASA ADS] [CrossRef] [Google Scholar]
 Casse, F., & Ferreira, J. 2000, A&A, 353, 1115 [NASA ADS] [Google Scholar]
 Clarke, C. J., & Alexander, R. D. 2016, MNRAS, 460, 3044 [NASA ADS] [CrossRef] [Google Scholar]
 Cunningham, A., Frank, A., & Hartmann, L. 2005, ApJ, 631, 1010 [NASA ADS] [CrossRef] [Google Scholar]
 de Valon, A., Dougados, C., Cabrit, S., et al. 2020, A&A, 634, A12 [Google Scholar]
 Downes, T. P., & Cabrit, S. 2007, A&A, 471, 873 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Eislöffel, J., & Mundt, R. 1998, AJ, 115, 1554 [CrossRef] [Google Scholar]
 Ellerbroek, L. E., Podio, L., Kaper, L., et al. 2013, A&A, 551, A5 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Facchini, S., Juhász, A., & Lodato, G. 2018, MNRAS, 473, 4459 [Google Scholar]
 FernándezLópez, M., Zapata, L. A., Rodríguez, L. F., et al. 2020, AJ, 159, 171 [Google Scholar]
 Ferreira, J., Dougados, C., & Cabrit, S. 2006, A&A, 453, 785 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Font, A. S., McCarthy, I. G., Johnstone, D., & Ballantyne, D. R. 2004, ApJ, 607, 890 [NASA ADS] [CrossRef] [Google Scholar]
 Garufi, A., Podio, L., Codella, C., et al. 2020, A&A, 636, A65 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Gaudel, M., Maury, A. J., Belloche, A., et al. 2020, A&A, 637, A92 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Guilloteau, S., Dutrey, A., Piétu, V., & Boehler, Y. 2011, A&A, 529, A105 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hartmann, L., & Kenyon, S. J. 1996, ARA&A, 34, 207 [NASA ADS] [CrossRef] [Google Scholar]
 Hennebelle, P., Lesur, G., & Fromang, S. 2017, A&A, 599, A86 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hirota, T., Machida, M. N., Matsushita, Y., et al. 2017, Nat. Astron., 1, 0146 [Google Scholar]
 JacqueminIde, J., Ferreira, J., & Lesur, G. 2019, MNRAS, 490, 3112 [Google Scholar]
 Launhardt, R., Pavlyuchenkov, Y., Gueth, F., et al. 2009, A&A, 494, 147 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lee, C.F., Mundy, L. G., Reipurth, B., Ostriker, E. C., & Stone, J. M. 2000, ApJ, 542, 925 [NASA ADS] [CrossRef] [Google Scholar]
 Lee, C.F., Stone, J. M., Ostriker, E. C., & Mundy, L. G. 2001, ApJ, 557, 429 [NASA ADS] [CrossRef] [Google Scholar]
 Lee, C.F., Li, Z.Y., Hirano, N., et al. 2018, ApJ, 863, 94 [NASA ADS] [CrossRef] [Google Scholar]
 Lee, C.F., Tabone, B., Cabrit, S., et al. 2021, ApJ, 907, L41 [NASA ADS] [CrossRef] [Google Scholar]
 Liang, L., Johnstone, D., Cabrit, S., & Kristensen, L. E. 2020, ApJ, 900, 15 [NASA ADS] [CrossRef] [Google Scholar]
 LópezVázquez, J. A., Cantó, J., & Lizano, S. 2019, ApJ, 879, 42 [CrossRef] [Google Scholar]
 Louvet, F., Dougados, C., Cabrit, S., et al. 2018, A&A, 618, A120 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Masciadri, E., & Raga, A. C. 2002, ApJ, 568, 733 [Google Scholar]
 Mitchell, G. F., Sargent, A. I., & Mannings, V. 1997, ApJ, 483, L127 [NASA ADS] [CrossRef] [Google Scholar]
 Mundt, R., & Fried, J. W. 1983, ApJ, 274, L83 [NASA ADS] [CrossRef] [Google Scholar]
 Mundt, R., Brugel, E. W., & Buehrke, T. 1987, ApJ, 319, 275 [NASA ADS] [CrossRef] [Google Scholar]
 Mundt, R., Ray, T. P., & Raga, A. C. 1991, A&A, 252, 740 [NASA ADS] [Google Scholar]
 Nixon, C., King, A., & Price, D. 2013, MNRAS, 434, 1946 [NASA ADS] [CrossRef] [Google Scholar]
 Nony, T., Motte, F., Louvet, F., et al. 2020, A&A, 636, A38 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ohashi, N., Hayashi, M., Ho, P. T. P., et al. 1997, ApJ, 488, 317 [NASA ADS] [CrossRef] [Google Scholar]
 Owen, J. E., Ercolano, B., & Clarke, C. J. 2011, MNRAS, 412, 13 [Google Scholar]
 Panoglou, D., Cabrit, S., Pineau Des Forêts, G., et al. 2012, A&A, 538, A2 [CrossRef] [EDP Sciences] [Google Scholar]
 Pety, J., Gueth, F., Guilloteau, S., & Dutrey, A. 2006, A&A, 458, 841 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Picogna, G., Ercolano, B., Owen, J. E., & Weber, M. L. 2019, MNRAS, 487, 691 [Google Scholar]
 Plunkett, A. L., Arce, H. G., Mardones, D., et al. 2015, Nature, 527, 70 [NASA ADS] [CrossRef] [Google Scholar]
 Podio, L., Eislöffel, J., Melnikov, S., Hodapp, K.W., & Bacciotti, F. 2011, A&A, 527, A13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pudritz, R. E., Ouyed, R., Fendt, C., & Brandenburg, A. 2007, in Protostars and Planets V, eds. B. Reipurth, D. Jewitt, & K. Keil, 277 [Google Scholar]
 Riols, A., Lesur, G., & Menard, F. 2020, A&A, 639, A95 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rodríguez, L. F., Dzib, S. A., Loinard, L., et al. 2012, Rev. Mexicana Astron. Astrofis., 48, 243 [Google Scholar]
 Shang, H., Shu, F. H., & Glassgold, A. E. 1998, ApJ, 493, L91 [NASA ADS] [CrossRef] [Google Scholar]
 Shang, H., Allen, A., Li, Z.Y., et al. 2006, ApJ, 649, 845 [NASA ADS] [CrossRef] [Google Scholar]
 Shang, H., Krasnopolsky, R., Liu, C.F., & Wang, L.Y. 2020, ApJ, 905, 116 [NASA ADS] [CrossRef] [Google Scholar]
 Shu, F. H., Ruden, S. P., Lada, C. J., & Lizano, S. 1991, ApJ, 370, L31 [NASA ADS] [CrossRef] [Google Scholar]
 Tabone, B., Cabrit, S., Bianchi, E., et al. 2017, A&A, 607, L6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Tabone, B., Raga, A., Cabrit, S., & Pineau des Forêts, G. 2018, A&A, 614, A119 [CrossRef] [EDP Sciences] [Google Scholar]
 Tabone, B., Cabrit, S., Pineau des Forêts, G., et al. 2020, A&A, 640, A82 [CrossRef] [EDP Sciences] [Google Scholar]
 Terquem, C., Papaloizou, J. C. B., Nelson, R. P., & Lin, D. N. C. 1998, ApJ, 502, 788 [Google Scholar]
 Ulrich, R. K. 1976, ApJ, 210, 377 [Google Scholar]
 Wang, L., & Goodman, J. 2017, ApJ, 847, 11 [Google Scholar]
 Wang, L.Y., Shang, H., Krasnopolsky, R., & Chiang, T.Y. 2015, ApJ, 815, 39 [NASA ADS] [CrossRef] [Google Scholar]
 Wang, L., Bai, X.N., & Goodman, J. 2019, ApJ, 874, 90 [Google Scholar]
 Wilgenbus, D., Cabrit, S., Pineau des Forêts, G., & Flower, D. R. 2000, A&A, 356, 1010 [NASA ADS] [Google Scholar]
 Zapata, L. A., Lizano, S., Rodríguez, L. F., et al. 2015, ApJ, 798, 131 [NASA ADS] [CrossRef] [Google Scholar]
 Zhang, Y., Higuchi, A. E., Sakai, N., et al. 2018, ApJ, 864, 76 [NASA ADS] [CrossRef] [Google Scholar]
 Zhang, Y., Arce, H. G., Mardones, D., et al. 2019, ApJ, 883, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Zhu, Z. 2019, MNRAS, 483, 4221 [Google Scholar]
It is obtained when a wideangle wind with velocity varying with angle as V_{w} ∝ cos θ and density varying as ∝ 1/(r^{2} sin^{2} θ) sweepsup a static, flattened isothermal core with density ∝ sin^{2} θ/r^{2}, and they mix instantly in the shell. We note that the radial shell expansion results from instant mixing, while the Hubble law derives from the identical radial falloff of wind and ambient density (both ∝1/r^{2}), which yields a shell speed that is constant over time (Shu et al. 1991). Finally, the parabolic shell shape derives from the combined θdependencies of the densities and wind speed (Lee et al. 2001).
We note that tan 14° ~ 1/4 hence V_{Z}(14°) ≃ 16/(τC), cf. Eq. (9).
A static medium ensures that, after full mixing, the shell expands in the same radial direction as the wind, while the 1/r^{2} ambient density decrease ensures that the expansion speed is constant over time (the density ratio between the wind and ambient medium being independent of radius); both properties then together yield the “Hubble law” V = r/τ (Shu et al. 1991; Lee et al. 2001).
Ram pressure equilibrium between the reverse shock in the wind and the forward shock in the static ambient medium imposes , where V_{w} and V_{s} are the wind and shell speeds, and ρ_{w} and ρ_{a} are the wind and ambient density. The massflux entering the shell from the wind side will then dominate over the sweptup mass if (V_{w}  V_{s}) < V_{s}.
All Tables
All Figures
Fig. 1 Substructures in the DG Tau B redshifted CO outflow. Left panels: ^{12}CO channel maps at selected lineofsight velocities. The white dashed line traces the θ = 17° outer limiting cone of the inner outflow as defined by DV20. Right panels: transverse PV diagrams across the flow averaged over a slice of ∆Z = 0.2″. Black dotted lines indicate the outer limits of the conical outflow at the specified height. The triangle (a) and the two circle symbols (c and d) highlight the onaxis height of one arch and two different cusps, respectively. The square symbol (b) is located in the extended outer flow. The symbols are represented both on the channel maps and PV diagrams. 

In the text 
Fig. 2 Principle of tomographic reconstruction method. Top panel: 3D representation of the two reference systems used: the outflow (X, Y, Z) in black and the observer (δx, δy, δz) in yellow. The colored circles illustrate the emitting rings defined by five parameters (see text). The colored dots trace the locations at ϕ = 0 and ϕ = π along the emitting rings. Bottom panels: transverse PV diagrams at δz = 2.5″ across the flow axis and averaged over a slice of ∆z = 0.2″. In the left panel is shown the schematic projection of the colored rings and corresponding colored dots in the PV diagram. Because of the flow inclination, the rings are in fact projected at slightly different heights. Their real projection is studied in Appendix D.1. The white dots in the right panel illustrate the outer limits of the PV diagram. The red curve shows the polynomial fitting of the two edges. 

In the text 
Fig. 3 Tomographic maps of V_{Z} (left panel) and specific angular momentum j (right panel) in the outflow referential. The black dashed line traces the conical fit of the region V_{Z} = 10–11 km s^{−1}. The white (resp. black) contours in the left (resp. right) panel show V_{Z} contours. The red dashed lines indicate the heights of the two extrema in specific angular momentum. 

In the text 
Fig. 4 Specific angular momentum j along curves of constant V_{Z}. The corresponding range of VZ is shown in the box. The uncertainty of specific angular momentum was obtained by propagating the V_{proj} uncertainty. The red dashed line corresponds to the median value of each curve. 

In the text 
Fig. 5 Comparison of the lowvelocity outer CO outflow with a classical WDS model. Left panel: ^{12}CO individual channel maps at different lineofsight velocities tracing the lowvelocity outer CO outflow. The white contours trace the model of a WDS of parabolic shape defined by C = 10^{−3} au^{−1}, dynamical age τ = 6000 yr and specific angular momentum j = 250 ± 50 au km s^{−1} with an inclination of i = 117°. Right panel: transverse PV cut at δz = 7.5″ averaged over a slice of ∆Z = 0.2″. The red ellipse traces the WDS model. (V − V_{sys}) units are km s^{−1}. 

In the text 
Fig. 6 Dynamical times (τ = Z/V_{Z}) derived from the tomographic map of V_{Z} are shown in color. The white lines correspond to a parabolic curve with 1/C varying from 30 to 135 au, in steps of 15 au. 

In the text 
Fig. 7 Comparison of ^{12}CO channel maps at three different lineofsight velocities (color maps) with predicted ellipses for parabolic WDS (white contours). The model used on each row is sketched in the leftmost panel. Top row: classical model with radial Hubble expansion (see Eq. (8) and Lee et al. 2000). Middle and bottom rows: modified model with “collimated” expansion (see Eq. (11)). Green, blue, and yellow contours highlight specific height ranges, indicated in the first column. White dotted lines outline the radial boundary defined by the shell ellipses with increasing velocity. Model parameters are listed in Table 2. δx and δz units are arcseconds. 

In the text 
Fig. 8 Kinematic evolution along the flow axis. Panels a: positionvelocity diagrams averaged over a slice of width ∆Z = 0.2″ at three different δz positions along the flow. The red contours trace solutions of collimated WDS used to fit A2, A3, and U5 (third line on Fig. 7). Panel b: in black is represented the maximal velocity of emission of the PV diagram. In red is shown the center velocity of the ellipses on the PV diagrams. Panel c: the gray region highlights the conical outflow domain. The black lines represent the parabolic morphology of the three solutions. The red dots show for each solution the region where the projected velocity reach V_{max}. 

In the text 
Fig. 9 Diskwind properties derived from the tomography along lines of constant V_{Z}. Red dashed curves show fits as a function of the poloidal velocity of the streamline V_{p}: the anchoring radius of the streamline r_{0,con} was fit by a thirdorder polynomial, the angle of the streamline with the flow axis θ by a power law and the product V_{P} × j by a constant value of 570 au km^{2} s^{−2}. Blue dotted lines represent the extrapolation used to model the highvelocity component that could not be mapped by tomography due to its low S/N. 

In the text 
Fig. 10 Comparison of observations with steady diskwind model. Left panels: individual ^{12}CO channel maps computed from the global diskwind model at selected lineofsight velocities (top row) are compared to observations (bottom row). The color scale is the same for all the channel maps. Right panels: transverse PV diagrams at two positions δz along the flow and averaged over a slice of width ∆Z = 0.2″. The background grayscale image shows the observations and the red contours trace the predictions from the diskwind model. (V − V_{sys}) units are km s^{−1}. 

In the text 
Fig. 11 Synthetic data cubes computed from the generic diskwind model with flow axis precession: with a precession angle of 0.2° (left panels) and 0.5° (right panels) and assuming a constant precession period τ_{p} = 400 yr (top row) or a constant precession spatial wavelength Λ_{p} = 800 au (bottom row). For each model are shown both a channel map at (V − V_{sys}) = 4.37 km s^{−1} and the corresponding tomography of j derived with the method described in Sect. 3. The color scale of j ranging from 0 to 140 au km s^{−1} is identical to the color scale of Fig. 3. 

In the text 
Fig. 12 Lineofsight velocity onaxis (at δx = 0) of the model diskwind flow surfaces as a function of their poloidal velocity. Orange circles show the front side (), and green circles show the back side () of the diskwind flow surfaces. In gray, we show the lineofsight velocity range where the arches and cusps are observed (see Appendix B). In light gray, is represented the low S/N domain where the cusps could not be characterized. 

In the text 
Fig. 13 Tomographic map of V_{Z} from Fig. 3 extrapolated from Z = 1200 au to Z = 3000 au using the conical wind model from Fig. 9. Gray hatched regions named SU05 represent the locus of emissivity enhancements required to explain the cusps onaxis positions at V − V_{sys} = 3.73–5.32 km s^{−1} (see Fig. B.1). White dots represent the solutions for the cusps at V − V_{sys} = 4.37 km s^{−1}. Red dotted lines indicate the height of the two extrema in specific angular momentum (from Fig. 3). The gray points located on the R = 0 axis represent the positions of jet knots observed by Podio et al. (2011). The red dots indicate the positions with associated errors of the ellipses identified in transverse PV diagrams (see Appendix C) and indicate that emissivity enhancements extend to the inner streamlines. 

In the text 
Fig. 14 Dynamical evolution of the enhancements. Left panel: dynamical age of the enhancements SU05 as a function of V_{P}. The red dotted line represents the mean value. Right panel: dynamical age versus timescale between adjacent cusps (τ_{i} − τ_{i+1}). The uncertainties are determined by propagating the errors on the cusp positions. 

In the text 
Fig. 15 ^{12}CO channel maps at selected lineofsight velocities for the steady conical diskwind model (bottom row) and observations (top row). White contours highlight the predicted emission from the modeled density enhancement SU0. The intensity color scale is the same for all channel maps. 

In the text 
Fig. 16 Specific angular momentum map (in color) for the model of ballistic rotating infall (Ulrich 1976) with R_{d} = 700 au and M_{*} = 1.1 M_{⊙}. The angular momentum varies from 25 au km s^{−1} to 900 au km s^{−1}. The green hatched and red filled regions represent respectively the limits of the conical outflow derived from the tomography in Sect. 3 and the parabolic shape of the classical Hubblelaw WDS model fit to the external CO outflow. The white and red contours outline the infalling streamlines with a specific angular momentum similar to the conical outflow j = 40–100 au km s^{−1} (white), and to the WDS model j = 200–300 au km s^{−1} (red). The black hatched region outline the streamlines with initial θ_{0} = 70 ± 5° reproducing the infall signatures seen in DG Tau B (DV20). 

In the text 
Fig. 17 Launching radius of the streamline as a function of its poloidal velocity assuming: conical extrapolation for the geometrical radius (green circles), steady thermal ejection conserving angular momentum (orange circles) and steady cold magnetocentrifugal ejection (blue circles). See text for more details. We do not derive the launching radius for the largest poloidal velocities due to the large uncertainties. 

In the text 
Fig. 18 Specific angular momentum ( j = R × V_{ϕ}) versus poloidal velocity V_{P} for steady and axisymmetric MHD disk winds. The black symbols represent values derived from the tomography and averaged along lines of constant V_{Z}, with their uncertainty. Red curves show the expected relation from selfsimilar cold magnetocentrifugal disk winds with r_{0} varying from 0.5 to 100 au and λ between 1.52 and 2.3 from Ferreira et al. (2006). The green dot corresponds to a warm MHD disk wind solution from Casse & Ferreira (2000) with r_{0} = 1 au and λ = 1.9. The green arrow shows the path followed by the solution with an increase of r_{0}. The gray box represents the estimated specific angular momentum for the outer low velocity layer on Sect. 4.1, taking a maximal poloidal velocity corresponding to a height of Z = 3300 au in the outflow referential. 

In the text 
Fig. 19 Schematic scenarios describing the different components of the DG Tau B redshifted outflow in the diskwind paradigm (see text). 

In the text 
Fig. 20 Dynamical timescales of the density enhancements responsible for the highvelocity cusps in the channel maps, as a function of the launching radius r_{0} of the streamline under the cold MHD diskwind hypothesis (assuming λ = 1.58). For each cusp, the red dashed line represents the best linear fit of τ(r_{0}) and its corresponding velocity. The one σ domain is shown in gray and the uncertainties of the fits are noted on the right. 

In the text 
Fig. A.1 Determination of the outflow axis PA. Top panels: Determination of the highvelocity component radial center as a function of projected heights above the disk for different outflow axis PA orientations. The black dotted lines show an angle of ±1° with respect to the central axis. Middle panel: ^{12}CO PV diagram perpendicular to the outflow axis. (V − V_{sys}) unit is km s^{−1}. Bottom panel: The black dashed profile shows the emissivity integrated from the PV diagram between (V − V_{sys}) = 6.64 km s^{−1} to (V − V_{sys}) = 9.77 km s^{−1}, this domain is indicated by the red dashed lines in the middle panel. The red curves show Gaussian fits used to derive the radial positions of the two edge peaks. The vertical black dotted lines in the bottom and middle panels indicate the positions of these peaks. The average of these two radial positions give the radial center at this height. 

In the text 
Fig. B.1 Characterization of the cusps. Left panels: Individual ^{12}CO channel maps at different lineofsight velocities illustrating the spatial evolution of the cusps. Six cusps, labeled U0U5, are clearly identified (U4 is only visible on three channel maps). We identify the projected heights of the bottom of the cusps (red circles) and their evolution in the channel maps (white dashed lines). Right panel: Longitudinal PV diagram at δx = 0″ averaged over a slice of ∆x = 0.2″. The velocity and height locations of the cusps derived in the channel maps are shown as red dots. The four orange crosses indicate the maximum onaxis projected velocity (back side) of the top of the fit four bright elliptical structures identified in transverse PV diagrams at the corresponding heights (see Fig. C.1). 

In the text 
Fig. C.1 Constraints on the radial velocity component V_{R}. Left panels: Transverse PV diagrams averaged over a slice of width ∆Z = 0.2″ at selected positions δz along the flow. In red are shown the ellipse fits. We show side by side the PV diagram with and without the fit for more visibility. Right panel: V_{Z} tomography of the outflow in the outflow referential, extrapolated beyond Z > 1200 au. The red arrows represent the velocity directions determined from the ellipse fits. The arrow colors trace the V_{Z} values determined from the ellipse fit. R and Z units are in au. 

In the text 
Fig. C.2 Illustration of the correspondences between ellipses and cusps. Left panels: transverse PV diagrams at δz = 2.9″ and 7.5″ (from Fig C.1), with red dots indicating the onaxis projected velocities of cusps U1 and U2 at δz = 7.5″. Each of these two cusps corresponds in position and velocity to the maximum velocity (back side) of an ellipse in the PV cut. Right panels: channel maps at the same lineofsight velocity as the orange crosses in the left panels, showing that elliptical structures seem to correspond to fainter cusps in positionposition space, but falling outside of the velocity domain where cusp characterization was possible. 

In the text 
Fig. D.1 Biases in the tomographic reconstruction for one single conical shell. Top panels: Black curves show computed transverse PV diagrams for one conical layer with constant V_{Z}, V_{R}, and j seen at different inclinations i and with θ_{V} = θ = 17°. The red dashed curve traces an elliptical fit. Bottom Panel: Computed relative bias in V_{Z} (green curves) and j (orange curves) for one conical layer as a function of the inclination and for different combinations of θ and θ_{v}: solid curves: θ = θ_{v} = 17°, dotted curves: θ = 17° & θ_{v} = 30°, and dashed curves: θ = θ_{v} = 30°. 

In the text 
Fig. D.2 Biases in the tomographic reconstruction for a stacking of conical shells. Left panels: The gray area shows the limits of the transverse PV diagram predicted for a stacking of conical shells with a shearing factor f_{sh} (see text). In black is shown the ellipse corresponding to a single layer of maximal radius R. The blue and red dots correspond respectively to the velocity of the ellipse at ±R and the maximal velocity of the PV diagram at ± R. The red dots are the ones used in our method, while the blue dots correspond to the true measurements. The differences in velocities will lead to a bias in our j and V_{Z} estimation. The red line shows the slope of the PV diagram , see text for more details. Right panel: Plot of the relative bias in V_{Z} (green curves) and j (orange curves) due to the shearing aspect of the PV diagram as a function of the x parameter (see text). The solid, dashed, and dotted lines correspond to different configurations of θ and θ_{v}. The black histogram, labeled C.O., shows the distribution of the x parameter derived at different heights Z and radii R in the conical outflow. The relative biases corresponding to the a and b ellipses shown in the left panels are represented. The orange (resp. green) colored contours trace the relative biases derived from applying the tomographic method to the synthetic data cube. Contours outline 30 to 90 % of the distributions in step of 15 %. 

In the text 
Fig. F.1 Schematic representation of the method used to create synthetic data cube. 

In the text 
Fig. G.1 ^{12}CO observed (first column) and synthetic diskwind channel maps at selected lineofsight velocities. The second and third columns correspond to the two diskwind models with axis precession presented in Sect. 5.2. The last column corresponds to a diskwind model with axisymmetric density enhancements, as presented in Sect. 5.3. δx and δz units are arcseconds. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.