Issue 
A&A
Volume 614, June 2018



Article Number  A119  
Number of page(s)  13  
Section  Interstellar and circumstellar matter  
DOI  https://doi.org/10.1051/00046361/201732031  
Published online  28 June 2018 
Interaction between a pulsating jet and a surrounding disk wind
A hydrodynamical perspective
^{1}
LERMA, Observatoire de Paris, PSL Research University, CNRS Sorbonne Universités UPMC Univ. Paris 06, 75014 Paris, France
email: benoit.tabone@obspm.fr
^{2}
Instituto de Ciencias Nucleares, Universidad Nacional Autónoma de México, Ap. 70543, 04510 Cd. Mx., Mexico
^{3}
Univ. Grenoble Alpes, CNRS, IPAG, 38000 Grenoble, France
^{4}
Institut d’Astrophysique Spatiale, CNRS UMR 8617, Université ParisSud, 91405 Orsay, France
Received:
3
October
2017
Accepted:
3
January
2017
Context. The molecular richness of fast protostellar jets within 20–100 au of their source, despite strong ultraviolet irradiation, remains a challenge for the models investigated so far. Aim.We aim to investigate the effect of interaction between a timevariable jet and a surrounding steady disk wind, to assess the possibility of jet chemical enrichement by the wind, and the characteristic signatures of such a configuration.
Methods. We have constructed an analytic model of a jet bow shock driven into a surrounding slower disk wind in the thin shell approximation. The refilling of the post bow shock cavity from below by the disk wind is also studied. An extension of the model to the case of two or more successive internal working surfaces (IWS) is made. We then compared this analytic model with numerical simulations with and without a surrounding disk wind.
Results. We find that at early times (of order the variability period), jet bow shocks travel in refilled pristine disk wind material, before interacting with the cocoon of older bow shocks. This opens the possibility of bow shock chemical enrichment (if the disk wind is molecular and dusty) and of probing the unperturbed disk wind structure near the jet base. Several distinctive signatures of the presence of a surrounding disk wind are identified, in the bow shock morphology and kinematics. Numerical simulations validate our analytical approach and further show that at large scale, the passage of many jet IWS inside a disk wind produces a stationary Vshaped cavity, closing down onto the axis at a finite distance from the source.
Key words: ISM: jets and outflows / HerbigHaro objects / shock waves / hydrodynamics / ISM: kinematics and dynamics
© ESO 2018
1. Introduction
Protostellar jets appear intimately linked to the process of mass accretion onto the growing star; their strikingly similar properties across protostellar age, mass, and accretion rate all point to universal ejection and collimation mechanisms (Cabrit 2002; Cabrit et al. 2007; Ellerbroek et al. 2013). Yet, jets from the youngest protostars – socalled Class 0 – are much brighter in molecules (e.g., Tafalla et al. 2000) than jets from more evolved protostars and premain sequence stars which are mainly atomic; Molecules have been traced as close as 20–100 au from the source (e.g., Lee et al. 2017; Hodapp & Chini 2014). The origin of this selective molecular richness remains an important issue for models of the jet origin. Three broad scenarios have been considered, with no fully validated answer so far.
In models of ejection from the stellar magnetosphere or the inner disk edge (e.g., Shu et al. 1994; Romanova et al. 2002; Matt & Pudritz 2005; Zanni & Ferreira 2013), the jet would be expected to be dustfree (the grain sublimation radius around a typical solarmass protostar is R_{sub} ~ 0.3 au, see for example Yvart et al. 2016). The lack of dust screening then makes the wind extremely sensitive to photodissociation by the accretion shock. Chemical models of dustfree winds by Glassgold et al. (1991) found that CO, SiO, and H_{2}O could no longer form at the wind base in the presence of a typical expected level of FUV excess^{1}. Raga et al. (2005) showed that H_{2} could form further out behind internal shocks. However, the key ions involved are also easily destroyed by FUV photons. Hence, molecule formation in a dustfree jet within 20–100 au of protostars remains an open issue.
A second proposed explanation is that the molecular component of jets may be tracing dusty MHD disk winds launched beyond R_{sub}, where dust can shield molecules against the FUV field and allow faster H_{2} reformation. Detailed models are successful at reproducing the higher molecule richness of Class 0 jets (Panoglou et al. 2012) the broad water line components revealed by Herschel/HIFI (Yvart et al. 2016), and the rotation signatures recently resolved by ALMA in the HH212 jet (Lee et al. 2017)and in the slow wider angle wind surrounding it (Tabone et al. 2017). However, the same disk wind models predict that the fastest, SiOrich streamlines in HH212 (flowing at ~ 100 km s^{−1}) would be launched from 0.05 – 0.2 au, within the dust sublimation radius (Tabone et al. 2017). Hence, this scenario still partly faces the unsolved question of molecule survival in a dustfree wind.
A third scenario is that molecules could be somehow “entrained” from the surroundings into the jet, assumed initially atomic. In a timedependent jet, travelling internal shocks will squeeze out highpressure jet material, which then sweeps up the surrounding gas into a curved bowshock. If the surrounding material is molecular, a partly molecular bowshock will result, with a more tenuous “wake” of shocked molecular gas trailing behind it (Raga & Cabrit 1993; Raga et al. 2005). As the next “internal working surface” (IWS) propagates into this wake, it may again produce a molecular jet bowshock. However, after the passage of many such IWS, the wake will be so shockprocessed and tenuous that not enough molecules may be left to produce molecular bowshocks close to the jet axis.
In the present paper, we revisit this last scenario in a new light by investigating whether a slower molecular “disk wind” surrounding the jet could help refill the wake and reinject fresh (unprocessed) molecules into the jet path. This new outlook is prompted by the discovery of a potential molecular disk wind wrapped around the dense axial jet in HH212 (Tabone et al. 2017). We explore this possibility by studying analytically and numerically the propagation of bow shocks driven by a timevariable, inner jet into a surrounding slower disk wind. This scenario may be seen as an extension of the recent modeling work of White et al. (2016) who studied the turbulent mixing layer between a jet and disk wind, with the novel addition of internal working surfaces (IWS) in the jet to produce a stronger coupling between the two outflow components.
Besides our main goal of exploring the impact of a DW on the chemical richness of Class 0 jets, our study has two other important motivations. First, we aim to identify specific signatures in the morphology and kinematics of jet bowshocks that could reveal the presence of a surrounding DW. Secondly, we aim to identify in which regions of space the pristine DW material would remain unperturbed, for comparison to theoretical DW models.
In the present exploratory study, we have limited ourselves to purely hydrodynamical and cylindrical flows, which allow us to develop an analytical model that greatly helps to capture the main effects of the twoflow interaction, and to understand the numerical results. Also, this is expected to be an optimal case for interaction between the two flows, as magnetic tension would tend to oppose mixing.
The paper is organized as follows. In Sect. 2, we build an analytical model (in the thin shell approximation) for the propagation of a bow shock driven by an IWS into a surrounding disk wind. The model is extended to the case of two or more successive IWS. In Sect. 3, we compare the analytic model with axisymmetric simulations of a variable jet + surrounding disk wind configuration, and compare the results with a “reference simulation” in which the same variable jet propagates into a stationary environment. Finally, the results are discussed in Sect. 4.
2. Analytical approach
2.1. Basic equations
We considered the “disk wind + jet” configuration shown in Fig. 1 where a cylindrical jet of radius r_{j} and timevariable velocity v_{j} directed along the zaxis is immersed in a planeparallel “disk wind” with uniform density ρ_{w} and timeindependent velocity v_{w} parallel to v_{j}.
Fig. 1
Schematic diagram of the flow around an internal working surface (IWS) in the frame of reference comoving with the IWS at velocity v_{z} = v_{j} (a similar configuration would apply for the leading working surface of a jet). The thick, horizontal line at the bottom of the graph is the jet (with a gap showing the position of the IWS). The working surface ejects jet material sideways at an initial velocity v_{0} into the slower disk wind, which in this frame of reference moves towards the outflow source at velocity v_{j} – v_{w}. The distance x is measured towards the outflow source.The shape of the thin shell bow shock is given by r_{b}(x) (see Eq. (5)), and terminates at the cylindrical radius r_{b; f}(t) with t the time elapsed since formation of the IWS (see Eq. (8)). 
The jet velocity variation is such that an internal working surface is produced within the jet beam. In the following derivations, we assume that the working surface is formed at t = 0 at the position of the source (i.e., z = 0), and that it then travels at a constant velocity v_{j} (for t > 0). Such a working surface could be produced, for example, by an outflow velocity with a constant value v_{1} < v_{j} for t < 0, jumping to a constant value v_{2} > v_{j} for t ≥ 0. Note that if the shock is produced at distance z_{s} > 0 at a time t_{s} > 0, the equations below remain valid with the transformation z → z – z_{s} and t → t – t_{s}.
In a frame of reference moving with the internal working surface (see Fig. 1) the overpressured shocked jet material which is ejected sideways from the working surface interacts with the slower moving, surrounding disk wind. In the strong radiative cooling limit, this sideways ejection leads to the formation of a thinshell bow shock, which sweeps up material of the surrounding disk wind, flowing towards the source at a relative velocity (v_{j}–v_{w}).
Assuming full mixing between jet and disk wind material, we can write the mass, r and xmomentum conservation equations at any point of radius r_{b} along this thinshell (r, x, and r_{b} being defined in Fig. 1) as:(1) (2) (3)
where ṁ is the mass rate, the rmomentum rate and the xmomentum rate of the mixed jet + disk wind material flowing along the thinshell bow shock, and ṁ_{o} and v_{0} are the mass rate and velocity (respectively) at which jet material is initially ejected sideways by the working surface. These equations have straightforward interpretations. As an illustration, we point out that Eq. (2) states that the radial momentum of the material flowing along the thin shell remains constant over time (which is due to the fact that the disk wind material adds no rmomentum), so that its radial velocity v_{r} decreases as ṁ increases, the rmomentum being shared with a larger amount of zero rmomentum material from the disk wind.
The mass rate ṁ_{o} and velocity v_{0} of the sideways ejected material are determined only by the properties of the working surface. For a highly radiative working surface, we would expect the postshock jet material to cool to ~ 10^{4} K before exiting sonically into the disk wind. Therefore, we would expect v_{0} ∼ 10 km s^{−1}. The mass rate ṁ_{0} will have values of the order of the (timedependent) mass loss rate in the jet beam.
We note that although our basic equations are similar to those of Ostriker et al. (2001) our approaches and derived equations will differ. They considered only the case of the leading jet bowshock propagating in a medium at rest (v_{w} = 0), so that the injected mass and momentum rates ṁ_{0} and ṁ_{0}v_{o} were expressed as a function of the velocity of the shock and the jet radius. Here, we keep ṁ_{0} and v_{0} as explicit parameters, so that we can consider a moving surrounding disk wind of arbitrary velocity v_{w}, and an arbitrarily small r_{j}.
2.2 Shape of the bow shock shell
For a disk wind with positionindependent density ρ_{w} and velocity v_{w}, the integrals in Eqs. (1) and (3) can be trivially performed, and from the ratio of Eqs. (2) and (3) one obtains the differential equation of r_{b}(x) :(4)
which can be integrated to obtain the shape of the thin shell bow shock as a function of x in the IWS reference frame:(5)
where we defined the characteristic scale^{2} (6)
Clearly, as the thin shell bow shock began to grow at t = 0, the solution given by Eq. (5) must terminate at a finite maximum radius r_{b,f} (see Fig. 1). The growth of this outer radius with time can be calculated combining Eqs. (1) and (2) to obtain:(7)
which can be integrated with the boundary condition r_{b,f} to obtain r_{b,f}(t) at the current time t:(8)
Now, in order to obtain the shape of the bowshock shell in the source frame (z, r) (see Fig. 2, cyan curve), when the IWS is located at distance z_{j} from the source, we simply insert the relation x = (z_{j} – z_{b }) into Eq. (5) and t = t ≡ z_{j}/v_{j} in Eq. (8). In the “narrow jet” limit where r_{j}→ 0, the thin shell bow shock has the simple cubic shape given by equation:(10)
Fig. 2
Schematic diagram showing the flow around a working surface of a jet (in this case the leading working surface, but the diagram also applies for an internal working surface). The jet is the horizontal, red rectangle at the bottom of the graph, with the source located at z = 0. The working surface in yellow is located at a distance z_{j} from the source and travels at a velocity v_{j}. It ejects material away from the axis at an initial velocity v_{0}. The jet is surrounded by a “disk wind”, which travels along the outflow axis at a velocity v_{w}. The shape of the thinshell bow shock (thick cyan line) is given by z_{b} as a function of r and ends at the edge of the bow wing (cyan point; z_{b f}, r_{b f}). The bow shock leaves behind a “cavity” (black region), which is partially refilled by the disk wind (brown region). The boundaries of the initial swept up cavity (in black dashed line) and of the refilled region (cyan dashdotted line) are given by z_{f} an z_{c} (respectively) as a function of cylindrical radius r (see Eqs. (14) and (15)). 
ending at the maximum “outer edge” radius r_{b, f} (see the cyan dot in Fig. 2) given by Eq. (8) evaluated at t = t_{j}:(11)
For a “wide jet” where is no longer negligible, the corresponding equations can also straightforwardly be obtained from Eqs. (5) and (8), and are given in Appendix 5. In the following, we consider the narrow jet regime, which leads to simpler equations.
2.3 The postbow shock cavity
Let us now consider the trajectory r_{f}(z_{f}) described in the z, r plane by the outer edge of the thin shell bow shock at earlier times, when the IWS travelled from its formation point z = 0 at t = 0 to its current location z_{j} at time t_{j}. This trajectory will define the shape of the volume swept out by the travelling and expanding bowshock into the slower disk wind (see Fig. 2, dashed black line).
At an earlier time t_{f} (0 ≤ t_{f} ≤ t_{j}), the bowshock terminated at an outer radius r_{f} ≤ r_{b, f} given by Eq. (11) with t _{j}= : t_{f} (12)
The distance z_{f} from the source where this radius r_{f} was reached is obtained from Eq. (10) by setting z_{b} = z_{f}, r_{b} = r_{f} and z_{j} = v_{j}t_{f:} (13)
Combining Eqs. (12) and (13) to eliminate t_{f}, and recalling that γ = (v_{j} – v_{w}/)v_{0} we then obtain the shape r_{f}(z_{f}) of the cavity swept by the (growing) edge of the bow shock wing associated with the travelling internal working surface (see dashed black curve in Fig. 2) :(14)
2.4 Refilling of the cavity by the disk wind
Of course, as soon as the bow shock wing has passed by, the disk wind (travelling in the zdirection at a velocity v_{w}, see Fig. 2) immediately starts to refill the sweptup cavity. For a given radius r_{f}(z_{f}) along the boundary of the sweptup volume, the refilling by the disk wind will thus start at the time t_{f} (given by Eq. (13)) when the bowshock edge reached this position; at the present time t_{j} the disk wind will have refilled a region of length (t_{j} – t_{f})v_{w} along the zaxis. The boundary between the windrefilled region and the emptied cavity thus has a locus z_{c}(r_{c}) (see the cyan dashdotted line in Fig. 2) given by:(15)
where for the second equality we have used Eqs. (12) and (13) and set r_{c} = r_{f}.
Therefore, the slower disk wind refills the cavity swept by the bow shock except for an inner, conical “hole” with halfopening angle α = arctan γ^{−1} = arctan[v_{0} / (v_{j}–v_{w})]. The conical cavity is attached to the wings of the bow shock at (z_{bf}, r_{bf}), and its vertex along the jet axis is located at a distance from the source z_{a} = v_{w}t_{j} = z_{j} (v_{w}/v_{j}) (see Eq. (15) with r_{c} = 0 and cyan asterisk in Fig. 2).
Figure 3 shows the analytical flow configurations obtained at three different evolutionary times (corresponding to z_{a} = v_{w}t_{j} = z_{j} (v_{w}/v_{j}), and), and for two choices of the wind velocity (v_{w} = 0 andv_{w} = 0.4v_{j}). In the two models, we have set . v_{o} = 0.2 The model with v_{w} = 0 (left frames of Fig. 3) produces a cavity which does not fill up. For v_{w} = 0.4v_{j} (right frames of Fig. 3), the bow shock has a more stubby shape compared to the v_{w} = 0 bow shock (L_{0} is larger) and the cavity which it leaves behind is partially refilled by the disk wind (brown region).
Fig. 3
Time evolution of the bow shock + cavity flow predicted by the analytic model for two choices of the wind velocity: v_{w} = 0 (leftframe) and v_{w} = 0.4v_{j} (right frame). The dark region is the empty part of the cavity (sweptup in the thin shell bow shock) and the brown region is the part of the cavity that has been refilled by the disk wind (this region being of course absent in the v_{w} = 0 model of the left frame). For both models, we show snapshots corresponding to t = 2L_{0}/V_{j}, 4L_{0}/V_{j}, and 8L_{0}/v_{j}), which result in positions z_{j} = 2L_{0}, 4L_{0}, and 8L_{0} for the working surface (see the labels on the top left of each frame). Both models have v_{0} = 0.2v_{j} 
2.5 Kinematics along the shell
From Eqs. (1)–(3), it is straightforward to show that for a narrow jet surrounded by a homogeneous disk wind, the radial and axial velocities (in the source rest frame) of the wellmixed thin shell material as a function of cylindrical radius r_{b} are :(16) (17)
where for the second equation we have also considered that v_{z} = v_{j} – v_{x} (see Figs. 1 and 2). In evaluating the radial velocities, one should keep in mind that the radius r_{b} is always smaller than the r_{b,f} value given by Eq. (11).
As expected, we find the following asymptotic limits:

