Shedding Light on the Ejection History of Molecular Outflows: Multiple Velocity Modes and Precession

Variable accretion has been well studied in the evolved stages of low-mass star formation. However, the accretion history in the initial phases of star formation is still a seldom studied topic. The outflows and jets emerging from protostellar objects could shed some light on their accretion history. We consider the recently studied case of W43-MM1, a protocluster containing 46 outflows driven by 27 protostellar cores. The outflow kinematics of the individual cores and associated knots in W43-MM1 indicate episodic protostellar ejection. We take the observed parameters of an individual core system (core #8) and perform 3D hydrodynamic simulations of such a system, including episodic changes in the velocity of the outflow. The simulations consist of a collimated jet emerging from a core, taking into account one- and two-velocity modes in the variation of the ejection velocity of the jet. In addition, we investigated the effect of including the precession of the jet in the one- and two-velocity-mode models. From the simulations, we constructed position–velocity diagrams and compared them with the observations. We find that including a second mode in the ejection velocity, as well as the precession, are required to explain the positions of the outflow knots and other position–velocity features observed in core #8 in W43-MM1.


INTRODUCTION
The time variability of gas accretion during star formation is an important and debated topic.While lowmass, pre-main sequence stars have been statistically studied on short time scales, the variability of accretion during the embedded protostellar phase is not well documented (see, e.g., reviews by Audard et al. 2014;Fischer et al. 2023).Yet, it is during this phase that most of the final stellar mass accumulates from the gas reservoir.Observations of few protostars have brought evidence of quasi-periodic ligh curves, interpreted as accretion variability (see e.g.V2492 Cyg Hillenbrand et al. 2013).
Jets and outflows develop as protostellar accretion occurs and they aid in the removal of angular momentum from the collapsing protostellar core.The morphology and kinematics of molecular outflows could thus provide fossil records of the accretion history of a protostar (Rohde et al. 2022).Molecular outflows from young proto-stars extend up to a few ∼0.1 pc, with collimation angles of a few degrees only, and velocities up to 100 km s −1 (see, e.g., the review by Bally 2016).They are often composed of a collimated jet at high velocity, plus a shell-like wind tracing gas at lower speeds.Collimated jets frequently appear as chains of knots interpreted as internal shocks produced by episodic variations in the protostar mass-loss rate and/or the ejection velocity.
In the massive protocluster W43-MM1, Nony et al. (2020) discovered a cluster of 46 outflows driven by 27 protostellar cores, with masses in the range from 1 to 100 M ⊙ .Their detailed investigation of outflow kinematics using PV diagrams reveled clear events of episodic ejection.A typical ejection timescale (period) of about 500 yr has been estimated from the positions of 86 knots along 38 outflow lobes.
Based on the observational evidence, Raga et al. (1990) proposed to include an intrinsic variability in the source of jets to model stellar outflows, and was soon applied to specific objects, as the case of HH34 (Raga & Kofman 1992).An analytical description of the formation of the shocks using the viscous Burgers equation was given in Kofman & Raga (1992).The solution of the resulting "internal working surfaces" for a sinusoidal variation can be found in Raga & Cantó (1998).This model has been widely used to explain Herbig-Haro (HH) objets, for instance, Castellanos-Ramírez et al. (2018) have modeled the optical HH1 jet as the result of a two-mode periodic velocity jet.The model reproduces remarkably well the chain of knots seen close to the source, as well as the time evolution which is resolved by multi-epoch HST observations.The motivation to include sinusoidal time dependencies in the ejection velocities is therefore phenomenological.The physical origin of the velocity variability is still unclear, but different explanations have been proposed.For example, the variability could be due to disk thermal instabilities caused by a runaway process.If the temperature in the accretion disk increases to reach the hydrogen ionization temperature, the opacity in the disk increases as κ ∝ T 10 (Audard et al. 2014) trapping the energy inside the disk and causing a rise in the temperature (thermal runaway).The viscosity is proportional to the temperature in the disk, and to the accretion rate, therefore an enhancement in temperature causes an increase in the accretion rate.Another way of inducing a thermal instability could be due to the presence of a companion (Lodato & Clarke 2004), or thermal instability could be triggered by the combination of gravitational and magneto rotational instabilities (Martin & Lubow 2011, 2013).Moreover, disk fragmentation could also cause mass accretion and luminosity outbursts, when the material migrate inwards onto the protostar.
The different types of instabilities mentioned before, could cause quasi-periodic variations in the accretion and ejection phenomena on timescales as long as 50 kyr (Audard et al. 2014).It would be expected that these enhancements are stronger than the underlying stochastic variability, which could also be present.Any particular form of periodic variation in a supersonic jet would form internal shocks due to the steepening of the waveform, and in the hyper-sonic limit can be described by the Burgers equation (Kofman & Raga 1992;Raga & Cantó 1998).Thus, the choice of periodic sinusoidal variations is a reasonable empirical description in jets and outflows with an observed periodicity in their knots.
Only few models have been compared with observations of molecular outflows using PV diagrams, and they strongly depend on the structure adopted for the protostellar envelope (e.g., Rohde et al. 2019).Two recent works use PV diagrams to compare their simulations with observations of molecular outflows.Rabenanahary et al. (2022) compared simulations of an outflow driven by a pulsed jet and observations in CARMA-7 by Plunkett et al. (2015).Rivera-Ortiz et al. (2023) carried out simulations of a jet aimed at reproducing the Cepheus-E outflow observed by de A. Schutzer et al. (2022).
Motivated by the observations of Nony et al. (2020), which show chains of knots arising from dense cores forming intermediate-to high-mass stars, we performed 3D hydrodynamic simulations of the outflows emanating from dense cores.In particular, we explore if a variability in velocity and precession can explain the kinematic structures emerging from the outflow of core MM1#8 in Nony et al. (2020).
The paper is organized as follows: In Section 2 we briefly present the observations of the W43-MM1 protocluser.In Section 3 we describe our 3D hydrodynamical numerical simulations.In Section 4 we report the results of the simulations, and the comparison with the observations.Finally, in Section 5 we summarize our main conclusions.

W43-MM1 PROTOCLUSTER DATA
This work makes use of ALMA Band 6 observations of the W43-MM1 protocluser.The continuum and lines images have an angular resolution of about 0.5 ′′ (∼ 2200 au at a distance of 5.5 kpc).A first catalogue of cores has been published by Motte et al. (2018), and the association between cores and molecular outflow identified in CO(2-1) has been presented in Nony et al. (2020).
In the following, we focus on the outflow driven by core #8.This core has a gas mass of 18±3 M ⊙ and a systemic velocity of V LSR = 97 km s −1 .(see Nony et al. 2020 andNony et al. 2023).Core #8 has been chosen for the unique morphology of its bipolar outflow.Indeed, the blue-shifted lobe is the longest of the region, with a total length of 0.44 pc, and therefore the more substructured (eight knots).Its red-shifted counterpart extends over 0.13 pc and has four knots.The relative positions and velocities of the knots reported in this study are taken from Table 3 of Nony et al. (2020), which we include in Table 1.Their velocity range is from 36 km s −1 to 79 km s −1 from the core V LSR .

The jet model
A common feature in low-mass protostellar outflows (i.e., HH objects) is a chain of aligned "knots" close to the source, which can be produced by a variability in the ejection.The most accepted model to produce the observed characteristics of this chain of knots is to consider a periodic velocity ejection variation (see for instance Biro & Raga 1994).In this scenario, high(er) velocity material catches up with previously ejected material, coalescing into a double shock structure also known as internal working surface.These internal working surfaces are precisely the knots, which move as a unit away from the source.The distance between successive knots is given by the variability of the ejection velocity.
In some cases, such as HH 1 and HH 34 (Raga et al. 2011(Raga et al. , 2012)), besides the chain of knots close to the source, a second series of one or more knots are observed at larger distances.In order to reproduce the morphology of HH 34, Raga & Noriega-Crespo (1998); Raga et al. (2011Raga et al. ( , 2012) ) included two or three modes in the ejection velocity variability.More recently, Castellanos-Ramírez et al. ( 2018) implemented a two-velocity variable mode to reproduce the morphology of HH 1.The variabilities in the jet velocity produce a series of twoshock internal working surfaces with maximum shock speeds comparable to the half-amplitude of the ejection velocity variability.
An alternative to a velocity variation that could produce a chain of knots is to consider a constant velocity and a time-dependent density.However, this scenario does not produce sufficiently strong shocks to account for the excitation state observed in the knots.In that case strong shocks are only produced when the material interacts with the environment, but there are no strong shocks in between knots.
Castellanos-Ramírez et al. ( 2018) point out that such variability could become important if the jet has a very substantial precession, in such a way that the denser regions of the jet can interact with the surrounding environment.
Interestingly, Nony et al. (2020) also proposed that the observed changes in the direction of the outflow of core #8 could be the consequence of periodic oscillations due to the precession of the molecular jet.The latter work prompted us to study core #8 in detail using 3D-hydrodynamical simulations of a jet with oneand two-mode velocity variability, which presents (or not) precession.In the following sections we describe the code used (Section 3.2), and summarize the set-up of the simulations (Section 3.3).

The code
To model the interaction of a jet from a core with its environment we use the guacho code (Esquivel et al. 2009;Esquivel & Raga 2013;Villarreal D'Angelo et al. 2018).The version of the code we used solves the ideal hydrodynamic equations (Equations 1-3) along with a rate equation for hydrogen (Equation 4): In the latter set of equations (Equations 1-4), ρ, u, and P are the (total) mass density, flow velocity, and m7 m10  2 for details on the simulations.Blue-green tinges show lower number densities (10 4 cm −3 ).Yellow-red tinges show higher number densities (10 5 -10 6 cm −3 ).The positions of the knots (in red) are clearly distinguished at the highest densities.
thermal pressure.E is the energy density, G and L are the energy gain and losses (respectively) due to radiative processes.An ideal gas equation of state is used to relate the total energy density with the kinetic energy density and thermal pressure as E = ρ|u 2 |/2+P/(γ−1), where γ = 5/3 is the ratio between the specific heat capacities.For the hydrogen rate equation (Equation 4), n HI is the neutral hydrogen (number) density, n e the electron density, and n HII the ionized hydrogen density.
We assume that all the material consists of hydrogen, that all the free electrons arise from ionization of hydrogen (n e = n HII ), and neglect the mass of such electrons to compute the mass density, i.e. ρ ≈ m H (n HI + n HII ), where m H is the hydrogen mass.In the source terms in Equation ( 4), α(T ) and c(T ) are the (case B) radiative recombination coefficient and the collisional ionization coefficient, respectively.The equations are solved in a regular Cartessian grid with a second order accurate Godunov type method that uses the HLLC approximate Riemann Solver (Toro 1999), and a linear reconstruction with the minmod slope limiter to ensure stability.We do not include radiative gains, but we compute the radiative cooling using the hydrogen ionization fraction following the prescription given by Biro et al. (1995), and add it to the energy equation with the implementation described in Villarreal D'Angelo et al. (2021).
In our simulations we do not consider magnetic fields.The magnetic field plays a crucial role in the launch-ing mechanism of the jets (Shu et al. 2000;Konigl & Pudritz 2000), and it is also regarded as important to maintain their collimation at large distances from the source (Blandford & Payne 1982).However, our models deal with scales larger than the injection of the jet, and since we impose a collimated flow, it is reasonable to disregard the magnetic effect as a first approximation.This can be further justified considering that the ram pressure in the jets produced in our models is on the order of 10 −4 dyn cm −3 (e.g. a number density of ∼ 10 6 cm −3 and a velocity of ∼ 100 km s −1 ), which is several orders of magnitude larger compared to the magnetic pressure of ∼ 4 × 10 −8 dyn cm −3 assuming a typical value of |B| = 1 mG.

Numerical Setup
We ran a set of 3D-hydrodynamical simulations of a collimated flow arising from a core.The extent of the computational domain is 0.17 × 0.17 × 0.5 pc in the (x, y, z) directions, with a uniform resolution of 8.33 × 10 −4 pc (i.e. a grid with 200 × 200 × 600 cells).
The jet is assumed to travel into a stationary environment with the characteristics of a clump with a (radial) density profile of the form:  2020) for core #8 (see Table 1).
where r 0 = 0.05 pc, and the value of ρ 0 is set to result in a mass of 2500 M ⊙ within a radius of 0.5 pc (yielding ρ 0 ≈ 2.1 × 10 −18 g cm −3 ).The temperature of the environment is set to T c = 1000 K.Even if this value is high, by running control simulations we found out that the medium temperature does not change our results; the important aspect being the pressure (nT ) contrast between the medium and the jet.
Therefore, we considered a jet also with a temperature of T j = 1000 K, and a (mean) mass loss rate of Ṁj = 2 × 10 −4 M ⊙ yr −1 .The density of the jet is constant, and its value is obtained for each model from the mass loss rate as: where v 1 is the mean velocity of the jet, and R j its cylindrical radius.The jet is imposed at every timestep in a cylindrical region of length L j = 1.6 × 10 16 cm, and radius R j = 1.6 × 10 16 cm.The injection velocity is along the cylinder symmetry axis.The interaction between the jet and the dense molecular envelope from the clump results in the formation of a layer of entrained gas which follows the motion of the central jet.The resulting ejected (jet) and entrained material is, at first order, representative of the observed molecular outflow.
For the velocity of the jet, we consider a sinusoidal variability in the ejection velocity.We first considered a one-mode velocity variability mode described of the form: Secondly, we considered a two-velocity variability mode described by: In the models, we explored three different values of the average velocity of the jet; ( vj = v 1 = 100, 120, and 150 km s −1 ).The velocity has a fast high-amplitude mode δv 1 ; we explore two different values δv 1 = 50 and 100 km s −1 , with a period of τ 1 = 0.15 kyr.For the twomode variability mode in the velocity, a second term is added (see Equation 8) corresponding to a slow lower-Lora et al. amplitude mode δv 2 ; we explored two different values δv 2 = 30 and 50 km s −1 , with a period of τ 2 = 0.6 kyr.
Besides the one-and two-velocity variability modes, we also investigated the effect of the precession of the jet.Therefore, we included models in which the jet precesses with a half-opening angle of α p = 5 • and a period of τ p = 2.4 kyr.In total we performed 12 simulations, the initial conditions and relevant parameters are summarized in Table 2.

RESULTS
We let the models described in Section 3.3 to run for 20 kyr.As an example, in Figure 1 we present two density render images of two different models (M7 and M10, see Table 2).The transparency in both panels of Figure 1 vary linearly with density, meaning that denser regions appear more opaque.Two different tinges represent two different density values; blue-green tinges show lower number densities (10 4 cm −3 ), and yellowred tinges show lower number densities (10 5 -10 6 cm −3 ).The rendered images were created using the Multi-code Analysis Toolkit for Astrophysical Simulation Data yt (Turk et al. 2011).

The one velocity mode case
In Figure 2 we show a density cut in the X-Z midplane of the simulation after an integration time of 15 kyr.The associated column density maps are shown in Figure A.1.The separation between successive knots in a single velocity mode case is given by ∆x ≈ v 1 τ 1 which for our models without precession; M1, M2 and M3 (left panels of Figure 2), are ∆x1 M1 = 0.015 pc, ∆x1 M2 = 0.018 pc, and ∆x1 M3 = 0.022 pc, respectively.
If we compare models M1 and M2, whose only difference is the mean jet velocity (by 20 km s −1 ), we can see that the chain of knots is evident at larger distances for the model with higher mean velocity (M2, middle row in the left column of Figure 2).The reason for this is a combination of the higher velocity of the flow itself, with the fact that the leading shock in faster models has the opportunity to clear more material ahead of the source.Thus, the subsequent internal working surfaces can travel more easily into the low-density canal carved by the leading shock.We can verify this in model M3, which has a higher mean velocity of 150 km s −1 .
For visual reference, in Figure 2 we over-plotted the positions of the knots of core #8 measured by Nony et al. (2020) and reported in Table 1.The one-velocity mode reproduces fairly well the knots closer to the source, but this model fails to reproduce the (blue) knots at farther distances (> 0.187 pc, B4 to B7).The change of Table 2. Initial conditions of the simulations; each of the columns correspond to (1) model number, (2) average outflow velocity, (3) fast-lowamplitude mode, (4) slow-high-amplitude mode, and (5) simulation including precession with half opening angle αp = 5 • and a precession period of τp = 2.4 kyr.All models include a velocity variability with period τ1 = 150 yr, while models 7 through 12 include a second variability mode with period τ2 = 600 yr.direction in the ejected material in the models with precession (right column in Figure 2) results in some of the material colliding with the "walls" of the low-density cavity carved be the leading shock (and in this case also by subsequent internal working surfaces).This oblique collision results in a density increase and a slight deceleration of the knots that collide with the wall.In the case of the one-mode variability models we see a density enhancement at a distance of ∼ 0.31 pc which can be associated with one of the knots farther away from the source (B5; top-right panel of Figure 2, and B7; middle-right panel of Figure 2).However, it fails to reproduce the other blue knots observed at greater distances (> 0.31 pc).

The two velocity mode case
We added to the simulations presented in Section 4.1 a second velocity mode with a 600 yr period, and a slightly lower amplitude than the fist mode (see Table 2).This second mode gives a separation among knots for model without precession of ∆x2 M7 = 0.06 pc, ∆x2 M8 = 0.07 pc, and ∆x2 M9 = 0.09 pc.The second velocity mode produces internal working surfaces that would start appearing at : from the source (Castellanos-Ramírez et al. 2018).This is strictly valid when δv 2 ≪ v 1 (Raga & Noriega-Crespo 1998)).In our models M7, M8, and M9 the starting position of the working surfaces of the second mode, are ∆x ws,7 = 0.07, ∆x ws,8 = 0.026, and ∆x ws,9 = 0.04, respectively.
For the two-mode velocity case we present Figure 3, which shows a density cut in the XZ-midplane at an integration time of 15 kyr.The associated column density maps are shown in Figure A.2.The simulations without and with precession are shown in the left and right panels, respectively (see all initial parameters used in Table 2).The impact of the second mode in the velocity is clear.In this case we obtain a series of knots periodically appearing at farther distances (∆x2), which now could account for the position of the four farthest blue knots observed in core #8.
In the simulations without precession (left panels of Figure 3), the second velocity mode results in a similar morphology to the case with one velocity mode, with a small variation in the distances among knots.Similarly to the one-velocity modes, the effect of increasing the average velocity of the jet (20 km s −1 between model M7 and M8, and 30 km s −1 between model M8 and M9; see left panels of Figure 3) results in the appearance of the farthest knots at greater distances.The same behaviour is obtained for the models including precession.However, when including the precession in the simulations (right panels of Figure 3) the farthest knots can be better distinguished.The latter is because the material collides with the surface of the cavity carved by previous high velocity material, thus increasing the local density resulting in better defined knots.

