<figure><img src="image_80.png" /><figcaption>**Fig. 3.** Influence of seven free parameters on the geometry and kinematics of shells driven by a pulsed conical jet in a singular flattened core, at an age of *t* = 700 yr. Each row corresponds to a different free parameter (as labeled in the left margin) and compares maps of number density, *n*<sub>H</sub> (left), radial velocity, *v*<sub>R</sub> (middle), and axial velocity, *v*<sub>z</sub> (right), for the reference model (left half of each map) and the modified model (right half of each map). The corresponding modified parameter values are marked at the bottom of the *n*<sub>H</sub> panel and model names at the bottom of the *v*<sub>z</sub> panel. Velocity colorbars are cropped to 0–50 km s<sup>-1</sup> for better visualization of the range detected in CO outflows. The reference model in this figure has a high mass-flux of 6 × 10<sup>-6</sup> M<sub>☉</sub> yr<sup>-1</sup> (see Table 2 for full list of model parameters). The main effect on the cavity shape is seen when varying the core density.</figcaption></figure>
<figure><img src="image_81.png" /><figcaption>**Fig. 3.** continued.</figcaption></figure>

changes, but with a 100 times smaller reference jet density. This
leads to $\dot{M} = 6 \times 10^{-8} \text{ M}_\odot \text{ yr}^{-1}$ in the corresponding reference
model (M_REF) that is typical of more evolved Class 1 jets. The
parameters of this second sequence of models (referred to as the
"medium-density" case) are summarized in Table A.1 and the
results are summarized graphically in Fig. A.1.

4.3. Results

Figure 3 shows the differences in shape and kinematics of the shells between the reference and modified models in the high-density case, by comparing their respective maps of $n_H$ (H nucleus number density), $v_R$, and $v_z$ at the same age of $t = 700$ yr. The $v_R$ maps show the lateral expansion induced by successive bowshocks, while the $v_z$ maps highlight the shear-like velocity gradient that develops along the main shell walls.

For an easier comparison between models, we use the density
maps in Fig. 3 to measure the full width $W(z)$ of the main shell
as a function of altitude $z$. We thus obtain Fig. 4 that overplots

the resulting shell shapes of each model. Finally, from the shell
width $W_{800}$ at an altitude $z = 800$ au, we derive the full open-
ing angle near the base, defined following Dutta et al. (2020) as
$\alpha_{800} = 2 \arctan(W_{800}/[2 \times 800 \text{ au}])$. Values of $\alpha_{800}$, $W_{800}$, and
$W_{12700}$ (the full width at the top of the computational box) are
summarized in the last three columns of Table 2.

The most impactful effect on both the morphology and
kinematics of the main shell is obtained here when decreas-
ing the core base density, $\rho_{a0}$, by a factor 100. In this case,
Figs. 3a and 4 clearly show that the leading shell opens twice
wider; $v_R$ reaches higher values, which cover broader areas
within the nested shells. The layer of strong $v_z$-gradient (between
0 and 50 km s$^{-1}$) along the main shell surface also becomes
thicker.

Increasing the jet opening angle (Figs. 3b and 4) has a more
moderate effect, with a maximum ~30% increase of the main
shell full width in comparison with the reference model in the
high-density case. Increasing the jet base density by a factor 100
also has a moderate influence (see Figs. 3c and 4).
<figure><img src="image_82.png" /><figcaption>**Fig. 4.** Full shell width $W(z)$ at $t = 700$ yr as a function of altitude, $z$, for the high-density jet models in Fig. 3. Colored curves have one parameter varied from the reference run, from among the ambient core base density $\rho_{a0}$, jet semi-opening angle $\theta_j$, initial jet base density $\rho_{j0}$, semi-amplitude of variability $\Delta V$, period $P$, variability profile (sawtooth instead of sinusoidal), or jet radius $R_j$. Model parameter values are listed in Table 2. The covered measurements range from $z \sim 800$ au to $z \sim 12$ 700 au.</figcaption></figure>

Finally, the shape and kinematics of the main shell are little affected when semi-amplitude $\Delta V$ (Fig. 3d), variability period $P$ (Fig. 3e) are increased, when a sawtooth profile $h(t)$ is superimposed (Fig. 3f), or when the jet radius $R_j$ (Fig. 3g) is decreased. However, we find that in comparison with the sinusoidal case, a sawtooth-like profile of variability leads to smoother shells borders and less unstable bowshocks and IWS, as well as reduced instabilities in general (cf. Fig. 3f).

We obtain the same behaviors for a reference jet density that is 100 times smaller (see Fig. A.1). Except for the ambient density, parameter changes only have modest effects on the shell width.

## 5. Millenia-long simulations of a conical pulsed jet in a stratified core

In this section, we investigate the long-term evolution of the shells driven by a pulsed conical jet through a flattened singular core. To do so, for the first time we present this type of simulation reaching up to 10 000 yr. Predicted sizes, position-velocity (PV) diagrams, and mass-velocity distributions are presented and qualitatively compared with the typical behavior in Atacama Large Millimeter/sub-millimeter Array (ALMA) observations of outflows.

### 5.1. Setup

In Sect. 4, we noted that a sawtooth jet velocity variability profile minimizes the development of instabilities without changing the overall shell shape and kinematics. We thus adopted such a sawtooth profile here, since it allows us to reach the desired long timescales at a more reasonable CPU cost. For consistency, we adopted the same jet and ambient parameters as in model M_SAWT in Table A.1, except for the variability period. We considered here a slightly longer value $P = 300$ yr, as inferred

from CO observations of both the HH46-47 outflow by Zhang et al. (2019) and the CARMA-7 outflow by Plunkett et al. (2015), before the inclination correction.

To keep a non-prohibitive computational time, we also use here a more diffusing Total Variation Diminishing Lax Friedrichs numerical scheme (TVDLF) at the shocks zone. This decreases the instabilities at the highest resolution, without changing the overall shell structure. Finally, to follow the outflow expansion, the computational domain is expanded to $3.0 \times 10^{17}$ cm = 20 000 au in $R$ and $z$ ($n_R \times n_z = 240 \times 240$ cells at AMR level 1), keeping the same resolution as in Sect. 4 and Table 2.

### 5.2. Long-time maps of density, mixing fraction, and velocities

Figure 5 shows the distribution of number density, $n_H$, fraction of core-originated material, $f_a$, and orthoradial and axial velocities, $v_R$ and $v_z$, in the outflow at $t = 10 000$ yr. Figure 5a shows that after 10 000 yr, the outflow structure differs from that seen at earlier ages (e.g., in Fig. 3). The initially single shell has split into two separate shells: an “outer shell” tracing the forward shock propagating at low speed into the ambient medium and an inner “main shell” tracing the jet-ambient interface (roughly delimited by the $f_a = 50\%$ border), where the wings of successive bowshocks pile-up. Both shells display a roughly parabolic or conical shape out to 20 000 au. Between these two shells is a slow cocoon of ambient material. Inside the main shell, we can discern several distinct “nested shells,” tracing the last bowshocks recently created by jet variability. This geometry is reminiscent of the “spider-like” structure observed at the base of the B5-IRS1 outflow, where an inner parabolic shell, with a jet shock at its apex, is nested inside the wide-angle low-velocity outflow cavity (Zapata et al. 2014).

Moreover, the mixing map (Fig. 5b) shows that a fraction of material from the surrounding core can go past the shock and mix with jet-originated material inside the main cavity up to the wings of the nested bowshocks. Some core-originated material can even reach particularly overdense areas bordering the jet walls, with $n_H$ between $10^2$ and $10^5$ cm$^{-3}$. In the following, we assume that this core-originated material remains molecular for the purposes of computing the predicted synthetic emission diagrams.

Finally, Figs. 5c and d shows the radial and axial velocity maps, respectively. The new bowshocks encounter less resistance than those at earlier ages and develop broader wings, as the outflow cavity has been cleared up by tens of older bowshocks. Nevertheless, the characteristic kinematic pattern remains similar to early times, with enhanced-$v_R$ in the bowshocks wings and a shear-like gradient of axial velocity $v_z$ along the main shell at the jet-ambient interface. This velocity structure creates a characteristic “bell-shaped” pattern in transverse position-velocity diagrams, presented in Sect. 5.4.2.

