© The Authors 2023

1 Introduction

Protostellar jets and outflows are an ubiquitous phenomenon of the star formation process from the early Class 0 to the late Class I phase, when the parental envelope is dissipated. They are an important agent of star formation feedback, which affects the gas’s physical and chemical properties from cloud scale down to the central parental cocoon, where mass accretion is occurring. On the one hand, they act as a source of energy and momentum into the parental cloud, and they disperse the material of the parental core, directly impacting the star formation efficiency and the final stellar mass (see Frank et al. 2014 for a review); on the other hand, they are thought to remove a significant fraction of the angular momentum from the star-disk system, enabling the gas in the accretion disk to reach the central protostar (Konigl & Pudritz 2000).

A natural consequence of the propagation of such high velocity outflows through the protostellar envelope and the ambient molecular medium are shock fronts (Reipurth & Raga 1999). These shocks both heat and compress the gas, which favors chemical reactions in the gas phase while they modify the dust grain properties and release a fraction of their material in the gas phase through, for example, sputtering and shattering, thereby leading to a different chemical composition than observed in the ambient, preshock gas. Commonly, the high relative abundance of CO and its low energy J transitions makes it a good tracer for outflows in the conditions of cold molecular clouds. Even more, the morphology, velocity, and sizes of a molecular outflow depend on the observed tracer, the luminosity, mass, and age of the outflow, which indicates its evolution and its inner structure (Bally 2016).

Over the years, there have been intense numerical efforts to simulate the chemical and dynamical structure of protostellar outflows and their evolution in the ambient interstellar gas. Smith (1994) and Smith et al. (1997) carried out 3D numerical simulations of dense molecular jets drilling through a molecular environment. They succeeded in reproducing the morphology of the “classical” CO bipolar outflows, which image the accumulated and accelerated cool gas. They showed that the simulations predict infrared shock structures remarkably similar to those found in highly collimated Class 0 outflows. They proved that in spite of the limitations in their description of the physical conditions (dust, magnetic field), these hydrodynamical models could provide a meaningful description of the dynamics of bipolar outflows from young stars, accurate enough images for comparison with observations, and the outflow mass and energy contribution to the interstellar medium.

Raga et al. (1990) showed that simulations with an arbitrarily imposed ejection velocity variability lead to the formation of chains of internal working surfaces traveling down the jet flow. Such features strikingly resemble Herbig-Haro (HH) flows, with a chain of aligned knots close to the outflow source, and a large “head”, resulting from the “turning on” of the jet flow at larger distances. Later on, the authors showed that a two-mode ejection velocity variability leads to the formation of chains of “short period knots”, which catch up with each other to form “long period knots”. This class of models has successfully accounted for the properties of HH flows and they nowadays dominate the literature on the theory of astrophysical jets (Canto 1985; Raga et al. 1990; Raga & Canto al. 2003; Noriega-Crespo et al. 2014; Frank et al. 2014; Rabenanahary et al. 2022). On the other hand, most of the hydrodynamical codes developed so far have not included a chemical network and thereby do not address the molecular gas composition, in contradiction with the observational evidence that outflows from Class 0 protostars are often chemically active (Bachiller et al. 2001; Arce et al. 2007; Lefloch et al. 2017; Ospina-Zamudio et al. 2019a; De Simone et al. 2020).

It is only with the advent of the first generation of chemo-hydrodynamical codes that it is now possible to accurately model the physical and chemical outflowing gas evolution, and to build spectroscopic diagnostics which can be tested against the observational constraints provided by the large (sub)millimeter arrays such as the Northern Extended Millimetre Array (NOEMA) and the Atacama Large Millimeter/submillimeter Array (ALMA), and to investigate in detail the processes that determine the jet-outflow properties in order to quantify the energetic and chemical feedback of newly born stars on their environment and eventually on galaxies.

In this work, we use WALKIMYA-2D, a new 2D hydrodynamical code coupled to a reduced gas phase chemical network allowing for the evolution of the CO chemistry in molecular outflows to be followed (Castellanos-Ramírez et al. 2018). The advantage of this code is that it computes the time evolution of chemical species in a full gas dynamic simulation. The new functionalities of WALKIMYA-2D allow us to investigate more thoroughly the dynamics and the physical structure of protostellar outflows.

We have selected the high-velocity outflow associated with the intermediate-mass Class 0 protostellar system Cep E-mm, located in the Cepheus OB3 association at a distance of 819 ± 16 pc (Karnath et al. 2019), whose CO emission observed at 1″ with the Institut de Radioastronomie Millimétrique (IRAM) interferometer is displayed in Fig. 1. This outflow has a luminosity of ≈100 L_⊙ (Lefloch et al. 1996) and a core mass of 35 M_⊙ (Crimier et al. 2010). The source has been the subject of several detailed studies at arcsec angular resolution, in particular in the CO rotational transitions in the (sub)millimeter domain (Lefloch et al. 1996, 2015; Ospina-Zamudio et al. 2019b), which have constrained the jet and outflow dynamical parameters (mass, mass-loss rate, momentum, density, temperature). The recent study by Schutzer et al. (2022) has brought a detailed view of the jet structure and its time variability, showing a complex interaction with the ambient protostellar material. It is also one of the few protostellar jets for which a full 3D picture of the gas kinematics is available. Overall, this outflow appears as an excellent testbed for chemo-hydrodynamical codes, with the prospect of getting more insight into the jet and outflow formation process.

The goal of this study is to obtain an accurate numerical model of the Cep E molecular outflow in agreement with the constraints provided by the large (sub)millimeter observatories. In particular, we aim to understand the dynamics of the jet and envelope interaction near the source. The synergy between the WALKIMYA-2D numerical simulations and the observations with the IRAM interferometer has allowed us to peer into the outflow formation mechanism with unprecedented detail and to retrieve the history of the mass-loss events that have shaped the outflow.

This paper is organized as follows: in Sect. 2, we describe the numerical setup and the boundary conditions, which were chosen in agreement with the observational constraints on Cep E and present the best fitting model. The following sections present our main results: the outflow formation in Sect. 3, the gas acceleration mechanism in Sect. 4, the mass-loss history in Sect. 5, and the CO emission in Sect. 6. Finally, we summarize our conclusions in Sect. 7.

