SciELO - Scientific Electronic Library Online

vol.57 número2On the effects of solar radiation pressure on the deviation of asteroidsLocation and stability of the triangular points in the triaxial elliptic restricted three-body problem índice de autoresíndice de materiabúsqueda de artículos
Home Pagelista alfabética de revistas  

Servicios Personalizados




Links relacionados

  • No hay artículos similaresSimilares en SciELO


Revista mexicana de astronomía y astrofísica

versión impresa ISSN 0185-1101

Rev. mex. astron. astrofis vol.57 no.2 Ciudad de México oct. 2021  Epub 15-Feb-2022 


Estimating the propagation of a uniformly accelerated jet

J. I. Castorena1 

A. C. Raga1 

A. Esquivel1 

A. Rodríguez-González1 

L. Hernández-Martínez2 

J. Cantó3 

F. Clever3 

1Instituto de Ciencias Nucleares, UNAM, México.

2Facultad de Ciencias, UNAM, México

3Instituto de Astronomía, UNAM, México.


We study the problem of a Herbig-Haro jet with a uniformly accelerating ejection velocity, travelling into a uniform environment. For the ejection density we consider two cases: a time-independent density, and a time-independent mass loss rate. For these two cases, we obtain analytic solutions for the motion of the jet head using a ram-pressure balance and a center of mass equation of motion. We also compute axisymmetric numerical simulations of the same flow, and compare the time-dependent positions of the leading working surface shocks with the predictions of the two analytic models. We find that if the jet is over-dense and over-pressured (with respect to the environment) during its evolution, a good agreement is obtained with the analytic models, with the flow initially following the center of mass analytic solution, and (for the constant ejection density case) at later times approaching the ram-pressure balance solution.

Key Words: Herbig-Haro objects; ISM; jets and outflows


Estudiamos el problema de un yet Herbig-Haro uniformemente acelerado, que se mueve a través de un medio homogéneo. Consideramos dos casos para la densidad de expulsión: una densidad independiente del tiempo, y una tasa de pérdida de masa constante. En ambos casos utilizamos los formalismos de equilibrio de presión hidrodinámica y de posición de centro de masa para obtener soluciones analíticas. También realizamos simulaciones numéricas de este flujo, y comparamos la dependencia temporal de la posición de los choques que conforman la superficie de trabajo con las predicciones analíticas. Encontramos que si la densidad y la presión del yet son mayores que los correspondientes valores del medio, durante toda la evolución, entonces se obtiene una buena coincidencia con los modelos analíticos: en un inicio el flujo sigue la solución analítica de centro de masa y en tiempos posteriores se acerca a la solución de equilibrio de presión hidrodinámica.


Jets from young stars and their associated HerbigHaro 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 typically occurred at times of ≈ 103 → 104 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 timevariability 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 H2 molecules in the resulting, accelerating jet head, presenting 1D hydrodynamic simulations including an H2 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 super-alfvénic leading head that evolves consistently with the computed timedependence 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. 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

u0(τ)={0,τ<0,aτ,τ0, (1)

where a is a constant acceleration and τ the ejection time. We assume a cylindrical flow, with a jet cross section σ which is circular and constant. We also assume that the jet moves into an environment of uniform density ρ a .

We consider two different forms for the ejection density:

  • 1. a constant mass loss rate per unit area m ̇, such that the ejection density is given by

ρ0τ=m˙u0τ, (2)

  • 2. a time-independent ejection density,

ρ0=const. (3)

We assume that the flow in the jet beam is freestreaming, satisfying the condition:

ux,t=xt-τ=u0τ, (4)

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

ρx,t=ρ0u0u0-t-τu˙0, (5)

where ρ 0(τ) is the ejection density (see equation 2 or 3) and ˙u 0 = du 0 /dτ 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 τ can be obtained from equations (1) and (4). Later, in § 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. 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. 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).

2.1. Center of Mass Equation of Motion

