Estimating the propagation of a uniformly
accelerated jet

Jets from young stars and their associated Herbig-Haro objects (HH) are now the most studied and best understood of all astrophysical jets. HH jets are collimated bipolar ejections from low, intermediate and sometimes high mass protostars or young stars (see, e.g., the review of Frank et al. 2014). The interaction between HH jets and the surrounding environment generates a main working surface known as the “head”. This structure (formed by a bow shock/jet shock pair) is a clear sign of the “turning on” of the jet, which tipically occurred at times of $\sim 10^{3}\to 10^{4}$ yr ago in observed HH jets.

This type of jets also have a structure of emitting knots (some of them resembling the jet head, and others forming aligned chains of more compact emission peaks) in the region between the outflow source and the jet head (see, e.g., Reipurth & Heathcote 1990 and the review of Reipurth & Bally 2001). Several possibilities of how to model these knots have been explored (see, e.g., Raga & Kofman, 1992, and Micono et al. 1998), but presently the favoured explanation is that they are the result of a time-variability in the ejection, which leads to the formation of “internal working surfaces” within the jet beam (see, e.g., Raga et al 1990). Some of the “classical HH objects” of the catalogue of Herbig (1974) are associated with either jet heads or knots.

A problem that has received little attention is the effect of a “slow turning on” of the outflow on the jet head. As far as we are aware, the only simulations studying the dynamics of the head of a “slow turning on” jet (as opposed to a jet which is instantaneously “switched on” at full velocity) are the ones of Lim et al. (2002). These authors focus their study on the survival of H₂ molecules in the resulting, accelerating jet head, presenting 1D hydrodynamic simulations including an H₂ formation/destruction chemical network.

Of course, there are a number of numerical studies of the ejection of MHD jets from star+accretion disk systems (e.g., Ouyed et al. 2003; Hiromitsu & Shibata 2005; Zanni et al. 2007; Ahmed & Shibata 2008; Romanova et al. 2018), in which a sudden switch on of the jet is definitely not imposed (the “switching on” of the outflow being a result of the simulated accretion+ejection flow). These models produce a highly superalfvénic leading head that evolves consistently with the computed time-dependence of the ejection.

In the present paper, we revisit the “slow turn on jet” problem. We study the dynamics of the head of a cylindrical jet with a linear ramp of increasing ejection velocity as a function of ejection time. We combine this ejection velocity with two forms for the ejection density variability:

• •

  a constant ejection density,

• •

  a constant mass loss rate (so that the ejection density is proportional to the inverse of the ejection velocity).

We obtain analytical models of this flow by assuming that the material in the jet beam is free streaming (before reaching the jet head), and that the motion of the working surface of the jet head can be described:

• •

  by using a ram-pressure balance condition,

• •

  by assuming that it coincides with the motion of the center of mass of the material that has entered the working surface.

The first of these assumptions (see, e.g., Raga & Cantó 1998) is correct for a “massless” working surface which instantaneously ejects most of the material sideways (into a jet cocoon), while the latter assumption (see Cantó et al. 2000) is appropriate for the “mass conserving” case in which the shocked material mostly stays within the working surface.

We also carry out axisymmetric numerical simulations of the flow, and compare the results with the ram-pressure balance and center of mass analytic models. This is the first time that a comparison between numerical simulations and the two analytic approximations (described above) has been attempted.

The paper is organized as follows. In §2 we present the center of mass and the ram-pressure balance solutions for a uniformly accelerated jet evolving in a homogeneous interstellar medium. In §3 we present the axisymmetric numerical simulations, and a comparison of the results with the analytic models is done in §4. A discussion of the results is held in §5, and concluding remarks are given in §6.

We consider a jet with an ejection velocity of the form

where $a$ is a constant acceleration and $\tau$ the ejection time. We assume a cylindrical flow, with a jet cross section $\sigma$ which is circular and constant. We also assume that the jet moves into an environment of uniform density $\rho_{a}$.

We consider two different forms for the ejection density:

1. 1.

   a constant mass loss rate per unit area ${\dot{m}}$, such that the ejection density is given by

   $$\rho_{0}(\tau)=\frac{\dot{m}}{u_{0}(\tau)},$$

   (2)