Fig. 1 CO 2–1 line emission in the Cep E-mm outflow as observed with the PdBI at 1″ resolution (Lefloch et al. 2015). Four main velocity components are detected: a) the outflow cavity walls (black) emitting at low velocities, in the range [−8;−3] km s−1 and [−19;−14] km s−1 in the northern and southern lobe, respectively; b) the jet, emitting at high velocities in the range [−135;−110] km s−1 in the southern lobe (blue) and in the range [+40;+80] km s−1 in the northern lobe (red); c) the southern terminal bow shock HH377, integrated in the range [−77,−64] km s−1 (green); d) the northern terminal bullet NB integrated in the velocity range [+84,+97] km s−1 (magenta). First contour and contour interval are 20 and 10% of the peak intensity in each component, respectively. The synthesized beam (1″.07 × 0″.87, HPBW) is shown in the bottom left corner. The main axis of the northern and southern lobes of the jet are shown with black arrows. (Schutzer et al. 2022).

2 Simulations

2.1 WALKIMYA-2D

In order to study the temporal evolution of the molecular composition of a gas, one needs to solve the rate of change of the abundances of the different species contained in it. It is required to construct a chemical network gathering the creation and destruction reactions for the different species. We have used the 2D chemo-hydrodynamical WALKIMYA-2D that solves the hydrodynamic equations and a chemical network, on an axisymmetric numerical adaptive mesh. A complete description of the code is presented in Castellanos-Ramírez et al. (2018). Initially, both the jet and the surrounding quiescent protostellar gas have a CO abundance of 1.6 × 10⁻⁴ relative to H₂. The chemical network is designed mainly to follow the evolution of the CO chemistry, as induced by the heating-cooling processes in the jet propagation and the shock interaction(s) with the ambient gas. For computational reasons (to keep computational times within reach), the network is reduced to 14 chemical species, including C, O, H₂O, OH, and CO (Castellanos-Ramírez et al. 2018). The reaction rates were obtained from the UMIST database (McElroy et al. 2013). WALKIMYA-2D and its chemical network were successfully benchmarked against laboratory experiments to study the NO formation by electrical discharges using both a zero-dimensional and a gas dynamic model approach Castellanos-Ramírez et al. (2018). Also, in the astrochemical context, the code was successfully benchmarked against the model of dark molecular cloud presented in McElroy et al. (2013), and has been used to explain the CO emission produced by the Orion Fingers in the Orion KL region (Rodríguez-González et al. 2022).

The energy loss rate is calculated adopting the same prescription as in Rivera-Ortiz et al. (2019, see references therein): for temperatures larger than 5800 K the cooling function considers the atomic contribution while in the lower temperature range, it uses a parametric molecular cooling function based on CO and H₂.

The adaptative mesh uses seven levels of refinement, yielding 4096 × 1024 cells, in a computational domain (5 × 1.25) × 10⁴ au, allowing a resolution as high as 12.2 au per cell side. We used a reflective boundary condition for the symmetry axis and a free outflow boundary condition for all the other borders. The size of the mesh is large enough that the outer boundaries do not affect the simulation. The jet is injected on the left side of the simulation box with the physical conditions indicated in Table 1.

2.2 Initial conditions: the protostellar core

In order to describe the protostellar core, we have adopted the physical conditions derived by Crimier et al. (2010) from a simple 1D modeling of the dust spectral energy distribution between 24 and 1300 µm, using the radiative transfer code 1D DUSTY (Ivezic & Elitzur 1997). The temperature is fitted by a broken radial power law T ∝r^β with β = −0.8 in the range 50–300 K and β = −0.4 in the range 7–50 K. In Crimier’s model, the density profile is assumed to start at a radius r_in = 70 au and to follow a single radial power law distribution n(r)= n₀ × (r₀/r)^α.

In order to avoid density singularities in the limit r → 0 in the numerical modeling, we have adopted the slightly modified density profile $$n\left(r\right)=\frac{n_0}{1+{\left(1/r_0\right)}^{\alpha }},$$(1)

where n₀ = 10⁹ cm⁻³, r₀ = 100 au and α = 1.9.

The simulation takes into account the role of gravity, in agreement with the source density stratification evidenced from Crimier’s analysis. Our model includes the resulting gravitational force as a source term, which reads: $$g=-\frac{2c_0^{2}r}{r_0^2\left(1+{\left(1/r_0\right)}^{\alpha }\right)}$$(2)

where c₀ is the local sound speed. The inclusion of this gravity term ensures hydrostatic equilibrium close to the source.

2.3 Initial conditions: the jet

The jet is bipolar and consists of a narrow (diameter < 400 au) central component surrounded by a collimated layer with a radial extent up to 1000 au, as shown by Schutzer et al. (2022). The jet appears young, with a dynamical age of 1400 yr, and compact, with a length of ~0.1–0.14 pc (18000–29000 au) for the southern and northern lobes, respectively. The physical conditions (temperature, H₂ gas density) in the outflow were estimated by Lefloch et al. (2015) and Ospina-Zamudio et al. (2019b) from a CO multi-rotational line study complemented with CO J = 2−1 observations at 1″ angular resolution with the IRAM interferometer. In the southern lobe, the jet appears to consist of a warm (T = 80–100 K) gas component of H₂ density n(H₂) = (0.5–1.0) × 10⁵ cm⁻³ component and a higher-excitation component of n(H) = (0.5–1.0) × 10⁶ cm⁻³ and temperature (T = 400–750 K), which the authors associated with the high-velocity knots. Similar physical conditions are found in the northern lobe, with a kinetic temperature T = 180–300 K and gas density n(H₂) = (0.6–2.0) × 10⁵ cm⁻³. Therefore, the jet is rather massive, with a total mass of 0.03 M_⊙ (0.09 M_⊙) in the southern (northern) lobe, after taking into account the revised distance to the source (Karnath et al. 2019).