We first consider a “mass conserving” working surface. For this case, we assume that the material going through the bow and jet shocks mostly stays within the working surface, with only a small fraction of this material being ejected sideways into the jet cocoon. The working surface position should then coincide with the center of mass of the free-streaming fluid parcels that have piled up within it. This center of mass position is given by:

xcm=xdmdm, (6)

with the differential element of mass given by:

dm=σρ0τu0τdτ+σρaxdx. (7)

In this equation, the first term on the right-hand-side corresponds to the ejected material, and the second term to the swept-up, stationary, environment entering the working surface. In equation (7), σ is the jet cross section, ρ 0(τ) is the ejection density, u 0(τ) the ejection velocity and ρ a (x) the (possibly positiondependent) environmental density.

Using equations (4-7) and considering that the ejection time τ 0 is integrated from 0 to τ, we obtain:

xcm0τρ0τ'u0τ'dτ'+0xcmρadx=0τxjρ0τ'u0τdτ'+0xcmxρadx, (8)

where we have set ρ a = const., and x j is the position that the fluid parcels would have if they were still free-streaming:

xj=t-τ'u0τ', (9)

see equation (4).

Also, from equations (1) and (4) we find:

t=xcmaτ+τ. (10)

Equation (10) can be inverted to obtain

τ=12t+t2-4xcma. (11)

Combining equations (8-10), we obtain:

xcm0τρ0τ'aτ'dτ'+ρaxcm2=0τt-τ'a2τ'2ρ0τ'dτ'. (12)

In order to carry out the remaining integrals, we have to specify ρ 0(τ’). We consider two forms for the ejection density (see equation 2).

  • (a) Constant mass loss rate.

Equation (12) takes the form:

xcmm˙τ+ρaxcm2=m˙aτ2t2-τ3. (13)

Using equations (10) and (13), we obtain

xcm2+m˙τρaxcm-m˙aτ33ρa=0, (14)

this is a quadratic equation with one positive solution, given by

xcmτ=m˙τ2ρa-1+1+4aρaτ3m˙. (15)

In the low density limit ρ a → 0, equation (14) takes a particularly simple form that leads to finding the solution

xcm=aτ23. (16)

This solution can also be obtained by carrying out a first-order Taylor series expansion of equation (15).

By substituting equation (11) into (13), we can also obtain x cm as an explicit function of t, giving

xcm=89xc1+34ttc3/2-1+98ttc, (17)

with t c m/aρ˙ a and x c m˙ 2 /aρ 2 a .

In the ttc limit, a second order Taylor series expansion of equation (17) gives:

xcm316at2. (18)

  • (b) Constant ejection density.

With a constant ρ 0, equation (12) takes the form:

xcmρ0aτ22+ρaxcm2=ρ0a2τ3t3-τ4, (19)

where t and τ are related through equation (10).

Substituting t as a function of τ we obtain:

xcm2+ρ0aτ23ρaxcm-ρ0a2τ46ρa=0, (20)

which has a positive solution

xcmτ=ρ0aτ26ρa-1+1+6ρaρ0. (21)

In the ρ a → 0 limit, equation (20) (or, alternatively, equation 21) takes the form

xcm=aτ22, (22)

with a quadratic dependency on the ejection time.

We can also find a solution to (19) as a function of evolutionary time t (using equation 11), giving:

xcm=a9β0β0β02-18+β02+63/2β02-22t2, (23)


β0=ρ0ρa, (24)

and therefore the position x cm has a quadratic dependency on the evolutionary time t.

Equation (23) takes the form of a constant acceleration condition, with the acceleration given by

g=2a9β0β0β02-18+β02+63/2β02-22, (25)

which in the ρ a → 0 low density limit becomes,

g231/2aβ0. (26)

2.2. Ram-Pressure Balance Equation of Motion

Figure 1 shows a schematic diagram of a jet head propagating through a stationary environment. This leading working surface has a structure consisting of two shocks: an upstream shock known as the jet shock (or Mach disk), which is slowing down the jet material, and a downstream bow shock which is accelerating the surrounding material. The two-shock working surface structure travels away from the outflow source at a velocity v ws .

