© A. de Valon et al. 2022

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ández-López et al. 2020). Despite their ubiquity, the exact origin of molecular outflows, their link to the high-velocity 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 swept-up material, tracing the interaction between an inner jet or a wide-angle 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₀ ≃ 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 wind-driven 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 high-angular 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 high-angular resolution observations also reveal striking signatures of multiple CO shells in a few sources (Zhang et al. 2019; Fernández-López et al. 2020). Under the classical paradigm where CO outflows trace swept-up 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 red-shifted 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 line-of-sight 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 swept-up 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 model-independent 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 wind-driven 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: 12CO channel maps at selected line-of-sight 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 on-axis 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⁻¹. Its opening angle decreases from 17° at (V − V_sys) = 2.8 km s⁻¹ until 12° at (V − V_sys) ≥ 5.0 km s⁻¹. The sheer-like velocity gradient across the conical layer is best seen in transverse position-velocity (hereafter PV) cuts (Fig. 1), where the flow width clearly narrows down at higher velocity, up (V − V_sys) > 5.0 km s⁻¹. The conical outflow is surrounded by a slower and wider outflowing component visible at (V − V_sys) < 2 km s⁻¹. 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:

Bow-shaped intensity enhancements are visible at low velocities (V − V_sys) = 0.88–1.51 km s⁻¹ (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⁻¹, 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⁻¹ 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⁻¹. We refer to these substructures as fingers.

At high velocities, from (V − V_sys) = 3.1 to almost 7 km s⁻¹, multiple U-shaped 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⁻¹ 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⁻¹ 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⁻¹ (see Fig. B.1). At higher velocities, the cusps could not be characterized because the outflow signal-to-noise 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 on-axis on the channel map at (V − V_sys) = 4.37 km s⁻¹. The apparent acceleration of each cusp, in (″)/km s⁻¹ 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 PV-diagrams as pseudo-ellipses (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 model-dependent study of these ellipses in Sect. 5.

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 model-independent 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 line-of-sight 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 edge-on 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 position-positionvelocity (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: $$\delta x=R\cos \varphi$$(1) $$\delta z=Z\sin i-R\sin \varphi \cos i$$(2) $$V_{\text{los}}=-V_{\text{z}}\cos i-V_{\varphi }\cos \varphi \sin i-V_{\text{R}}\sin \varphi \sin i.$$(3)

A transverse PV diagram corresponds to a pseudo-slit 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 edge-on 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: $$Z=\frac{\delta z}{\sin i}$$(4) $$R=\delta x_{\max }$$(5) $$V_Z\simeq \frac{V_{\text{los}}\left(\delta x_{\max }\right)+V_{\text{los}}\left(-\delta x_{\max }\right)}{-2\cos i}$$(6) $$V_{\varphi }\simeq \frac{V_{\text{los}}\left(\delta x_{\max }\right)-V_{\text{los}}\left(-\delta x_{\max }\right)}{2\sin i}.$$(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⁻¹. We used in the following a mean value of 0.1 km s⁻¹. 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⁻¹. 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 VZ (left panel) and specific angular momentum j (right panel) in the outflow referential. The black dashed line traces the conical fit of the region VZ = 10–11 km s−1. The white (resp. black) contours in the left (resp. right) panel show VZ 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 semi-opening angles varying from 12° for the highest velocities to 17° for the lowest velocities. V_Z radially decreases from ≈ 14 km s⁻¹ to 5 km s⁻¹. 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⁻¹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⁻¹ in the inner radius to ≈ 100 au km s⁻¹ 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⁻¹ is consistent with the previous estimate of DV20 of ≈65 au km s⁻¹.

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⁻¹. 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⁻¹ 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 VZ. The corresponding range of VZ is shown in the box. The uncertainty of specific angular momentum was obtained by propagating the Vproj 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 model-dependent 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 low-amplitude 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⁻¹ 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⁻¹, absorbed by the medium. In that case, our reconstruction method is not applicable. By consequence, we did not apply our method for line-of-sight velocities (V − V_sys) < 2 km s⁻¹. In the following, we investigate to which extent wind-driven 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 Wind-driven 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 wide-angle 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 wind-driven 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¹ (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 Hubble-law kinematics of the shell can be empirically described by two parameters, C and τ, through: $$\begin{array}{ccc}Z=C\times R^2& V_Z=\frac{Z}{\tau }& V_{\text{R}}=\frac{R}{\tau },\end{array}$$(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: $$V_Z\left(\theta \right)=\frac{1}{\tau C}{\left(\frac{Z}{R}\right)}^2={\left(\tau C{\tan }^2\theta \right)}^{-1}$$(9) $$V_{\text{R}}\left(\theta \right)=\frac{1}{\tau C}\left(\frac{Z}{R}\right)={\left(\tau C\tan \theta \right)}^{-1}.$$(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 line-of-sight 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 low-velocity outer flow component, and then the bright conical outflow and its discrete structures.

Fig. 5 Comparison of the low-velocity outer CO outflow with a classical WDS model. Left panel: 12CO individual channel maps at different line-of-sight velocities tracing the low-velocity 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 − Vsys) units are km s−1.

4.1 Low-velocity outer flow

The wide and low-velocity outflow at (V − V_sys) < 2 km s⁻¹ 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 line-of-sight velocities (see Fig. 1 left panels). Although such WDS models do not usually consider rotation, we included rotation to properly fit the large left-right 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 swept-up 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 left-right asymmetry.

Figure 5 (left panels) shows that the outer contour of the low-velocity outflow and its increased spatial offset with velocity are well fit by a WDS obeying Eq. (8) with parameters C = 10⁻³ au⁻¹, τ = 6000 yr and j = 250 ± 50 au km s⁻¹. 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 blue-shifted 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 wind-driven 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 quasi-parabolic curves above Z ≃ 400 au (see Fig. 6). Hence we investigate below whether the conical outflow could be made of successive nested parabolic wind-driven 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 left-right 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 large-scale 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 blue-shifted lobe (Eislöffel & Mundt 1998), hence i ≤ 115° for the redshifted lobe. This limit agrees within 2^° with the large-scale 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 aspect-ratio 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 on-axis (δx = 0) where the projected velocity is greatly affected by the radial velocity component V_R (see Eq. (3) with $\varphi =\pm \frac{\pi }{2}$). To keep a small number of model parameters, we thus chose to add an ad hoc free parameter η that modified the radial velocity as: $$\begin{array}{ccc}Z=C{R}^2& V_Z=\frac{Z}{\tau }& V_{\text{R}}=\eta \frac{R}{\tau }.\end{array}$$(11)

In Appendix E, we show that the ellipse aspect ratio in this modified model is set at A = (η tan² 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 (forward-directed) 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⁻¹ 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 best-fit 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 HH46-47. The fit velocity values V_Z at θ = 14° (in the region of the conical flow) are also listed². Not surprisingly, their range of ≃7.5–16 km s⁻¹ is similar to our “model-independent” tomographic results for V_Z in the conical flow region.

The Hubble law in the WDS model predicts an ever-increasing shell speed at higher altitudes, until it reaches the polar wind speed, which is ≃125 km s⁻¹ according to the red-shifted jet speed in DG Tau B (Eislöffel & Mundt 1998). In contrast, the maximum line-of-sight 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⁻¹ (see Fig. 8b). Therefore, all successive wind-driven 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⁻¹. This velocity limit is close to the molecule dissociation limit ≃20 km s⁻¹ 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 shock-dissociation of ambient CO. Figure 8c shows that the corresponding truncation region for the best-fit 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 low-velocities (arches), ellipse flanks at mid-velocities (fingers), and ellipse bottom at high-velocities (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, WDS-A0, 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 Hubble-law dynamics.

Another serious issue is that the best-fitting 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, low-velocity 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 Hubble-law dynamics. We discuss in Sect. 6.1 the implications of these results in the context of more general wind or jet-driven shell scenarios.

Fig. 6 Dynamical times (τ = Z/VZ) derived from the tomographic map of VZ 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 12CO channel maps at three different line-of-sight velocities (color maps) with predicted ellipses for parabolic WDS (white contours). The model used on each row is sketched in the left-most 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 Disk-wind 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: position-velocity 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 Vmax.

5.1 Steady disk-wind 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 Z-axis 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² s⁻² for all the streamlines. The fits were achieved using nonlinear least squares, the equations are as follows: $$r_0\simeq 0.18-5.33V_{\text{p}}+44.29V_{\text{p}}^2-73.23V_{\text{p}}^3$$(12) $$\theta \simeq 12.63+348.57V_{\text{p}}^{-2.55}$$(13) $$j\times V_{\text{p}}\simeq 570.$$(14)

Here r₀ is in au, V_P in km s⁻¹, θ in degrees, and j in au km s⁻¹. 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⁻¹, 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⁻¹. The parameters in the velocity range V_p = 14–20 km s⁻¹, 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 almost-vertical high-velocity 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² s⁻² (see Fig. 9). We assumed optically thin emission throughout the outflow, which is justified by the observed ratios of ¹³CO/¹²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⁻¹ as well as its variation with line-of-sight 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⁻¹ (see Sect. 2, Fig. 10). To describe completely this extended low-velocity component, we would need to extrapolate the disk-wind model at larger ejection radii. However, as this component falls at absorbed cloud velocities, we are not able to derive model-independent constraints on the dynamics of this component. Proper modeling would require time-consuming 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 disk-wind 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 disk-wind models which could explain the various substructures observed.

Fig. 9 Disk-wind properties derived from the tomography along lines of constant VZ. Red dashed curves show fits as a function of the poloidal velocity of the streamline Vp: the anchoring radius of the streamline r0,con was fit by a third-order polynomial, the angle of the streamline with the flow axis θ by a power law and the product VP × j by a constant value of 570 au km2 s−2. Blue dotted lines represent the extrapolation used to model the high-velocity component that could not be mapped by tomography due to its low S/N.

Fig. 10 Comparison of observations with steady disk-wind model. Left panels: individual 12CO channel maps computed from the global disk-wind model at selected line-of-sight 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 disk-wind model. (V − Vsys) 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 disk-wind 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 ${\tau }_{\text{p}}\propto V_{\text{p}}^{-1}\propto r_0^{0.5}$. 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 high-velocity 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 $r_0^{0.5}$ 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 disk-wind 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 − Vsys) = 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 on-axis at δx = 0 for the front side and the back side ($\varphi =\pm \frac{\pi }{2}$) of each conical wind layer in our model, predicted from Eq. (3). We also represent the domain of line-of-sight 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⁻¹ 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 ($\varphi \approx -\frac{\pi }{2}$) and the arches from the front side ($\varphi \approx +\frac{\pi }{2}$) of the enhanced ring. Figure 12 also shows that wind layers with poloidal velocities V_p > 9 km s⁻¹ 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 line-of-sight velocity V_los. The projected height of the cusp corresponds to an equation Z(R) derived from Eq. (2), with $\varphi =-\frac{\pi }{2}$. 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⁻¹. 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 line-of-sight 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 disk-wind 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 Hubble-law dynamics in the underlying flow.

Interestingly enough, the density enhancements SU3-4 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⁻¹, due to the low S/N of the cusps at the corresponding projected velocities (V − V_sys) > 5.32 km s⁻¹ (see Fig. 12).

We computed the dynamical age of these emissivity enhancements with ${\tau }_{\text{w}}=\frac{Z}{V_Z}$. 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 on-axis (at δx = 0), corresponding to $\varphi =-\frac{\pi }{2}$. 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 disk-wind 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⁻¹), since we did not have constraints on their location at higher velocities. Figure 15 shows computed and observed channel maps at different line-of-sight 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⁻¹ 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 disk-wind 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⁻¹, we do not have model-independent 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 Line-of-sight velocity on-axis (at δx = 0) of the model disk-wind flow surfaces as a function of their poloidal velocity. Orange circles show the front side ($\varphi =\frac{\pi }{2}$), and green circles show the back side ($\varphi =-\frac{\pi }{2}$) of the disk-wind flow surfaces. In gray, we show the line-of-sight 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 VZ from Fig. 3 extrapolated from Z = 1200 au to Z = 3000 au using the conical wind model from Fig. 9. Gray hatched regions named SU0-5 represent the locus of emissivity enhancements required to explain the cusps on-axis positions at V − Vsys = 3.73–5.32 km s−1 (see Fig. B.1). White dots represent the solutions for the cusps at V − Vsys = 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 $\tau =\frac{Z}{V_Z}$ of the enhancements SU0-5 as a function of VP. 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 12CO channel maps at selected line-of-sight velocities for the steady conical disk-wind 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 small-scale redshifted CO outflow in DG Tau B: 1) a stacking of multiple shells swept-up 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 wind-driven shells

We discuss here a scenario where the redshifted CO outflow in DG Tau B can be accounted for by a stacking of multiple swept-up shells resulting from the interaction of an episodic wide-angle 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 low-velocity 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⁻¹, in the same sense as the disk (see Sect. 4.1).

Studies of rotation signatures in pre-stellar and proto-stellar 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⁻¹ (0.2–2 × 10⁻³ km s⁻¹ pc, Gaudel et al. 2020). Our estimate for the outer flow in DG Tau B, j_outer ≃ 250 ± 50 au km s⁻¹ 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² radial density decrease³. A rotating infalling envelope, in contrast, has a non-radial motion and a flatter radial density law ∝ 1/r^1.5 (Ulrich 1976). However, the calculations of López-Vázquez et al. (2019) for such an ambient medium show that the WDS expansion remains quasi radial, except close to the mid-plane, 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ópez-Vá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, free-falling 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 mixing-layer. The static shell shape and mixing-layer properties were recently computed in the case of a free-falling 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 mixing-layer is predicted to be $j\approx 0.15\sqrt{G{M}_{*}{R}_{\text{d}}}\simeq 120$ au km s⁻¹, 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 Rd = 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 Hubble-law 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 Hubble-law velocity field, which requires an expanding shell (Shu et al. 1991). (2) The shell models fitting the arches at low-velocity and cusps at high velocity do not agree with the observed emission morphology in channel maps at mid-velocities, 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 Hubble-law velocity field, regardless of the shell detailed morphology. Therefore, any model where the shell is expanding quasi-radially and at nearly constant speed over time will fail to reproduce our observations. This, in particular, discards all models where the wide-angle wind and the ambient medium share a similar power-law in r, and where they instantly mix in the shell (e.g., the models of López-Vá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 swept-up 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⁻¹ (see Fig. 3), raises additional issues for the swept-up shell scenario. We first note that a wide-angle “X-wind” 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⁻¹, a factor two to three times less than observed. In addition, the wind cannot dominate the swept-up shell mass unless it is slower than twice the shell speed⁴, which in the present case would require V_w ≤ 15 km s⁻¹ (see Table 2). This is inconsistent with a wind originating from close to the protostar. The “X-wind” model, for example, has V_w ≃ 150 km s⁻¹ (Shang et al. 1998).

The low observed expansion speeds ≃6–14 km s⁻¹ 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 swept-up. 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 swept-up 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 swept-up material originates from the rotating disk atmosphere, at radii R₀ ≃ j²/(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⁻¹ (see last column of Table 2 for WDS-A2, WDS-A3 and WDS-U5). 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 large-scale 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 swept-up shells with enough angular momentum and appropriate kinematics could be generated. However, at this stage and taking into account all the above-mentioned difficulties, we do not favor this scenario as the origin of the small-scale 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 disk-wind 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 quasi-steady 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 photo-evaporated disk winds (PDW; Alexander et al. 2014), cold magneto-centrifugal ejection (Pudritz et al. 2007) or a combination of the two processes in the so-called warm or magneto-thermal 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 Photo-evaporated disk winds

In photo-evaporated disk winds, the high energy radiation (UV and X-rays) of the central accreting protostar heats the surface layers of the disk to high temperatures (10³–10⁴ K) up to significant radial distances. Beyond the gravitational radius $r_{\text{g}}=\left(G{M}_{*}\right)\text{/}{c}_{\text{s}}^2$, 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 high-energy irradiation. Extreme-UV (EUV) heating creates an isothermal ionized layer on the disk surface with temperature T ≃ 10⁴ K (c_s = 10 km s⁻¹) which drives a fast wind (~35 km s⁻¹) with mass-loss rates 10⁻⁹–10⁻¹⁰ M_⊙ yr⁻¹ (Font et al. 2004; Wang & Goodman 2017). X-ray irradiation results in cooler and slower flows (c_s ≃ 3–5 km s⁻¹, v ≃ 15–20 km s⁻¹) but penetrates at higher densities and therefore can drive mass-loss rates up to a few 10⁻⁸ M_⊙ yr⁻¹ (Picogna et al. 2019). Non ionizing far-UV (FUV) heating mostly drives slow mass-loss (v ≃ 18 km s⁻¹) from the outer disk regions.

The conical morphology of the DG Tau B CO outflow matches expectations from both self-similar PDW models by Clarke & Alexander (2016) and hydrodynamical simulations by Owen et al. (2011) and Wang & Goodman (2017). In such disk-wind solutions, the specific angular momentum is conserved along the streamlines and is equal to the Keplerian value at the foot-point radius. The observed values of j in the CO conical outflow imply footpoint radii in the range r₀ = 1.6 to 8.2 au. Figure 17 shows that these foot-point 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 self-similar models of Clarke & Alexander (2016).

The values of r₀ are in line with expectations from EUV dominated PDW models. Indeed, for a 1M_⊙ star and c_s ≃ 10 km s⁻¹, 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⁻¹ require sound speeds ≤2–8 km s⁻¹ at radii r₀ ≤ 10 au, excluding EUV-driven models.

The large derived mass flux of 1.7–2.9 × 10⁻⁷ M_⊙ yr⁻¹ for the CO conical flow (DV20) within r₀ ≤ 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 X-ray driven models. Indeed, recent X-ray driven photo-evaporation models of Picogna et al. (2019) predict mass loss rates up to 10⁻⁷ M_⊙ yr⁻¹ for stellar X-ray luminosities L_X ≥ 10³¹ erg s⁻¹. 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 X-ray photons. Figure 9 in Picogna et al. (2019) shows that only 10% of the total mass flux originates from disk foot-point radii below 10 au. In summary, current PDW models fail to account for the combination of large mass flux and small foot-point radii of r₀ = 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 thermo-chemical 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 magneto-centrifugal 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 large-scale 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₀)², where r_A is the poloidal Alfven radius and r₀ 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 magneto-thermal winds have low to moderate λ values but can extract significant mass and angular momentum from the disk.

Terminal velocities depend on both the foot-point radii r₀, λ 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): $$R{V}_{\varphi }=\lambda R_0{V}_{\varphi }\left(R_0\right)$$(15) $$V_{\text{p}}=V_{\varphi }\left(R_0\right)\sqrt{2\lambda -3+\beta },$$(16)

where β encompasses all pressure effects, including thermal and turbulent Alfvén waves (see Ferreira et al. 2006). The streamline foot-point 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₀,λ) predicted by cold magneto-centrifugal disk-wind models (Ferreira et al. 2006). The mean poloidal velocities and specific angular momentum coincide with a line of constant λ_cold ≃ 1.58 with foot-point 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₀ 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₀ = 1 au. We see that using the cold MHD curves would lead to overestimate λ ≈ 2. 5 and underestimate r₀ ≈ 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 (Jacquemin-Ide et al. 2019).

The derived minimum foot-point radius of r_in ≥ 0.7 au for the CO wind streamlines is in good agreement with the thermo-chemical predictions of Panoglou et al. (2012), who demonstrate that CO molecules magnetically launched from foot-point radii ≥1 au survive in the case of accretion rates in the disk ≥10⁻⁷ M_⊙ yr⁻¹. Similar results are obtained for warm magneto-thermal wind solutions (Wang et al. 2019).

The streamline foot-point 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 B-field. 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⁻⁸–10⁻⁷ M_⊙ yr⁻¹ are predicted by the magneto-thermal 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 ${\dot{M}}_{\text{acc}}\left(r\right)\propto r^{\xi }$ (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 (${\dot{M}}_{\text{acc}}\left(r_{\text{in}}\right)={\dot{M}}_{\text{acc}}\left(r_{\text{out}}\right)-{\dot{M}}_{\text{DW}}$) and the definition of ξ, we derived the following expression: $$\frac{{\dot{M}}_{\text{DW}}}{{\dot{M}}_{\text{acc}}\left(r_{\text{in}}\right)}={\left(\frac{r_{\text{out}}}{r_{\text{in}}}\right)}^{\xi }-1.$$(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⁻⁹ M_⊙ yr⁻¹, giving a mass accretion rate onto the star of ≈6 × 10⁻⁸ M_⊙ yr⁻¹. We took this value as a lower limit to ${\dot{M}}_{\text{acc}}\left(r_{\text{in}}\right)$. From the measured mass flux in the conical CO flow Ṁ_DW = 2.3 × 10⁻⁷ M_⊙ yr⁻¹ and the disk-wind 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 X-wind. However, the outer outflow could be tracing the interaction between the envelope and outer disk-wind streamlines located outside of the conical outflow (r₀ > 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 disk-wind 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₀, 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₀ > 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 VP for steady and axisymmetric MHD disk winds. The black symbols represent values derived from the tomography and averaged along lines of constant VZ, with their uncertainty. Red curves show the expected relation from self-similar cold magneto-centrifugal disk winds with r0 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 r0 = 1 au and λ = 1.9. The green arrow shows the path followed by the solution with an increase of r0. 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 Disk-wind 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 over-densities). Indeed, in the jet-wind interaction scenario, the over-densities 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⁻¹, 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 low-amplitude 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 disk-wind 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 over-densities at the positions of individual jet bow shocks, illustrated by the red regions in Fig. 10 (right panels) in Tabone et al. (2018). These over-densities globally reproduce the observed shapes of the perturbations in the DG Tau B conical layer: for the perturbations closer-in, 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 over-densities 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 semi-amplitude 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 mis-aligned 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 mis-aligned 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 ${\tau }_o\text{/}{\tau }_{\text{p}}\simeq 0.37\mu \text{/}\sqrt{1-\mu }$, where µ = m₂/(m₁ + m₂) 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 semi-amplitude wiggling angle α is dominated by the precession motion, which translates into the condition V₀ ≤ V_Z tan(α), that is V₀ ≤ 0.1 km s⁻¹ where V₀ is the orbital velocity of the flow source. Combining these two constraints, we derive a companion mass ratio µ ≤ 2.5 × 10⁻³ and binary separation: a ≤ 0.5 au. This separation is smaller than the launching radius of the CO disk wind (r₀ ≥ 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 mis-aligned binary system. In such circumstances the inner rim of the circumbinary disk can break from the outer disk and precess. However, large mis-alignments 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.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⁻² 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, $a^3=M_{\text{tot}}{\tau }_o^2$, 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 semi-amplitude of the wiggling (α ≃ 0.5°) the orbital velocity of the CO outflow source is constrained at: V₀ = V_CO × tan(α) = 0.1 km s⁻¹, using an average velocity of ≃10 km s⁻¹ for the CO outflow. This in turn gives the orbital radius of the outflow source around the center of mass of the binary r₀ = 1.3 au and the ratio µ between the mass of the companion and the total mass of the system: µ = r₀/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 W43-MM1 in Nony et al. (2020), CARMA-7 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 disk-wind 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 mid-plane (Z ≈ 0), the expression of the travel time for the density launched from a fixed radius r₀ is: $$\tau \left(r_0\right)=t-t_{\text{eject}}\left(0\right)$$(18) $$t_{\text{eject}}\left(r_0\right)=t_{\text{eject}}\left(0\right)+\frac{r_0}{V_{\text{prop}}}$$(19) $$\tau \left(r_0\right)=\tau \left(0\right)-\frac{r_0}{V_{\text{prop}}},$$(20)

where t is the current time and t_eject(r₀) 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 disk-wind 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₀ (≃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₀) 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⁻¹. The same study could be achieved with PDW models. In that case, the radius of ejection were multiplied by λ² (=2.5) and therefore V_prop values were also multiplied by 2.5, giving an average radial velocity V_prop ≈ −0.5 km s⁻¹. 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⁻¹. The accreting density must also propagate toward the inner regions of the conical flow (r₀ < 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 wind-driven 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⁴–10⁵ 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 mid-plane, where the radial accretion velocity is ≃V_Kep/10³ (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 disk-wind paradigm (see text).

Fig. 20 Dynamical timescales of the density enhancements responsible for the high-velocity cusps in the channel maps, as a function of the launching radius r0 of the streamline under the cold MHD disk-wind hypothesis (assuming λ = 1.58). For each cusp, the red dashed line represents the best linear fit of τ(r0) 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 ¹²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 ¹²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⁻¹.

• 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⁻¹. This suggests that it is tracing either a swept-up 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 wind-driven 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 wide-angle 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 disk-wind foot-point radii were derived at r₀ = 1.6–8.2 au (resp. 0.7–3.4 au) in the limiting cases where thermal (resp. magneto-centrifugal) processes dominate. However, none of the current photo-evaporated wind models can reproduce the large observed mass flux (2.3 × 10⁻⁷ M_⊙ yr⁻¹) ejected from r₀ < 10 au. In contrast, an MHD disk-wind 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; Jacquemin-Ide 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⁻³ 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⁻¹, consistent with the surface accretion flow predicted in MHD wind-driven 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 shear-like 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 co-funded 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 CONACyT-280775 and UNAM-PAPIIT 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⁻¹ to (V − V_sys) = 9.77 km s⁻¹. This gave the radial emissivity profile of the high-velocity 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 high-velocity 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 miss-estimation 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 non-nominal PAs) with height. Fig. A.1 also shows no clear signature of wiggling in the high-velocity 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 high-velocity emission.

Fig. A.1 Determination of the outflow axis PA. Top panels: Determination of the high-velocity 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: 12CO PV diagram perpendicular to the outflow axis. (V − Vsys) unit is km s−1. Bottom panel: The black dashed profile shows the emissivity integrated from the PV diagram between (V − Vsys) = 6.64 km s−1 to (V − Vsys) = 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 12CO channel maps at different line-of-sight velocities illustrating the spatial evolution of the cusps. Six cusps, labeled U0-U5, 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 on-axis 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 on-axis maximal velocity of the fit ellipses was compared to the cusp location in the on-axis 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 line-of-sight 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⁻¹, 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⁻¹ 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 on-axis high-velocity portion of some ellipses in transverse PV cuts, and vice-versa, 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 VR. 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: VZ 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 VZ 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 on-axis 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 line-of-sight velocity as the orange crosses in the left panels, showing that elliptical structures seem to correspond to fainter cusps in position-position 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₀ + 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: $$R\left(\varphi \right)=\frac{\delta z\tan \theta +R_0\sin i}{\sin i-\sin \varphi \cos i\tan \theta }.$$(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: $\frac{\partial \delta x}{\partial \varphi }=0$. For a conical layer, the solution is obtained for: $$\sin \varphi =\frac{\tan \theta }{\tan i}.$$(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 VZ, VR, 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 VZ (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,ϕ: $$\frac{R_{\text{est}}-R_{\text{real}}}{R_{\text{real}}}=\sqrt{1-\frac{{\tan }^2\theta }{{\tan }^{2}i}}-1$$(D.3) $$\frac{Z_{\text{est}}-Z_{\text{real}}}{Z_{\text{real}}}=-\frac{{\tan }^2\theta }{{\tan }^{2}i}$$(D.4) $$\frac{V_{\text{est,Z}}-V_{\text{real,Z}}}{V_{\text{real,Z}}}=\tan \theta \times \tan {\theta }_V$$(D.5) $$\frac{V_{\text{est,}\varphi }-V_{\text{real,}\varphi }}{V_{\text{real,}\varphi }}=\sqrt{1-\frac{{\tan }^2\theta }{{\tan }^{2}i}}-1.$$(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 shear-like 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₀ with the same opening angle θ. The velocity shear in V_Z of the conical layers was defined by a shearing parameter $f_{\text{sh}}=\frac{\partial V_Z}{\partial R}$. The specific angular momentum of each layer, j, was assumed constant with Z and vary between 20 and 90 au km s⁻¹, 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⁻¹ au⁻¹ and with the three (θ, θ_V) configurations studied in the previous subsection.

In order to study efficiently this bias, we defined the a-dimensional parameter $x=\frac{\left[V_{\text{los}}\left(\pm R\right)\right]}{R}\times {\left(\frac{\partial \left[V_{\text{los}}\left(\pm R\right)\right]}{\partial R}\right)}^{-1}$, where we defined [V_los(±R)] = 1/2(V_los(R) + V_los(−R)). This a-dimensional 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 fsh (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 VZ estimation. The red line shows the slope of the PV diagram $\frac{\partial V_{\text{los}}\left(\pm R\right)}{\partial R}$, see text for more details. Right panel: Plot of the relative bias in VZ (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⁻¹), 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 Wind-driven 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: $$\zeta =\frac{4V_{\text{CM}}\tau C\cos i}{{\eta }^2{\sin }^{2}i}$$(E.1) $$R=-\eta \frac{\tan i}{2C}\left(\sin \varphi \pm \sqrt{{\sin }^2\varphi -\zeta }\right)$$(E.2) $$\delta x=R\cos \varphi$$(E.2) $$\delta z=C{R}^2\sin i-R\sin \varphi \cos i.$$(E.4)

Equation E.2 has no solution in the case ζ > 1 ($V_{\text{CM}}\cos i>\frac{{\eta }^2{\sin }^{2}i}{4\tau C}$). 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²ϕ > ζ. Developing Eq. E.4, we obtained the following: $$\delta z=\eta \frac{{\sin }^2\varphi \sin i}{2C}\left(\eta {\tan }^{2}i+1\right)-V_{\text{CM}}\tau \tan i$$(E.5) $$\pm \eta \frac{\sin \varphi \sin i}{2C}\left(\eta {\tan }^{2}i+1\right)\sqrt{{\sin }^2\varphi -\zeta }.$$(E.6)

We reformulated δz as z₀ + ∆Z sin β with: $$z_0=\eta \frac{\sin i}{2C}\left(\eta {\tan }^{2}i+1\right)-V_{\text{CM}}\tau \tan i$$(E.7) $$\text{Δ}Z=\eta \frac{\sin i}{2C}\left(\eta {\tan }^{2}i+1\right)\sqrt{1-\zeta }$$(E.8) $$\sin \beta =\frac{-1+{\sin }^2\varphi \pm \sin \varphi \sqrt{{\sin }^2\varphi -\zeta }}{\sqrt{1-\zeta }}.$$(E.9)

Similarly, manipulating Eq. E.3, δx could be reformulated as δx = ∆X cosβ with the following: $$\text{Δ}X=\mp \eta \frac{\tan i}{2C}\sqrt{1-\zeta }.$$(E.10)

Hence (δx, δz) trace an ellipse of center (0, z₀) and aspect ratio: $$\left|\frac{\text{Δ}Z}{\text{Δ}X}\right|=\left(\eta {\tan }^{2}i+1\right)\left|\cos i\right|.$$(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 power-law 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⁵.

Appendix G Channel maps of disk-wind models

Fig. G.1 12CO observed (first column) and synthetic disk-wind channel maps at selected line-of-sight velocities. The second and third columns correspond to the two disk-wind models with axis precession presented in Sect. 5.2. The last column corresponds to a disk-wind model with axisymmetric density enhancements, as presented in Sect. 5.3. δx and δz units are arcseconds.

References

All Tables

All Figures

Fig. 1 Substructures in the DG Tau B redshifted CO outflow. Left panels: 12CO channel maps at selected line-of-sight 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 on-axis 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 VZ (left panel) and specific angular momentum j (right panel) in the outflow referential. The black dashed line traces the conical fit of the region VZ = 10–11 km s−1. The white (resp. black) contours in the left (resp. right) panel show VZ 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 VZ. The corresponding range of VZ is shown in the box. The uncertainty of specific angular momentum was obtained by propagating the Vproj uncertainty. The red dashed line corresponds to the median value of each curve.
In the text

Fig. 5 Comparison of the low-velocity outer CO outflow with a classical WDS model. Left panel: 12CO individual channel maps at different line-of-sight velocities tracing the low-velocity 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 − Vsys) units are km s−1.

In the text

	Fig. 6 Dynamical times (τ = Z/VZ) derived from the tomographic map of VZ 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 12CO channel maps at three different line-of-sight velocities (color maps) with predicted ellipses for parabolic WDS (white contours). The model used on each row is sketched in the left-most 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: position-velocity 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 Vmax.

In the text

Fig. 9 Disk-wind properties derived from the tomography along lines of constant VZ. Red dashed curves show fits as a function of the poloidal velocity of the streamline Vp: the anchoring radius of the streamline r0,con was fit by a third-order polynomial, the angle of the streamline with the flow axis θ by a power law and the product VP × j by a constant value of 570 au km2 s−2. Blue dotted lines represent the extrapolation used to model the high-velocity component that could not be mapped by tomography due to its low S/N.

In the text

Fig. 10 Comparison of observations with steady disk-wind model. Left panels: individual 12CO channel maps computed from the global disk-wind model at selected line-of-sight 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 disk-wind model. (V − Vsys) units are km s−1.

In the text

Fig. 11 Synthetic data cubes computed from the generic disk-wind 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 − Vsys) = 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 Line-of-sight velocity on-axis (at δx = 0) of the model disk-wind flow surfaces as a function of their poloidal velocity. Orange circles show the front side ($\varphi =\frac{\pi }{2}$), and green circles show the back side ($\varphi =-\frac{\pi }{2}$) of the disk-wind flow surfaces. In gray, we show the line-of-sight 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 VZ from Fig. 3 extrapolated from Z = 1200 au to Z = 3000 au using the conical wind model from Fig. 9. Gray hatched regions named SU0-5 represent the locus of emissivity enhancements required to explain the cusps on-axis positions at V − Vsys = 3.73–5.32 km s−1 (see Fig. B.1). White dots represent the solutions for the cusps at V − Vsys = 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 $\tau =\frac{Z}{V_Z}$ of the enhancements SU0-5 as a function of VP. 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 12CO channel maps at selected line-of-sight velocities for the steady conical disk-wind 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 Rd = 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 Hubble-law 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 magneto-centrifugal 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 VP for steady and axisymmetric MHD disk winds. The black symbols represent values derived from the tomography and averaged along lines of constant VZ, with their uncertainty. Red curves show the expected relation from self-similar cold magneto-centrifugal disk winds with r0 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 r0 = 1 au and λ = 1.9. The green arrow shows the path followed by the solution with an increase of r0. 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 disk-wind paradigm (see text).
In the text

Fig. 20 Dynamical timescales of the density enhancements responsible for the high-velocity cusps in the channel maps, as a function of the launching radius r0 of the streamline under the cold MHD disk-wind hypothesis (assuming λ = 1.58). For each cusp, the red dashed line represents the best linear fit of τ(r0) 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 high-velocity 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: 12CO PV diagram perpendicular to the outflow axis. (V − Vsys) unit is km s−1. Bottom panel: The black dashed profile shows the emissivity integrated from the PV diagram between (V − Vsys) = 6.64 km s−1 to (V − Vsys) = 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 12CO channel maps at different line-of-sight velocities illustrating the spatial evolution of the cusps. Six cusps, labeled U0-U5, 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 on-axis 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 VR. 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: VZ 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 VZ 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 on-axis 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 line-of-sight velocity as the orange crosses in the left panels, showing that elliptical structures seem to correspond to fainter cusps in position-position 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 VZ, VR, 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 VZ (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 fsh (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 VZ estimation. The red line shows the slope of the PV diagram $\frac{\partial V_{\text{los}}\left(\pm R\right)}{\partial R}$, see text for more details. Right panel: Plot of the relative bias in VZ (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 12CO observed (first column) and synthetic disk-wind channel maps at selected line-of-sight velocities. The second and third columns correspond to the two disk-wind models with axis precession presented in Sect. 5.2. The last column corresponds to a disk-wind model with axisymmetric density enhancements, as presented in Sect. 5.3. δx and δz units are arcseconds.
In the text