PV diagrams
The upper panel of Figure 5 (a) of Nony et al. (2020) shows the PV diagrams along the molecular outflow of core #8.In order to compare our results to the observations, we also created from our simulations the synthetic PV diagrams of line-of-sight (LoS) velocities as a function of projected distance (see Figures 4 and 5).The PV diagrams were computed placing a virtual slit along the jet axis with a velocity resolution of 2 km s −1 .We averaged the PV from 20 cells in the x-direction (perpendicular to the jet axis), resulting in an spatial resolution of ∼ 0.015 pc, corresponding approximately to the smoothed spatial resolution in the observed PVs.
An additional free parameter when computing PV diagrams is the relative orientation between the jet and Lora et al. the LoS.The PVs shown here correspond to an angle of 60 • between the jet axis and the LoS (30 • with respect to the plane of the sky).This angle is close to the most probable inclination for an outflow, and it is also consistent with the estimations for the observed outflow (see Figure 12a of Nony et al. 2020).With this angle, LoS velocities are reduced by a factor of two compared to the intrinsic jet velocities, while projected distances are less impacted.We over-plotted in Figures 4 and 5 the CO contours used by Nony et al. (2020).The PV diagrams in the one-velocity mode case without precession (Models 1, 2, and 3; left panels of Figure 4) show the continuous chain of knots driven by δv 1 .For models 1, 2, and 3, at distances > 0.18 pc, the PV diagrams show a clear plateau at constant velocity of ∆V ≃ −60km s −1 .Such plateau is observed even when including the precession (Models 4 and 5; see upper and middle right panels of Figure 4).This is not in agreement with the variation of velocity observed in the data, with a velocity of knots B3 and B4 about 20 km s −1 lower than B5 and B6 and even 40 km s −1 lower than B1 and B2.
The PV diagrams in the two-velocity mode case without precession (Models 7, 8, and 9; left panels of Figure 5) show clearly the knots at distances greater than 0.25 pc, whereas the plateau observed in the analogous onevelocity mode (Models 1, 2 and 3) is not observed.The velocity variations observed in the models with precession (10 and 11), increasing rapidly up to ≃ 0.12 pc and decreasing at larger distances, qualitatively reproduce the observations.In order to quantitatively determine which of the twovelocity mode models best reproduces the observations, we measured the positions and velocities of the knots in the PV diagrams of our simulations, following the method described in Nony et al. (2020).The average velocity of the knots is ∆V knot ≃ 63, 71 and 97 km s −1 for models M7-M10, M8-M11 and M9-M12, respectively.These values are about 60% of the average injection velocity.Models M7 and M10 have an average knotvelocity similar to that reported by Nony et al. (2020) for core #8 (∆V knot = 50 km s −1 ), and therefore our best model candidates.
Interestingly, based on the model predictions for the knots velocity and position, we found between knots B1 and B2 of the observed PV diagram (Figure 5 (a) of Nony et al. 2020) an over-density which was not previously reported.Although this structure does not correspond to all our criteria for defining a knot -an overdensity at a local velocity maximum -, it could arguably correspond to an independent ejecta.We labeled this structure B1a and report its coordinates in Table 1.
In Figure 6 we compare the relative positions and velocities of the observed knots of core #8 outflow (B0 to B7) with that of the knots found in models M7 and M10.We associated each of the observed knots to its closest counterpart in the models, and measured the relative velocity and position offsets.We first note that the inclusion of precession in model M10 compared to model M7 introduces only a small difference in the positions and velocities of the knots (up to 8 km s −1 and 0.04 pc, respectively).Model M10, which includes precession, reproduces slightly better the observations than model M7.The average position and velocity differences are indeed smaller: 0.018 pc vs 0.025 pc, and 11.5 vs 13.1 km s −1 , for models M10 and M7 respectively.Model M10 performs specially better in reproducing the knots at large distances (B5 and B6).However, neither model reproduces the velocity of knots B3 and B4 and, conversely, both models predict a knot at high velocity (86 km s −1 at 0.06 pc) which does not have an observational counterpart.