### 5.3. Deceleration of the jet-driven shells

In contrast to a wind-driven shell in a flattened singular core, which expands at constant speed over time (Shu et al. 1991), a jet-driven shell is expected to decelerate. However, no analytical estimate of that deceleration exists in the case of a non-uniform ambient medium, making numerical simulations necessary.

Figure 6 shows the time evolution of the shell widths during the 10 000 yr-long simulation leading to the snapshot of
<figure><img src="image_83.png" /><figcaption>**Fig. 5.** Maps at $t = 10\,000$ yr of (a) hydrogen nuclei density $n_H$, (b) fraction of core-originated material $f_a$, (c) orthoradial velocity $v_R$, and (d) axial velocity $v_z$ of the shell driven by a conical pulsed jet in a flattened singular core. Parameters are identical to model M_SAWT in Table A.1 except for a longer variability period $P = 300$ yr. Velocities in the colorbars are cropped to the range typically detected in CO outflow observations, namely between 0 and 50 km s$^{-1}$. The white contour in the left-half of panels of (b), (c), and (d) shows $f_a = 50\%$.</figcaption></figure>

Fig. 5. We plot the full width as function of time at two different heights: $z = 800$ au (for comparison with Dutta et al. 2020), and $z = 20\,000$ au (top of the computational domain, after the jet head reaches it). We denote the corresponding shell widths as $W_{800}$ and $W_{20\,000}$.

Figure 6 shows that the formation of two separate shells (denoted as the outer and main shells in Fig. 5) occurs around $t = 1000$ yr. Both shells are seen to decelerate at late times. Deceleration is stronger at lower altitudes, where the ambient core is denser. At $z = 800$ au, the main shell stops expanding after 6000 yr and reaches a final width $W_{800} \approx 2 \times 10^{16}$ cm = 1500 au. The corresponding final opening angle is $\alpha_{800} = 86^\circ$. This behavior is consistent with observations suggesting that the base opening angle of CO outflows stops increasing after $t \approx 8000$ yr, with a final value (uncorrected for inclination) $\alpha_{\text{obs}} \approx 90^\circ-100^\circ$ (Velusamy et al. 2014).

The widths $W_{800}$ of the main shell at $z = 800$ au also fit very well within the observed range of flow widths at the same projected height (indicated by grey bands in Fig. 6), measured by Dutta et al. (2020) in a sample of 22 CO outflows in Orion. On large core scales of $z = 20\,000$ au = 0.1 pc, the main shell reaches a width $W_{20\,000} = 15\,000$ au at an age of $10^4$ yr. This is similar to the observed CO outflow width at the same (deprojected) height in HH46-47 (Zhang et al. 2016), indicated in blue in Fig. 6. On intermediate scales of $z = 8000$ au, the main shell width at an age of $10^4$ yr in our simulation is 9000 au (see Fig. 5). This is also in good agreement with cavity widths observed at the same (projected) distance in scattered light, lying in the range 1100–8500 au in 75% of cases (cf. semi-opening angles reported in Habel et al. 2021). Therefore, a jet driven into a flattened singular core seems able to reproduce typical observed outflow widths on both small and large scales for realistic long ages of $\ge 10\,000$ yr.

<figure><img src="image_84.png" /><figcaption>**Fig. 6.** Temporal evolution in the simulation of Fig. 5 of the full width of the main shell (jet-ambient interface, crosses) and outer shell (forward shock, circles) measured at altitudes $z = 800$ au (black) and $z = 20\,000$ au (red). Analytical predictions for the equivalent “modified X-wind” model of Lee et al. (2001) are shown as dotted black and red lines of slope $t^{1/2}$ (see text). The full range of CO outflow widths observed at $z = 800$ au by Dutta et al. (2020) is indicated by the light grey band (with second and third quartiles in darker grey, and median as a thin green line). The full width of the HH46-47 outflow at $z = 20\,000$ au, from Zhang et al. (2016), is shown in blue.</figcaption></figure>
For comparison purposes, the dotted lines in Fig. 6 show the predicted evolution of shell widths at $z = 800$ and $20\,000$ au for the equivalent wide-angle “modified X-wind” model of Lee et al. (2001), which has the same mean axial velocity and total mass-flux as our jet and propagates in the same flattened singular core. In this model, the wind velocity drops away from the axis as $V_w = V_{w0} \cos \theta$, and the base of the shell is a parabola, $z = R^2/R_0$, expanding self-similarly over time as $R_0 = V_0 t$, where $V_0 = V_{w0} \sqrt{\eta}$ with $\eta$ the (fixed) ratio of wind-to-ambient density in the equator. The shell width at any fixed height $z$ is then given by $W_z = 2R = 2\sqrt{z}V_0 t$. In the equivalent wide-angle wind model considered here, $V_{w0} = 120 \text{ km s}^{-1}$ and $\eta = 2.2 \times 10^{-3}$ (Lee et al. 2001), yielding $V_0 = 5.6 \text{ km s}^{-1}$. Observed outflow shells fitted with the same analytical model<fnref n="3" /> have similar values of $V_0 \approx 1.5–9 \text{ km s}^{-1}$ (Lee et al. 2000; Zhang et al. 2019).

The black dotted line in Fig. 6 shows that while the wide-angle wind-driven shell predicts comparable widths to the jet-driven shell at early times $\le 300$ yr (as noted in Sect. 3), its self-similar expansion without deceleration exceeds observed CO outflow widths at $z = 800$ au (grey band in Fig. 6) after only 2000 yr, and the discrepancy increases over time as $t^{1/2}$. Recent simulations of X-wind-driven shells including magnetic fields in both the wind and the ambient medium predict similar shell widths at $500–1000$ yr as this simple analytical model and confirm that the shell expands at constant rate in a self-similar fashion, $R \propto t$ (Shang et al. 2020). We note that this ongoing expansion is tightly linked to the adopted $1/r^2$ ambient density distribution (Shu et al. 1991). Steady wind-driven shells could form if the core has a shallower density gradient and if mixing is inefficient at the wind-core interface (Smith 1986; Liang et al. 2020).

## 5.4. Synthetic predictions

### 5.4.1. Longitudinal position-velocity diagrams

Figure 7 presents longitudinal position-velocity (PV) diagrams cut along the jet axis from our long-term simulation at $t = 10\,000$ yr, assuming various inclination angles, $i$, of the blueshifted lobe from the line of sight.

Most importantly, this figure shows that at ages typical of Class 0 outflows, negative velocities all vanish below an inclination angle of $i < 70^\circ$ from the line of sight. In particular, there is no more observable blue/red overlap at $i = 60^\circ$, unlike what has been predicted for jet-driven shells in uniform media (Lee et al. 2001).

This behavior is induced by the transverse deceleration of the main shell on late timescales, as discussed in Sect. 5.3. The slow expansion of the main shell restricts the sideways expansion of the bowshock wings propagating inside it, forcing them to adopt more forward-directed velocities. This may be seen in Fig. 8, where we plot the direction of velocity vectors in the outflow. Along the main shell, where bowshock wings pile up, velocity vectors are almost parallel to the shell walls. This produces much less blue-red overlap in the projection than is otherwise expected for a bowshock in a uniform medium. In the outer shell, which traces the forward shock expanding into the ambient medium, the velocity vectors are widely open (perpendicular to the shell) but the local expansion velocity is so low that the emission falls near rest velocity once it is projected.

Challenges to the notion of jet-driven shells producing too much blue-red overlap over a broad range of view angles thus

no longer appears justified when realistic long ages and ambient core stratifications are considered. A longitudinal PV cut extending up to $1.5 \times 10^{17}$ cm from the source was recently obtained in the CARMA-7 outflow with ALMA (Plunkett et al. 2015), with both high resolution and high sensitivity.