v_{r} has an initial value v_{0} for r_{b} → 0 (i.e., as it leaves the jet working surface) and goes to 0 at large radii (as the radial momentum of the thin shell bow shock is shared with an increasing mass of disk wind),

v_{z} has an initial value v_{j} when it leaves the working surface (r_{b} → 0), and for large radii tends to the disk wind velocity v_{m}.
Equations (16) and (17) give the velocity of the wellmixed material within the thin shell bow shock. These velocities correspond to the Doppler velocities observed in an astronomical observation provided that the emission does indeed come from fullymixed material.
Another extreme limit is if the emission is actually dominated by the gas that has just gone through the bow shock and which is not yet mixed with the thin shell flow material. In this case, the axial and radial velocities of the emitting material would correspond to the velocity directly behind a highly compressive radiative shock. For such a shock, the velocity of the post shock flow (measured in the reference system moving with the bow shock) is basically equal to the projection of the incoming flow velocity parallel to the shock front. It is straightforward to show that in this case the immediate post shock radial and axial velocities of the emitting material in the source rest frame are given – in the narrow jet limit – by:(18) (19)
We note that while v_{z,ps} has the same asymptotic limits as v_{z} in the full mixing case (see Eqs. (17) and (19)), the radial postshock velocity v_{r},_{ps} tends to zero both for r_{b} → 0 and for r_{b} → ∞ (see Eq. (18)), reaching a maximum value of (v_{j} – v_{w})/2 for a radius equal to This peak value for the radial velocity is a general result of bow shock kinematics, valid regardless of the shape of the bow shock, which was first derived by Hartigan et al. (1987).
By combining Eqs. (16)– (19) with Eq. (10) it is straightforward to obtain the axial and radial velocities as a function of the distance z along the symmetry axis. Examples of these dependencies are shown in the following section.
2.6 A dimensional example
We now consider a particular model of an internal working surface moving at a velocity v_{j} = 100 km s^{−1}, located at z_{j} = 10^{16} cm along the zaxis and ejecting side way material at a rate ṁ_{0} = 10^{−8} M_{ʘ} yr^{−1} with a lateral ejection velocity v_{0} = 10 km s^{−1}.
For the surrounding disk wind, we assume a number density of atomic nuclei n_{w} = ρ_{w}/1.4 m_{H} = 10^{4} cm^{−3} and velocities v_{w} = 0 and v_{w} = 40 km s^{−1}. With these parameters, we obtain L_{0} = 5.2 × 10^{14} cm (for v_{w} = 0) and L_{0} = 8.7 × 10^{14} cm (for v_{w} = 40 km s^{−1}). Note that Jṁ_{0} and ρ_{w} are only involved in the shape and kinematic equations through L_{0} ∞ (ṁ_{0}/ρ_{w})^{1/2}, so that only their ratio actually matters in defining the flow properties.
For these two working surfaces, we obtain the shapes, and the radial and axial velocities (as a function of z) shown in Fig. 4. From this figure, it is clear that for the vw = 40 km s^{−1} model we obtain a flatter working surface than for the vw = 0 case, because L_{0} is larger.
Fig. 4
Shape of the bow shock and the cavity (top), radial velocities v_{r} (center) and axial velocities v_{z} (bottom) for the two models discussed in the text. The solid curves show the velocities of the wellmixed material within the thin shell flow, and the dashed curves show the immediate postbow shock velocities. The dotted red line shows v_{z} = v_{w}. 
The velocities of the fully mixed thin shell material (shown with solid lines in Fig. 4) have the following behaviors:

v_{r} has a value of v_{0} = 10 km s^{−1} at z = z_{j}, and monotonically decreases toward (but not reaching) zero for decreasing values of z;

v_{z} (lower plots of Fig. 4) has a value of v_{j} = 100 km s^{−1} at z_{j}, and monotonically decreases for lower values of z, towards (but not reaching) an asymptotic limit of v_{w}.
The immediate postbow shock velocities (shown with dashed lines in Fig. 4) have the following behaviors:

v_{r},_{ps} (central plots of Fig. 4) is zero at the apex of the bow shock surface at z = zj, and rapidly grows to a maximum value of 50 km s^{−1} (for v_{w} = 0) and 30 km s^{−1} (for v_{w} = 40 km s^{−1}), this value corresponding to (v_{j} – v_{w})/2, as discussed in the previous section is reached at . The radial velocity then decreases again at smaller z until the end of the bowshock wings;

the axial velocity v_{zps} has the same qualitative behavior as the wellmixed v_{z} (see above), but with a different functional form that approaches faster its limit v_{w} in the bowshock wings.
We expect that in reality, due to incomplete mixing, the emitting material will have axial and radial velocities between the fullymixed layer and immediate postbow shock velocities shown in Fig. 4. The difference between these two velocities is particularly important for the radial component of the velocity of the emitting material.
2.7. Successive bow shocks
In the previous section, we assumed that the bow shock associated with an internal working surface travels through undisturbed disk wind material. However, we saw that the cavity formed behind it is only partially refilled by fresh disk wind. Therefore, a second bow shock will travel into a disk wind structure containing an empty, conical cavity left behind by the first bow shock.
We now assume that the variable ejection velocity of the jet produces a second working surface at z = 0 at a time τ_{j}, which also travels along the jet axis with the same velocity v_{j} as the first working surface. Figure 5 illustrates three steps of the propagation of this second working surface.
Fig. 5
Time evolution of two successive bow shocks and the cavities predicted by the analytical model. The first bow shock is ejected at t = 0, the second shock is formed at time t = τ_{j} when the first bow shock is at z = v_{j}τ_{j} taken here as 6L_{0} (top). At a time t = τ_{j} + 2L_{0}/v_{j}, the second bow shock is still propagating in the undisturbed disk wind material (center). At a time t_{c} the second bow shock catches up the emptied cavity of the first bow shock at its vertex (bottom). 
At t = τ_{j} (Fig. 5, first panel) the first working surface is at a distance z_{j1} = v_{j}τ_{j} from the outflow source and its cavity is partially filled with fresh disk wind material, while the second working surface has not yet expanded. At a time τ_{j} < t < t_{c} (Fig. 5, center) the second bow shock travels in unperturbed, pristine disk wind material that refilled the cavity behind the first bowshock; hence its shape and kinematics are still given by the same equations derived above for the leading internal working surface, and any molecules present in the disk wind can enter the second bowshock. At a time t = t_{c} (Fig. 5, bottom panel) the apex of the second bowshock just catches up with the vertex of the conical cavity emptied by the first bow shock and not refilled by the disk wind.
To obtain the time t_{c}, we note that at any time t (τ_{j} < t < t_{c}), the position of the apex of the second bow shock is z_{j2} = v_{j}(t – τ_{j}) and the position of the vertex of the empty cavity behind the first bow shock is z_{a1} = v_{w}t. By equating these two quantities we get(20)
This interaction occurs at a distance l_{c} from the source:(21)
where Δz = T_{j}v_{j} denotes the distance between two successive IWS. Unless v_{w} is very close to v_{j}, we find that l_{c} is of the order of a few times the typical IWS spacing.
Our model thus predicts that no more pristine unperturbed disk wind material can remain close to the jet axis beyond z = l_{c}. When the second IWS reaches z > l_{c}, the central region of its bow shock shell propagates into the emptied cavity left behind by the previous IWS. This second bow shock shell will in general become less curved than the first one, because its central region travels into a low density cavity instead of unperturbed disk wind.
3. Numerical simulations
In the previous section, we proposed a simple analytical “thinshell” model that describes the morphology and the kinematics of a bow shock produced by a pulsating jet travelling in a surrounding disk wind. We especially show that the disk wind refills up part of the cavity carved by the bow shock, allowing us to observe pristine disk wind close the source. For successive bow shocks, bow shocks travels in an undisturbed disk wind up to a critical distance l _{c}. Above this altitude, bow shock shells interact with each other and analytical models can only be heuristic.
In this section, we present numerical simulations that start with the simple configuration adopted above, first to determine to what extent the analytical model can be used to describe a realistic situation – e.g., with a partial mixing – and secondly to study briefly the long term evolution of the interacting bow shock shells.
3.1. Numerical method and setup
We carry out numerical simulations of a variable ejection jet surrounded by a wide “disk wind” outflow. We have implemented the new HD numerical code Coyotl which solves the “2.5D” Euler ideal fluid equations in cylindrical coordinates:(22)
where U is the vector of conserved quantities(23)
with fluxes in the r and z directions given, respectively, by(24) (25)
n_{i} are passive scalars used to separate the jet from the diskwind material in the flow. Assuming an ideal equation of state, the total energy density e is(26)
where the cooling function Λ(T) is the parametrized atomic/ionic cooling term of Raga & Canto (1989), which approximates the cooling curve of Raymond et al. (1976) for temperatures above 10^{4} K.
The numerical scheme is based on a second order Godunov method with an HLLC Riemann solver (Toro 1999). The calculation of the fluxes and data reconstruction uses the second order scheme described by Falle (1991). This algorithm solves Euler equations in a true cylindrical coordinate system as written in Eq. (22) and calculates the cell gradients through the center of gravity of the cylindrical cells.
We ran two simulations: a reference simulation called noDW model, with v_{w }= 0 (i.e., a jet in a stationary ambient medium) and a simulation with v_{w }= 0.4v_{j} called DW model. To follow the refilling of the cavity close to the source and the interaction between various shells, we integrate equations on a 2000 au × 350 au domain, with a resolution of 1 au per cell. All jet and wind parameters except v_{w} are kept equal between the two simulations, and are summarized in Table 1.
Model parameters.
Our initial conditions have an inner, constant velocity jet filling the r < r_{j} region at all z, and the disk wind (or external stationary medium) filling the rest of the computational domain. This setup differs from the standard jet initialization in which the jet is introduced only in a small region around z ≃ 0 and then propagates through the domain, producing a transient with a leading jet bow shock that sweeps aside the ambient medium (Blondin et al. 1990; de Gouveia dal Pino & Benz 1993). By initializing our simulations with a jet that already extends across the whole range of z, we do not perturb the surrounding medium with the transient leading jet bow shock and we can directly follow the interaction between an IWS and the unperturbed disk wind close to the outflow source.
In order to form IWS in the jet, we impose a sawtooth ejection variability with a mean velocity v_{j}, a velocity jump ∂v_{j} and a period Δτ_{j}.The ejection velocity is assumed to rise linearly from v_{j} – ∂v_{j}/2 to v_{j} + ∂v_{j}/2 during a timelapse ƞΔτ_{j} and to linearly decrease over a time (1 – η)Δτ_{j} back to a velocity v_{j} – δv_{j}/2. Using the small amplitude velocity variability approximation of (Raga et al. 1990), we estimate that for our chosen parameters in Table 1 this variability will produce an internal working surface in the jet at time t_{S }= 5 yr and distance z_{S }= (v_{j} – δv_{j}/2)^{2}ηΔτ_{j}/δv_{j }= 75 au from the central source,and that the working surface will travel in the jet at a velocity v_{j }= 96 km s^{−1}.
We adopt a jet radius r_{j }= 20 au consistent with the width of the HH 212 molecular jet obtained from VLBI measurements
by Claussen et al. (1998) and with the widths of atomic jets estimated by Dougados et al. (2000) close to the source. The density contrast between the jet and the wind is chosen to be sufficiently high (ρ_{j} /ρ_{w }= 29) to produce wide bow shock shell. Temperatures are chosen to insure a transverse pressure equilibrium between the jet and the wind. Note that simulations are not very sensitive to the jet temperature since strong internal shocks cooled down by atomic lines set the temperature of the IWS at T ∼ 10^{4} K.
The boundary conditions are reflecting on the symmetry axis (r = 0) and outflowing in the outer radial and axial cells. On the z = 0 boundary, we introduce the jet by imposing fixed constant physical conditions for r < r_{j}. For r > r_{j} we impose either the disk wind physical conditions, or a reflecting condition (for the reference simulation with v_{w} = 0). In order to avoid numerical problems due to the zvelocity shear between the jet and the surrounding disk wind we put a velocity gradient on three cells (i.e., 3 au) at the outer edge of the jet inflow.
3.2. Single bow shock propagation
Figure 6 shows snapshots of the noDW simulation (two frames on the left) and of the DW simulation (two frames on the right) after a t = 48 yr time integration, which is larger than the ejection variability period of 33 yr (see Table 1). The first internal working surface (IWS) has travelled to a distance of 995 au from the source, and a second IWS to 355 au. In this subsection, we study successively the shape of the first bow shock shell, the refilling of the cavity behind it, and the kinematics of the shell, comparing each of them with our analytical predictions.
Fig. 6
Maps of density and pressure for the reference noDW simulation with v_{w} = 0 (left) and the DW simulation with v_{w} = 0.4v_{j}(right) at a time t = 48 yr. Color scales on top are in g cm^{−3} for density and in dyn cm^{−2} for pressure. White contours show the locus of 50% mixing ratio between jet and diskwind/ambient material. The cyan curve shows a fit (to the numerical results) by the analytic shell shape in Eq. (11), with L_{0} = 65 au (left) and 108 au (right). The cyan dot indicates the maximum radius of the shell, the cyan asterix indicates the predicted vertex of the empty conical cavity left behind the shell, and the cyan dashdotted line is the analytical predicted boundary between the emptied cavity and the region refilled from below by fresh disk wind (see Fig. 2 and Eq. (15)). 
3.2.1. Shape of the bow shock shell
The cyan curves in Fig. 6 show that the bow shock shells in the two simulations can be well fitted with the cubic analytic solution for the thinshell (Eq. (5)), with values for the characteristic scale L_{0} = 65 au for the reference noDW simulation and L_{0} = 108 au for the DW simulation. In the simulations, the sideways ejection velocity v_{0} and massflux ṁ_{0} (see Sect. 2) are a result of the IWS shock configuration, which compresses the jet material and ejects it sideways. Since the two simulations only differ in the presence or lack of a surrounding disk wind, the jet IWS in the two simulations have similar characteristics. We then expect that ṁ_{0}v_{0} is the same, and that L_{0} should vary with the wind velocity as L_{0}k (v_{j} – v_{W})^{−1} (see Eq. (6)). The values of L_{0} found above by fitting the shell shape are indeed consistent with this expectation.
The cyan dots on the leading bow shock wings in Fig. 6 indicate the maximum radius of the bow shock shell as observed in the numerical simulations, r_{max} = 133 au in the noDW and r_{max} = 137 au in the DW cases. Assuming that it corresponds to the current position of the edge of the thin shell (r_{bf}, z_{bf}), as defined in Fig. 2, our analytic model predicts that r_{bf} depends on L_{0} and v_{0} through Eq. (11) or (A.2). With L_{0} = 65, 108 au, we would deduce v_{0} = 27, 19 km s^{−1} for the noDW and DW simulations, respectively.
To obtain a direct measurement of v_{0}, we plot in Fig. 7 a transverse cut of the radial velocity v_{r} at the position of the leading working surface (z = 993 au) in the DW simulation. Inside the IWS, because of both the adiabatic expansion and the mass flux across the IWS, v_{r} increases from zero to a maximum velocity of 14 km s^{−1}. This direct measurement of v_{0} is smaller than the values of 27, 19 km s^{−1} inferred from the maximum radius of the shell in our simulations using Eq. (11) (see above). However, we note that taking v_{0} = 14 km s^{−1} (its real value), the predicted r_{bf} would be r_{bf} = 111 au in the noDW case and r_{bf} = 121 au in the DW case, only slightly smaller than the r_{max} found in our simulations.
Fig. 7
Transverse cut across the flow at the IWS location (z = 993 au) in the noDW timeframe shown in Fig. 6. This cut shows the radial velocity as a function of distance from the jet axis in solid line. We also plot the radial velocity weighted by the abundance of the jet tracer with a dashed line. The radial velocity first grows outwards, reaches a maximum velocity of ≈ 14 km s^{−1} at a radius of ∼ 25 au (somewhat larger than the 20 au initial jet radius), and then remains with values > 10 km s^{−1} until it drops sharply to 0 at r ∼ 50 au. The velocity maximum at r ∼ 25 au corresponds to the shock against the jet material. The second maximum at r ∼ 50 au is the shock that propagates in the diskwind and the zero radial velocity material at larger radii is the undisturbed disk wind. 
This slight difference in outer radius between the analytic model and the numerical simulations could be a result of several effects:

the analytic model assumes a working surface with a timeindependent, sideways ejection, while the numerical simulation has an IWS with time–dependent sideways ejection that depends on the evolution of the IWS shocks. The IWS in the simulations produces a higher sideways velocity at early times (v_{0} ∼ 18 km s^{−1})^{3}, closer to the values deduced from the analytic cavity shapes;

in the numerical simulation, the sideways ejection from the IWS is not highly supersonic. The thermal gas pressure is therefore expected to be an additional source of sideways momentum for the shell (an effect not included in our momentum conserving analytic model); this will act to produce a higher effective v_{0};

similarly, the thermal pressure in the head of the bow shock driven into the surrounding environment will result in a sideways push which is not present in the momentum conserving analytic model;

the numerical simulations do not have instant mixing between the sideways ejected jet material and the shocked environment (or disk wind), as assumed in the analytic model. Since the immediate post shock velocity in the radial direction is generally greater that the radial mean shell velocity (see example in Fig. 4), the growth rate of the bow shock can be enhanced. In the reference frame of the IWS (see Fig. 1), the nonmixed material will “slide” along the shell surface, extending r_{b},_{f} to larger values.
Lee et al. (2001) found in their simulations similar disagreements between direct measurements of the sideways momentum ejected by the IWS and the momentum estimated from the fitted shape of their analytic shell model.
3.2.2. Cavity refilling
The asterisk in cyan in each panel of Fig. 6 indicates the location of the vertex of the emptied cavity as predicted from the analytic model (see Fig. 2). For the noDW simulation, this point is located at the shock formation position (z_{s} = 75 au) whereas for the DW simulation this point is located at . We also plot in cyan dash– dotted the line connecting this vertex to r _{max}, which traces the boundary predicted by the analytical model between the emptied sweptout conical cavity and the unperturbed surrounding medium/refilled disk wind (see black conical region in Fig. 2). Three important features can be seen.
First, in both numerical simulations, the emptied cavity predicted by the analytical thinshell model, i.e. the conical volume inside the dashdotted cyan line, is not really empty, but partially filled with a cocoon of low density and pressure material. No unperturbed ambient gas or disk wind can be left inside this volume (in black in Fig. 3), which was entirely sweptout by the growing shell during the IWS propagation. Hence this cocoon is made of shocked material that did not fully mix in the shell, and reexpanded in the lowpressure cavity behind it, refilling it “from above”. The white contour, which denotes the surface of 50% jet/environment mixing fraction (obtained following a passive scalar) shows that the cocoon is mainly filled with jet material close to the axis, where the shell mass is dominated by gas ejected from the IWS. Further from the axis and closer to the theoretical boundary (cyan dashdot line) it is filled by ambient material that was swept up by the bowshock and reexpanded behind it.
Secondly, the boundary with unperturbed ambient or disk wind material closes back to the axis at the predicted vertex position (see cyan asterisk Fig. 6), but is delimited by a weak shock that extends slightly outside from the predicted analytical boundary (dashdotted cyan line Fig. 6) . In the noDW model, the analytical boundary represents the trajectory of the edge of the bow shock (z_{c} = z_{f} in the case v_{w} = 0). Hence, this weak shock is produced by the supersonic motion of the highpressure edge of the bow shock (r_{b},_{f}, z_{b},_{f}) in the static surrounding medium. This launches a weak outward shock that repels the boundary of the unperturbed ambient material slightly outside the predicted cavity boundary (in cyan dashdotted line). In the presence of a supersonic disk wind, the weak shock front is advected away from the source so that it still closes back onaxis at the predicted vertex position z_{a}. Hence, Eq. (15) gives a strong limit on the boundary between perturbed and unperturbed material.
Thirdly, in the presence of a diskwind, the region between the predicted cavity boundary (cyan dashdotted line) and the weak shock front outside it is refilled by fresh disk wind material coming from below. To analyse this process we show in Fig. 8 density and velocity maps of the region around the leading bow shock of the DW simulation. The dashed white contours show 10 %, 0.1 %, and 0.001 % jet material mixing fractions. Material that went through the bowshock and reexpanded in the cocoon has also been partially mixed with jet material. As a consequence, regions where no jet material is observed are regions where the disk wind has refilled the cocoon from below. The location of the last, outer contour (corresponding to a 10^{−5} jet material mixing fraction) shows that the weak outer shock front propagates into unmixed, fresh DW material. This material manages to cross the weak shock to refill from below the bottom part of the sweptout cavity.
Fig. 8
Zoom of the leading IWS of the simulation with a surrounding disk wind at time t = 48 yr. Left: density stratification (with the logarithmic color scale given by the top bar in g cm^{−3}), center: radial velocity (with the linear scale of top bar in km s͔–1) and right: axial velocity structure (with the linear scale of the top bar in km s^{−1}). The white contours show the surfaces of 50 % (solid line), 10 %, 0.1 % and 0.001 % (outer contour) jet material mixing fractions. The cyan asterisk is the predicted vertex of the cavity, the cyan dashdotted line in the predicted boundary of the cavity, and the cyan curve is the fitted shape of the bow shock. 
The weak shock provides a slight push outwards to the refilling DW, with radial velocities that vary from + 6 km s^{ =1} to + 3 km s^{−1} along the shock front (middle panel of Fig. 8), similar to the adiabatic sound speed of the disk wind (c_{sw} ≈ 2.8 km s^{−1}). The weak shock also reduces the DW inflow velocity v_{z} to values slightly below v_{w} = 0.4v_{j} = 38.4 km s^{−1} (right panel of Fig. 8). However, refilling remains efficient up to the locus predicted by our analytical model (dashdotted cyan line), as the jet mixing fraction there remains very small (≃ 0.1 %). The presence of the weak shock does not appear to significantly modify the extent of DW refilling compared to analytical expectations.
In summary, we can therefore distinguish in our simulations three refilling regions behind the bowshock:

a low density cocoon trailing the bow shock, that is refilled from above by shell material reexpanding into the emptied cavity. This region is mainly composed of jet material close to the apex of the bow shock, and of shocked sweptup disk wind material behind the wings of the bowshock;

an intermediate region (outside the cyan dashdotted line and inside the weak shock closing the cavity) refilled from below by weakly shocked disk wind material;

a region upstream of the weak shock closing the cavity, that is refilled by unperturbed fresh disk wind keeping its initial physical conditions.
3.2.3. Kinematics
We now compare the kinematics in both simulations with our analytical predictions. Figure 9 shows “positionvelocity” (PV) diagrams for v_{r} and v_{z} as a function of distance z along the flow axis. In order to enhance the contribution from the material that has just been shocked, each pixel in a snap shot has been weighted by the cube of the pressure p^{3} times the elementary volume 2πrΔrΔz. Using this weighting, the maximum intensity (in yellow and orange shades) at each position in the PV diagrams then traces the velocity in the shell. The separation between material originating mainly from the jet or mainly from the surroundings/disk wind is done using a passive scalar. The predicted mixed shell velocities (Eqs. (16) and (17)) are shown in blue, and the predicted immediate postbowshock velocities (Eqs. (18) and (19)) are shown in magenta. Following the discussion of Fig. 7, we take v_{0} = 14 km s^{−1}.
Fig. 9
Longitudinal positionvelocity (PV) diagrams for the no DW simulation (top) and the DW simulation (bottom). From left to right: v_{r} for the jet material only, v_{r} for the surrounding material only, v_{z} for the jet material only, vz for the surrounding material only, and density stratification. The ordinate of all frames is position along the outflow axis (in au). The color scale in the PVs is scaled by volume × cube of pressure so as to be maximum (in red and yellow shades) for shocked material in the shell, while the color scale for density is in g cm^{−3}. Blue curves are predicted velocities in the full mixing hypothesis (Eqs. (16) and (17)), while magenta curves are the predicted immediatepost shock velocities (Eqs. (18) and (19)). 
In the v_{r} PV diagrams of the surrounding material, the expansion velocities of shocked material in the shell (orange shading) are always larger than predicted by the (blue) fullmixing curve (except very close to the bow shock apex where the shear is maximum). The simulation more closely follows the immediate postbow shock velocity curve (magenta), indicating that highpressure shocked material in the shell is not fully mixed in our simulations. Conversely, the v_{r} PV diagram of the jet material decreases monotonically along the bow shock wing with velocities always slightly smaller than the full mixing velocity curve (in blue).
Concerning the velocity along the jet zaxis, the v_{z} values for jet material lie close to, or slightly above the full mixing curve in blue. The v_{z} PV diagrams for the surrounding material generally show smaller v_{z} than predicted by the full mixing curves. The highpressure sweptup shell material (in orange) lies close to the immediate postshock velocities (magenta curve).
The relatively small v_{r} velocities and large v_{z} velocities observed in the jet dominated material indicate that even if the full mixing hypothesis does not hold, the momentum is still conserved: if the velocities of the sweptup surrounding material are greater than expected from the full mixing hypothesis, then the velocities of the jet material (in the IWS restframe) must be smaller than the predicted full mixing velocities (and viceversa).
As predicted, the most striking difference between the disk wind model and the reference noDW model is the saturation of the v_{z} velocity in the bowshock wings to a nonzero value of v_{z} ≈ v_{w}. Even if this asymptotic limit does not depend on any mixing, the incomplete mixing obtained in the simulations produces a more rapid convergence to v_{w} than predicted in the case of full mixing (blue curve).
3.3. Longterm evolution
Figure 10 of the reference noDW simulation (three left frames) and of the diskwind simulation (three right frames) at times t = 71, 119, and 167 yr. From this figure, we see that the morphologies of the regions perturbed by the jet after the passage of several IWS are very different in the two cases.
Fig. 10
Density maps for the noDW reference simulation (three frames on the left) and the DW simulation (three frames on the right) at t = 71, 119, and 167 yr. The white contours indicate the surface of 50% (solid line) and 90% (dashed line) jet material mixing fractions. The black lines in the diskwind simulation show a cone of α = 11° opening halfangle, which circumscribes the boundary the region perturbed by the jet and its IWS. The black dashed lines show the predicted trajectory of the edge of the bow shock (see Eq. (14)). The density color scale is given by the right bar (in g cm^{−3}). 
In the noDW simulation, the region perturbed by the jet behind the leading bowshock expands into a roughly cylindrical shape, which tapers off close to the position of the outflow source (where it becomes a weak, radially expanding shock). This is the standard shape of the perturbed region in simulations of variable, radiative jets propagating into a uniform static medium, seen since the early work of Stone & Norman (1993) and Biro & Raga (1994).
In the disk wind simulation, in contrast, the region perturbed by the jet behind the leading bowshock takes a conical shape, tapering off at large distances from the outflow source. For the parameters of our DW simulation, the halfopening angle of the perturbed region is α ≈ 11° (see Fig. 10). This cone is located outside the predicted trajectory of the edge of the bow shock (drawn in black dashed line) given by Eq. (14). This broadening occurs because, as discussed in item 3 of Sect. 3.2.2, the edge of the bow shock drives a weak outer shock into the undisturbed DW, which propagates away at a speed close to c_{s} ∼ .8 km s^{−1}. Taking into account the advection of the weak shock by the DW, one predicts that this will broaden the disturbed region by an angle ß = arctan c_{s}/v_{w} = 4°, in agreement with the observed cone opening. Obviously, in the noDW simulations, this weak outer shock travels laterally without being advected, and no limiting cone forms.
In this surprisingly simple configuration adopted by our jet + disk wind simulation, the overall longterm effect of the disk wind is to stop the perturbations from travelling beyond this “opening cone” of the sideways ejection from the IWS. Another important effect of the DW is to push the locus of 50 % ambient material (white contour) closer to the jet axis than in the noDW case, due to the disk wind partial refilling behind each bowshock. Hence, the first few IWS close to the source can still sweep up (possibly molecular) DW material. The internal IWS are also more curved than in the noDW case, where material ejected sideways meets a very lowdensity cocoon, producing flattopped internal bowshocks (see Fig. 10).
4. Summary
In this paper we have presented a first exploration of an hydrodynamical flow composed by an inner, variable jet surrounded by a slower, steady, cylindrical disk wind. The jet variability produces IWS which drive bow shocks into the disk wind, producing a strong coupling between the two components of the flow.
We have developed a standard thin shell model for the bow shock driven into the disk wind by a single IWS, for a jet of arbitrarily small radius (see Sect. 2), deriving the shape of the bow shock and the refilling by the continuing disk wind of the cavity left behind by the bow shock. The model was extended
to give a qualitative description of the flow resulting from two or more successive IWS bow shocks plowing through the disk wind (Sect. 3).
The appropriateness and limitations of the predictions of bow shock shapes and kinematics from this analytic model have been checked with axisymmetric numerical simulations: one of a variable jet + disk wind configuration, and a second reference simulation with the same variable jet surrounded by a stationary environment. We compared the analytic model with the numerical simulations, and we found a relatively good agreement, giving us an understanding of the main features of the simulated flows. These features are:

the bow shocks of the numerical IWS have cubic morphologies which can be reproduced quite convincingly with the thinshell, momentum conserving analytic model (see Eqs. (10) and (11) and Fig. 6);

the kinematics in the simulated bow shocks has a behavior which approximately follows the kinematics predicted from the analytic model for the fullymixed layer (for jetdominated material) or the immediate postbow shock gas (for highpressure sweptup ambient gas) (see Figs. 4 and 9);

these bow shocks leave behind cavities which are partially refilled by the slower disk wind (see Figs. 3 and 8);

thanks to this refilling, subsequent IWS will propagate into fresh disk wind material up to a distance from the source l_{c}= Δz/(v_{j}/v_{w} –1) (see Figs. 5 and 10).
The main contribution of this paper is thus to provide a numerically validated, simple analytic model which can be used to model bowlike shapes of knots observed close to the outflow sources in high velocity, collimated optical and molecular outflows (LavalleyFouquet et al. 2000; Podio et al. 2015). As shown by our simulations, this shape modeling (in the narrow jet limit) allows one to estimate the sideways ejection velocity from the IWS and the length scale of the bowshock. From this, constraints could be inferred on the massloss rate from the IWS and the surrounding flow properties (see Eq. (6)).
Another important contribution of this paper is to predict the regions where a surrounding disk wind would remain unperturbed. A quite dramatic result of our jet + disk wind simulation is that the perturbations of the disk wind by the IWS bow shocks are confined inside a cone. Therefore, all of the gas outside this confinement cone is unperturbed disk wind material. Also, there are pockets of undisturbed disk wind material within this cone, in the refilled region between the source and the last IWS, and also ahead of the latest IWS when it is at z < l_{c} (see the three right hand frames of Fig. 10). These are the regions in which one still finds a record of the undisturbed characteristics of the disk wind, which could be useful for comparisons with disk wind models.
Finally, another result of observational interest is that we identify several distinctive signs of a cylindrical DW around a timevariable jet: (i) bow shocks that close upon the axis at a finite distance from the source (at a fraction v_{w}/v_{j} of the distance to the bow apex), (ii) a nonzero (= v_{w}) asymptotic value of longitudinal velocity in the far bowshock wings, (iii) internal bowshocks that are curved rather than flattopped, (iv) a predominance of DW material ahead of the first few IWS, which (if the DW is chemically richer and/or dustier than the jet) should produce different emission signatures compared to the more distant IWS.
Extensions of the analytic model to more complex jet + disk wind flows do not appear very attractive (as, e.g., relaxing the assumption of a cylindrical uniform disk wind) as quite complex expressions are obtained, and are therefore not straightforwardly applicable to model observed structures. On the other hand, the numerical simulations presented here can be extended in many directions, including a more realistic disk wind model (e.g., with a radial dependence of the density and velocity, and a velocity not aligned with the outflow axis); studying the effect of a nontop hat jet cross section; going from the HD to the MHD equations; including a chemical/ionic network and the associated cooling functions. If future comparisons between jet + disk wind models and observations are sufficiently promising, the extensions listed above (as well as other easily imagined possibilities) will become worthy of exploration.
Acknowledgments
We thank the anonymous referee for useful comments. This work has been done within the LABEX PLAS@PAR project, and received financial state aid managed by the Agence Nationale de la Recherche, as part of the “Programme d’Investissements d’Avenir” under the reference ANR11IDEX 000402. AR acknowledges support from the DGAPA (UNAM) grant IG100218. This research has made use of NASA’s Astrophysics Data System.
Appendix A: General equations
In the analytic part of this work, equations ruling the geometry and the kinematics of a bow shock travelling in a disk wind are given for simplicity in the narrow jet limit r_{j} → 0. In this appendix, we give equations valid for an arbitrary jet width that we have used to fit numerical simulations. For the definition of the quantities, we refer to Figs. 1 and 2.
A.1. Shape of the bowshock and of the cavity
Eqs. (1) to (8) are valid for an arbitrary width of a jet. Inserting z_{b} = z_{j} – x into Eq. (5) we get the shape of the bow shock (z_{b}, r_{b}) (see thick cyan line Fig. 1):(A.1)
Eq. (8) gives straightforwardly the radius r_{bf} of the edge of the bow shock shell(A.2)
Combining Eqs. (A.1) and (A.2) we get the trajectory of the outer edge of the cavity (black dotted line Fig. 1):(A.3)
The boundary of the partially refilled cavity (cyan dashdotted line Fig. 1) obtained from Eqs. (A.2) and (A.3) is given by:(A.4)
Hence, in the wide jet case, the boundary between the refilled region and the empty cavity has a conical shape.
A.2. Kinematics
Integration of Eqs. (1)–(3) gives the fully mixed radial and axial velocities:(A.5) (A.6)
Immediate postshock velocities obtained by considering the velocity component tangential to the shock surface are:(A.7) (A.8)
The flat UV flux in their UV1–UV2 models is ~30–500 that in BP Tau (Bergin et al. 2003), for a windmass flux corresponding to a 1000 times larger accretion rate (~3 × 10^{−5} M_{ʘ} yr^{−1}).
Noting that L_{0} is the radius where the sweptup xmomentum is equal to 3 times the injected rmomentum ṁ_{0} v_{0}, Eq. (5) is equivalent to Eq. (22) in Ostriker et al. (2001).
References
 Bergin, E., Calvet, N., D’Alessio, P., & Herczeg, G. J. 2003, ApJ, 591, L159 [NASA ADS] [CrossRef] [Google Scholar]
 Biro, S., & Raga, A. C. 1994, ApJ, 434, 221 [NASA ADS] [CrossRef] [Google Scholar]
 Blondin, J. M., Fryxell, B. A., & Konigl, A. 1990, ApJ, 360, 370 [NASA ADS] [CrossRef] [Google Scholar]
 Cabrit, S. 2002, in Proceedings of “Star Formation and the Physics of Young Stars”, eds. J. Bouvier & J.P. Zahn (Les Ulis, Fr: EDP Sciences), EAS Pub. Ser, 3, 147 [Google Scholar]
 Cabrit, S., Codella, C., Gueth, F., et al. 2007, A&A, 468, L29 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Claussen, M. J., Marvel, K. B., Wootten, A., & Wilking, B. A. 1998, ApJ, 507, L79 [NASA ADS] [CrossRef] [Google Scholar]
 de Gouveia dal Pino, E. M., & Benz, W. 1993, ApJ, 410, 686 [Google Scholar]
 Dougados, C., Cabrit, S., Lavalley, C., & Ménard, F. 2000, A&A, 357, L61 [NASA ADS] [Google Scholar]
 Ellerbroek, L. E., Podio, L., Kaper, L., et al. 2013, A&A, 551, A5 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Falle, S. A. E. G. 1991, MNRAS, 250, 581 [NASA ADS] [CrossRef] [Google Scholar]
 Glassgold, A. E., Mamon, G. A., & Huggins, P. J. 1991, ApJ, 373, 254 [NASA ADS] [CrossRef] [Google Scholar]
 Hartigan, P., Raymond, J., & Hartmann, L. 1987, ApJ, 316, 323 [NASA ADS] [CrossRef] [Google Scholar]
 Hodapp, K. W., & Chini, R. 2014, ApJ, 794, 169 [NASA ADS] [CrossRef] [Google Scholar]
 LavalleyFouquet, C., Cabrit, S., & Dougados, C. 2000, A&A, 356, L41 [NASA ADS] [Google Scholar]
 Lee, C.F., Stone, J. M., Ostriker, E. C., & Mundy, L. G. 2001, ApJ, 557, 429 [NASA ADS] [CrossRef] [Google Scholar]
 Lee, C.F., Ho, P. T. P., Li, Z.Y., et al. 2017, Nat. Astron., 1, 0152 [NASA ADS] [CrossRef] [Google Scholar]
 Matt, S., & Pudritz, R. E. 2005, ApJ, 632, L135 [NASA ADS] [CrossRef] [Google Scholar]
 Ostriker, E. C., Lee, C.F., Stone, J. M., & Mundy, L. G. 2001, ApJ, 557, 443 [NASA ADS] [CrossRef] [Google Scholar]
 Panoglou, D., Cabrit, S., Pineau Des Forêts, G., et al. 2012, A&A, 538, A2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Podio, L., Codella, C., Gueth, F., et al. 2015, A&A, 581, A85 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Raga, A., & Cabrit, S. 1993, A&A, 278, 267 [Google Scholar]
 Raga, A. C., & Canto, J. 1989, ApJ, 344, 404 [NASA ADS] [CrossRef] [Google Scholar]
 Raga, A. C., Binette, L., Canto, J., & Calvet, N. 1990, ApJ, 364, 601 [NASA ADS] [CrossRef] [Google Scholar]
 Raga, A. C., Williams, D. A., & Lim, A. J. 2005, Rev. Mex. Astron. Astrofis., 41, 137 [Google Scholar]
 Raymond, J. C., Cox, D. P., & Smith, B. W. 1976, ApJ, 204, 290 [NASA ADS] [CrossRef] [Google Scholar]
 Romanova, M. M., Ustyugova, G. V., Koldoba, A. V., & Lovelace, R. V. E. 2002, ApJ, 578, 420 [NASA ADS] [CrossRef] [Google Scholar]
 Shu, F., Najita, J., Ostriker, E., et al. 1994, ApJ, 429, 781 [NASA ADS] [CrossRef] [Google Scholar]
 Stone, J. M., & Norman, M. L. 1993, ApJ, 413, 198 [NASA ADS] [CrossRef] [Google Scholar]
 Tabone, B., Cabrit, S., Bianchi, E., et al. 2017, A&A, 607, L6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Tafalla, M., Myers, P. C., Mardones, D., & Bachiller, R. 2000, A&A, 359, 967 [NASA ADS] [Google Scholar]
 Toro, E. F. 1999, Riemann solvers and numerical methods for fluid dynamics: a practical introduction (Berlin: Springer) [CrossRef] [Google Scholar]
 White, M. C., Bicknell, G. V., Sutherland, R. S., Salmeron, R., & McGregor, P. J. 2016, MNRAS, 455, 2042 [NASA ADS] [CrossRef] [Google Scholar]
 Yvart, W., Cabrit, S., Pineau des Forêts, G., & Ferreira, J. 2016, A&A, 585, A74 [Google Scholar]
 Zanni, C., & Ferreira, J. 2013, A&A, 550, A99 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