CONCLUSIONS AND DISCUSSION
Jets and outflows in star forming regions often present substructure in the form of knots, which indicates that the material ejection could be episodic (e.g., Qiu et al. 2019;Li et al. 2020).Significant ejection variability has been observed in the high-mass protocluster W43-MM1 (Nony et al. 2020), with the outflow of core #8 presenting the most numerous knots and therefore constituting the target of our study.
We ran 3D hydro-dynamical simulations of a protostellar jet and studied the effect of including a periodic velocity ejection variation, with one and two velocity modes.The values of the velocity-modes were selected to match the observed periodicity of core #8 reported in Lora et al.Nony et al. (2020).Additionally, we explored the effect of the precession of the jet.
We ended up with 12 different models (see Table 2), for which we studied the knots formed in the simulations making use of density (and column density) maps, as well and PV diagrams (see Figures 2-5).
We summarize our results as follows: • Models with one-velocity ejection modes (M1 to M6) do not match the observed velocity variations of core #8.Instead, models with two-velocity modes (M7 to M12) have a better agreement with the observed velocity and positions of the knots observed in core #8.
• The observed outflow shows evidence of precession.Precession leads to periodic oscillations of the outflow direction, and it also has an impact on the velocity of the knots (see, e.g., the comparison between M7 and M10).
• Model M10, with two-velocity modes, precession, and an average jet injection velocity v 1 = 100 km s −1 reproduces fairly well the observations of core #8 even without fine tuning.Very interestingly, our simulations predict the position and velocity of a new knot (B1a, see Figure 6) that was not previously reported in Nony et al. (2020).
Both the variability modes of the ejection and the precession of the jet are likely associated to mechanisms occurring at the protostellar disk scale (∼ 100 au).Perturbations in a binary system have been proposed as a possible phenomena at play (Terquem et al. 1999;Raga et al. 2009).
The physical processes in the molecular outflow of core MM1#8 could be different from those of the intermediate-mass Cepheus E core, which was recently studied with CO observations and numerical simulations (de A. Schutzer et al. 2022;Rivera-Ortiz et al. 2023).Rivera-Ortiz et al. (2023) showed that one-velocity ejection variations successfully reproduce the observed features, and especially the bimodal distribution of the dynamical timescales of the knots.On the contrary, we showed in this work that a two-velocity ejection mode is necessary to reproduce the observed velocity variations among the knots.The differences might come from the different environment where these objects are, or the existence of a companion.
In general, the morphological differences in distinct outflows might come from the different physical mechanisms acting in the disk.For example, thermal instabilities and the existence of a companion, might give rise to quasi-periodic variations in the accretion rate, which might translate in a periodic difference in the positions of the knots in an outflow, depending on the outburst periodicity (Audard et al. 2014).Morphological differences might also come from the evolutionary stage of the outflows, their inclination, and precession, which might translate into shorter or longer outflows.In addition, the environment density where the outflow propagates, could also shape its morphology.
Our modelling of the outflow from core #8 is a step forward to a better understanding of intermediate-mass outflows.This intermediate-mass object (18M ⊙ ) indeed has the potential to form a star of 6-9M ⊙ (when assuming 1/3 to 1/2 efficiency).Moreover, it lies within a cluster of outflow driven by low-to high-mass cores (Nony et al. 2020).Then, we are studying an object in the border between an intermediate and a high-mass protostars.Our work thus demonstrates that the twomode velocity variability, initially proposed for optical jets from low-mass protostars, could also explain the morphology of a molecular outflow from an intermediate mass core.Further simulations of protostellar outflows should be carried to explore the outflows temporal variation and the merging of knots.