A qualitative comparison with our predictions is presented in Fig. 9. Since the PV of CARMA-7 presents significant blue-red overlap, it is believed to be close to the plane of the sky, hence, we considered an inclination of $i = 85^\circ$ to the line of sight. The observed PV cut, reproduced in Fig. 9a, shows a striking quasi-periodic series of velocity peaks. The differences in dynamical timescales between successive identified velocity peaks yield an apparent period of variability in CARMA-7 of $\Delta\tau_{dyn} \approx 300$ yr (see Plunkett et al. 2015). However, for a quasi edge-on inclination, our synthetic PV diagrams with $P = 300$ yr in Fig. 7 predict a much wider knot spacing than observed in Fig. 9a. As noted by Plunkett et al. (2015), the value of $\Delta\tau_{dyn}$ may need to be corrected for projection effects by a typical factor 10. We thus present in Fig. 9 a model with jet variability period that is ten times shorter, namely, $P \approx 30$ yr.

Two synthetic PV diagrams are presented in Fig. 9: one including all material from jet and core (Fig. 9b) and the other including only core-originated material (Fig. 9c). Figure 9b is dominated by periodic sawtooth structures tracing the time variable jet and the sideways ejected material from its internal working surfaces (IWS). Figure 9c is dominated by structures with apparent “Hubble-law” acceleration, tracing ambient gas swept-up in the successive nested bowshock wings. We may remark that the jet and IWS remain faintly visible, through ambient gas dynamically entrained along the jet borders (cf. mixing map in Fig. 5b).

We find that the predicted structures in Fig. 9c are qualitatively similar to what is observed in CO emission in CARMA-7 (Fig. 9a). Since our model is very simplified (e.g., it does not include any jet precession), this qualitative agreement can be considered as promising. It also confirms that a suitable inclination-correction is essential to estimate the true period of velocity variations in a quasi edge-on outflow.

### 5.4.2. Transverse PV diagrams

Figures 10a, c show synthetic transverse PV diagrams for our $10\,000$ yr-old simulation of an outflow with $P = 300$ yr (same simulation as in Fig. 5). We adopted an inclination angle $i = 55^\circ$ from the line of sight, the mean inclination of the HH46-47 jet determined by Hartigan et al. (2005). First, each diagram forms a characteristic “bell-like” shape, which peaks at high-velocity (the jet) and broadens smoothly down to rest velocity. This shape is a direct consequence of the deceleration of bowshock wings as they expand and interact inside the shell, which produces a velocity decreasing away from the jet axis. Second, while most of the mass is piled-up near rest velocity, and in the $v_z$ “shear layer” along the shell walls (responsible for the two “horns” along the edges of the bell), bright ellipses are also present at intermediate velocities. Those ellipses trace intersections of the line of sight with individual bowshocks expanding inside the main shell.

These predicted characteristics (bell-shape with nested ellipses) bear striking qualitative resemblance with observed transverse PV cuts at high resolution and sensitivity, recently obtained with ALMA across the HH46-47 outflow by Zhang et al. (2019) and shown for comparison in Figs. 10b, d. The agreement appears even better than with models of wide-angle wind-driven shells (Zhang et al. 2019), which tend to overpredict

<footnote n="3">³ Lee et al. (2000) use different notations, $C = 1/R_0$ and $v_0 = 1/t$.</footnote>
<figure><img src="image_85.png" /><figcaption>**Fig. 7.** Longitudinal position-velocity diagrams along the flow axis, inferred from the simulation in Fig. 5 at an age $t = 10\,000$ yr, for different inclination angles, $i$, from the line of sight, ranging from 90° (edge-on) to 30°.</figcaption></figure>

<figure><img src="image_86.png" /><figcaption>**Fig. 8.** Orientation of the velocity vectors for the $t = 10\,000$ yr-long simulation in Fig. 5, in the reference frame of the ambient core. Arrows are shown for velocities above 0.45 km s⁻¹, and are not scaled with the velocity modulus. Density is displayed in the background color image.</figcaption></figure>
<figure><img src="image_87.png" /><figcaption>**Fig. 9.** Observed and modeled longitudinal position-velocity (PV) diagrams in CARMA-7. (a) Longitudinal PV diagram of CARMA-7 inside a 374 au-wide cut along the outflow axis, as observed in $^{12}$CO(2–1) by Plunkett et al. (2015). (b) Synthetic longitudinal PV diagram along a 400 au-wide cut, at $t = 500$ yr, and inclination $i = 85^\circ$ from the line of sight. Parameters are the same as for the model shown in Fig. 5 except for a shorter jet variability period $P = 30$ yr. (c) Same as (b) but including only the core-originated material.</figcaption></figure>

the outflow width at high projected velocities (see red ellipses on top of the observed PVs in Figs. 10b,d).

### 5.4.3. Mass-velocity distribution

As an additional diagnostic tool, Fig. 11 shows the mass-velocity (MV) distribution of entrained core material derived from the long term simulation in Fig. 5 at an age $t = 10\,000$ yr, assuming an inclination of $i = 55^\circ$ from the line of sight. This MV distribution is confronted to that derived by Zhang et al. (2016) from CO observations of the HH46-47 outflow, which is viewed at the same inclination angle.

Two observed MV distributions were computed by Zhang et al. (2016) and are plotted in Fig. 11: the first one (open symbols) was derived from $^{12}$CO without correction for

optical depth. The shallow power-law slope $\gamma \approx -2$ is very typical of the slopes reported previously in other CO outflows, before optical-depth correction. The second observed MV distribution (filled symbols) was obtained after applying a velocity-dependent optical-depth correction to $^{12}$CO and $^{13}$CO data, and adding (optically thin) $^{18}$O data at low velocities. Because optical depth increases at lower flow velocity, the corrected distribution exhibits a much steeper slope $\gamma \approx -3.4$ (in log-log). Similar steeper slopes have been obtained in other outflows after applying a velocity-dependent optical-depth correction, such as in B5-IRS1 (Yu et al. 1999).

Figure 11 shows a remarkably good qualitative agreement between the simulated and the observed MV distribution in HH46-47 after correction for optical depth, which should be closest to the true mass distribution in the outflow. This
<figure><img src="image_88.png" /><figcaption>**Fig. 10.** Modeled and observed transverse position-velocity (PV) diagrams in HH46-47. Left panels: transverse position-velocity diagrams perpendicular to the outflow axis, inferred from the jet-driven shells simulation in Fig. 5 and Sect. 5.2, at *t* = 10 000 yr inside 400 au-wide cuts at projected heights of (a) 6.6 × 10¹⁶ cm and (c) 8.9 × 10¹⁶ cm, with an outflow inclined by an angle *i* = 55° from the line of sight. Here only the core-originated material mass contribution is shown. Right panels: ¹²CO(2–1) emission of HH46/47 along two 450 au-wide cuts at projected heights of (b) 8.1 × 10¹⁶ cm and (d) 1.0 × 10¹⁷ cm from the central source, with red ellipses showing best-fit models by wide-angle wind-driven shells. Adapted from Zhang et al. (2019).</figcaption></figure>

agreement further reinforces the attractiveness of the jet-driven
shell model at realistic long ages and with ambient stratification
as a possible origin for CO outflows.

5.5. *Model approximations*

Here, we briefly discuss the main physics left out in our simula-
tions and how it might affect the results. First, we considered
an equilibrium atomic cooling curve, without following the
out-of-equilibrium ionization and cooling behind shock fronts.
Therefore, we cannot compute realistic synthetic emission maps
in atomic and ionic lines, whose flux is proportional to the elec-
tron density. This approximation is sufficient, however, for our
purpose of determining the overall shape, dynamics, and mass
distribution in the dense shells of cooled post-shock gas. We also
neglected chemistry, noting that CO dissociation only affects the
observed mass distribution above 20 km s⁻¹ (Downes & Cabrit
2003; Moraghan et al. 2008) and our comparisons with obser-
vations are made at lower velocities. Such approximations allow
us to carry out for the first time, in a cost-effective way, long-
term simulations up to 10⁴ yr that can be compared with actual
outflows observed with ALMA.