Fig. 1 Jet head working surface in a reference frame at rest with respect to the outflow source (left) and in a reference frame moving along with the working surface (right). The bow shock is shown in red and the jet shock in blue. The color figure can be viewed online. 

The right panel of Figure 1 shows the situation seen in a reference frame at motion with the jet head. For a hypersonic flow, the two working surface shocks are strong, so that the post-shock gas pressures are given by:

Pbs=2γ+1ρavws2, (27)

for the bow shock, and

Pjs=2γ+1ρjvj-vws2, (28)

for the jet shock (assuming that the jet and the environment have the same specific heat ratio γ). In this equation, v j and ρ j are the velocity (with respect to the outflow source) and the density of the material that is presently entering the jet shock.

If the working surface is moving with a constant velocity the condition

Pbs=Pjsρavws2=ρjvj-vws2, (29)

is satisfied. This is called the ram-pressure balance condition. This condition is also valid for a working surface moving with a variable velocity as long as the inertia of the material between the bow shock and the jet shock is negligible, which is the case if most of the material is ejected sideways from the working surface into a jet cocoon.

From equation (29) we find that

vws=βvj1+β,withβρjρa. (30)

Clearly, the shock velocity associated with the bow shock is v bs = v ws and the shock velocity of the jet shock is:

vjs=vj1+β. (31)

We now consider the equation of motion v ws = dx ws /dt for the working surface. Setting v j = u 0(τ) (i.e., the velocity of the material entering the working surface is equal to the ejection velocity at the corresponding ejection time τ), we then have:

1+ρaρjdxwsdt=u0τ. (32)

The ejection time τ of the material entering the working surface is obtained from the free-streaming relation


where t is the evolutionary time (see equation 4). Solving for t and differentiating with respect to τ we obtain

dtdτ=1+1u0dxwsdτ-xwsu02du0dτ. (33)

Combining equations (32) and (33), we obtain:

dxwsdτρaρj=u0-xwsdlnu0dτ. (34)

The density of the jet material entering the working surface is ρ j = ρ(x ws ,t) (see equation 5), which for our chosen ejection velocity variability (see equation 1) becomes:

ρj=ρ01-xwsaτ2. (35)

Finally, from equations (34) and (35) we obtain:

dxwsdτ=ρ0ρaaτ1-xwsaτ2. (36)

In order to proceed it is now necessary to specify the form of ρ 0(τ). We therefore consider the two cases of constant mass loss rate and constant ejection density.

  • (a) Constant mass loss rate.

For m ̇ = const. equation (36) takes the form:

dxwsdτ=m˙ρaaτ-xwsτ. (37)

It is convenient to use the dimensionless variables:

η=aρa2m˙2xws,y=aρam˙τ. (38)

In terms of these variables, equation (37) is

dηdy=y-ηy. (39)

In Appendix A, it is shown that an approximate analytic solution of equation (39) can be constructed as a non-linear average of a “near” and a “far field” analytic solution (see equations 51 and 52), which in terms of the respective dimensional variables is

xws=m˙2aρa2aρam˙τ-52+3254aρam˙τ-158-45. (40)

  • (b) Constant ejection density.

Setting ρ 0 = const. in (36) and taking the square of the equation, we obtain:

ρaaρ0dxwsdτ2+xws-aτ2=0. (41)

Proposing a power law solution, we find that

xws=aρ08ρa-1+1+16ρaρ0τ2. (42)

We therefore find a quadratic dependency on the ejection time.

Using equations (11) and (42) we then obtain:

xws=8aβ02-1+1+16β028+β02-1+1+16β022t2, (43)

where β 0 is given by equation (24). This implies a quadratic dependence on the evolutionary time.

In the β01 limit, equation (43) becomes

xwsa2β0t2, (44)

and the β0≫1 limit leads to

xwsa4t2. (45)


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 (∝ 1/r) terms are included after each time step 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 parallelized with the Message Passing Interface. We used a computational grid with a size of 0.05 and 0.2pc along the r and x directions, respectively. The spatial resolution was ≈ 6.89au, corresponding to 600 × 6000 cells along the radial and axial directions.