2. 2.

   a time-independent ejection density,

   $$\rho_{0}=\mathrm{const.}$$

   (3)

We assume that the flow in the jet beam is free-streaming, satisfying the condition:

where $u(x,t)$ is the velocity as a function of distance $x$ from the source at an evolutionary time $t$ and $u_{0}(\tau)$ is the ejection velocity (at the ejection time $\tau<t$). For an ejection time $\tau<0$, $x=0$, i.e., no material is ejected. Also, for a cylindrical, free-streaming flow, the density is given by:

where $\rho_{0}(\tau)$ is the ejection density (see equation 2 or 3) and ${\dot{u}}_{0}=du_{0}/d\tau$ is the derivative of the ejection velocity with respect to the ejection time (for a derivation of this equation, see Raga & Kofman, 1992). The relationship between the evolutionary time $t$ and the ejection time $\tau$ can be obtained from equations (1) and (4). Later, in section 2.1, we explicitly show this relationship.

When the flow starts (at $t=0$), a working surface is formed at $x=0$. This “jet head” then travels away from the outflow source for increasing times. In order to describe the time evolution of this working surface, we use two analytic approximations:

1. 1.

   center of mass- we assume that the position of the working surface coincides with the center of mass of the ejected, free-streaming material,

2. 2.

   ram-pressure balance- we assume that the motion of the working surface is determined by the jet/environment ram-pressure balance condition.

The first of these approximations is appropriate for a working surface in which the shocked gas mostly remains within the jet head, and the latter approximation is valid for a working surface that ejects most of the matter sideways (into a jet cocoon).

We therefore develop four analytic models, with the two ejection densities and the two analytic approximations described above. The models all share the linearly accelerating ejection velocity given by equation (1).

In order to check our analytical solutions we computed a set of 2D axisymmetric simulations for both the constant mass loss rate and constant ejection density cases. For the simulations we use a code that solves the ideal gas-dynamic (Euler) equations in a fixed two dimensional grid. The code uses a second order Godunov type method with the HLLC (Toro et. al., 1994) approximate Riemann solver, including a linear reconstruction of the primitive variables with the minmod slope limiter to avoid spurious oscillations. In order to use cylindrical coordinates, the appropriate geometrical source ($\propto 1/r$) terms are included after each timestep in an operator splitting fashion. Additionally, we incorporated the cooling function dependent on density, metallicity and temperature in the energy equation (also as source term) as proposed by Wang et al. (2014), which we compute assuming a solar metallicity.

In order to stabilize the method an artificial diffusion with a (dimensionless) value of $0.01$ was added to all of the equations, and a Courant number of $0.2$ was used. The code is written in fortran90 and is paralelized with the Message Passing Interface. We used a computational grid with a size of 0.05 and 0.2 pc along the $r$ and $x$ directions, respectively. The spatial resolution was $\sim 6.89\,\mathrm{au}$, corresponding to $600\times 6000$ cells along the radial and axial directions.

We computed four numerical simulations, two with a constant ejection density (models $n_{100}$ and $n_{5000}$) and two with a constant mass loss rate (models $\dot{m}_{100}$ and $\dot{m}_{5000}$). The top (purple background) and bottom (black background) subpanels of each subfigure presented in Figure 2 display the numerical density and temperature stratifications at an evolutionary time of 2500 yr for these four models.

Generally, the highest densities and temperatures are found in the on-axis regions of the jet heads. The jet heads of the models with constant ejection density (models $n_{100}$ and $n_{5000}$) have traveled to larger distances from the outflow source than the ones of constant mass loss rate ($\dot{m}_{100}$ and $\dot{m}_{5000}$). This qualitative difference is due to the fact that in the constant mass loss rate models the injected momentum rate ${\dot{M}}_{j}v_{j}$ scales linearly with $v_{j}$ (and hence also increases linearly with time, see equation 1), while this scaling is quadratic (with ejection velocity and with time) in the constant ejection density models. This means that the constant ejection density jets have a higher momentum content.

As a result of the lower jet beam density and pressure, the constant mass loss rate models develop an incident/reflected crossing shock structure, which (in the incident shock region) leads to the production of faster off-axis motions in the jet head (see the two bottom panels of Figure 2). This kind of structure is always found in jet simulations with appropriate parameters, i.e. simulations where the jet is under-dense and under-pressured with respect to the environment during its evolution (see, e.g., Raga 1988; Downes & Ray 1999).

