© The Authors 2026

1 Introduction

Hypersonic collimated jets, surrounded by a lower velocity outflow, are commonly observed in Class 0/I objects, corresponding to the earliest evolutionary phases of protostars (Ray & Ferreira 2021). In the past decades, understanding the formation of protostellar outflows, the role of jets, and their interaction with the ambient medium has been a real challenge. Progress in 2D and 3D hydrodynamical simulations (Rabenanahary et al. 2022; Rivera-Ortiz et al. 2023) has clearly supported, both morphologically and dynamically, the jet-driven bow-shock mechanism initially proposed by Raga & Cabrit (1993). These models consider isolated forming protostars, whose jets may be morphologically modified by the interaction of a binary component that adds a velocity-component tangent to the orbital motion (de A. Schutzer et al. 2022) or precession (Rivera-Ortiz et al. 2023) as a possible consequence of the interaction of binary systems. However, source multiplicity in protostellar systems appears to be very common (see, e.g., Reipurth et al. 2014 for a review). For instance, the Small Millimeter Array (SMA) study by Chen et al. (2013), a sample of 33 Class 0 protostars, all located less than 400 pc away, showed that more than two-thirds of the source sample displayed evidence of a companion at a distance of 50–5000 au. Similar results were obtained consistently toward BHR71, Cep E, or IRAS4A, for example (Zapata et al. 2018; Ospina-Zamudio et al. 2018; Taquet et al. 2020; Toledano-Juárez et al. 2023). Hence, the possibility that several outflows could propagate inside a single parental protostellar envelope and interact with the ambient gas or other outflows can no longer be discarded or considered a low-probability phenomenon.

As intermediate-to-high-mass stars form in dense clusters (e.g., Orion with ∼106 stars pc⁻³; Rivilla et al. 2014), dynamical interactions among cluster members may be very common (Bonnell & Bate 2005). In the case of massive protostars with circumstellar disks and envelopes of a few thousand astronomical units in size, their impact areas should increase dramatically, producing frequent interactions at the center of the cluster. The rich cluster of outflows detected in the CO 2–1 and SiO 5–4 lines in the W43-MM1 high-mass star forming region by Nony et al. (2020, 2023) (see their Fig. 1), as part of the ALMA-IMF Large Program, suggests that multiplicity is present in the 27 protostellar cores of the cluster; there, multiple mass-ejection phenomena interact on a 5 × 10³ au scale. Recently, observations on a 10³ au scale have become commonly accessible with the large (sub)millimeter facilities ALMA and NOEMA, and their new instruments reveal faint structures within the outflow; so, only now can observational constraints be put on flow–jet interactions from multiple sources in a common protostellar core (De Simone et al. 2022; Cortes-Rangel et al. 2023; Durán et al. 2025; Arazi et al. 2026).

The present study aims to shed light on the physical processes involved in the collision of outflow–jet envelopes in protostellar binary systems and to propose observational signatures that could provide empirical evidence of such events. However, numerous free parameters – such as the physical properties of the outflows, variability, and the angle between interacting outflows – depend solely on the specific object under study. Therefore, we focused on the interaction between the molecular outflows ejected by the Cep-A and Cep-B components in the Cep E-mm system. Cep-A and Cep-B form a binary protostellar system, separated by approximately 1000 au, with supersonic molecular outflows and velocities on the order of ∼100 km/s (Lefloch et al. 1996; Ospina-Zamudio et al. 2018; de A. Schutzer et al. 2022). The collision of these molecular outflows could be reflected in the morphology of the system structure.

In this paper, we have used a uniform-mesh 3D hydrodynamical code (Rodríguez-González et al. 2019), that accounts for the cooling rate at low temperatures, to model the collision of the jets powered by these forming protostars within a dense and stratified protostellar envelope. We adopted initial conditions for the jet and the parental core based on the observational studies of Class 0 protostellar core Cep E-mm (de A. Schutzer et al. 2022; Lefloch et al. 2015; Crimier et al. 2010), which can be considered as a prototype for such observational and numerical studies (Rivera-Ortiz et al. 2023).

Our article is organized as follows. Section 2 briefly describes the astrophysical source used to determine the initial jet and protostellar conditions in the numerical simulations. Section 3 describes the hydrodynamical code used to carry out the numerical simulations and the physical constraints we adopted for the modeling based on the CepE observations. The results of our simulations are discussed in Section 4; we first show how the dynamics of both jets and their outflows can be strongly perturbed under specific jet dynamical conditions, both on small and large scales. Interestingly, the numerical signatures of the interaction are consistent with the observational study of de A. Schutzer et al. (2022). We also present the various observational signatures to be expected based on preliminary exploration of the outflow’s corresponding parameter space. Finally, our conclusions are presented in Section 5.

2 Cepheus E

Cepheus E-mm is the prototype of an intermediate-mass star forming core located in the Cepheus molecular cloud at a distance of d = 820 pc (de A. Schutzer et al. 2022), with a bolometric luminosity of L_b ∼ 100 L_⊙ (Lefloch et al. 1996). Its millimeter emission peak is centered at α(J2000)= 23^h03^m12.8^s, δ(J2000)= +61 ^o42′26″. The Cep E-mm core is a binary system of two protostars, Cep E-A and Cep E-B. According to Ospina-Zamudio et al. (2018) they have a separation of 1.35″ (1000 au). They both power molecular jets with velocities up to 80–120 km s⁻¹ detected in the millimeter lines of various molecular species, such as CO and SiO, as can be seen in Fig. 1 of Ospina-Zamudio et al. (2018). The authors also determined the physical properties of components A and B. Based on 1.3 mm continuum observations, they showed that CepE-A has a size of 875 au, a mass of 0.56 M_⊙, a 1.3 mm total flux of 169 mJy, and a dust temperature of T_d = 60 K, while CepE-B has a size of 670 AU, a mass of 0.13 M_⊙, a 1.3 mm total flux of 40 mJy, and the same dust temperature of T_d = 60 K.