3.1. Initial and Boundary Conditions

We model jets propagating into a uniform, quiescent environment, and consider the values n a = 100 cm−3 and n a = 5000 cm−3 for the environmental density. For convenience, from now on we will use number density instead of mass density. We consider a mean molecular weight of µ = 1.3. We also consider the cases of jets with a constant mass loss rate and a constant ejection density.

In all simulations, we consider the ejection velocity variability given by equation (1) with a = 100 kms−1 per millenium = 3.17 × 10−4 cms−2. An initial jet radius r j = 300AU (the cross section being σ = πr j 2) and a temperature T = 100K (for both the jet and the environment) were imposed in all models.

We therefore compute four models, two with constant ejection density (which we label n 100 and n 5000, with the subscript giving the ambient density in cm−3) and two with constant mass loss rate (which we label m ̇100 and m ̇5000). For the models m˙ 100 and m˙ 5000 we choose a total mass loss rate M˙=2.24×10-8Myr-1, which corresponds to a mass loss rate at the jet source per unit area m˙ = 3.36 × 10−13 gcm−2 s−1. For the models n 100 and n 5000, we choose a density n j = 1.66×104 cm−3 (for a gas with 90% H and 10% He). The model parameters are summarized in Table 1.

In the constant mass flux cases, we find a maximum jet density of n j = 1.56×107 cm−3, and a minimum jet density of n j = 6.24×103 cm−3, corresponding to an ejection time τ = 1yr, and τ = 2500yr, respectively.

We should note that the parameters we have chosen are inspired in the characteristics of “classical” HH jets such as HH 34 and HH 111, which have:

  • radii ≈ 100 AU (Reipurth et al. 2002) and ≈ 300 AU (Reipurth et al. 1997), for HH 34 and HH 111, respectively. We chose the larger of these two radii in order to have a better resolution of the jet in our numerical simulations. These jet radii are measured in the “jet knot chain” at distances of ≈ 10'' (corresponding to ≈ 6 × 1016 cm) from the outflow sources, and not in the considerably broader HH 34S and HH 111V “heads”;

  • full spatial velocities (i.e., combining radial velocities and proper motions) between 200 and 300 km s−1 (Hartigan et al. 2001; Reipurth et al. 2002) and lengths (out to HH 34S and HH 111V) of ≈ 8 × 1017 cm. We have then made a choice of a = 100kms−1 per millenium (see above), which after an evolutionary time of ≈ 2500 yr produces a jet head with a position and velocity similar to the ones of HH 34S and HH 111V (see the following section). The simulations are stopped at this time (i.e., at 2500 yr);

  • mass loss rates in the 10 10-810-7M range (depending on the emission lines used for deriving this value and on the position along the jets, see Podio et al, 2006). We have therefore chosen densities that produce mass loss rates of this order of magnitude;

  • even though the ambient densities in HH 34 and HH 111 are of ≈ 10 cm−3 (19), we have chosen higher environmental densities (of 100 and 5000 cm−3, see above) in order to have a substantial braking effect and to obtain larger differences between the center of mass and rampressure balance solutions.


Mj˙ nj na
Myr-1 cm-3 cm-3
n100 1.30 x 10-7+ 1.66×104 100
n5000 5000
m˙100 3.35×10-7 1.55×104+ 100
m˙5000 5000

*These values correspond to an ejection time τ = 1000yr.


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 m ̇100 and m ̇5000). The top two (purple background) and bottom two (black background) panels presented in Figure 2 display the numerical density and temperature stratifications at an evolutionary time of 2500 yr for these four models.

Fig. 2 Numerical density and temperature stratifications obtained from our four models (see the text and Table 1), at an evolutionary time of 2500 years. For each model, we show the number density (top half) and temperature (bottom half) in logarithmic color scales. The axial and radial axes are given in units of 1017 cm. In the constant mass loss rate models the jet density and the pressure drop quite considerably, leading to the formation of internal “crossing shocks”. The color figure can be viewed online. 

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 distances from the outflow source larger than the ones of constant mass loss rate (m ̇100 and m ̇5000). This qualitative difference is due to the fact that in the constant mass loss rate models the injected momentum rate 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 larger 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 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 the bow shock lies close to the prediction from the center of mass analytic model, and the jet shock position lies somewhat below.