The upper left panel of Figure 3 shows the position of the bow shock (dash-dotted blue line) and jet shock (dotted blue line) as a function of time for the ${\dot{m}}_{100}$ constant mass loss rate model. These positions were obtained by searching for the jumps along the symmetry axis of the pressure stratifications for both the bow shock and jet shock. We compare these shock positions with the analytical solution of the center of mass (dashed green line) and ram-pressure balance (solid red line) analytic models. We find that the position of both the bow shock lies close to the prediction from the center of mass analytic model, and the jet shock position lies somewhat below.

The upper right panel of Figure 3 shows the time-dependent position of the jet head for the $n_{100}$ constant ejection density model. For times $<1500$ yr, the jet and bow shock positions (obtained from the numerical simulation) quite closely follow the center of mass analytic solution, At larger times, the jet and bow shock positions lie above the center of mass solution, and begin to approach the ram-pressure balance solution. At all times, the positions of the jet and bow shock lie in between the working surface positions predicted by the two analytic models.

The lower left panel of Figure 3 shows the evolution of the jet head for the ${\dot{m}}_{5000}$ constant mass loss rate model. For this model, the positions of the jet and bow shock initially follow the center of mass analytic solution, but for times $t>1500$ yr begin to show large deviations, moving slower than the predictions of the two analytic models. These strong deviations are not surprising given the rather extreme departures from a constant cross section, cylindrical jet beam (assumed in the analytic models) found in the numerical simulation (see the bottom panel of Figure 2).

The strong departures from the simple, cylindrical structure assumed in the analytic models occurs because at increasing times the jet density (and pressure) drop quite considerably in this model, leading to the formation of internal “crossing shocks”. These shocks then lead to the formation of complex off-axis structures in the jet head (see the two bottom panels of Figure 2).

The lower right panel of Figure 3 shows the dynamics of the jet head for the $n_{5000}$ constant ejection density model. This model shows jet and bow shock positions that initially follow the center of mass analytic model, and at times $t>1500$ yr tend to approximate the ram-pressure balance model.

Finally, we calculate the relative position differences $\Delta x_{a}/x_{cm}=x/x_{cm}-1$, with $x_{cm}$ being the position of the “center of mass jet head” (see equations 17 and 21). The distance from the source $x_{a}$ corresponds either to the bow shock or jet shock positions (obtained from the numerical simulations), or to the position of the “ram-pressure balance jet head” (see equations 40 and 43). Figure 4 shows the relative position differences (with respect to the center of mass position for the jet head) for all of our computed models. The first thing to note (see the red lines in the two right panels of Figure 4) is that the relative position difference between the center of mass and ram pressure balance solutions does not strongly depend on time for the constant ejection density models (a result that can be seen by comparing equations 21 and 42). This result does not strictly hold for the constant mass loss rate models (red lines in the left panels of Figure 4).

For the constant mass loss rate models we see that for $t>550$ yr the bow and jet shock positions remain below the center of mass working surface position, with negative $\Delta x$ values (see the left panels of Figure 4). For the constant ejection density models (right panels of Figure 4) we see an initial convergence to the center of mass solution of the bow and jet shock positions, and for $t>1300$ yr we see progressively larger, positive values for $\Delta x/x_{cm}$ (starting to approach the ram-pressure balance working surface position).

We have presented a first attempt to compare the analytic ram-pressure balance (Raga & Kofman 1992) and center of mass (Cantó et al. 2000) analytic formalisms (for obtaining the motions of working surfaces in variable ejection jets) with axisymmetric numerical simulations. To this effect, we choose the relatively simple problem of the leading head of a jet produced with an ejection velocity that increases linearly with time. We explore the cases of a time-independent mass loss rate (in which the ejection density scales as the inverse of the ejection velocity) and of a constant ejection density. From the simulations we compute the on-axis positions of the jet shock and bow shock, and compare them with the predictions of the analytic models (in which the jet and bow shock are assumed to be spatially coincident).