The main axis of jets Cep E-A and Cep E-B are close to 90^◦, almost perpendicular to one another. As can be seen in Fig. 1, the high-velocity jet powered by source A propagates at (radial) velocities of about +65 km s⁻¹ and −110 km s⁻¹ in the northern (red contours in Fig. 1) and southern (blue contours in Fig. 1) lobes, respectively, regarding the ambient cloud gas (V_lsr= –10.9 km s⁻¹). In comparison, source B and its jet have been much less studied. Nevertheless, the jet’s radial velocities are on the same order of magnitude as those observed in the jet from source A, with typical values of ±70 km s⁻¹ in both red-shifted and blueshifted lobes regarding the ambient cloud (see right panel in Fig. 1).

Close to the driving protostar, A, a spatial overlap is observed in the high-velocity emission of the A and B jets. The blueshifted component of the high-velocity B jet has propagated over a distance of about 5000 au westward of the protostar. In contrast, the redshifted component seems to have propagated about 3300 au eastward. The spatial distributions of the outflows appear to overlap, at least in projection, in the plane of the sky. Near to source B in Fig. 1 (left panel), there is a decrease in CO emission (by 4″) from the jet at velocities of >50 km s⁻¹ (blue contours), which is not seen for the northern jet. We suggest that this reduction in CO emission occurs due to an interaction between the primary jet and a second collimated flow, where the pressure gradients are the forces changing the individual outflow dynamics. Observationally, this effect manifests in projection as a decrease in the velocity in the original plane of motion or a smaller projected size for this secondary outflow. In the right panel of Fig. 1, one can see the flow ejected by source B, and it could be because this flow moves in one direction (to the right in this image), but it does not evolve in the same manner in the opposite direction. It can be because of an interaction of the jet (launched by source B) with the medium through which it propagates.

In addition to the overlap of the CO and SiO signatures of the jets and their outflows, de A. Schutzer et al. (2022) showed that the redshifted component of the B jet overlaps with a minimum of emission in the southern low-velocity cavity walls of the A jet. The detection of SiO, CO, CS, and H₂CO (Fig. 2) attests to the presence of a somewhat “active” shock, and the most likely hypothesis is that the B jet is impacting the southern low-velocity cavity driven by jet A. Furthermore, there is a 7^◦ misalignment between the northern jet axis and the southern jet axis that can be seen in Figs. 1 and 2 (left panel), where we note an asymmetric structure in the cavity (de A. Schutzer et al. 2022). Figure 2 also shows the flow of source B with the cyan arrow. This is one of the most convincing cases of jet interaction in a protobinary system; that is, where there is evidence of the force produced by the ram pressure gradient of the envelope of one jet onto the other.

Next, we estimated the probability of the two outflows coinciding spatially. First, we approximated that both outflows, produced from sources A and B, had a biconic shape. An observer from source B would see two spherical triangles projected in the sky from their point of view, with a height of L cos i and base o L cos i tan α, where i is the outflow-A inclination angle with respect to the plane of the sky, and α is the A-cone opening angle. Then, assuming source A as the pole of a spherical system and at a distance, r, from the origin of this system at source, the solid angle covered by a triangular section could be approximated as Ω = ϕ(1 − cos θ), where $\tan \theta =\frac{L\cos i}{D(1+\sin i)}$ and $\tan \varphi =\frac{L\cos i\tan \alpha }{D(1+\sin i)}=\tan \theta \tan \alpha$. Additionally, the outflow produced by B would make a solid angle of ω = 2π(1−cos β), where β is the opening angle of the B cone. Finally, the probability can be estimated as the ratio of the solid angle subtended by the outflows, Ω₁ + Ω₂ + ω, and the solid angle of a semisphere, 2π: $$P=\frac{{\text{Ω}}_1+{\text{Ω}}_2+\omega }{2\pi }.$$(1)

In Fig. 3, we estimate the probability of an encounter using opening angles of 90^◦ and 15^◦ for the Cep E-A cavity and jet, respectively, and opening angles of 60^◦ and 15^◦ for the cavity and jet produced by source B, respectively. We used an extension of L = 8000 au, where these opening angles seem reasonable, and the separation of the sources is a free parameter, since they have an apparent separation of 2136 au (dashed line), which could be larger due to projection effects. The lines represent the probability of the encounter between cavities (black) and jets (blue). This model still underestimates the collision probability, P ∼ 80%, which is already high in the case of the cavity–cavity encounter since we did not take into account their 3D structure. Therefore, we expect the cavities associated with the Cep E binary system are spatially overlapping, while their driving jets are not necessarily colliding (P ∼ 20%).