The jet asymptotic radial velocities were determined from the CO 2−1 line profiles as +65 km s⁻¹ and −125 km s⁻¹ in the northern and southern lobe, respectively (Lefloch et al. 2015). The jet proper motions in both lobes were measured by Noriega-Crespo et al. (2014) from combining multiple mid-infrared IRAC 4.5µm observations and H₂ 2.12 µm images at a difference of 16 years. We note that the revised distance of Cep E-mm (820 pc instead of 730 pc previously) does not significantly affect the tangential velocities, nor the inclination angle of the jet with respect to the plane of the sky derived by Lefloch et al. (2015), which is now 47° (40° previously estimated) and the estimated dynamical age of the jet (about 10³ yr). Based on these observations, the molecular jet velocity is estimated ~100 km s⁻¹ (150 km s⁻¹) in the northern (southern) lobe, respectively. Kinematical analysis of the molecular gas knots in the jet shows evidence for time and velocity variability of the mass-ejection process in Cep E (Schutzer et al. 2022).

In the simulation, the jet is launched at z = 0, with a radius r_j, a gas density n_j(H) = 10⁶ cm⁻³, and a temperature T_j = 300 K. In order to account for knot formation inside the jet, we introduced variability in the gas injection, which we modeled following the equation $$V_j=V_{j,0}\left[1+{\delta }_{\upsilon }\cos \left(2\pi t/\tau \right)\right]$$(3)

where τ = 130 yr is the injection mass period and t is the evolutionary time. The jet injection velocity V_j,0 and the relative amplitude variability δ_υ are free parameters of the simulation. We adopted δ_υ = 0.05–0.08, which implies velocity variations of 10–15 km s⁻¹, consistent with the typical variations reported by Schutzer et al. (2022) in Cep E and in other outflows such as IRAS 04166+2706 (Santiago-García et al. 2009).

2.4 Comparison with observations

We have run four models M1–M4, whose initial jet parameters are listed in Table 1. For the sake of simplicity, our WALKIMYA-2D simulations do not include the effect of jet precession. Comparison between numerical simulations and observations are made with the northern outflow lobe of Cep E, whose entrained gas dynamics is better revealed in CO interferometric observations (Schutzer et al. 2022): morphology, gas acceleration, time-variability, and knot formation.

The simulations give us access to the physical structure of the outflow and its chemical properties. Filtering the emission of the high- and low-velocity components allows direct comparison between the simulations and the observational signatures of the distinct outflow components and therefore it is possible to quantify the physical processes involved in the outflow formation. In order to compare our numerical simulations with the CO observations, the emissivity has been computed directly from the hydrodynamical simulations assuming local thermodynamic equilibrium (LTE). We have constructed synthetic CO J = 2−1 maps of the outflow cavity and the jet, by integrating the emission in the velocity intervals [3;7 km s⁻¹] and [50;150 km s⁻¹], respectively. Those maps were subsequently convolved with a Gaussian profile whose size (FWHM) corresponds to the synthetic beamsize of the interferometer (1″ or 820 au).

We first explored the parameter space and searched for the model which best reproduces the qualitative and quantitative properties of the Cep E northern outflow lobe. Overvall, we found that all four models M1–M4 show similar results. Therefore, in what follows we present the results of model M1, that is, the simulation that best accounts for the molecular gas observations of Cep E, whose initial jet parameters are summarized in Table 1. The initial jet velocity is 200 km s⁻¹, with a velocity fluctuation amplitude of 0.08, and an ejection radius of 100 au.

We compare the numerical results to the main observational features of the northern outflow lobe, identified in the recent study by Schutzer et al. (2022): morphology, gas acceleration, time-variability, and knot formation.

3 Outflow formation

Figure 2 displays the molecular gas distribution at five timesteps of the outflow formation, from t = 0 to t = 2000 yr by steps of 500 yr. (A video is available at Fig. 2).

The jet is launched at t = 0 into gas of density 10⁹ cm⁻³. The density contrast between the jet (10⁶ cm⁻³) and the ambient gas (10⁹ cm⁻³) makes the shock propagate into the ambient gas at a velocity much lower than the injection velocity υ_j,0, an effect that has been studied in detail already by Raga et al. (1990). Turning on the jet creates a high-velocity, collimated narrow component which propagates in the ejection direction. It is surrounded by a cocoon of gas that expands laterally, creating a low-velocity (0–10 km s⁻¹) “cavity” component, with dense walls surrounding the jet. The first working surface interacts with the very dense gas close to the jet source, which makes it decelerate, causing the second ejected knot to catch up with the first bow-shock. As shown in Fig. 2, this process repeats a few times until the bowshock leaves the inner envelope and propagates into the protostellar gas, reaching its terminal velocity. It is only after ~200 yr and the launch of 3 knots that the jet manages to drill out the inner 1000 au of the protostellar envelope and to propagate into the surrounding gas (second panel from top). As a consequence the overpressure resulting from the accumulation of knots pushes the protostellar envelope laterally away from the jet, causing the formation of a wide-angle outflow cavity, and a density increase in the low-velocity cavity walls at the base of the jet. This implies that a kinematic age computed from the bowshock position and velocity is actually a lower limit to the real age of the outflow.

As time proceeds, the model shows a leading jet head and a series of traveling “internal working surfaces” (IWS). These structures form as the result of the ejection time-variability and are observed as “knots”, or overdensity regions with a small size of a few ~ 100 au. At t = 1000 yr, 3 knots are easily detected along the jet, while at t = 2000 yr, 8 eight knots are easily identified. It is worth noticing they tend to expand radially as they propagate along the jet, tracing the wings of inner bowshocks, since the slower material interacts with the faster material ejected at later times. These wings appear to drive the formation of complex structures inside the outflow cavity as can be seen in the bottom panel of Fig. 2. Fast moving knots tend to accumulate at the head of the bow, increasing the density at the tip. After 1500 yr, the knots located at distance >15 000 au display wings with a typical size of 1000 au.