For three of our four simulations (models $n_{100}$, $n_{5000}$ and ${\dot{m}}_{100}$), a good agreement is found between the numerical simulations and the analytic models. In these three models the jet and bow shock positions initially follow the center of mass analytic solution, and (for two of these models) at later times deviate towards the faster moving, ram-pressure balance analytic model. This result is completely satisfying, as one would expect that at early times the shocked jet and environment material will mostly remain within the working surface (as assumed in the center of mass formalism), and that a substantial sideways leakage of material (as assumed in the ram-pressure balance model) should occur at later times. This is discussed in more detail below.

It is clear that right after a working surface is formed, the two associated shocks are very close to each other (with a separation $d\ll r_{j}$, where $r_{j}$ is the local jet beam radius). Because the surface $S\approx 2\pi r_{j}l$ through which the gas exits the working surface is then small ($S\ll\pi r_{j}^{2}$), most of the mass going through the shocks remains within the working surface. At later times, the separation between the two shocks increases, and part of the mass processed by the working surface shocks is ejected sideways into a jet cocoon. For a jet with a time-indepenent ejection velocity, the separation $d$ between the two shocks attains a steady, maximum value determined by the balance between the incoming and outflowing mass rates. Such a balance between inflow and outflow from the working surface is also approximately attained in the case of a variable velocity jet, but the resulting shock separation (and also mass within the working surface) still has a (slower) time-dependence.

It is possible to obtain analytic estimates of the separation between the two shocks (and therefore the mass of the gas within the working surface) with different degrees of complexity (and accuracy), following the approaches developed for determining the shock standoff distance in blunt body flows (see, e.g., the book of Hayes & Probstein 2003 and Cantó & Raga 1998). Here, we present a very simple estimate of the mass within a working surface.

In the stationary state (i.e., after the initial regime of growing sepration between the two working surface shocks), the mass $M_{ws}$ within the two working surface shocks is:

where ${\dot{M}}_{in}$ is the mass being fed through the two shocks and $t_{s}$ is the timescale for sideways leakage of the material. Assuming that the two shocks are highly radiative (with a cooling distance¹¹1The cooling distance is usually defined as the distance from a plane-parallel shock to the point where the temperature has dropped to a certain value (see for example Hartigan et al., 1987). to $\sim 10^{3}$ K much smaller than the jet radius $r_{j}$), this leakage will occur at the post-cooling sound speed $c_{0}=\sqrt{5kT/3\mu m_{H}}\approx 3$ km s^-1 (for $T=10^{3}$ K and $\mu=1.3$), and we then have

Also, the total mass that has been fed into the working surface is

where $t_{d}=x_{ws}/v_{ws}$ (with $x_{ws}$ being the position and $v_{ws}$ the velocity of the working surface) is the dynamical timescale.

Combining equations (46-48) we obtain:

In this equation, we have considered that the $M_{ws}/M_{in}\leq 1$ condition has to be satisfied (with $M_{ws}=M_{in}$ for the initial, mass conserving regime, and $M_{ws}<M_{in}$ when the sideways ejection from the working surface becomes important), while our analytic estimate can give values $>1$ because of its implicit assumption of a balance between the mass rates into and out of the working surface (which is not satisfied in the early evolution of the working surface, as discussed above). Clearly, for the case of an accelerating working surface (as the one that we obtain from our variable ejection models), the derivation of equation (49) should be done evaluating the appropriate integrals (see section 4 of the paper of Lim et al. 2002), but the simple derivation shown above still gives the correct scaling properties.

We now consider that the transition between the mass conserving and efficient mass loosing working surface regimes (which can be described with the “center of mass” and “ram-pressure balance” equations of motion, respectively) occurs when $M_{ws}/M_{in}\sim 1/2$ (i.e., that the working surface has lost half of the processed jet material). Therefore, the position $x_{t}$ of the working surface at which the transition between the center of mass and ram pressure balance regimes occurs, can be estimated setting $M_{ws}/M_{in}=1/2$ in equation (49), which gives:

Setting $c_{0}=3$ km s^-1 (see above), $r_{j}=300$ au (see section 3) and a velocity $v_{ws}\approx 200$ km s^-1, we then obtain $x_{t}\approx 3\times 10^{17}$ cm, This is the correct order of magnitude for the location of the transition between the “mass conserving” center of mass and the “efficient mass loosing” ram-pressure balance regimes for the working surface found in our $n_{100}$ and $n_{5000}$ numerical simulations (see Figure 3).