Fig. 3 Position of the jet head as a function of time for the constant mass loss rate models (left panels) and the constant ejection density models (right panels). The red solid line shows the ram-pressure balance solution and the green dashed line shows the center of mass solution. The dash-dotted line shows the bow shock and the dotted line the jet shock positions obtained from the numerical simulations. The color figure can be viewed online. 

The upper right panel of Figure 3 shows the timedependent 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 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 more slowly 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 ∆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).

Fig. 4 The relative differences of the ram-pressure balance solution (solid red line), the numerical bow shock (dashdotted blue line) and jet shock positions (dotted blue line) with respect to the center of mass working surface solution (the dashed, green line showing the ∆x/x cm = 0 line) as a function of time t. The results obtained for our four models (see Table 1) are shown (constant mass loss rate models on the left, and constant ejection density models on the right). The color figure can be viewed online.  

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 ∆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 ∆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) formalisms (for obtaining the motions of working surfaces in variable ejection jets) with axisymmetric numerical simulations. To this effect, we chose the relatively simple problem of the leading head of a jet produced with an ejection velocity that increases linearly with time. We explored 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 computed the on-axis positions of the jet shock and bow shock, and compared 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 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 satisfactory, 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 drj, where r j is the local jet beam radius). Because the surface S ≈ 2πr j l through which the gas exits the working surface is then small (Sπrj2), 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-independent 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) timedependence.

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 separation between the two working surface shocks), the mass M ws within the two working surface shocks is:

MwsM˙ints, (46)

where 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 distance4 to ≈ 103 K much smaller than the jet radius r j ), this leakage will occur at the post-cooling sound speed c0=5kT/3μmHkm s-1 (for T=103 K and μ=1.3), and we then have

ts=rj2c0. (47)

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

MinM˙intd, (48)

where t d = x ws /v ws is the dynamical timescale (with x ws being the position and v ws the velocity of the working surface).

Combining equations (46-48) we obtain

MwsMinmin1,rj2xwsvwsc0 (49)

In this equation, we have considered that the M ws /M in ≤ 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 § 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 losing 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 ≈ 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:

xt=vwsc0rj. (50)

Setting c 0 = 3 km s−1 (see above), r j = 300 au (see § 3) and a velocity v ws ≈ 200 km s−1, we then obtain x t ≈ 3 × 1017 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 losing” 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. Since 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 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 rampressure 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 rampressure 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 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 § 5.


Ahmed, I. & Shibata, K. 2008, PASJ, 60, 871 [ Links ]

Cantó, J., Raga, A.C., & D’Alessio, P. 2000, MNRAS, 313, 656 [ Links ]

Cantó, J. & Raga, A. C. 1998, MNRAS , 297, 383 [ Links ]

______. 2003, RMxAA, 39, 261 [ Links ]

Downes, T. P. & Ray, T. P. 1999, A&A, 345, 977 [ Links ]

Frank, A., Ray, T. P., Cabrit, S., et al . 2014, in Protostars and Planets VI, ed. H. Beuther, R. S. Klessen, C. P. Dullemond, and T. Henning (Tucson, AZ: UAP), 451 [ Links ]

Hartigan, P., Morse, J. A., Reipurth, B., Heathcote, S., & Bally, J. 2001, ApJ, 559, 157 [ Links ]

Hartigan, P., Raymond, J., & Hartmann, L. 1987, ApJ , 316, 323 [ Links ]

Hayes, W. D. & Probstein, R. F. 2003, Hypersonic Inviscid Flow (Mineola, NY: Dover) [ Links ]

Heathcote, S. & Reipurth, B. 1990, AJ, 104, 2193 [ Links ]

Herbig, G. 1974, LicOB, 658, 1 [ Links ]