One expects the gas density to decrease as a result of its radial expansion in the course of the jet propagation. However, looking in detail at the jet density field at e.g. t = 1000 yr (second panel in Fig. 2), it appears that the gas density does not decrease monotonically along the narrow jet. The density distribution of fast material along the propagation direction is analyzed in Sect. 5.2.

At t = 1500 yr, a second outflow cavity has now formed. The first outflow cavity has reached a size of ≈15 000 au and expanded up to a radius of about 5000 au (6″ at the distance of Cep E-mm), and the head of the outflow has reached a distance of 2.8 × 10⁴ au (34″). At later times, the second cavity expands to reach a similar radius of ≈5000 au.

Inspection of Figure 2 reveals that a complex network of relatively dense structures (n ~ 10⁵–10⁶ cm⁻³) forms inside the outflow cavity very early. One can note the presence of “filaments” along the jet in the first 1000 yr. As can be seen in the panel at t = 1000 yr, a filament has formed along the jet. Located at R = 500 au, it is detected up to z = 2000 au and appears as a thin, yellow, wiggling structure parallel to the jet. We also note the presence of a “shell”, or a bow, which formed in the dense protostellar envelope and moves slowly (~20 km s⁻¹) as it propagates over 8000 au in 2000 yr, whereas the jet reached 40 × 10³ au. This dense shell (or bow) connects the jet envelope of entrained gas to the cavity walls.

Our numerical results on the outflow morphology appear in good agreement with the NOEMA observations of the Cep E northern lobe (Fig. 1). As can be seen in Table 2, both the length and the radius of the outflow cavity are correctly reproduced; this good match between observations and simulations is obtained for a computational timescale of 1500 yr, similar to the estimated CepE dynamical age (1400yr). Even more, the jet momentum obtained from the simulation is 3.78 M_⊙ km s⁻¹, comparable to the momentum obtained by Lefloch et al. (2015) of 2.5 M_⊙ km s⁻¹, taking into account the inclination angle.

In order to better understand the role of the core initial conditions in shaping the outflow morphology, we also carried out simulations without the inclusion of the high-density core given by Crimier nor the related gravitational term. In these cases, it turned out that the low-velocity outflows formed in the ambient medium of the jet propagation were actually far too collimated with respect to what is observed, in particular close to the protostar and the jet launch region. Therefore, it appears that the wide opening of the low-velocity outflow cavity depends on the density of the dense inner region and the contribution of the latter to the gravitational field, two parameters that are related to the evolutionary stage of the system. The importance of including a density distribution and a gravitational term is in agreement with Raga & Cabrit (1993) and Cabrit et al. (1997) who proposed that a steep radial density decrease would produce a wider opening angle for jet-driven outflows.

Fig. 2 Outflow formation and impact of the jet on the gas density distribution as a function of time. Snapshots for 0, 500, 1000, 1500, and 2000 yr. The injection velocity variability forms internal working surfaces, or knots, and a structured cavity. The leading bowshock accumulates the mass of several knots. A movie is available online.

4 Gas acceleration

Cep E displays evidence of gas acceleration along the jet over a distance of 1000 au away from the protostar (Fig. 3). This effect was reported and observationally characterized by Schutzer et al. (2022). We show in the bottom panel of Fig. 4, the synthetic Position-Velocity diagram of the CO gas along the main axis of the jet. Two kinematical components are clearly identified: (i) a low velocity component, between +0 and +10 km s⁻¹ up to 8000 au (10″) from the source, associated with the outflow cavity walls; (ii) a high velocity component, that displays a periodical behavior, reflecting the mass injection variability and the internal shocks, where the fast material catches up with the slower material. The jet is detected at a radial velocity close to V_j = 65 km s⁻¹ at z > 2000 au from the source.

Our simulations show that the effective velocity of the outflow is significantly lower than the injection jet velocity. This appears as a result of the jet interaction with the dense central envelope. In this case, the Position-Velocity diagram in Fig. 4 shows that the CO gas emission accelerates from ambient velocity (V = 0 km s⁻¹) to the terminal jet velocity on a short scale <5 × 10³ au. Close to the source, gas structures accelerated up to V ~ 50 km s⁻¹ are also detected along the jet axis up to 1000–2000 au (bottom panel of Fig. 4).

An observational signature of envelope gas can also be found in the Position-Velocity diagram of the CO J = 2−1 emission across the jet main axis, as can be seen in Fig. 5. The CO emission contours show how the ambient material initially at rest at V_lsr = −11 km s⁻¹ is gradually accelerated as one gets closer to the location of the protostar at α = 0.0″. In addition to the jet mainly detected up to V ~ 60 km s⁻¹, signatures of high-velocity knots (up +90 km s⁻¹) are also identified close to the protostar (Fig. 5).

It is worth remembering that the synthetic Position-Velocity diagrams use a wide slit in the direction of the projected jet axis and the contribution of all the possible velocity projections along the line of sight, using an angle of 47° with respect to the plane of the sky, whereas the Cep E observational plots (Fig. 3) display the radial component of the jet velocity. For the sake of comparison with the Cep E jet, we have checked whether the solution found by Schutzer et al. (2022), (V − V_lsr)/V_j = exp(−δ₀/δ) could provide a reasonable fit to the synthetic jet velocity profile along the main axis. The yellow curve drawn in the bottom panel of Fig. 4 traces the best fitting solution of Schutzer et al. obtained for a length scale δ₀ = 690 au, applied to the radial jet velocity V_r = 65 km s⁻¹. Taking into account the inclination of the jet with respect to the line of sight, the agreement between numerical and observational jet velocity profiles is very satisfying both qualitatively and quantitatively.

To conclude, our simulation satisfyingly accounts for the acceleration of material from the protostellar envelope by the Cep E jet, from ambient velocity up to reaching the radial terminal jet velocity V_j ~ 90 km s⁻¹. This process occurs over a scale of 5000 au.

Fig. 3 Cep E northern outflow lobe. Position-Velocity diagram of the CO J = 2−1 line emission along the jet main axis. Positions are in arcsec offset relative to the location of the driving protostar Cep E-A. First contour and contour interval are 10% of the peak intensity. We have superimposed in yellow the best fitting solution (V − Vlsr)/V0 = exp(−δ0/δ) where δ0 = 690 au (Schutzer et al. 2022).