Fig. 1 Left: CO J= 2–1 emission detected in the Cep E protostellar core with the IRAM interferometer at 1″ (820 au) resolution. The black contours trace the emission from the low-velocity gas (|V −Vlsr| < 8 km s−1), and the blue/red contours trace the high-velocity jet (V − Vlsr > 50 km s−1, powered by component A of the Cep E protostellar core. The data are from Lefloch et al. (2015) and de A. Schutzer et al. (2022). Right: SiO J= 5–4 emissions detected in the Cep E protostellar core with the IRAM interferometer at 1.4″ (1150 au) resolution in the jets powered by protostars A (red and white contours) and B (blue and red contours). The velocity intervals of integrated flux are reported (in km s−1) next to the lobes of the jets from A and B. The first contour and contour interval are 10% and 20% of the maximum peak intensity. The data are from Ospina-Zamudio et al. (2018).

Fig. 2 Interaction between southern low-velocity outflow lobe of A, at [−8; −6] km s−1 (black contours) and the eastern lobe of the jet from B (red contours), as observed with the IRAM interferometer in the molecular lines of CO J=2–1 and CS J=2–1 (left), and SiO J=2–1 and H2CO (right).

Fig. 3 Probability of two outflows overlapping for two wide cavities (continuous black line) and for two collimated jets (blue line). The dashed line represents the observed separation of Cep E-A and Cep E-B, which could be larger because of projection effects.

3 Numerical model

To describe the collision between the envelopes and the main flow and characterize the observational signatures produced by this interaction, we selected Cep E to obtain a significant number of physical constraints required to produce a consistent model based on the observed physical conditions of Cep E–South, which may lead to a strong interaction between molecular outflows capable of altering the evolution of this object. We also searched for the expected numerical signatures at length scales accessible to the (sub)millimeter interferometers ALMA and NOEMA in the largest outflows. Cep E is a protobinary system, where each component produces a high-velocity outflow (see Ospina-Zamudio et al. (2018)). As discussed before, there is evidence of interaction between the powerful outflow produced by the component Cep E-A and the jet produced by the component Cep E-B. Also, we considered that each material injection zone is embedded in a density stratification associated with a protostellar core, as modeled by Crimier et al. (2010). Then, we modeled the collision of the outflows produced by two high-velocity jets moving in this stratified environment. This section describes the numerical model. First, we explain the code used to perform our simulations. Then, we give the physical constraints obtained from the Cep E protostellar source.

3.1 Numerical setup

The presented grid of numerical simulations (Table 1) were carried inside a physical domain [x, y, z] ∈ [(0, 1), (0, 4), (0, 1)] × 10⁴au with a spatial resolution of 4.88 au in each direction. The numerical box was a rectangular prism, and, considering that the origin of the Cartesian plane lies at the center of a square face of the box, with coordinates [0, 0, 0], we injected high-velocity gas at two positions into the grid to reproduce the propagation of two protostellar jets, which are similar to those observed in Cep E.

The first jet was the main one. It was launched from position [0, 0, 0]. The other injection zone represented the secondary jet, which we injected at the position [1500, 1500, 250] au.

We performed our 3D hydrodynamic simulations using the GUACHO code¹ (Esquivel et al. 2009). This code solves the near-conservation laws governing the gas dynamics and the neutral hydrogen rate in a Cartesian grid. The complete set of equations is $$\frac{\partial \rho }{\partial t}+\nabla \cdot (\rho v)=0,$$(2) $$\frac{\partial \rho v}{\partial t}+\nabla \cdot (\rho v\cdot v+p)=0,$$(3) $$\frac{\partial e}{\partial t}+\nabla \cdot [v(e+p)]=G_{\text{rad}}-L_{\text{rad}},$$(4)

where, ρ, v, p, and e are the density, the velocity, the thermal pressure, and the total energy density of the gas, respectively. G_rad and L_rad are the energy gains and losses due to radiation. We assumed a standard ideal gas law for the closure of the gas total energy density, and the gas thermal pressure is related by e = ρv²/2 + P/(γ − 1), where P is the total pressure and the ratio between specific capacities is γ = 1.4. The code is described in detail in Rodríguez-González et al. (2019). In our simulations, energy was lost via radiative cooling, which was included as described by Raga & Reipurth (2004) for the atomic gas, which reached temperatures of thousands of kelvins. For lower temperatures, we included the parametric molecular cooling function presented in Kosiński & Hanasz (2007) and Rivera-Ortiz et al. (2019). We also considered heating resulting from the hydrogen’s cosmic-ray ionization rate proposed in Henney et al. (2009). The two protostars are separated by a distance of 2136 AU, which makes the Cep E outflow age a fraction of the orbital period of the system, and the interaction between their supersonic jets occurs at large distances from both protostars, where the escape velocity is much lower than the outflow velocities. With this configuration, gravity and self-gravity were not considered in the models investigated in this work. Additionally, in low-mass star-forming regions, magnetic pressure is significantly weaker than thermal and ram pressure in the supersonic outflow at a distance more than thousands of astronomical units from the protostellar disk, since a large faction of magnetic energy is already converted to kinetic energy and the magnetic field is on the order of 10⁻⁵–10⁻³ G for similar jets in the pre- and post-shock gas, respectively, in other HH outflows (Pudritz & Ray 2019). The effects of gravity, magnetic fields, and even the radiation on the accretion disk and outflow launching have been extensively studied in previous works (e.g., Meliani et al. 2006; Mignon-Risse et al. 2021b,a; Oliva & Kuiper 2023). While on large scales the magnetic pressure gradient is small, and protostellar jets at early stages are magnetic in origin, we focused our analysis on the effects of the ram pressure gradient of one jet onto the other as a possible cause of the asymmetries observed.

Fig. 4 Two protostellar sources inject high-velocity material. The gas density around each source was stratified with a power law. In our simulations, the jet injections were not meant to collide if they preserved their width, but, since they expand thermically and by entrainment, the evolved jets – that is, high-velocity material – could collide.

3.2 Physical constraints

In Fig. 4, we used one jet (the primary jet), which was propagating in the y-direction and a secondary jet that propagated in the x-direction. We decided to avoid a direct collision of the two jets, that is, the crossing of the two jet axes, since the direct collision has a lower probability of occurrence; however, as they evolved, they could still interact. Both propagation directions were perpendicular in the xy-plane, but they had an impact parameter of 250 au in the zy-plane, so the material injected by the secondary jet came out of the primary cavity with more momentum and had a strong impact on this cavity.

Both of the jets were considered as cylinders with a radius of r_J=63 au and an injection length of r_L=42 au, and both jets had the same thermodynamic properties; that is, a numerical density of n_j = 1 × 10⁶ cm⁻³ and a temperature of T_j = 250 K. The primary jet was injected using a time–velocity variation, as was proposed in Rivera-Ortiz et al. (2023): $$V_j=V_{j,m}[1+\delta \upsilon \cos (2\pi t\text{}/\text{}\tau )],$$(5)

where V_j,m is the mean jet-injection velocity = 200 km s⁻¹, and the relative amplitude variability is δv_j = 0.08, which implies velocity variations of 16 km s⁻¹ and τ = 130 yr, which is twice the period determined by de A. Schutzer et al. (2022). The mass injection can be calculated using (Rivera-Ortiz et al. 2023) $$\frac{\dot{M}}{{10}^{-6}{\text{M}}_{\odot }{\text{yr}}^{-1}}=1.5{\left[\frac{r_j}{50\text{\hspace{0.17em}au}}\right]}^2\left[\frac{n_j}{{10}^6{\text{cm}}^{-3}}\right]\left[\frac{v_j}{165\text{\hspace{0.17em}km}{\text{s}}^{-1}}\right]=\mathrm{2.88.}$$(6)

The secondary jet was turned on t_s = 700 yr after the primary jet was initially ejected. This timing was determined based on the estimated age of the flow and the location of material dispersion in the southeastern low-velocity cavity, which is located about a half distance considering the jet length and a dynamical age of 1500 yr (de A. Schutzer et al. 2022; Rivera-Ortiz et al. 2023). We changed this value to inject the secondary jet at different distances from the main jet’s wall cavity. The velocity of the secondary jet was V_j,s = 100 km s⁻¹. We define $$\begin{array}{l}r={[x^2+y^2+z^2]}^{1\text{}/\text{}2}\hfill \\ r_b={[{(x-x_b)}^2+{(y-y_b)}^2+{(z-z_b)}^2]}^{1\text{}/\text{}2},\hfill \end{array}$$

as the radial distances to the primary source, r, and to the secondary source, r_b. Then, considering both envelopes, the initial gas density and temperature distributions were imposed as follows: $$n(r)=\frac{n_0}{1+{(r\text{}/\text{}{r}_0)}^{{\alpha }_1}}+\frac{n_0}{1+{(r_b\text{}/\text{}{r}_0)}^{{\alpha }_1}},$$(7)

which is the profile proposed by Crimier et al. (2010). Considering that each term on the right-hand side of this equation represents the density profile corresponding to each source, with a parameter α₁ = 1.9, this profile corresponds to a gas mass of 0.3 M_⊙ up to 500 au from the central source and a gas mass of 0.14 M_⊙ up to 300 au from the secondary source. This value is in reasonable agreement with the mass reported by Ospina-Zamudio et al. (2019) and the mass–luminosity relationship for young stellar objects (Motte & André 2001).

The region also had a temperature profile given by $$T(r)=\frac{1}{2}\left[\frac{T_a}{1+{(r\text{}/\text{}{r}_0)}^{{\alpha }_2}}+\frac{T_b}{1+{(r_b\text{}/\text{}{r}_0)}^{{\alpha }_2}}\right],$$(8)

where T_a = T_b = 300 K was the central temperature of each core and α₂ is −0.8 in the range between 50–300 K and −0.4 in the range 7–50 K, according to Crimier et al. (2010).

Fig. 5 Upper panel: snapshot of gas density of a single jet (model M1) at t = 1875 yr. One can see the work-surface structure due to the consideration of the variable jet. The envelope shows a higher density at the “base” of the jet and a lower one toward the front. Lower panels: snapshots of gas density of the jet-interaction model (model M2) at three evolutionary times: t=1375, 1625, and 1875 yr (top, middle, and bottom panels). Considering the estimated age of the Cep E protostellar outflow of ∼1500 yr (de A. Schutzer et al. 2022), different stages of evolution around this age are shown. All the snapshots are present on the xy-plane.

4 Results

We carried out two hydrodynamical simulations to analyze the effect of the interactions of the two jets in the dynamics and morphological evolution of the gas injected by the protostellar objects. The first model (M1) had a single jet where the material is injected with a velocity variability, described by Eq. (5), into a stratified interstellar medium, allowing its dynamic evolution. We allowed the material to evolve for 1875 yr. That is, there are more than 14 periods of maximum velocity in the injection of the jet. In the second model (M2), we considered the evolution of a single jet for 660 yr, and we turned on the second jet, which propagated perpendicularly regarding the first of them. We have followed the gas evolution since the secondary jet reached the limit of the simulation box (also at 1875 yr).

4.1 Morphology evolution

The morphology of the single-jet (model M1) evolution is presented in Fig. 5 (top). There is a snapshot of the central slice of the density distribution of the M1 model, where only the evolution of a single jet has been considered. As expected, in the density, the working surfaces resulting from the velocity variation imposed on the jet’s injection zone and the material environment’s envelope are shocked by shock waves expanding laterally with subsonic velocities and low Mach numbers. In the region closest to the material injection point, a very dense envelope is visible. From this figure, one can see how the density on the outer wall is greater at the base of the jet; this is because there is a much denser environment in this region, and we have seven knots, or visible working surfaces, that structure the gas in the cavity from the center. For this model, one can see that the envelope reaches a width of 6000 au, while the length is approximately 30 000 au.

We also ran a simulation considering the mass and energy input rates of two different jets. We are interested in evaluating the effect of a secondary jet (which we considered continuous) on the evolution of a variable jet that we consider the primary one. We added the secondary jet, assuming that it originated within the envelope formed by the supersonic movement of the envelope towards the interstellar medium, where it was propagating. Figure 5 shows the central-slice density evolution of the gas in model M2 (color maps), with two winds colliding. The upper, middle, and bottom panels present the evolutionary time 1375, 1625, and 1875 yr, and one can see the secondary jet position in [(3500), (5250)] au in this image. One can also see the evolution of the gas injected by the secondary jet at 1375 yr (top panel) in the disturbance of the main envelope where the gas from the secondary jet impacted the wall (bottom wall in this figure). Moreover, at 1375 yr, the frontal part of the envelope does not seem to have significant changes; however, as the interaction continues, the envelope deforms, and this perturbation grows, creating a considerable asymmetry on the same side where the secondary jet was injected. The cavity reaches a width of about 8500 au at 1875 yr, which is about 30% wider than the cavity formed in single-jet evolution.

Note that the disturbance in the primary jet envelope is related to the first moment at which the two jets interact, a relic of the collision. The interaction between the material coming out of the secondary jet that impacts and intervenes in the flow of gas that makes up the primary jet changes the general direction of the flow. The size of this region of impact increases, forming a bubble; the evolutionary times 1625 and 1875 yr are presented in the middle and bottom panels, respectively. At 1875 yr (bottom panel), a bubble of about 2000 au forms, containing material injected by the primary and secondary jet. The densest part of the shell of the bubble forms because of the supersonic collision of the jet and original cloud gas. Therefore, the gas envelope formed by the thermal expansion of the envelope has a faster evolving region due to the kinetic energy injected by the secondary jet into this cavity. The primary jet outflow suffers a deviation in its trajectory of approximately five degrees. The deviation of the flow has its most significant effect in the interaction with a part of the cavity that moves at much slower velocities than the front part of the shell. In this region of impact, the envelope has a lateral movement of a few km s⁻¹. The collision between the main jet and the inner part of the envelope has an essential effect on the overall evolution of the bubble. The surface where the jet hits partially reflects it, making its movement conditioned by the shape of the wall where it slides.

The cavity formed by a single jet appears to have a lower density because much of the injected material remains on the working surfaces, where two shocks are formed: one toward the source and the other away from the source; these are created by the variability of the jet. The lateral growth induced by the secondary jet can be observed as jet envelopes much larger than those expected in the single-jet models. The jet beam also appears almost three times wider, and the knots or working surfaces are broader than the primary jet flow. In the scenario with two jets, the working surfaces appear much less dense or compact and, in some cases, seem to have dissipated (see Fig. 5, bottom panel), creating a region at the front of the cavity that appears to expand due to internal pressure.

4.2 Force analysis

To obtain the effect of the interaction between molecular outflows and to compare it with other forces that may play an important role in the evolution of young outflows, we computed the material derivative of the velocity matrices under the assumption of a stationary state. This procedure allowed us to obtain acceleration matrices along each axis, from which the force per unit volume can be derived using Eq. (3).

Figure 8 presents the resulting force field, where arrows indicate the 3D direction of the force, and a color-scale encodes its magnitude. In both panels, the figure illustrates that the force due to the gas flow is strongest in the jet-injection region, at the knots or working surfaces formed by the temporal variability of the primary jet velocity, and at the leading edge where the jet encounters the ambient gas cloud. Both the upper and lower panels also reveal the expansion of the envelope, albeit with lower forces. Notably, this expansion is more pronounced up to two knots behind the leading knot, highlighting that the future evolution of this object will be highly asymmetric in its frontal region.

From the yx-plane (upper panel), the interaction between the flows further reveals (a) a change in the force of the secondary flow along the direction of propagation of the primary jet; (b) the formation of a protrusion or secondary bubble generated by the momentum injected by the secondary source; and (c) broadening forces that affect the evolution of the primary jet. From the yz-plane (lower panel), one can observe that the secondary flow develops an expansion triggered both by its interaction with the primary jet and the resulting force from the combined flows.

We note that the deformation of the primary jet is caused by the ram pressure gradient associated with the secondary jet, which introduces large transverse velocity gradients, producing a net transverse force – that is, an advection – described by the term of Eq. (3) ∇ · (ρv · v), breaking the symmetry of the primary jet. This is consistent with Fig. 8, where the vectors show that the deformation of the primary jet correlates with the jets closest distance region.

To quantify the effect of the magnetic field with respect to the supersonic motion of the two flows present in Cep E, we can calculate the magnetic pressure, P_B, with the ram pressure, P_ram; one can be calculated as $P_B=\frac{B^2}{4\pi }$, where B is the magnetic field, while the ram pressure is P_ram = µm_hnv². The force produced by each pressure is f = ∇P, which we approximated to $f=\frac{P}{\text{Δ}x}$, assuming the pressure changes in a scale, ∆x, that is, the cell size in our simulations. For values of B = 1 mG, v = 100 km s⁻¹, and n ∼ 10⁶ cm⁻³, the ratio $\frac{f_{\text{ram}}}{f_B}\simeq {10}^{-4}$, which helped us to scale these forces. For lower velocities of around a few km s⁻¹, both forces are comparable, where B is as large as 1 mG, but typically different values, away from the protostellar core, of 10⁻⁶–10⁻⁴ G (Pudritz & Ray 2019), make the magnetic force very small compared with the ram pressure. This value for B is derived from the analysis by Lefloch et al. (2015), in which a grid of models was created to reproduce the CO emission in the Cep E jet shocks.

4.3 Jet deviation and pseudo-precesion

The collision between the outflows appears different when observed from another reference axis. For instance, the yz-plane (Fig. 6) shows the evolved interaction between the outflows. In the region of the first knot close to the source, the secondary jet can be noticed as a dense circular dot above the principal jet. On the other hand, the flow diversion caused by the interaction between the outflows is noticeable because the nodes move away from the center of the box where the injection is generated, creating the impression of descending.

On the other hand, to characterize the deviation of the main flow, we placed some boxes centered on the barycenter of each knot produced by the working surfaces. As seen in this plane, that is, the yz-plane, one can observe how the main flow is deflected due to the moment injected into the secondary flow. The secondary jet was injected at 250 au over the plane of the primary jet. Therefore, the main flow deviated in the z-direction, producing an asymmetric envelope evolution. This asymmetry in the evolution of the envelope is driven by the flows of the main jet, which interact with the inner part of the envelope, thereby changing the flow direction. Furthermore, one can evaluate the asymmetric evolution of the nodes and observe variations among the different planes. Figure 9 shows the position of the knots via green and pink points in the yx- and yz-planes, respectively, of the M2 simulation. The yz-plane shows that the global deviation of the jet is 3.5 degrees in the opposite direction to the injection point of the secondary jet. However, in the xy-plane, the deviation concerning the x-axis is negative for the first five nodes (away from the direction of the secondary source) and positive for the last two. This false precession is caused by the interaction between the jet and the internal wall of the cavity on the opposite side of the injection due to the secondary jet. The maximum deflection due to this process is about 2.3 degrees in the opposite direction to the secondary jet and 5.11 degrees in the other direction. In addition to the model M2 jet being wider than the model M1 jet by approximately 15% (see Fig. 5), it was observed with a precession effect that makes it appear much more expansive than in the case of single jet evolution.

At t = 1875 yr, the total mass injected by the primary jet in models M1 and M2 is ∼0.005 M_⊙ (see Eq. (6)). The redistribution of this material within the computational domain is illustrated in Fig. 7. More than half of the injected mass is concentrated in the working surfaces (blue crosses), while a dense bow shock forms at the head of the cavity, corresponding to the swept-up ambient medium that receives most of the jet momentum. The mass distribution along the jet axis shows that the internal knots have comparable masses, typically ∼4 × 10⁻⁴ M_⊙, and together they account for about 70% of the total injected mass. Approximately 10% of the jet material is found in the surrounding cavity, which has been carved out by the jet and filled with material originating from the jet itself. This matter has been expelled laterally from the beam due to the action of internal shocks propagating within the jet. The remaining ∼20% of the mass is distributed along the jet beam between the knots.

This contains material from the environment swept up during the entire evolution of the jet. However, the total injected mass (cumulative mass) does not agree with the injected one. It is likely that a large amount of injected material, about 30% in M1, is not in the knots but in the cavity dug by the jet. The last knot, the one furthest from the injection zone, is much more massive than the internal one because a large amount of the environment has been swept away by the bow shock formed by the supersonic evolution produced by this jet.

Fig. 6 Snapshots of gas density of two-jet collisional model (model M2) in the Y–Z plane. The color represents the density, and the white boxes show the position of the knots. Using the boxes, we can find the mass of the knots.

Fig. 7 Mass of the knots of models M1 (blue circles) M2 (red circles) and the accumulated mass are shown by crosses. The figures show that the masses in the knots of the model with two-jet interaction were systematically greater for all knots than those of the M1 model, except for the furthest knot, which was difficult to determine.

4.4 Gas dynamics

To analyze the motion of the gas within the envelope in detail, we plot velocity fields to the density maps for the numerical simulation at 1875 yr; see Fig. 10. However, to analyze the dynamic evolution of this gas in detail, we divided maps with a velocity field into three velocity channels (in Fig. 10): (a) 5–15 km s⁻¹ (top panel), (b) 15–40 km s⁻¹, and (c) >40 km s⁻¹. The lower velocity channel mainly shows the lateral expansion of the envelope, and the higher velocity channel shows mostly the motion of the gas injected by the jets. Therefore, the channel of velocities between 15 and 40 km s⁻¹ shows the gas dynamics inside the cavity.

In Fig. 10a, in the lower velocity channel, at <7 km s⁻¹ we can see how the envelope expands laterally, moving against the surrounding environment. In the frontal cavity, and in particular where the largest disturbance of the envelope is found, it can be seen that the expansion grows quasi-spherically with a velocity greater than 10 km s⁻¹, and the gas with higher temperatures, in this velocity range, is pushing the cavity. For the high-velocity gas (see Fig. 10), the jet material moving with a high Mach number propagates in the jet-flow line, but for the collision with the secondary jet curves, this happens within the propagation of the main flow. This deviation causes the main flow to interact with the inner wall of the envelope, deflecting the jet flow due to its collision with the wall. Changing the direction of the main flow produces an asymmetry in the envelope that moves in the direction of the injection region opposite the injection of the secondary flow and is a relic of the initial interaction of the two flows. At the rear of this flow, one can observe the creation of a small bubble caused by the interaction between the secondary flow and the inner wall of the envelope. This protuberance reduces the density of the wall, creating a hole in the envelope. It is more visible, particularly in the channel of velocities between 15 and 40 km s⁻¹, where you can see the supersonic expansion of the frontal cavity. This expansion appears spherical, producing a significant increase in the frontal bubble, which produces the great asymmetry observed. In Fig. 1, a decrease in CO emission is also observed near source B. One can observe similar behavior in Fig. 10c, where, near the source, no material is observed at higher speeds (>40 km s⁻¹), which means that this effect is due to the collision.

Fig. 8 Resulting force field, where arrows indicate the 3D direction of the force and a color scale encodes its magnitude in the y–x plane (top panel) and the y–z plane (lower panel).

Fig. 9 Position of the knots with respect to the xy- (green and square) and zy- (pink and circular) planes. The dashed line represents the angle between the first and the last knot. For the xy-plane, an angle of 0.84 degrees is observed, and for the zy-plane, we see an angle of 3.51 degrees. Furthermore, the angles for the knots are presented. These angles were calculated using two consecutive knots, representing the inclination of the line drawn between the two knots.

4.5 Observational predictions

To predict the observation of this kind of jet collision, we made emission maps and integrated them over the z-axis, which we considered as the line of sight. That is, we assuming that the main flow of the primary jet propagates along the y-axis and the main flow of the secondary jet propagates along the x-axis. Therefore, by integrating over the z-axis, we can analyze the formation of the hole in the envelope of our object. To obtain the molecular emission, we consider that all the gas has a fraction of CO concerning the total density of χ_CO = 1.67×10⁻⁴. To find the emission, we used Eq. (3) of Rodríguez-González et al. (2023): $$j_{\text{CO}(2\to 1)}=\frac{1}{4\pi }\frac{g_1}{Z(T)}{n}_{\text{CO}}\cdot e^{\frac{-E_{2\to 1}}{k{T}_{2\to 1}}}{A}_{2\to 1}{E}_{2\to 1},$$(9)

where, g₁ = 3 is a degeneracy factor, $Z(T)={\displaystyle \sum }_i^N{e}^{-T_{\text{levCO}}\text{}/\text{}T(i,j)}$ is the partition function, with T_levCO as the temperature of the corresponding energy levels, A_2→1 = 7.16 × 10⁻⁷ s⁻¹ is the spontaneous emission coefficient, E_2→1 = hν_2→1 is the energy of the transition with ν_2→1 = 230.538 GHz, and h is the Planck constant (Rybicki & Lightman 1985). We assumed that n_CO = χ_CO ∗ n_g, where n_g is the numerical gas density.

Figure 11 shows a normalized CO emission map for velocities of 1 < v < 20 km s⁻¹ in the direction of the observer, considering the jet is travelling along the plane of the sky. In the upper panel (between 1 and 4 km/s), we show that the emission delineates the outer cavity walls, including the bubble generated by the secondary jet in the bottom left corner. A second bubble can be noticed by the lack of emission at a distance of 2 × 10⁴ au, corresponding to the broadest part of the outflow. Also, the leading bow shock shows no emission, which can be explained by the fact that it is expanding at larger velocities. The second panel (between 4 and 7 km s⁻¹) shows the inner cavity walls, while the emission is spread more evenly across the outflow. Again, the leading shock, as the secondary bubble, shows no emission in this velocity range. The third panel (between 7 and 10 km s⁻¹) shows more emission inside both bubbles, and under x = 10⁴ au the main jet outer layer can trace the deviation due to the collision between the momenta of the flows. In the fourth and fifth panels (between 10 and 13 km s⁻¹ and between 13 and 20 km s⁻¹), the emission is stronger in the leading part, which means there is a significant expansion, probably due to the bullets catching up with each other in a previous simulation time. That is comparable to Cep E, where the outflow southern lobe appears filled with material connecting the jet and the cavity. This could explain the hole observed in the southern cavity of Cep E in the low-velocity CO (2 → 1) emission. While other mechanisms intended to change the CO emission remain to be explored, this scenario is a plausible explanation for the observed maps.

Fig. 10 Velocity field for low (from 5 to 15 km s−1; panel a), intermediate (from 15 to 40 km s−1; panel b), and high velocities (>40 km s−1; panel c).

5 Conclusions

We analyzed the collision of a system of two outflows produced by the collimated jets ejected from two protostars, taking some of the physical constraints imposed by the protostellar binary system Cep E, and we carried out a 3D hydrodynamic numerical model analysis. The presence of free parameters prevents the creation of a fine-tuned model, but we reproduced the physical signatures associated with this complex collision in Cepheus E, where the momentum in all three spatial directions and the ram pressure can significantly influence the structure morphology and dynamical evolution of the primary envelope; this is driven by the protostar Cep-A. To quantify the interaction between Cep E outflows, we produced a fiducial model consisting of a single source ejecting a collimated hypersonic jet that propagates in a density-stratified environment on the y-axis, which we named a the primary jet. Then, we produced a second model that includes the fiducial parameters. After 660 yr, a secondary source ejects a collimated jet that moves parallel to the x-axis from x = 1500 au y = 1500 and z = 250 au to avoid direct collision between both simulated jets. Nevertheless, since they expand perpendicularly to their main axis, they affect each other’s dynamics and the morphology of their cavities. We summarize our conclusions in the following points:

• We estimate the probability of an encounter between different outflow components: ∼80% for a cavity–cavity encounter and ∼20% for a jet–jet encounter;

• The fiducial model shows that the primary high-velocity jet produces a low-velocity cavity comparable to axisymmetric simulations such as Rivera-Ortiz et al. (2023);

• The propagation of a secondary high-velocity jet inside the cavity generated by the primary jet induces perturbations inside the primary cavity, showing a morphology of density structure associated with the first interaction. This collision creates an intermediate velocity bubble that travels behind the primary jet’s leading bow shock, where it tends to expand and feed back to the primary outflow cavity;

• The primary jet’s morphology is also affected. Compared with the fiducial model and axisymmetric simulations, the primary jet becomes wider faster by interacting with the secondary jet. Then, the internal knots produced by the variability in the primary jet reach a larger size that is approximately twice that of the unperturbed knots; therefore, its density is reduced considerably. This can also explain why M2 is shorter than M1 at 1875 yr, since the dynamical friction depends on the surface of the propagating knots;

• The secondary jet produces a shocked region inside the primary cavity wall, deeply affecting the initial morphology. This can be observed as a stationary bubble aligned to the secondary jet’s propagating axis, where the density and temperature increase significantly. This could explain the hole observed in the southern cavity of Cep E in the low-velocity CO (2 → 1) emission;

• The collision region of the jets, that is, the zone where both jets are at a minimum distance from each other, is not significantly affected compared to the dynamics and morphology of the main cavity. This likely depends on the impact parameter between both jets, which is a parameter that remains to be explored in more detail in future simulations;

• For our simulation, the primary jet exhibits a pseudo-precession produced by the outflows interaction. In the projected xy-plane where both jets are propagating, the primary jet obtains a maximum deviation of ∼3.5^◦ from its initial direction. This could explain why the northern jet in Cep E exhibits precession around its main propagating axis. The southern jet’s main axis is not aligned with the northern lobe, and its precession amplitude is also different from that of its northern counterpart;

• Given its density, radius, and velocity, the main jet injects mass at a rate of 4.3 × 10⁻⁸ M_⊙ yr⁻¹. We calculate the mass of each of the moving knots in models M1 and M2 at t = 1875 yr, which are, in both cases, around 5 × 10⁻⁵ M_⊙. That corresponds to a mass-injection rate of ∼1.3 × 10⁻⁷, which is around three times the previous jet injection rate. This factor agrees with the entrainment effect discussed by Rivera-Ortiz et al. (2023).

Since the interaction between the two supersonic outflows occurs at thousands of astronomical units from the protostars, we do not consider the following elements that are likely to be important for the evolution of the protostellar flows: the disks surrounding the protostellar objects, radiative transfer, magnetic fields, or gravitational forces of the binary system. In the particular case of Cepheus E, the collimated flows (or 690 jets) from the north and south exhibit similarities in phase and period of the jet injection, but not in spatial frequency. Our analysis thus focuses on the purely hydrodynamic effects of the interaction between two outflows. Therefore, there is a large range of possibilities for future studies of the molecular outflows ejected by protostellar Class 0 objects.

Fig. 11 Emission of M2 model considering channels of velocities. For these maps, the velocities in the direction of the observer were considered. The upper panel shows the emission for the velocities between 1 and 4 km s−1. The emission delineates the outer cavity walls, including the bubble generated by the secondary jet in the bottom left corner. On the other hand, a second bubble can be noticed showing the lack of emission at a distance of 2 × 104 au. The second panel, 4 and 7 km s−1, shows the inner cavity walls and the leading shock, like the secondary bubble, shows no emission in this velocity range. The third panel, showing the emission between 7 and 10 km s−1, allows us to observe more emission inside both bubbles, and near the base it shows that the main jet’s outer layer can trace the deviation due to the collision between the flows. The fourth and fifth panels present the higher velocities, between 10 and 13 km s−1 and 13–20 km s−1, respectively. The emission is stronger in the leading part, meaning a significant expansion exists.

Acknowledgements

The authors are grateful to B. Lefloch for his participation and valuable contributions to the development of this article, and appreciate the comments of the anonymous referee, who helped improve the content of this work. We acknowledge the support of the UNAM-PAPIIT grants IN110722, IN102724, and IG101125, as well as the Miztli-UNAM supercomputer project LANCAD-UNAM-DGTIC-123 (2022-1) and LANCAD-UNAM-DGTIC-128 (2023-1). E.A. and A.D. are grateful for the grants, CONAHCYT numbers 4050143 and 825308, respectively. RO is grateful to the postdoctoral program SECIHTI EPM 519149, and the project SECIHTI ACADEMIAS-2025-A-54 managed by Sociedad Astronomica de Zacatecas.

References

All Tables

All Figures

Fig. 1 Left: CO J= 2–1 emission detected in the Cep E protostellar core with the IRAM interferometer at 1″ (820 au) resolution. The black contours trace the emission from the low-velocity gas (|V −Vlsr| < 8 km s−1), and the blue/red contours trace the high-velocity jet (V − Vlsr > 50 km s−1, powered by component A of the Cep E protostellar core. The data are from Lefloch et al. (2015) and de A. Schutzer et al. (2022). Right: SiO J= 5–4 emissions detected in the Cep E protostellar core with the IRAM interferometer at 1.4″ (1150 au) resolution in the jets powered by protostars A (red and white contours) and B (blue and red contours). The velocity intervals of integrated flux are reported (in km s−1) next to the lobes of the jets from A and B. The first contour and contour interval are 10% and 20% of the maximum peak intensity. The data are from Ospina-Zamudio et al. (2018).

In the text

	Fig. 2 Interaction between southern low-velocity outflow lobe of A, at [−8; −6] km s−1 (black contours) and the eastern lobe of the jet from B (red contours), as observed with the IRAM interferometer in the molecular lines of CO J=2–1 and CS J=2–1 (left), and SiO J=2–1 and H2CO (right).
In the text

	Fig. 3 Probability of two outflows overlapping for two wide cavities (continuous black line) and for two collimated jets (blue line). The dashed line represents the observed separation of Cep E-A and Cep E-B, which could be larger because of projection effects.
In the text

	Fig. 4 Two protostellar sources inject high-velocity material. The gas density around each source was stratified with a power law. In our simulations, the jet injections were not meant to collide if they preserved their width, but, since they expand thermically and by entrainment, the evolved jets – that is, high-velocity material – could collide.
In the text

Fig. 5 Upper panel: snapshot of gas density of a single jet (model M1) at t = 1875 yr. One can see the work-surface structure due to the consideration of the variable jet. The envelope shows a higher density at the “base” of the jet and a lower one toward the front. Lower panels: snapshots of gas density of the jet-interaction model (model M2) at three evolutionary times: t=1375, 1625, and 1875 yr (top, middle, and bottom panels). Considering the estimated age of the Cep E protostellar outflow of ∼1500 yr (de A. Schutzer et al. 2022), different stages of evolution around this age are shown. All the snapshots are present on the xy-plane.

In the text

	Fig. 6 Snapshots of gas density of two-jet collisional model (model M2) in the Y–Z plane. The color represents the density, and the white boxes show the position of the knots. Using the boxes, we can find the mass of the knots.
In the text

	Fig. 7 Mass of the knots of models M1 (blue circles) M2 (red circles) and the accumulated mass are shown by crosses. The figures show that the masses in the knots of the model with two-jet interaction were systematically greater for all knots than those of the M1 model, except for the furthest knot, which was difficult to determine.
In the text

	Fig. 8 Resulting force field, where arrows indicate the 3D direction of the force and a color scale encodes its magnitude in the y–x plane (top panel) and the y–z plane (lower panel).
In the text

Fig. 9 Position of the knots with respect to the xy- (green and square) and zy- (pink and circular) planes. The dashed line represents the angle between the first and the last knot. For the xy-plane, an angle of 0.84 degrees is observed, and for the zy-plane, we see an angle of 3.51 degrees. Furthermore, the angles for the knots are presented. These angles were calculated using two consecutive knots, representing the inclination of the line drawn between the two knots.

In the text

	Fig. 10 Velocity field for low (from 5 to 15 km s−1; panel a), intermediate (from 15 to 40 km s−1; panel b), and high velocities (>40 km s−1; panel c).
In the text

Fig. 11 Emission of M2 model considering channels of velocities. For these maps, the velocities in the direction of the observer were considered. The upper panel shows the emission for the velocities between 1 and 4 km s−1. The emission delineates the outer cavity walls, including the bubble generated by the secondary jet in the bottom left corner. On the other hand, a second bubble can be noticed showing the lack of emission at a distance of 2 × 104 au. The second panel, 4 and 7 km s−1, shows the inner cavity walls and the leading shock, like the secondary bubble, shows no emission in this velocity range. The third panel, showing the emission between 7 and 10 km s−1, allows us to observe more emission inside both bubbles, and near the base it shows that the main jet’s outer layer can trace the deviation due to the collision between the flows. The fourth and fifth panels present the higher velocities, between 10 and 13 km s−1 and 13–20 km s−1, respectively. The emission is stronger in the leading part, meaning a significant expansion exists.

In the text