All Tables
All Figures
Fig. 1
Schematic diagram of the flow around an internal working surface (IWS) in the frame of reference comoving with the IWS at velocity v_{z} = v_{j} (a similar configuration would apply for the leading working surface of a jet). The thick, horizontal line at the bottom of the graph is the jet (with a gap showing the position of the IWS). The working surface ejects jet material sideways at an initial velocity v_{0} into the slower disk wind, which in this frame of reference moves towards the outflow source at velocity v_{j} – v_{w}. The distance x is measured towards the outflow source.The shape of the thin shell bow shock is given by r_{b}(x) (see Eq. (5)), and terminates at the cylindrical radius r_{b; f}(t) with t the time elapsed since formation of the IWS (see Eq. (8)). 

In the text 
Fig. 2
Schematic diagram showing the flow around a working surface of a jet (in this case the leading working surface, but the diagram also applies for an internal working surface). The jet is the horizontal, red rectangle at the bottom of the graph, with the source located at z = 0. The working surface in yellow is located at a distance z_{j} from the source and travels at a velocity v_{j}. It ejects material away from the axis at an initial velocity v_{0}. The jet is surrounded by a “disk wind”, which travels along the outflow axis at a velocity v_{w}. The shape of the thinshell bow shock (thick cyan line) is given by z_{b} as a function of r and ends at the edge of the bow wing (cyan point; z_{b f}, r_{b f}). The bow shock leaves behind a “cavity” (black region), which is partially refilled by the disk wind (brown region). The boundaries of the initial swept up cavity (in black dashed line) and of the refilled region (cyan dashdotted line) are given by z_{f} an z_{c} (respectively) as a function of cylindrical radius r (see Eqs. (14) and (15)). 