Fig. 4 Best fitting model (Ml) to the Cep E northern outflow, (top) Density distribution at 1250 yr. The jet has propagated over a distance of 2.4 × 104 au and formed a low-velocity cavity of 3000 au radius. Five knots and their associated wings are detected along the jet.(bottom) Position-Velocity diagram of the CO emission along the jet main axis obtained using the 47° outflow inclination angle with respect to the plane of the sky.

Fig. 5 Cep E northern outflow lobe. Position-Velocity diagram of the CO J = 2−1 line emission across the jet main axis at δ = 0″.Positions are in arcsec offset relative to the location of the driving protostar CepE-A. First contour and contour interval are 10% of the peak intensity (Schutzer et al. 2022).

5 Mass-loss history

5.1 The jet

We first measured the mass of outflowing material as a function of time. The mass integration was performed over the whole velocity range. The result is displayed in Fig. 6 and shows that the outflow mass increases linearly with time at a rate of 1.6 × 10⁻³ M_⊙yr⁻¹. At a time of 1500 yr in the simulation (similar to the dynamical age of Cep E) the mass of outflowing gas amounts to 2.3 M_⊙, a value in good agreement with the observational determination in the northern lobe (see Table 2).

In a second step, we have studied the variations of the jet mass as a function of time. We took into account all the gas at velocity V > 50 km s⁻¹ to estimate the jet mass and follow its variations. The variations of the jet mass as a function of time are displayed in Fig. 7. After an initial delay of ~200 yr, which corresponds to the time needed for the knots to drill the envelope, the jet mass increases linearly with time, at a rate Ṁ = 2.3 × 10⁻⁵ M_⊙ yr⁻¹ in very good agreement with the observational determination by Schutzer et al. (2022), (2.7 × 10″⁵ M_⊙ yr⁻¹). As discussed previously in Sect. 4, the protostellar material appears to feed the high-velocity jet through an entrainment process. An observational signature of this effect is found in the Position-Velocity diagram of the CO J = 2−1 emission across the jet main axis, displayed in Fig. 5. The CO emission contours show how the ambient material initially at rest at V_lsr = −11 km s⁻¹ is gradually accelerated as one gets closer to the location of the protostar at α = 0.0″. In addition to the jet mainly detected up to V ~ 60 km s⁻¹, signatures of higher-velocity knots (up +90 km s⁻¹) are also identified close to the protostar (Fig. 5).

Third, in order to quantify the magnitude of the entrainment effect to the jet mass, we have disentangled the contributions of the injected material from the entrained material. Theoretically, the jet mass injection rate per unit of time can be estimated as a function of the injected material density n_j,·, the jet velocity υ_j and the jet injection radius r_j $$\left[\frac{\dot{M}}{M_{\odot }{\text{yr}}^{-1}}\right]=1.5\times {10}^{-6}{\left[\frac{r_j}{50\text{\hspace{0.17em}au}}\right]}^2\left[\frac{n_j}{{10}^6{\text{cm}}^{-3}}\right]\left[\frac{{\upsilon }_j}{165\text{\hspace{0.17em}km\hspace{0.17em}}{\text{s}}^{-1}}\right].$$(4)

Numerically, the jet mass injection rate in the computational domain is in reasonable agreement with the simple, theoretical description from Eq. (4), (M = 7.3 × 10⁻⁶ M_⊙ yr⁻¹ with the conditions for model M1). In both numerical and theoretical cases, the jet mass injection rate is lower than the jet mass rate Ṁ of the simulation by a factor of ≈3. We consider the difference between both values as very significant as it highlights the importance of the entrainment process of circumstellar material in the formation of molecular jets.

Fig. 6 Model M1. Outflow mass evolution as a function of time.

Fig. 7 Model M1. Jet mass evolution as a function of time. The points are obtained from the simulation every 100 yr. The black continous line is a linear fit of the points starting from 300 yr, with a slope of 2.3 ×10−5M⊙ yr−1 and the dashed line corresponds to the material injection rate of Ṁ = 7.3 × 10−6 used in the simulation.

5.2 The knots

We now examine the physical properties of the knots which form along the jet in the numerical simulation and their evolution with time. We then compare them with their observational counterparts in the Cep E jet.

5.2.1 Formation and evolution

We have studied the gas density distribution along the jet as a function of time in the simulation up to 2000 yr. We found that the evolution is characterized by very short timescales on the order of 25 yr. In order to illustrate and better explain the process at work, we present here in Fig. 8 a montage of the jet density profiles at intervals of 25 yr over a time of 125 yr, between 1450 and 1575 yr. The mean density distribution appears steady and constant at distances further than 10⁴ au, while at shorter distances, there is a decreasing density profile. This is a fossil signature of the steep initial density profile of the envelope, since it is expected to be difficult to remove the dense envelope material close to the source, even when it has been processed by the countinuous injection of the jet. The peaks, which appear on top of the mean density distribution, result from the ejection velocity variability in the jet source (Sect. 2.3), which causes fast material to shock with lower velocity gas. Each of these peaks corresponds to a local overdensity structure in the jet, whose properties were derived from a Gaussian profile fitting, leading to a typical size of 1000 au, a density of about 10⁵ cm⁻³, and a typical mass of 10⁻³ M_⊙. The physical properties of these peaks are actually very similar to those of the CO knots reported by Schutzer et al. (2022) in the Cep E jet, and, for that reason, they are probably their counterparts in the numerical simulation. Based on Fig. 8, it appears that at distances larger than 10⁴ au the density and the knot distributions are almost steady, which reflects the fact that the knots are propagating into the low-density material of the jet at a speed of about 525 au in 25 yr. The high density in the inner 10⁴ au makes the situation drastically different. Rivera-Ortiz et al. (2019) have shown that gas flowing into an environment of similar density can create reverse shocks that propagates at a significant fraction of the original velocity relatively to the shock velocity. This is precisely the present situation with a jet of density 10⁶ cm⁻³ propagating into gas of density 2–4 × 10⁵ cm⁻³. This effect is seen in Fig. 8, where the number of peaks varies in the inner 10⁴au. As an example we follow the formation of the knot labeled Kbc: at t = 1450 yr Kc looks as a short and wide peak, that gets taller and thiner up to 1500 yr, when it appears to have merged with Kc. Immediatly after this, a small peak Kd forms behind it at 1550 yr, which moves forward at a slightly lower velocity, creating a very wide peak and giving a similar density distribution at 1450 yr and 1575 yr, an interval of 125 yr, which is approximately the period of variability used in the simulation. This process gets repeated several times, accelerating the fossil envelope material, as described in Sect. 4, relating the process of gas acceleration with the formation of knots.