ACKNOWLEDGMENTS
We would like to thank the referee for giving very constructive comments and suggestions, which widely help to improve this paper.VL acknowledges the support of CONAHCyT.VL wishes to thank A. Raga for discussions and suggestions throughout this work, and the life sailed through the years.VL and AE acknowledge support from PAPIIT-UNAM grant IN113522.TN and RGM acknowledge support from UNAM-PAPIIT project IN108822 and from CONACyT Ciencia de Frontera project ID 86372.TN also acknowledges support from the postdoctoral fellowship program of the UNAM.

Figure 1 .
Figure 1.Example of two synthetic nebular images (renders) of the jet material of simulations M7 (left panel) and M10 (right panel).See Table2for details on the simulations.Blue-green tinges show lower number densities (10 4 cm −3 ).Yellow-red tinges show higher number densities (10 5 -10 6 cm −3 ).The positions of the knots (in red) are clearly distinguished at the highest densities.

Figure 2 .
Figure2.Density cuts in the XZ-midlplane for the single-velocity mode models.In the left column we present the models without precession, and in the right column the models with precession (with a half opening angle of 5 • and a precession period of 2400 yr).In different rows we show models with different velocity variability.The average velocity of the jet and the amplitude of the fluctuation are displayed at the right of each row.Blue and red crosses mark the position of the blue-and red-shifted knots, respectively, reported byNony et al. (2020) for core #8 (see Table1).