In the text 
Fig. 3
Time evolution of the bow shock + cavity flow predicted by the analytic model for two choices of the wind velocity: v_{w} = 0 (leftframe) and v_{w} = 0.4v_{j} (right frame). The dark region is the empty part of the cavity (sweptup in the thin shell bow shock) and the brown region is the part of the cavity that has been refilled by the disk wind (this region being of course absent in the v_{w} = 0 model of the left frame). For both models, we show snapshots corresponding to t = 2L_{0}/V_{j}, 4L_{0}/V_{j}, and 8L_{0}/v_{j}), which result in positions z_{j} = 2L_{0}, 4L_{0}, and 8L_{0} for the working surface (see the labels on the top left of each frame). Both models have v_{0} = 0.2v_{j} 

In the text 
Fig. 4
Shape of the bow shock and the cavity (top), radial velocities v_{r} (center) and axial velocities v_{z} (bottom) for the two models discussed in the text. The solid curves show the velocities of the wellmixed material within the thin shell flow, and the dashed curves show the immediate postbow shock velocities. The dotted red line shows v_{z} = v_{w}. 

In the text 
Fig. 5
Time evolution of two successive bow shocks and the cavities predicted by the analytical model. The first bow shock is ejected at t = 0, the second shock is formed at time t = τ_{j} when the first bow shock is at z = v_{j}τ_{j} taken here as 6L_{0} (top). At a time t = τ_{j} + 2L_{0}/v_{j}, the second bow shock is still propagating in the undisturbed disk wind material (center). At a time t_{c} the second bow shock catches up the emptied cavity of the first bow shock at its vertex (bottom). 

In the text 
Fig. 6
Maps of density and pressure for the reference noDW simulation with v_{w} = 0 (left) and the DW simulation with v_{w} = 0.4v_{j}(right) at a time t = 48 yr. Color scales on top are in g cm^{−3} for density and in dyn cm^{−2} for pressure. White contours show the locus of 50% mixing ratio between jet and diskwind/ambient material. The cyan curve shows a fit (to the numerical results) by the analytic shell shape in Eq. (11), with L_{0} = 65 au (left) and 108 au (right). The cyan dot indicates the maximum radius of the shell, the cyan asterix indicates the predicted vertex of the empty conical cavity left behind the shell, and the cyan dashdotted line is the analytical predicted boundary between the emptied cavity and the region refilled from below by fresh disk wind (see Fig. 2 and Eq. (15)). 

In the text 
Fig. 7
Transverse cut across the flow at the IWS location (z = 993 au) in the noDW timeframe shown in Fig. 6. This cut shows the radial velocity as a function of distance from the jet axis in solid line. We also plot the radial velocity weighted by the abundance of the jet tracer with a dashed line. The radial velocity first grows outwards, reaches a maximum velocity of ≈ 14 km s^{−1} at a radius of ∼ 25 au (somewhat larger than the 20 au initial jet radius), and then remains with values > 10 km s^{−1} until it drops sharply to 0 at r ∼ 50 au. The velocity maximum at r ∼ 25 au corresponds to the shock against the jet material. The second maximum at r ∼ 50 au is the shock that propagates in the diskwind and the zero radial velocity material at larger radii is the undisturbed disk wind. 

In the text 
Fig. 8
Zoom of the leading IWS of the simulation with a surrounding disk wind at time t = 48 yr. Left: density stratification (with the logarithmic color scale given by the top bar in g cm^{−3}), center: radial velocity (with the linear scale of top bar in km s͔–1) and right: axial velocity structure (with the linear scale of the top bar in km s^{−1}). The white contours show the surfaces of 50 % (solid line), 10 %, 0.1 % and 0.001 % (outer contour) jet material mixing fractions. The cyan asterisk is the predicted vertex of the cavity, the cyan dashdotted line in the predicted boundary of the cavity, and the cyan curve is the fitted shape of the bow shock. 

In the text 
Fig. 9
Longitudinal positionvelocity (PV) diagrams for the no DW simulation (top) and the DW simulation (bottom). From left to right: v_{r} for the jet material only, v_{r} for the surrounding material only, v_{z} for the jet material only, vz for the surrounding material only, and density stratification. The ordinate of all frames is position along the outflow axis (in au). The color scale in the PVs is scaled by volume × cube of pressure so as to be maximum (in red and yellow shades) for shocked material in the shell, while the color scale for density is in g cm^{−3}. Blue curves are predicted velocities in the full mixing hypothesis (Eqs. (16) and (17)), while magenta curves are the predicted immediatepost shock velocities (Eqs. (18) and (19)). 

In the text 
Fig. 10
Density maps for the noDW reference simulation (three frames on the left) and the DW simulation (three frames on the right) at t = 71, 119, and 167 yr. The white contours indicate the surface of 50% (solid line) and 90% (dashed line) jet material mixing fractions. The black lines in the diskwind simulation show a cone of α = 11° opening halfangle, which circumscribes the boundary the region perturbed by the jet and its IWS. The black dashed lines show the predicted trajectory of the edge of the bow shock (see Eq. (14)). The density color scale is given by the right bar (in g cm^{−3}). 

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