We have reported in Fig. 9 the mass of the knots as a function of their location along the jet at t = 1500 yr. It appears that the knots tend to increase their mass by a factor of 2 as they move away from the injection site. The simulation shows that the mass increase occurs on a very short length scale as the second knot at z ~ 3000 au has already a mass of 1.7 × 10⁻³ M_⊙. Beyond 8000 au, all knots display rather similar and steady masses, ~l.7 × 10⁻³ M_⊙. The terminal knot is an exception due to its higher mass (>3 × 10⁻³ M_⊙). This high value is consistent with the fact that several knots have already reached the tip of the jet and are accumulating there, as can be seen in panels at t = 1500 and 2000 yr in Fig. 2.

Fig. 8 Gas density distribution along the jet main axis for different snapshots from 1450 to 1575 yr (dashed lines). Knots can be identified as overdensities in the density distribution (thick lines) in the order of 105 c−3. To follow their behavior some knots have been labeled as Ka. Kb, Kc and Kd.

Fig. 9 Mass of the knots detected at 1500 yr in Model Ml as a function of their position z along the jet, projected by 47°, which corresponds to the angle between the Cep E jet main axis and the plane of the sky. The dashed line corresponds to the mass injected to the simulation in a single period.

5.2.2 Observational comparison

In the later times of the simulation, it appears that the spatial separation between the knots located close to the source is actually shorter by about a factor 0.5 than the separation between those located in the older region of the jet. This is consistent with the bimodal distribution of dynamical ages reported by Schutzer et al. (2022).

The physical parameters of the knots (size, mass, density) such as derived from the numerical simulation, are therefore in good agreement with the observations. Schutzer et al. (2022) noticed that the distribution of knot dynamical ages appeared to be bimodal, with the presence of two timescales: a short time interval of 50–80 yr close to the protostar (8000 au), corresponding to the first 4–5 knots, and a longer timescale (150–200 yr) at larger distances from the protostar.

Interestingly, the numerical simulation reports the same trend, as can be seen in Fig. 10, in which the ages of all the knots identified along the jet have been reported. The first four younger knots are separated by a timescale of about 100 yr, whereas the following knots are separated by a timescale of about 200 yr, or twice the value measured close to the jet launch region. Our simulations show that a high-velocity knot can overcome the spatial separation to the previous éjecta, which was slown down as a result of the interaction with the dense gas of the protostellar envelope. The typical timescale for the collision is 500 yr (Fig. 10).

We propose that the bimodal distribution is actually the signature of the interaction of the knots with the high-density gas inside the shell. As soon as the knots run out of the shell and propagate into lower density gas, they encounter free motion conditions.

6 CO emission

Thanks to the detailed modeling of the CO chemistry by WALKIMYA-2D, we could produce synthetic CO integrated emissivity maps of the outflow jet and cavity, assuming LTE and optically thin emission. As an example, Fig. 11 displays the CO emissivity of the low-velocity outflow cavity integrated between 3 and 7 km s⁻¹ (black contours) and of the high-velocity jet, integrated between 50 and 150 km s⁻¹ (blue contours), as calculated by Model Ml at 1500 yr. The CO image reveals the previously identified highly collimated structure of the jet with a typical width of ~1000 au, which slightly increases with distance from the source. An important result is that the CO emission both from the high-velocity (50–150 km s⁻¹) jet (blue contours) and the low-velocity (3–7 km s⁻¹) cavity (black contours) drops ~2000 au before the terminal bowshock, whose location is marked by a red spot in the plot. The lack of low-velocity CO emission over a few 1000 au near the apex of the outflow cavity is consistent with the terminal bowshock being very efficient at dissociating the ambient molecular gas.

In order to gain more insight into the time variations of the CO abundance (relative to H₂) in the outflow, we have computed the CO abundance distribution along the high-velocity jet at 4 different times of the simulation in model M1: 0, 500, 1000, and 1500 yr. The distributions are shown in Fig. 12. For the sake of simplicity, we have reported the ratio X_CO of the CO abundance to its initial value, taken equal to 10⁻⁴, a standard value in dense interstellar clouds (Lacy et al. 1994).

In the very early stages of the simulation (t < 500 yr), when the jet has run only a few 1000 au, X_CO is well below its canonical value of 1.0. For instance, we measure X_CO ≈ 0.5 at 5000 au from the protostar at t = 500 yr. At later times in the simulation (1000 and 1500 yr in model Ml, see Fig. 12), two regimes are identified in the X_CO distribution observed along the jet: The first regime is close to the jet launch region X_CO = 1 The length of this region increases with time (2000 au at t = 1000 yr, 10 000 au at t = 1500 yr) as the jet propagates at velocity V_j = 140 km s⁻¹. The second regime, at a larger distance, X_CO decreases from 1 to a few 10⁻¹ with increasing distance to the protostar. The decrease is rather smooth and does not display abrupt variations along the jet. The minimum value is reached at the terminal bowshock before X_CO jumps abruptly again to 1, which stands as the canonical value in the ambient, quiescent gas in which the jet propagates.