Figure 3 .
Figure 3. Same as Figure 2 but for the models with two-velocity modes.The parameters of these are indicated at the right of each row of plots.

Figure 4 .
Figure 4. Color map of the PV diagrams of the models with one velocity mode (M1-M6, see Table2).We over-plotted the CO iso-contours of the PV diagram of core#8 reported by Nony et al. (2020) at 5, 10, 22 to 67 times the units of σCO = 2.5 mJy beam −1 .With blue crosses, we show the positions of the observed blue knots.As a reference, we draw over the PV diagram of model M5 three yellow lines which represent the "Hubble laws" connecting knots with low velocity structures, reported in Figure 5a of Nony et al. (2020).

Figure 5 .Figure 6 .
Figure 5. Color map of the PV diagrams of the models with two-velocity modes (models M7-M12).Symbols have the same meaning as in Figure 4.

APPENDIXA.
COLUMN DENSITY MAPSIn this appendix we present the column density maps of the simulations presented in the main body of this article.
Figure A.1 shows the column density with the same projection as the PV diagrams for the models with one velocity mode (see Figure4).
Figure A.2  shows the column density with the same projection as the PV diagrams for the models with two velocity modes (see Figure5).In both Figures A.1 and A.2 the left panels show the models without precession, and the right panels show the models with precession (with a half opening angle of 5 • and a precession period of 2400 yr).The rows of Figures A.1 and A.2 show the models with different velocity variability, different average velocity of the jet, and different amplitude of the fluctuation of the velocity.The blue and the red crosses mark the position of the blue-and red-shifted knots for core #8 (see

Figure A. 1 .
Figure A.1.Column density in the same projection as the PV figures for the models with one-velocity mode.The white dashed line denotes the position of the virtual slits used to produce the PV diagrams.

Figure A. 2 .
Figure A.2. Same as Fig. A.1, for the models with two-velocity modes.

Table 1 .
Parameters of the blueand red-shifted knots of core #8.Taken from Table3ofNony et al.