Overall, we neglected the infall motions. This allowed
us to properly compare with simulations of outflows driven
by X-winds into the same stratified singular cores (Lee et al.
2001; Shang et al. 2020), where infall is likewise neglected.
In such cores, infall propagates inside-out at the ambient
sound speed, *a*, and generates a shallower density distribution,
1/*r*¹.⁵, inside the sonic radius, *r*<sub>inf</sub> = *at* (Terebey et al. 1984).
Contrary to wind-driven flows where it can lead to a steady
shell (Liang et al. 2020), a 1/*r*¹.⁵ slope does not greatly reduce
the width of jet-driven flows, as compared to 1/*r*² (Moraghan
et al. 2008). The main effect of infall would thus be to add
shear, entrainment, and extra ram pressure at the base of the
outer shell (cf. Liang et al. 2020, for the wind-driven case).
For our simulation parameters, infall would occur inside
r<sub>inf</sub> ≈ 1000 au (a/0.5 km s<sup>-1</sup>)×(t/10<sup>4</sup> yr) and, hence, it would
affect only a small fraction of our full computational domain,
extending up to 20 000 au.

Our simulation parameters do not explore the “long-period”
regime where the ambient core would have time to partly refill
the cavity in between jet outbursts. This interesting situation
would occur if major outbursts happen only every few 10⁴ yr, for
instance, due to tidal interaction in wide eccentric binaries. The
<figure><img src="image_89.png" /><figcaption>**Fig. 11.** Simulated mass-velocity (MV) distribution of core-originated material at *t* = 10 000 yr from the jet-driven shell simulation in Figs. 5 and 10, inclined by *i* = 55° from the line of sight (connected black symbols). Blue and red symbols show observed MV distributions in the blueshifted and redshifted lobes of the HH46-47 outflow, as derived by Zhang et al. (2016). Open squares assume optically thin CO emission, while filled circles include a velocity-dependent correction for optical depth (see text). The dashed lines show power-laws of slope *γ* = −3.4 (red) and *γ* = −2.7 (blue).</figcaption></figure>

observed spacing of jet knots in Class 0/1 sources indicates much shorter variability timescales, however, with multiple periodic modes of decades, centuries, and thousands of years (Raga et al. 2012; Lee 2020). In addition, a recent study of infrared variability towards embedded protostars with *Spitzer/IRAC* presented the conclusion that Class 0 protostars undergo a major burst on average every 438 yr, with a 95% confidence interval of 161–1884 yr (Zakri et al. 2022). The long-period regime thus seems quite rare among the youngest, embedded Class 0/1 objects that drive observable CO outflows.

Finally, we have neglected magnetic fields, both in the jet and in the ambient medium. Magnetic pressure in the jet would act to reduce the postshock compression and increase the cooling length, by typically an order of magnitude (Hartigan et al. 1994). This effect is observed directly in resolved internal shocks of stellar jets (Hartigan & Wright 2015). High-resolution radiative numerical simulations of pulsed magnetized jets show that the reduced cooling tends to broaden the nested bowshock shells, but the effect appears rather modest (de Colle & Raga 2006), justifying its omission here.

Conversely, a magnetic field in the ambient medium would tend to resist against the sideways shell expansion. A significant reduction in shell width compared to the purely hydrodynamic case requires, however, strong fields near equipartition (i.e., a ratio of thermal to magnetic pressure $\beta \approx 1$). This is demonstrated, for example, in the simulations of Shang et al. (2020) of X-wind-driven shells into cores of varying degrees of flattening and magnetization<fnref n="4" />. For the standard moderate core flattening,

$n=2$, adopted in the present work, the maximum reduction in shell width due to ambient magnetic field is only 25%, therefore, our predicted cavity shapes and dynamics should remain valid. We note that an added complexity at higher core magnetization, in the general case of non-zero rotation, would be the probable launching of a massive slow MHD wind from the Keplerian disk formed around the protostar (see e.g., Lesur 2021, and references therein). The interaction of an inner pulsed jet (or wide-angle X-wind) with an outer disk wind could significantly affect the formation of outflow cavities, but the long-term evolution has only been examined in the hydrodynamical case so far, to our knowledge (Tabone et al. 2018).

Concerning observational predictions, we focused here on signatures of the shell shapes, kinematics, and mass distribution at velocities below 20 km s⁻¹, which are well traced by low-excitation CO emission observable with ALMA. Our simulations may also be used to assess the detectability of warmer jet-originated material inside the cavity volume. From Fig. 5 we estimate that the nested shells (driven by each of the jet pulses) generate typical shock speeds of 30 km s⁻¹ inside the cavity. The pre-shock density inside the cavity is ≈ 100 cm⁻³ at *z* ≤ 3000 au, and drops off at higher altitudes. Using the atomic shock model grid of Hartigan et al. (1994), the maximum predicted surface brightness in [S II] for solar abundances is then ≈5 × 10⁻⁶ erg s⁻¹ cm⁻² sr⁻¹. Such extended [S II] emission inside outflow cavities would be worthwhile to search for, but it might be difficult to isolate against the much brighter axial jet and scattered light from the cavity walls, especially since young protostars with powerful CO outflows are often located in regions of high optical extinction.

When the jet is dense enough to be partly molecular, another possible tracer of nested shells inside the main outflow cavity is H₂ ro-vibrational emission, the most spectacular example so far being the Class 0 outflow HH 212. Each large H₂ bowshock, produced by a major jet pulse, is seen to connect to a separate CO shell at the flow base, nested inside the main cavity (see Fig. 2 in Lee et al. 2015), a morphology consistent with our predictions for a pulsed jet-driven outflow. Modeling the H₂ line emission could be helpful to further discriminate between this scenario and the pulsed wide-angle wind model proposed by Zhang et al. (2019).

## 6. Conclusions

We confirm in this paper that the swept-up shell driven by a cylindrical pulsed jet opens much wider inside an ambient core with steeply decreasing density, than within a homogeneous core (which was the configuration commonly adopted until now for analytical models and most simulations of jet-driven bowshocks). At early times of a few hundred years, the jet-driven shell can open as wide as with a wide-angle X-wind, when considering the same flattened singular ambient core and the same injected mass-flux and velocity variability.

Then, we investigated the impact that several parameters in our model have on the general morphology, opening angle, and kinematics of jet-driven shells. The parameter which is by far the most impactful is the ambient core density, followed by the jet density and jet opening angle in a less impactful fashion. Within the range of values explored in this paper, the other parameters do not significantly affect these diagnostics.

<footnote>is seen only in simulations with $n \ge 4$ and $\alpha_b = 1$, corresponding to midplane values of $\beta \le 3.5$ that approach equipartition.</footnote>

<footnote n="4">⁴ The ratio of thermal to magnetic pressure in the core models of Shang et al. (2020) can be recovered as $\beta = 2\alpha_b^{-2} (v_A/a)^{-2}$ where $\alpha_b = [0, 0.1, 1]$ is their scaling parameter and $(v_A/a)$ is the ratio of Alfvénic speed to sound speed in the magnetostatic solution of Li & Shu 1996 (Eq. (69)), which increases with the core flattening parameter $n$. A significant reduction in shell width compared to the hydrodynamical case</footnote>
Finally, running a representative simulation up to $t = 10\,000$ yr reveals drastic changes in the long-term. After several hundred years, the initial shell splits into a slow parabolic outer shell fully made of core-original material (tracing the forward shock) and an inner, faster main shell tracing the jet-ambient interface, which stops expanding at the base after 8000 yr, unlike wide-angle wind-driven shells. This main shell encompasses a mixed-material cavity inside which successive bowshocks driven by the pulsed jet expand and slow down by interacting with previous ones, producing a strongly sheared velocity field parallel to the (roughly conical) main shell walls. Both the morphology and the velocity fields are very different from analytical predictions of ballistic jet bowshock models in uniform media (Ostriker et al. 2001).

The long-term simulation of our basic model shows none of the caveats of steady jet bowshocks in uniform media (excessive length to width ratio, excessive blue-red overlap, and overly low speeds at large shell widths). On the contrary, it shows very promising similarities with the most recent observations of Class 0 outflows observed at high resolution with ALMA, in terms of predicted shell widths, full opening angle (≈90° after $10^4$ yr), longitudinal and transverse position-velocity cuts, and mass-velocity distribution (Dutta et al. 2020; Plunkett et al. 2015; Zhang et al. 2016, 2019). Some comparisons even show more resemblance to observations than the widely used wide angle “modified X-wind” model of Lee et al. (2001). This is the case for the moderate outflow widths at 800 au, the characteristic “bell-shape” of transverse PV cuts, and the steep mass-velocity relation after CO opacity correction.

More generally, this paper shows that a realistic modeling of the surrounding core density stratification, as well as long integration times of at least $10^4$ yr, are both essential to reliably predict the properties of outflows driven by a pulsed jet, and to confront them with observations. In the future, we plan to extend our simulations to include chemistry, magnetic field, rotation, and precession to model specific observed protostellar outflows in more detail.

*Acknowledgements.* We are grateful to the anonymous referee for insightful suggestions that improved the paper, and to Doug Johnstone for supportive discussions and thoughtful comments. We acknowledge support by the Programme National Physique et Chimie du Milieu Interstellaire (PCMI) of CNRS/INSU with INC/INP co-funded by CEA and CNES. Part of the computations were carried out on the OCCIGEN cluster at CINES (https://www.cines.fr/) in Montpellier (project named lut6216, allocation A0090406842 and A0100412483) and on the MesoPSL cluster (http://www.mesopsl.fr/) of PSL University at Observatoire de Paris.

## References

Arce, H. G., Shepherd, D., Gueth, F., et al. 2007, in Protostars and Planets V, eds. B. Reipurth, D. Jewitt, & K. Keil, 245

Arce, H. G., Borkin, M. A., Goodman, A. A., Pineda, J. E., & Halle, M. W. 2010, ApJ, 715, 1170

Béthune, W., Lesur, G., & Ferreira, J. 2017, A&A, 600, A75

Cabrit, S., Raga, A., & Gueth, F. 1997, in Herbig-Haro Flows and the Birth of Stars, eds. B. Reipurth, & C. Bertout, 182, 163

Dalgarno, A., & McCray, R. A. 1972, ARA&A, 10, 375

de Colle, F., & Raga, A. C. 2006, A&A, 449, 1061

Dougados, C., Cabrit, S., Lavalley, C., & Ménard, F. 2000, A&A, 357, A61

Downes, T. P., & Cabrit, S. 2003, A&A, 403, 135

Dutta, S., Lee, C.-F., Liu, T., et al. 2020, ApJS, 251, 20

Federrath, C., Schrön, M., Banerjee, R., & Klessen, R. S. 2014, ApJ, 790, 128

Ferreira, J., & Deguiran, R. 2013, High Energy Density Phys., 9, 67

Frank, A., Ray, T. P., Cabrit, S., et al. 2014, in Protostars and Planets VI, eds. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning, 451

Habel, N. M., Megeath, S. T., Booker, J. J., et al. 2021, ApJ, 911, 153

Hartigan, P., & Wright, A. 2015, ApJ, 811, 12

Hartigan, P., Morse, J. A., & Raymond, J. 1994, ApJ, 436, 125

Hartigan, P., Heathcote, S., Morse, J. A., Reipurth, B., & Bally, J. 2005, AJ, 130, 2197

Ireland, L. G., Zanni, C., Matt, S. P., & Pantolmos, G. 2021, ApJ, 906, 4

Kaastra, J. S., & Mewe, R. 2000, in Atomic Data Needs for X-ray Astronomy, eds. M. A. Bautista, T. R. Kallman, & A. K. Pradhan, 161

Keppens, R., Teunissen, J., Xia, C., & Porth, O. 2021, Comput. Math. Applic., 81, 316

Krumholz, M. R., & Federrath, C. 2019, Front. Astron. Space Sci., 6, 7

Lee, C.-F. 2020, A&ARv, 28, 1

Lee, C.-F., Mundy, L. G., Reipurth, B., Ostriker, E. C., & Stone, J. M. 2000, ApJ, 542, 925

Lee, C.-F., Stone, J. M., Ostriker, E. C., & Mundy, L. G. 2001, ApJ, 557, 429

Lee, C.-F., Hirano, N., Zhang, Q., et al. 2015, ApJ, 805, 186

Lesur, G. R. J. 2021, A&A, 650, A35

Li, S. 2005, J. Comput. Phys., 203, 344

Li, Z.-Y., & Shu, F. H. 1996, ApJ, 472, 211

Liang, L., Johnstone, D., Cabrit, S., & Kristensen, L. E. 2020, ApJ, 900, 15

MacDonald, J., & Bailey, M. E. 1981, MNRAS, 197, 995

Masson, C. R., & Chernin, L. M. 1993, ApJ, 414, 230

Matzner, C. D., & McKee, C. F. 1999, ApJ, 526, L109

Moraghan, A., Smith, M. D., & Rosen, A. 2008, MNRAS, 386, 2091

Motte, F., & André, P. 2001, A&A, 365, 440

Ostriker, E. C., Lee, C.-F., Stone, J. M., & Mundy, L. G. 2001, ApJ, 557, 443

Plunkett, A. L., Arce, H. G., Mardones, D., et al. 2015, Nature, 527, 70

Porth, O., Xia, C., Hendrix, T., Moschou, S. P., & Keppens, R. 2014, ApJS, 214, 4

Raga, A., & Cabrit, S. 1993, A&A, 278, 267

Raga, A. C., Canto, J., Binette, L., & Calvet, N. 1990, ApJ, 364, 601

Raga, A. C., Rodríguez-González, A., Noriega-Crespo, A., & Esquivel, A. 2012, ApJ, 744, L12

Rohde, P. F., Walch, S., Seifried, D., et al. 2019, MNRAS, 483, 2563

Santiago-García, J., Tafalla, M., Johnstone, D., & Bachiller, R. 2009, A&A, 495, 169

Schure, K. M., Kosenko, D., Kaastra, J. S., Keppens, R., & Vink, J. 2009, A&A, 508, 751

Shang, H., Allen, A., Li, Z.-Y., et al. 2006, ApJ, 649, 845

Shang, H., Krasnopolsky, R., Liu, C.-F., & Wang, L.-Y. 2020, ApJ, 905, 116

Shu, F. H., Ruden, S. P., Lada, C. J., & Lizano, S. 1991, ApJ, 370, L31

Shu, F. H., Najita, J., Ostriker, E. C., & Shang, H. 1995, ApJ, 455, L155

Smith, M. D. 1986, MNRAS, 223, 57

Suttner, G., Smith, M. D., Yorke, H. W., & Zinnecker, H. 1997, A&A, 318, 595

Tabone, B., Raga, A., Cabrit, S., & Pineau des Forêts, G. 2018, A&A, 614, A119

Terebey, S., Shu, F. H., & Cassen, P. 1984, ApJ, 286, 529

van der Marel, R. P., Kallivayalil, N., & Besla, G. 2009, in The Magellanic System: Stars, Gas, and Galaxies, eds. J. T. Van Loon, & J. M. Oliveira, 256, 81

van Marle, A. J., & Keppens, R. 2011, Computers and Fluids, 42, 44

Velusamy, T., & Langer, W. D. 1998, Nature, 392, 685

Velusamy, T., Langer, W. D., & Thompson, T. 2014, ApJ, 783, 6

Völker, R., Smith, M. D., Suttner, G., & Yorke, H. W. 1999, A&A, 343, 953

Yu, K. C., Billawala, Y., & Bally, J. 1999, AJ, 118, 2940

Zakri, W., Megeath, S. T., Fischer, W. J., et al. 2022, ApJ, 924, L23

Zapata, L. A., Arce, H. G., Brassfield, E., et al. 2014, MNRAS, 441, 3696

Zhang, Y., Arce, H. G., Mardones, D., et al. 2016, ApJ, 832, 158

Zhang, Y., Arce, H. G., Mardones, D., et al. 2019, ApJ, 883, 1
# Appendix A: Effect of variables with lower jet density

In this appendix, we redo the process of Sect. 4, but starting from a reference model with a jet base density that is two lower by orders of magnitude: $ρ_{j0} = 1.8 \times 10^{-19} \text{ g cm}^{-3}$. In this way, we can further probe the effects of the base jet density on the morphology and kinematics of the driven shells. We also test the robustness of the conclusions we got from Sect. 4 at other jet density values.

Thus, additionally to the reference model, we launch seven modified models, each having the same parameters as the reference model except for a single one. Those models are:

(i) M_DENSA: with an ambient density scaling $ρ_{a,0}$ of $1.6 \times 10^{-20} \text{ g cm}^{-3}$ instead of $1.6 \times 10^{-18} \text{ g cm}^{-3}$;

(ii) M_THETA: with a jet semi-opening angle $θ_j$ of 7° instead of 3°;

(iii) M_DENSJ: with a jet base density $ρ_{j0}$ 100 times larger than the reference case;

(iv) M_VARAMP: with a variability semi-amplitude $ΔV$ of 90 km/s instead of 60 km/s;

(v) M_PER: with a variability period $P = 300$ yr instead of 115 yr;

(vi) M_SAWT: with a sawtooth instead of sinusoidal velocity variability profile;

(vii) M_RAD: with an initial jet radius $R_j$ of 20 au instead of 50 au.

Figure A.1 plots and compares the inferred shell shapes at $t = 700$ yr for the reference and modified models. Table A.1 summarizes the parameters of the simulations, and lists the full opening angle at 800 au, $α_{800}$, and full shell widths at $z = 800$ au and 12700 au ($W_{800}$ and $W_{12700}$) measured at $t = 700$ yr in each model. The effect of parameters changes on the outflow structure are discussed in Sect. 4.3.

<figure><img src="image_90.png" /><figcaption>Fig. A.1: Full shell width $W(z)$ at $t = 700$ yr as a function of altitude $z$ for the models in Table A.1, with jet densities 100 times smaller than in Fig. 4. Colored curves have one parameter varied from the reference run, among the ambient core base density, $ρ_{a0}$, jet semi-opening angle, $θ_j$, initial jet base density, $ρ_{j0}$, semi-amplitude of variability, $ΔV$, period, $P$, type of variability profile (sawtooth instead of sinusoidal), or jet radius, $R_j$.</figcaption></figure>
<table><caption>Table A.1: Parameters of conical pulsed medium-density jet simulations with the resulting opening angles and full widths</caption><thead><tr><th colspan="2">Fixed parameters</th></tr></thead><tbody><tr><td>mean jet velocity</td><td>v<sub>0</sub> = 120 km/s</td></tr><tr><td>jet density variation</td><td>ρ<sub>j</sub>(t) = ρ<sub>j0</sub>[v<sub>0</sub>/v<sub>j</sub>(t)] × (R<sub>j</sub><sup>2</sup> + z<sub>0</sub><sup>2</sup>)(R<sub>j</sub><sup>2</sup> + [z + z<sub>0</sub>]<sup>2</sup>)<sup>-1</sup>, with z<sub>0</sub> = R<sub>j</sub>/ tan θ<sub>j</sub> (constant mass-flux)</td></tr><tr><td>jet temperature</td><td>T<sub>j</sub> = 100 K</td></tr><tr><td>ambient core temperature</td><td>T<sub>a</sub> = 100 K</td></tr><tr><td>core density profile</td><td>Flattened singular core ρ<sub>a</sub>(r) = ρ<sub>a0</sub> sin<sup>2</sup> θ r<sub>0</sub><sup>2</sup> (r<sup>2</sup>)<sup>-1</sup>, with r<sub>0</sub> = 2.5 × 10<sup>15</sup> cm</td></tr><tr><td>radiative cooling function</td><td>Λ(100 K ≤ T < 10<sup>4</sup> K) from Dalgarno & McCray (1972)<br>Λ(T ≥ 10<sup>4</sup> K) from Schure et al. (2009)</td></tr><tr><td>simulation domain</td><td>(R, z) = (7.0 × 10<sup>16</sup> cm, 1.9 × 10<sup>17</sup> cm) = (4679 au, 12 700 au)</td></tr><tr><td>number of cells</td><td>n<sub>R</sub> × n<sub>z</sub> = 56 × 152 at AMR level 1</td></tr><tr><td>maximal AMR level</td><td>5</td></tr><tr><td>maximum resolution</td><td>ΔR = Δz = 7.8 × 10<sup>13</sup> cm = 5.2 au at AMR level 5</td></tr><tr><td>snapshot age</td><td>700 yr</td></tr></tbody></table><table><thead><tr><th>Parameter</th><th>Reference Model<br>M_REF</th><th>Modified<br>parameter<sup>a</sup></th><th>Modified<br>Model name</th><th>α<sub>800</sub><sup>b</sup><br>(°)</th><th>W<sub>800</sub><sup>c</sup><br>(10<sup>16</sup> cm)</th><th>W<sub>12 700</sub><sup>d</sup><br>(10<sup>16</sup> cm)</th></tr></thead><tbody><tr><td>core base density ρ<sub>a0</sub></td><td>1.6 × 10<sup>-18</sup> g cm<sup>-3</sup></td><td>1.6 × 10<sup>-20</sup> g cm<sup>-3</sup></td><td>M_DENSA</td><td>112</td><td>3.5</td><td>12.4</td></tr><tr><td>jet semi opening angle θ<sub>j</sub></td><td>3°</td><td>7°</td><td>M_THETA</td><td>81</td><td>2.0</td><td>8.4</td></tr><tr><td>jet base initial density<sup>e</sup> ρ<sub>j0</sub></td><td>1.8 × 10<sup>-19</sup> g cm<sup>-3</sup></td><td>1.8 × 10<sup>-17</sup> g cm<sup>-3</sup></td><td>M_DENSJ<sup>f</sup></td><td>88</td><td>2.3</td><td>7.6</td></tr><tr><td>semi-amplitude ΔV</td><td>60 km/s</td><td>90 km/s</td><td>M_VAMP</td><td>77</td><td>1.9</td><td>8.8</td></tr><tr><td>jet variability period P</td><td>115 yr</td><td>300 yr</td><td>M_PER</td><td>74</td><td>1.8</td><td>7.4</td></tr><tr><td>jet variability profile h(t)</td><td>1 + ΔV/v<sub>0</sub> sin 2πt/P</td><td>1 + ΔV/v<sub>0</sub> {1 - 2 · mod(t/P, 1)}</td><td>M_SAWT</td><td>74</td><td>1.8</td><td>7.0</td></tr><tr><td>jet radius R<sub>j</sub></td><td>7.5 × 10<sup>14</sup> cm</td><td>3.0 × 10<sup>14</sup> cm</td><td>M_RAD</td><td>65</td><td>1.5</td><td>6.4</td></tr></tbody></table>

<footnote>a In each modified model, only one parameter at a time is changed with respect to the reference model.</footnote>

<footnote>b Shell full opening angle at z = 800 au. We obtain α<sub>800</sub> = 72° for the reference model.</footnote>

<footnote>c Full shell width at z = 800 au. We obtain W<sub>800</sub> = 1.7 × 10<sup>16</sup> cm for the reference model.</footnote>

<footnote>d Full shell width at z = 12700 au (top of the grid). We obtain W<sub>12 700</sub> = 8.0 × 10<sup>16</sup> cm for the reference model.</footnote>

<footnote>e Medium jet density case, yielding a one-sided mass-flux M̄ = 6 × 10<sup>-8</sup> M<sub>☉</sub>/yr for the reference values of θ<sub>j</sub> and R<sub>j</sub></footnote>

<footnote>f This model is equivalent to H_REF.</footnote>
15 Forme et cinématique du choc d'étrave créé dans un environnement homogène statique par un jet variable; comparaison entre simulations numériques et prédicitions analytiques (Tabone et al. 2018). 33

16 Distributions intensité-vitesse $I_{CO}(v)$ observées dans les flots L1448, Orion A, NGC2071, L1551, et Mon R2 comparées à celles simulées pour un flot poussé par un jet, en prenant en compte la dissociation et la température d'excitation (lignes pleines). Les lignes tiretées représentent les distributions masse-vitesse $m(v)$ simulées. Figure tirée de Downes & Cabrit (2003). 34

17 A gauche: Diagramme Position-Vitesse (PV) en CO le long de l'axe du flot pour la simulation de flot poussé par un jet variable de Völker et al. (1999). A droite: PV du modèle ballistique de choc d'étrave de Ostriker et al. (2001) (courbes pleines), superposé sur le PV du flot HH212 observé par Lee et al. (2000) en CO (niveaux de gris). Figure tirée de Lee et al. (2001). 35

18 Trois panneaux supérieurs: Images de l'émission $H_2$ issues de simulations 3D de flots poussés par 3 configurations de jets variables de Völker et al. (1999). Panneau inférieur: image $H_2$ de HH212 (nuances de gris, Zinnecker et al. 1998) + contours d'émission CO (Lee et al. 2000) + modèles ballistiques de chocs d'étrave d'Ostriker et al. (2001) en courbes noires. Tiré de Lee et al. (2001). 36

19 Exemple d'une partie d'une grille rectangulaire adaptative en deux dimensions dans MPI-AMRVAC, raffinée ici avec des blocs allant du niveau 1 à 3 ($2^2 = 4$ cellules par niveau d'un bloc). Ici chaque bloc peut contenir $block_n x^D=2$ cellules de même niveau suivant chaque dimension. Le bloc de niveau 1 en haut à droite est d'abord raffiné une fois en 4 blocs de niveau 2. Puis, l'un deux (le bloc inférieur gauche) est ensuite raffiné en 4 blocs de niveau 3. 40

20 Schématisation de l'approche Volume Finis, avec des variables initialisées au centre de chaque cellule et des flux calculés au centre des faces 44

21 Eventail de Riemann: Solutions d'ondes caractéristiques au problème de Riemann hydrodynamique (Riemann fan) en 1 dimension. 45

22 Solutions au problème de Riemann hydrodynamique (Riemann fan) avec en plus les sous-structures résolues par les solveurs HLL,HLLC,TVDLF, et HLLCD 46

23 Grille "totale" $\mathbb{D}$ de calcul (jaune+gris, indice (G)) avec des rangées de 2 cellules fantômes sur chaque bord: Grille de calcul intégrée/interne $M = \mathbb{D} \setminus \mathbb{B}$ (jaune, indice (M)) + rangées de cellules fantômes $\mathbb{B}$ (gris, indicé (B)). 57

24 Schématisation de deux conditions aux bords (miroir et figées) sur un exemple en $z = 0$ et avec deux rangées de cellules fantômes par bord. Les couleurs indiquent la valeur en densité $\rho$ dans chaque cellule du domaine ($\rho$ augmente en allant successivement du cyan au bleu foncé, puis au violet et au rouge). Les flèches jaunes et cyan montrent les vecteurs des vitesses totales $\vec{v} = \vec{v}_R + \vec{v}_z$ dans ces mêmes cellules et dans le référentiel au repos de l'enveloppe. La première et la seconde ligne confrontent l'état des cellules dans les deux conditions aux bords, respectivement à l'instant $t = 0$ et $t > 0$. 59
<table>
  <tbody>
    <tr>
      <td>41</td>
      <td>Même principe que la Fig. 40, mais avec un modèle de référence ayant un flux de masse du jet 100 fois plus faible : M&#775; = 6 &times; 10<sup>-8</sup> M<sub>&odot;</sub> yr<sup>-1</sup>.</td>
      <td>101</td>
    </tr>
    <tr>
      <td>41</td>
      <td>Suite.</td>
      <td>102</td>
    </tr>
    <tr>
      <td>42</td>
      <td>Largeur totale (W(z)) de la couche principale à t = 700 ans, en fonction de l'altitude z, pour les modèles de jet "haute densité" (haut) et de "densité moyenne" (bas) des Figs. 40 &amp; 41. Les courbes colorées sont obtenues en variant un seul paramètre par rapport à la simulation de référence (en pointillés et croix noires) : la densité de base du cœur &rho;<sub>a0</sub>; la demi-ouverture angulaire du jet &theta;<sub>j</sub>; la densité de base initiale du jet &rho;<sub>j0</sub>; la demi-amplitude de variabilité &Delta;V; la période P; le profil de variabilité (en dent de scie au lieu de sinusoidal); ou le rayon du jet R<sub>j</sub>. Les valeurs du paramètre que nous faisons varier dans chaque modèle sont indiquées dans la Table 5.</td>
      <td>104</td>
    </tr>
    <tr>
      <td>43</td>
      <td>Cartes à t = 10 000 ans de la (a) densité de noyaux d'hydrogène n<sub>H</sub> = n(H + H<sup>+</sup>) + 2n(H<sub>2</sub>), (b) fraction de gaz originaire du cœur f<sub>a</sub>, (c) vitesse transverse v<sub>R</sub>, et (d) vitesse axiale dans notre simulation de jet conique pulsé à travers un cœur singulier aplati. Les paramètres sont ceux de la Sect. V.3.1. L'échelle de couleur des vitesses est tronquée sur l'intervalle typiquement observé dans les flots CO : entre 0 et 50 km/s. Le contour blanc sur la moitié gauche de chaque panneau trace l'iso-contour f<sub>a</sub> = 50 %.</td>
      <td>106</td>
    </tr>
    <tr>
      <td>44</td>
      <td>Evolution dans le temps de la largeur totale de la couche principale (interface jet-cœur, croix) et externe (choc en avant, cercles) pour la simulation de la Fig. 43i. Les largeurs sont mesurées aux altitudes z = 800 au (noir) et z = 20 000 au (rouge). Les prédicitions analytiques issues du modèle équivalent de "vent X grand angle modifié" de Lee et al. (2001) (voir Sect. V.3.4.2) sont indiquées en traits noirs/rouges pointillés avec des pentes t<sup>1/2</sup> (voir texte). Est aussi indiqué en gris clair l'intervalle de largeurs observées à une hauteur projetée de 800 au dans 22 flots par Dutta et al. (2020) (deuxième et troisième quartiles en gris foncé, médiane en vert). La largeur totale à z = 20 000 au du flot HH46/47 (observé en CO par Zhang et al. 2016) est représentée en bleu.</td>
      <td>108</td>
    </tr>
    <tr>
      <td>45</td>
      <td>Diagramme PV longitudinal calculé le long de l'axe du flot simulé dans la Fig. 43i à un âge de t = 10 000 ans, pour différentes inclinaisons i du flot par rapport à la ligne de visée : 30&deg; &le; i &le; 90&deg; (edge-on). Les deux lobes du flots sont représentés (le lobe rouge aux offsets positifs, et le lobe bleu aux offsets négatifs).</td>
      <td>111</td>
    </tr>
    <tr>
      <td>46</td>
      <td>Champ d'orientation des vecteurs vitesses pour la simulation à t = 10 000 ans de la Fig. 43i, dans le référentiel au repos du cœur environnant. Pour plus de lisibilité, les flèches sont montrées pour des vitesses absolues au dessus de 0.1 km/s, et ont toutes la même longueur quel que soit le module de la vitesse. En arrière-plan couleur: carte de densité n<sub>H</sub>.</td>
      <td>112</td>
    </tr>
    <tr>
      <td>47</td>
      <td>Système de coordonnées employé dans l'outil Python Aniriniil pour l'établissement des carte synthétiques. Le plan du ciel correspond à (&delta;, &alpha;); z est confondu avec l'axe du jet reconstruit en 3 dimensions et incliné d'un angle i par rapport à &Gamma; qui est ici compté positivement en pointant vers nous.</td>
      <td>117</td>
    </tr>
    <tr>
      <td>48</td>
      <td>Représentation du cube de données synthétiques déduit par l'outil Python Aniriniil.</td>
      <td>118</td>
    </tr>
  </tbody>
</table>
49 Comparaison des cartes observées et simulées du flot HH46/47. *Panneau inférieur*: Carte d'intensité intégrée en $^{12}$CO(1 − 0) du flot HH46/47 obtenue avec ALMA par Zhang et al. (2016). Les intervalles d'intégration sont indiqués dans l'encadré. *Panneau supérieur*: Carte de masse intégrée sur les mêmes intervalles de vitesses 2 km/s ≤ |vobs| ≤ 10 km/s, issue de la simulation de la Fig. 43i à t = 10⁴ ans. Ici, seule la matière originaire du coeur est représentée. Le flot est incliné de i = 55° par rapport à la ligne de visée, à une distance de 450 pc. Le domaine fini de simulation ne restitue qu'une partie de tout le flot incliné simulé et tronque ainsi la carte au-delà de 20", après projection sur le plan du ciel. 121

50 Comparaison des diagrammes PV transverses simulés et observés dans HH46/47. *Panneaux de droite*: Émission en $^{12}$CO(2 − 1) de HH46/47 le long de deux coupes larges de 450 au à des distances projetées de 8.1 × 10$^{16}$ cm (b) et 1.00 × 10$^{17}$ cm (d) par rapport à la source centrale. Les ellipses rouges montrent les modèles analytiques de deux couches paraboliques produites par un vent X modifié et ayant le meilleur ajustement aux observations. Figures adaptées de Zhang et al. (2019). *Panneaux de gauche*: Diagrammes calculés à partir de la simulation de la Fig. 43i à t ≃ 10 000 ans avec une inclinaison de i = 55° par rapport à la ligne de visée, à une distance de 450 pc. Seule la matière originaire du coeur est représentée. Les coupes ont une largeur de 400 au et sont effectuées aux mêmes altitudes projetées que les observations. 122

51 Comparaison des diagrammes PV longitudinaux observés et simulés pour HH46/47: (a) PV observé avec ALMA en $^{12}$CO(2 − 1) par Zhang et al. (2019), dans une coupe large de 450 au. (b) PV calculé pour la simulation de la Sect. V.3.2, le long d'une coupe large de 400 au à un âge de t = 9900 ans (33 périodes), et une inclinaison de i = 55° sur la ligne de visée. (c) Comme (b) mais incluant seulement la contribution du gaz issu du coeur ambiant. 124

52 Comme la Fig. 51 mais cette fois les deux diagrammes PV longitudinaux synthétiques (b) et (c) incluent uniquement la matière originaire du coeur, et diffèrent seulement dans la phase du cycle de variabilité où le jet est observé: 9900 ans (exactement 33 périodes) pour (b) et 10⁴ ans (33+1/3) pour (c). 126

53 Comparaison de la distribution masse-vitesse (MV) observée dans HH46/47 et simulée. *Symboles noirs connectés*: MV de la matière originaire du coeur, évalué à t ≃ 10$^{4}$ans à partir de la simulation de la Fig. 51, incliné de i = 55° par rapport à la ligne de visée. *Symboles bleus et rouges*: distributions MV dans les lobes bleu et rouge du flot HH46/47, dérivées par Zhang et al. (2016). Les carrés vides supposent une émission en CO optiquement mince, tandis que les cercles pleins incluent une correction d'épaisseur optique dépendant de la vitesse (voir texte). Les lignes tiretées montrent les pentes des lois de puissance γ = −3.4 (rouge) et γ = −2.7 (bleu). 127

54 Comparaison des diagrammes PV longitudinaux observé dans CARMA-7 et simulés: (a) PV observé avec ALMA en $^{12}$CO(2-1) par Plunkett et al. (2015) dans une coupe large de 374 au. (b) PV calculé à un âge de t = 500 ans pour une simulation de période P = 30 ans (voir Sect. VI.3.1) le long d'une coupe large de 400 au, et pour un angle d'inclinaison de i = 85° par rapport à la ligne de visée. (c) Comme (b) mais incluant seulement la contribution de la masse issue du coeur. . . . 130
<table>
  <tr>
    <td>55</td>
    <td>Comparaison du PV longitudinal observé de CARMA-7 (a) avec deux diagrammes PV synthétiques incluant uniquement la matière originaire du coeur, et différant par l'âge de la simulation: (b) t = 500 ans, (c) t = 10 000 ans.</td>
    <td>131</td>
  </tr>
</table>
LISTE DES TABLES

<table><tr><td>1</td><td>Liste des paramètres libres principaux et initiaux ($t = 0$) de l'enveloppe environnante utilisée dans les Chapitres V et VI</td><td>56</td></tr><tr><td>2</td><td>Liste des paramètres principaux du jet pulsé utilisés dans les Chapitres V et VI</td><td>62</td></tr><tr><td>3</td><td>Résolution et temps de calcul en fonction du niveau maximal de raffinement</td><td>77</td></tr><tr><td>4</td><td>Paramètres des simulations des jets pulsés cylindriques de la Sect. V.1</td><td>83</td></tr><tr><td>5</td><td>Paramètres des simulations de jet pulsé conique à hautes densités, avec les ouvertures angulaires résultantes et les largeurs totales associées (Fig. 40)</td><td>96</td></tr><tr><td>6</td><td>Paramètres des simulations de jet conique à "densité moyenne", avec les largeurs totales et angles d'ouverture résultants</td><td>97</td></tr></table>
## RÉSUMÉ

L'efficacité observée dans la formation des étoiles jusqu'au parsec ne dépasse pas quelques pourcents, sans oublier le décalage de la fonction initiale de masse (IMF) à seulement ~30% de la masse du cœur pré-stellaire initial. Comprendre en détail le processus derrière d' aussi faibles efficacités reste encore à ce jour une question ouverte. En outre, les simulations numériques les plus récentes ont démontré que la turbulence et les champs magnétiques à eux seuls ne peuvent suffire à reproduire de telles valeurs. Elles montrent que la rétroaction des flots protostellaires joue un rôle crucial en perturbant les écoulements d'accrétion, en évacuant la matière des coeurs, et/ou en maintenant la turbulence. Malheureusement, que ce soit en termes de volume de nuage affecté, d'impulsion injectée, de masse entraînée, ou d'impact sur le disque et l'enveloppe en effondrement : l'importance de cette rétroaction dépend fortement de la géométrie sous-jacente du vent protostellaire. Cette dernière reste encore débattue : "vent X grand angle" rapide, vent de disque MHD plus lent, ou jet collimaté ? De toute évidence, afin d'évaluer fiablement l'impact de la rétroaction des flots sur la formation stellaire, il est d'une importance cruciale de déterminer la géométrie de vent la plus réaliste (et/ou les géométries que nous pouvons exclure). Pour apporter une nouvelle contribution quant à cette question, nous présentons des simulations numériques de flots poussés par un jet pulsé collimaté, lancé à travers un cœur pré-stellaire stratifié. Nous comparons nos simulations avec les observations ALMA récentes, ainsi qu'avec les prédictions analogues pour un vent X grand angle. Nos simulations sont les premières à combiner sur une échelle de 0.1 pc la variabilité du jet, la stratification en densité de l'enveloppe et des échelles de temps de 10 000 ans comparables aux flots jeunes observés. Les prédictions de nos simulations en termes de largeur de flot, de diagrammes position-vitesse, et de distribution masse-vitesse, montrent une ressemblance frappante avec les observations ALMA de flots CO tels que HH46/47 et CARMA-7. L'accord est même plus prometteur qu'avec les modèles de flots poussés par un "vent X grand angle". Ces résultats pourraient avoir une implication majeure sur le rôle des flots dans la régulation de la formation stellaire.

## MOTS CLÉS

Flots protostellaires, formation stellaire, simulations hydrodynamiques, jets, protoétoiles

## ABSTRACT

A long-standing open question in star formation is the process responsible for its low efficiency on parsec scales (a few %), and for shifting down the Initial Mass Function (IMF) to only ~30% of the prestellar core mass distribution. The most recent numerical simulations show that neither turbulence nor magnetic fields can, alone, reproduce these low efficiencies, and that feedback by protostellar outflows must play a crucial role by disrupting accretion streams, expelling material from cores, and/or sustaining turbulence. Unfortunately, the magnitude of outflow feedback (affected cloud volume, injected momentum, entrained mass, impact on the disk and infalling envelope) depends strongly on the underlying protostellar wind geometry, which remains uncertain and heavily debated: a fast wide-angle "X-wind", a slower MHD disk wind, a narrow jet ? Clearly, if we want to reliably assess the role of outflow feedback in star formation, it is of utmost importance to determine which wind geometry is the most realistic (and/or which one can be excluded). As a new contribution towards this goal, we present, for the first time, numerical predictions for outflows driven by a narrow pulsed jet in a stratified prestellar core. We compare our simulations against recent ALMA observations and analogous predictions for a wide-angle X-wind. Our simulations are the first to combine jet variability, ambient density-stratification, and long timescales up to 10 000 yrs (typical of young outflows) on scales up to 0.1 pc. We find that the predicted widths, position-velocity diagrams, and mass-velocity distribution, show striking resemblance with ALMA observations of CO outflows such as HH46/47 and CARMA-7, and in closer agreement than models based on a wide-angle "X-wind". The results obtained in this work could have major implications for the feedback of protostellar outflows on star formation.

## KEYWORDS

Keywords : Protostellar outflows, star formation, hydrodynamical simulations, jets, protostars