To summarize, our simulation shows that CO is partly destroyed in the early times of the jet launch, close to the protostar (X_CO < 10⁻¹). The destruction process appears to stop ≈900 yr after the beginning of the simulation and the CO abundance of the ejected material remains steady, equal to its initial value. We propose that CO dissociation occurs as a consequence of the violent shocks caused by the first high-velocity knots (and the jet) when they are launched and impact the dense, protostellar envelope at the beginning of the simulation. It is only when the knots finally drill out the envelope and escape out of the inner protostellar region, entraining a fraction of the ambient protostellar envelope, that the local value of X_CO tends to return to its initial value. This process can also be influenced by the gas entrainment process.

Interestingly, we note that the slope of the X_CO distribution along the jet at t = 1000 yr and t = 1500 yr is similar, suggesting that it does not vary with time. Also, the velocity variations between subsequent knots (~ 10 km s⁻¹) do not appear to significantly affect the X_CO distribution. Therefore it seems that the jet kept the memory of the initial launch process in the X_CO axial distribution.

This numerical result has several implications on the observational side. First, lower values of X_CO may actually reflect the conditions of the jet formation process. Of course, precession and knot interaction with the ambient gas (or the cavity) may alter this conclusion. We speculate that in the latter case, the X_CO variations occur on a much smaller length scale, corresponding to the size of the knots, i.e. 1000 au typically. Our model suggests that the distribution of X_CO along the jet could probably provide constraints robust enough to discriminate between the processes at work.

On the observational side, Gusdorf et al. (2017) studied the emission of the OI 63 µm line at ~6″ resolution with the Stratospheric Observatory For Infrared Astronomy (SOFIA) and detected both the signature of the Cep E southern jet and the terminal bowshock HH377. A column density ratio N(OI)/N(CO) ~ 2.7 was measured in the jet, indicating that that the jet is essentially atomic in the region of HH377. The authors proposed that the OI emission could arise from dissociative J-type shocks or with a radiative precursor, caused by the knots propagating at a different velocity in the jet (Lehmann et al. 2020). Interestingly, no signature of OI was detected in the jet toward shock position BI, located halfway between the protostar and HH377.

As discussed above, our numerical simulations propose an alternative explanation for the origin of OI in the Cep E jet. A detailed map of the OI emission along the southern jet with SOFIA, would help to confirm whether the CO dissociation is localized only toward HH377 or if it is also present along the jet, hence could trace the history of the early stages of the mass-ejection process.

Fig. 10 Dynamical age of the knots identified in model Ml at t = 1500 yr in terms of a multiple number и¡ of the shorter dynamical age.

Fig. 11 Synthetic map of Model Ml at t = 1500 yr considering a projection angle with respect to the plane of the sky of 47°, as observed in Cep E-mm. CO J = 2-1 integrated emissivity contour maps of the cavity (Black contours, 3–7 km s−1) and the jet (Blue contours, 50–150 km s−1) components. Scaling is logarithmic. The response of the IRAM interferometer was modeled by a Gaussian of 830 au diameter (FWHM), corresponding to a beam size of 1″ (FWHP) at the distance of Cep E-mm, which is drawn with a gray disk (bottom left). The location of the frontal shock is marked with a red point.

Fig. 12 Model M1. Variations of the CO/H2 abundance relative to its initial value (10−4) along the high-velocity jet at 4 different times of the simulation: 0, 500 (dotted), 1000 (solid black), 1500 yr (solid red).

7 Conclusions

Using the reactive hydrodynamical code Walkimya-2D, which includes a chemical network based on CO (Castellanos-Ramírez et al. 2018), we have carried out a set of simulations in order to reproduce the morphology and the physical properties of the molecular jet-driven outflow of the intermediate-mass protostellar source Cep E, as derived from previous observational studies (Lefloch et al. 2015; Gusdorf et al. 2017; Ospina-Zamudio et al. 2019b; Schutzer et al. 2022). Outflow precession was not considered in this work. We have obtained a very satisfying agreement, both qualitatively and quantitatively, with the observations when modeling a time-variable jet with initial density n(H)= 10⁶ cm⁻³, radius r_j = 100 au, temperature T_j = 300 K, velocity V_j = 200 km s⁻¹ propagating into the density stratified protostellar envelope as modeled by Crimier et al. (2010). Our main results are as follows:

• The jet and cavity morphologies (width and length) are consistent with the observations and the expected kinematics using a variable jet ejection δV/V_j = 0.08. The best fitting solution is obtained at a numerical timescale t_dyn ≃ 1500 yr, consistent with the jet dynamical timescale estimated observationally (Schutzer et al. 2022, ~1400 yr).

• The jet terminal velocity (≃90 km s⁻¹) is different from the injection velocity V_j (200 km s⁻¹), as a result of the first ejections being decelerated by the dense envelope and the entrainment of a layer of ambient protostellar material. It implies that the jet dynamical timescale is actually lower than the duration of the ejection phase.

• The jet acceleration reported observationally by Schutzer et al. (2022) is consistently reproduced in the simulations. It appears to be the result of ambient material entrainment by the knots on a length scale of ~700 au, in agreement with the observational data.

• We reproduce the properties of the knots along the jet assuming a periodic ejection. We found evidence for knot interactions in the dense inner protostellar region, where the densities of the local gas and the jet are comparable, which lead to the formation of secondary shocks in the close protostellar environment. At larger distances from the protostar, the lower ambient gas density allows the knots to propagate freely, without significant interaction. We propose that this process could account for the bimodal distribution of knots observed along the Cep E jet. Knots have a typical size of 1000–2000 au, with a mass of ~1.5 × 10⁻³ M_⊙, and density of 10⁵ cm⁻³. The mass carried away by the knots in the jet translates into a steady ejection mass rate of 2.3 × 10⁻⁵ M_⊙ yr⁻¹, a factor of 3 higher than the mass injection rate into the jet (7.3 × 10⁻⁶ M_⊙ yr⁻¹). This difference is the signature of the entrainment and the subsequent acceleration of ambient material by the jet.

• The shock interaction of the jet knots with the protostellar envelope in the early times of the simulation lead to the dissociation of CO. The destruction process appears to stop after ≈900 yr, time from which the CO abundance of the ejected material remains steady and equal to its initial (canonical) value. As a consequence, the older part of the jet is characterized by a lower CO gas abundance, which gradually decreases as one moves closer to the jet head.