To summarize, we find that numerical simulations of a jet with a linearly accelerating ejection velocity vs. time produce a leading head which at early times can be analytically modeled with the center of mass formalism, and (for appropriate flow parameters) at later evolutionary times approaches the ram-pressure balance solutions (the distance from the source at which the transition between these two regimes occurs can be estimated with equation (50)). This result holds provided that during its time-evolution the jet does not become under dense/under pressured, in which case strong departures from the analytic solutions are found.

This paper presents a first attempt to compare analytic solutions based on both the “ram-pressure balance” and “center of mass” formalisms with hydrodynamical (axisymmetric) simulations. For this comparison, we have studied the case of a linearly increasing ejection velocity as a function of time (and two different forms for the ejection density, see §2).

This simple ejection variability could apply to the early evolution of an HH outflow, and is relevant for the problem of entrainment of environmental molecular material into the jet head (see Lim et al. 2002) and the potential destruction of molecules by shocks, particularly in the jet head. As a direct application to the observed, evolved HH jets is not straightforward, the interest of our work is primarily theoretical.

While full analytic solutions for different forms of the ejection variability have been found in the past with the “center of mass” formalism (see Cantó et al, 2000 and Cantó & Raga 2003), no solutions (except the trivial one for the constant ejection case) have been previously found with the “ram-pressure balance” equation of motion for the working surfaces. We were surprised to find that the “linearly accelerating” ejection velocity case does lead to a full, analytic ram-pressure balance solution for the case of a constant ejection density (though we managed to find only an approximate analytic solution for the constant mass loss rate case). This result implies that ram-pressure balance analytic solutions might also exist for other ejection velocity/density variability combinations.

Comparing the center of mass and the ram-pressure balance solutions with axisymmetric numerical simulations (with the same parameters) we find that for the case of over-dense jets (with higher densities in the jet side of the leading working surface) the numerical results initially follow the center of mass solution, and at later times approach the ram-pressure balance analytic solution. This result shows a satisfying consistency between the analytic approaches and the numerical solutions with the full (axisymmetric) hydrodynamic equations. Not entirely unexpectedly, we do not find a good agreement between the analytic and numerical results for cases in which the jet density drops to values lower than the environmental density, as the jet beam then develops internal shocks (in the numerical solutions) that affect the motion of the leading working surface.

Clearly, it would be interesting to extend this comparison between numerical simulations and analytic solutions of variable jets to the case of jets with “internal working surfaces” within their beams. This problem, which would clearly have much more direct applications to observed HH jets, is left for a future study.

ARa and JC acknowledge support from the DGAPA-UNAM grant IG100218. We thank an anonymous referee for helpful comments which lead (among other things) to many of the issues discussed in section 5.

In §2.1 we established that for a constant mass loss rate the ram-pressure balance assumption leads to equation (37). In terms of dimensionless variables (see equation 38) the equation of motion takes the form presented in (39). By integrating it numerically we obtain the $\eta(y)$ dependency shown the top frame of Figure 5.

This figure also shows the“near” and “far field” solutions

which correspond to the $y\rightarrow 0$ (near) and $y\gg 1$ (far) limits, respectively. These two solutions can be obtained straightforwardly taking the appropriate limits in equation (39).

Since the solutions in (51) represent two regimes with analytical integrals, a non-linear average between this approximations is computed:

to obtain an approximation to the full $\eta(y)$ solution. These average has a relative error $\epsilon_{rel}=\eta^{(a)}/\eta-1<0.04$ with respect to the “exact” solution $\eta(y)$ (obtained by numerical integration of equation 39), and is therefore a reasonable analytical approximation of the solution. The dependency of $\epsilon_{rel}$ on $y$ is shown in the bottom panel of Figure 5.

Another approximation to the numerical $\eta(y)$ solution is

which coincides with the exact solution in the $y\rightarrow 0$ and $y\rightarrow\infty$ limits, and has a maximum relative deviation of $\approx 0.09$. This approximation has a functional form similar to the solution for the constant mass loss rate case obtained with the “center of mass formalism” (see equation 15).