Hiromitsu, K. & Shibata, K. 2005, ApJ , 634, 879 [ Links ]

Lim, A. J., Raga, A. C., Rawlings, J. M. C., & Williams, D. A. 2002, MNRAS , 335, 817 [ Links ]

Micono, M., Massaglia, S., Bodo, G., Rossi, P., & Ferrari, A. 1998, A&A , 333, 1001 [ Links ]

Ouyed, R., Clarke, D. A., & Pudritz, R. E. 2003, ApJ , 582, 292 [ Links ]

Podio, L., Bacciotti, F., Nisini, B., et al. A&A , 456, 189 [ Links ]

Raga, A. C. 1988, ApJ , 335, 820 [ Links ]

Raga, A. C., Cantó, J., Binette, L., & Calvet, N. 1990, ApJ , 364, 601 [ Links ]

Raga, A. C. & Binette, L. 1991, RMxAA , 22, 265 [ Links ]

Raga, A. C. & Cantó, J. 1998, RMxAA , 34, 73 [ Links ]

Raga, A. C. & Kofman, L. 1992, ApJ , 386, 222 [ Links ]

Reipurth, B. 1991, ESOC, 33, 247 [ Links ]

Reipurth, B. & Bally, J. 2001, ARA&A , 39, 403 [ Links ]

Reipurth, B. & Heathcote, S. 1990, A&A , 229, 527 [ Links ]

______. 1991, A&A, 246, 511 [ Links ]

Reipurth, B., Hartigan, P., Heathcote, S., Morse, J. A., & Bally, J. 1997, AJ, 114, 757 [ Links ]

Reipurth, B., Heathcote, S., Morse, J., Hartigan, P., & Bally, J. 2002, AJ, 123, 362 [ Links ]

Romanova, M. M., Blinova, A. A., Ustyugova, G. V., Koldoba, A. V., & Lovelace, R. V. E. 2018, NewA, 62, 94 [ Links ]

Toro, E. F., Spruce, M., & Speares, W. 1994, ShWav, 4, 25 [ Links ]

Wang, Y., Ferland, G. J., Lykins, M. L., et al. 2014, MNRAS, 440, 3100 [ Links ]

Zanni, C., Ferrari, A., Rosner, R., Bodo, G., & Massaglia, S. 2007, A&A , 469, 811 [ Links ]

4 The 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).


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 η(y) dependency shown the top frame of Figure 5.

Fig. 5 Top frame: the exact (i.e., numerical) solution of equation (39) (solid curve), the “near” (short dashes) and the “far field” (long dashes) analytic solutions given by equation (51). Bottom frame: relative deviation of the approximate solution of equation (52) with respect to the “exact” (numerical) solution of equation (39). 

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

ηny=y2,ηfy=23y3/2, (51)

which correspond to the y → 0 (near) and y ≫ 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:

ηa=ηn-5/4+ηf-5/4-4/5, (52)

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

Another approximation to the numerical η(y) solution is

ηb=2y9-1+1+9y, (53)

which coincides with the exact solution in the y → 0 and y →∞ limits, and has a maximum relative deviation of ≈ 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).

Received: October 13, 2020; Accepted: March 29, 2021

J. Cantó and F. Clever: Instituto de Astronomía, Universidad Nacional Autónoma de México, Apdo. Postal 70-264, 04510, CDMX, México (

J. I. Castorena: Instituto de Ciencias Nucleares, Universidad Nacional Autónoma de México, Apdo. Postal 70-543, 04510, CDMX, México (

A. Esquivel, A. C. Raga, and A. Rodríguez-González: Instituto de Ciencias Nucleares, Universidad Nacional Autónoma de México, Apdo. Postal 70-543, 04510, CDMX, México (raga, ary,

L. Hernández-Martínez: Facultad de Ciencias, Universidad Nacional Autónoma de México, Av. Universidad 3000, Circuito Exterior S/N, Delegación Coyoacán, C.P. 04510, Ciudad Universitaria, CDMX, México (

Creative Commons License This is an open-access article distributed under the terms of the Creative Commons Attribution License