Our simulations underline the importance of mass-ejection time variability in the molecular outflow formation process and its interaction with the protostellar envelope. More work should be done, both observationally and numerically, in order to investigate the role of knots and their importance in the dynamical evolution of other young protostellar systems. Interferometric observations at subarcsec angular resolution of CepE and other young protostellar jets should be undertaken as they would allow to make a significant step forward by bringing extremely useful constraints on the knot internal structure and the dynamical processes at work in the jet. In parallel, the high spatial resolution accessible in the WALKIMYA-2D numerical simulations (~ 10 au) gives access to a novel view on the structure of protostellar jets and their dynamics, which will help interpret a new harvest of observations.

Movie

Movie 1 associated with Fig. 2 (CepE-2x) Access here

Acknowledgments

P.R.-R.O., A.S., and B.L. acknowledge support from (a) the European Union’s Horizon 2020 research and innovation program under the Marie Skłodowska-Curie grant agreement No 811312 for the project “Astro-“Chemical Origins” (ACO); (b) the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program for the Project “The Dawn of Organic Chemistry” (DOC) grant agreement No 741002; (c) the UNAMPAPIIT grant IN110722. Some of the computations presented in this paper were performed using the GRICAD infrastructure (https://gricad.univ-grenoble-alpes.fr), which is partly supported by the Equip@Meso project (reference ANR-10-EQPX-29-01) of the programme Investissements d’Avenir supervised by the Agence Nationale de la Recherche and the Miztli-UNAM supercomputer, project LANCADUNAM-DGTIC-123 2022-1.

References

All Tables

All Figures

Fig. 1 CO 2–1 line emission in the Cep E-mm outflow as observed with the PdBI at 1″ resolution (Lefloch et al. 2015). Four main velocity components are detected: a) the outflow cavity walls (black) emitting at low velocities, in the range [−8;−3] km s−1 and [−19;−14] km s−1 in the northern and southern lobe, respectively; b) the jet, emitting at high velocities in the range [−135;−110] km s−1 in the southern lobe (blue) and in the range [+40;+80] km s−1 in the northern lobe (red); c) the southern terminal bow shock HH377, integrated in the range [−77,−64] km s−1 (green); d) the northern terminal bullet NB integrated in the velocity range [+84,+97] km s−1 (magenta). First contour and contour interval are 20 and 10% of the peak intensity in each component, respectively. The synthesized beam (1″.07 × 0″.87, HPBW) is shown in the bottom left corner. The main axis of the northern and southern lobes of the jet are shown with black arrows. (Schutzer et al. 2022).

In the text

	Fig. 2 Outflow formation and impact of the jet on the gas density distribution as a function of time. Snapshots for 0, 500, 1000, 1500, and 2000 yr. The injection velocity variability forms internal working surfaces, or knots, and a structured cavity. The leading bowshock accumulates the mass of several knots. A movie is available online.
In the text

Fig. 3 Cep E northern outflow lobe. Position-Velocity diagram of the CO J = 2−1 line emission along the jet main axis. Positions are in arcsec offset relative to the location of the driving protostar Cep E-A. First contour and contour interval are 10% of the peak intensity. We have superimposed in yellow the best fitting solution (V − Vlsr)/V0 = exp(−δ0/δ) where δ0 = 690 au (Schutzer et al. 2022).

In the text

Fig. 4 Best fitting model (Ml) to the Cep E northern outflow, (top) Density distribution at 1250 yr. The jet has propagated over a distance of 2.4 × 104 au and formed a low-velocity cavity of 3000 au radius. Five knots and their associated wings are detected along the jet.(bottom) Position-Velocity diagram of the CO emission along the jet main axis obtained using the 47° outflow inclination angle with respect to the plane of the sky.

In the text

	Fig. 5 Cep E northern outflow lobe. Position-Velocity diagram of the CO J = 2−1 line emission across the jet main axis at δ = 0″.Positions are in arcsec offset relative to the location of the driving protostar CepE-A. First contour and contour interval are 10% of the peak intensity (Schutzer et al. 2022).
In the text

	Fig. 6 Model M1. Outflow mass evolution as a function of time.
In the text

	Fig. 7 Model M1. Jet mass evolution as a function of time. The points are obtained from the simulation every 100 yr. The black continous line is a linear fit of the points starting from 300 yr, with a slope of 2.3 ×10−5M⊙ yr−1 and the dashed line corresponds to the material injection rate of Ṁ = 7.3 × 10−6 used in the simulation.
In the text

	Fig. 8 Gas density distribution along the jet main axis for different snapshots from 1450 to 1575 yr (dashed lines). Knots can be identified as overdensities in the density distribution (thick lines) in the order of 105 c−3. To follow their behavior some knots have been labeled as Ka. Kb, Kc and Kd.
In the text

	Fig. 9 Mass of the knots detected at 1500 yr in Model Ml as a function of their position z along the jet, projected by 47°, which corresponds to the angle between the Cep E jet main axis and the plane of the sky. The dashed line corresponds to the mass injected to the simulation in a single period.
In the text

	Fig. 10 Dynamical age of the knots identified in model Ml at t = 1500 yr in terms of a multiple number и¡ of the shorter dynamical age.
In the text

Fig. 11 Synthetic map of Model Ml at t = 1500 yr considering a projection angle with respect to the plane of the sky of 47°, as observed in Cep E-mm. CO J = 2-1 integrated emissivity contour maps of the cavity (Black contours, 3–7 km s−1) and the jet (Blue contours, 50–150 km s−1) components. Scaling is logarithmic. The response of the IRAM interferometer was modeled by a Gaussian of 830 au diameter (FWHM), corresponding to a beam size of 1″ (FWHP) at the distance of Cep E-mm, which is drawn with a gray disk (bottom left). The location of the frontal shock is marked with a red point.

In the text

	Fig. 12 Model M1. Variations of the CO/H2 abundance relative to its initial value (10−4) along the high-velocity jet at 4 different times of the simulation: 0, 500 (dotted), 1000 (solid black), 1500 yr (solid red).
In the text