Äîêóìåíò âçÿò èç êýøà ïîèñêîâîé ìàøèíû. Àäðåñ îðèãèíàëüíîãî äîêóìåíòà : http://www.mrao.cam.ac.uk/~krause/krause04.aa.small.pdf Äàòà èçìåíåíèÿ: Mon Feb 14 22:05:32 2005 Äàòà èíäåêñèðîâàíèÿ: Mon Oct 1 22:34:58 2012 Êîäèðîâêà: Ïîèñêîâûå ñëîâà: êàñêàä êåìáëà

Astronomy & Astrophysics manuscript no. krause (DOI: will be inserted by hand later)

25th October 2004

Very Light Jets II: Bipolar large scale simulations in King atmospheres
Martin Krause
1 2

1,2

Landessternwarte K¨ onigstuhl, D-69117 Heidelberg, Germany Astrophysics Group, Cavendish Laboratory, Madingley Road, Cambridge CB3 0HE, United Kingdom

Received 25th October 2004/ Accepted ¡date¿ Abstract Hydrodynamic jets, underdense with respect to their environment by a factor of up to 10 4 , were computed in axisymmetry as well as in 3D. They finally reached a size of up to 220 jet radii, corresponding to a 100 kpc sized radio galaxy. The simulations are "bipolar", involving both jets. These are injected into a King type density profile with small stochastic density variations. The back-reaction of the cocoons on the beams in the center produces armlength asymmetries of a few percent, with the longer jets on the side with the higher average density. Two distinguishable bow shock phases were observed: an inner elliptical part, and a later cylindrical, cigar-like phase, which is known from previous simulations. The sideways motion of the inner elliptical bow shock part is shown to follow the law of motion for spherical blast waves also in the late phase, where the aspect ratio is high, with good accuracy. X-ray emission maps are calculated and the two bow shock phases are shown to appear as rings and elongated or elliptical regions, depending on the viewing angle. Such structures are observed in the X-ray data of several radio galaxies (e.g. in Abell 2052 and Hercules A), the best example being Cygnus A. In this case, an elliptical bow shock is infered from the observations, a jet power of 1047 erg/s is derived, and the Lorentz factor can be limited to > 10. Based on the simulation results and the comparison to the observations, the emission line gas producing the alignment effect in radio galaxies at high redshift is suggested to be cooled gas entrained over the cocoon boundary. Key words. Hydrodynamics ­ Instabilities ­ Shock waves ­ Galaxies: active­Radio continuum: galaxies­X-rays: galaxies: clusters

1. Intro duction
The research on simulations of extragalactic jets has been extended recently into the regime of very light jets (compare e.g. Krause 2002; Saxton et al. 2002; Carvalho & O'Dea 2002a,b; Zanni et al. 2003; Bicknell et al. 2003), down to density contrasts (jet over ambient density) of = 10-5 (Krause 2003). The interest has been stimulated not only by the difficulty to fix that parameter for real sources, but also by a first success in explaining some parameters of observed extragalactic jets. Low jet densities are needed to get large radio cocoon and bow shock widths. Based on the results from a grid of simulations and comparison to the radio and X-ray data of the radio galaxy Cygnus A (Carilli & Barthel 1996; Smith et al. 2002), a density contrast of roughly < 10-3 has been claimed for this source (Rosen et al. 1999; Krause 2003). All these simulations have been carried out in axisymmetry and a scale of a few dozen jet radii. The aim of this paper is
Send offprint requests to : M.Krause@lsw.uni-heidelberg.de M.Krause, e-mail:

to verify and extend the earlier simulation results on the 100 kpc (=200 jet radii, Rj ) scale and in three dimensions. Krause (2003) found that very light jets start their life in a blastwave-like phase. During that phase, terminating when the bow shock reaches R1 Rj /2 1/4 , the bow shock is spherical because of pressure predominance, obeying the blastwave equation of motion:
r 0 t t

M(r )r dr = 2

dt
0 0

E (t )dt

.

(1)

Here, M(r) is an arbitrary spherically symmetric ambient gas mass distribution and E (t) is the energy injection law. For stationary energy injection E = Lt and constant ambient density 0 , one gets: r= 5Lt3 4 0
1/5

.

(2)

This equation relates the total luminosity, L, of the jet to observables like the bow shock radius and its velocity (via time derivation of (2)). With information on the bow shock propagation at later phases, it can be hoped that


2

M. Krause: Very Light Jets II
vp ext ext ext R grid boundary
initial condition: grid boundary jet inflow

the jet power could be reconstructed from the bow shock shape, which is probably already observed in the case of Cygnus A by Chandra (see sect. 5). It is therefore important to check the evolution of the bow shock at later phases. This was difficult for most simulations so far because of the assumed reflection symmetry in the equatorial plane of the system. Interaction of the backflow with this boundary disturbs the sideways evolution of the bow shock. Therefore, the symmetry assumption has been dropped in the present work, and the simulations are bipolar, evolving back-to-back jets in opposite directions, a technique that has emerged only recently (Reynolds et al. 2001; Krause & Camenzind 2002, 2003; Basson & Alexander 2003). A fully 3D simulation is presented in sect. 3, and a large scale axisymmetric simulation is presented in sect. 4. Comparisons to observations are presented in sect. 5. Part of this is an update of previously published results, using a more involved analysis.

Z

v p jet jet jet

Figure 1. Sketch of the cylindrical grid.

ities, the following boundary conditions were applied: f (iZ , iR , i ) := f (iZ , 6 - iR , i + i ) , iR = 0, 1, 2 f (iZ , 3, i ) := < f (iZ , 3, i ) > |iZ =const , where i denotes the number of grid cells corresponding to the angle . The second equation indicates that directly on the axis , the scalars were averaged over i for constant iZ , since all the different i refer to the same physical point. Due to the staggered mesh, the Z and components of vectors are shifted by half a grid cell away from the axis. Therefore for these components, there arise no special problems, and the boundary condition is: v
Z,

2. Numerics
The magneto-hydrodynamic (MHD) code NIRVANA was employed (Ziegler & Yorke 1997). It solves the hydrodynamic equations in three dimensions (3D) for density , velocity v , and internal energy e: + · (v ) = 0 t v + · (v v ) = - p - t e + · (ev ) = -p · v - 2 , t (3) (4) (5)

(iZ , iR , i ) := v

Z,

(iZ , 5 - iR , i + i )

, iR = 0, 1, 2

The radial vector components are not shifted away from the axis. So in principal, they all denote the same physical point. However, the flow should be allowed to cross the axis from one side to the other. This is only possible if the radial velocity takes a reasonable value there. Therefore, two possibilities arise for the boundary conditions: vR (iZ , iR , i ) := -vR (iZ , 6 - iR , i + i ) vR (iZ , 3, i ) := , iR = 0, 1, 2 I II

where denotes an external gravitational potential and is cooling function. Here, bremsstrahlung ( = 7.5 â a 1020 T (1. + 4.4 â 10-10 T ) cm3 erg/s) has been used for the 3D computation. The code was vectorised and parallelised by OpenMP like methods, and successfully tested on the NEC SX-5 (Krause & Camenzind 2002) at the high performance computing center in Stuttgart (Germany), where the computations have been carried out.

0: (vR (iZ , 4, i ) - vR (iZ , 2, i ))/2 :

2.1. Boundary conditions for the 3D cylindrical grid
For the 3D simulation, a cylindrical grid was employed. The disadvantage of the cylindrical coordinates (Z, R, ) compared to the Cartesian ones is the appearance of internal boundaries. The grid somehow has to be connected across these boundaries. In the direction, periodic boundary conditions were applied. For the boundary on the axis, no analogue could be found in the literature. The boundary condition here is similar to the periodic case: One side of the grid should know about the other side. Therefore, three cells were used below the axis, to which consequently the index iR = 3 was assigned. With this choice and the staggered mesh, equations (3-5) are well defined everywhere on the grid. For the scalar quant-

In order to check the influence of these two different boundary conditions on the simulation, the 3D jet was simulated twice, the first time with the case I, and the second time with the case II boundary condition. For both runs, the integrated quantities (directionally split energy and momentum, and mass) and also the timesteps were equal. Also, from comparison of the contour plots, no difference could be found. Hence, the flow seems to have enough possibility to flow past the axis, and the detailed behaviour at that line does not influence the result by much. The simulation result also shows that this approach in general works in regions of undisturbed jet flow. Where the jet is dominated by instabilities and wants to bend, the axis looks like an obstacle. Therefore, the method seems to be problematic if details about the jet beam are of special interest. Nevertheless, the results confirm that the representation of the beam is acceptable, and for the other regions the output was fine. The big advantage of the cylindrical grid is the reduction of the number of cells to compute. The reduction of the physical computational volume is 22%. Since the -direction can be covered with


M. Krause: Very Light Jets II

3

fewer cells than a third dimension in Cartesian coordinates, this approach saves a significant amount of memory and CPU time. The computation would not have been possible on a 3D Cartesian grid of similar central resolution.
R [kpc]

8 7 6 5 4 3 2 1 0 0 0.5 radial bow shock position global fit: 1.02 + 4.44 t0.58 1 1.5 2

3. 3D simulation 3.1. Simulation setup
A cylindrical grid was used for the jet simulation (compare section 2.1). The size of the computational domain was: Z [-69 kpc, 69 kpc], R [0, 57 kpc] and [0, 2 ]. 2042, 805, and 57 grid points were used in the Z,R and directions, respectively. With a jet radius of rj = 0.55 kpc, this gives a resolution of 8 points per beam radius (ppb). The resolution in direction scales linearly down from 16 ppb on the jet boundary to 0.2 ppb on the edge of the grid. The grid was initialised with an isothermal King cluster atmosphere: e (r) =
e,0

Figure 3. Radius where the bow shock hits the line ( = 0, Z = 0) versus time, with a globally fitting function.

1+

r2 a2

-3 /2

,

(6)

where r = R2 + Z 2 denotes the spherical radius, e,0 = 1.2 â 10-25 g/cm3 is the characteristic density, = 0.75 and a = 35 kpc is the core radius. In order to break the bipolar and axial symmetry, this density profile was modified by random perturbations of the following kind: 1. With factor 2. With factor in the 10% probability the density was increased by a between 1 and 1.4. 5% probability the density was increased by a between 1 and 2, only if the cell was unmodified first process and the Z coordinate was positive.

The jet is injected in the middle of the grid in the region Z [-0.55, 0.55 kpc], R [0, 0.55 kpc], and [0, 2 ]. This region has the constant values: jet = 6.68â10-28 g/cm3 , vZ = ±0.4c, where c denotes the speed of light. The plus sign applies for the positive Z region and the minus sign for the negative one. The kinetic jet luminosity is Lkin = 1.04 â 1046 erg/s for both jets together. The pressure was set in order to match the external pressure at that position. This gives a slightly varying density contrast across the grid of = jet /ext 7 â 10-3 and an internal Mach number M = 10. The jacket of the injection cylinder is a further boundary. This boundary is open, so the material can leave the grid here into the central kiloparsec of the simulated radio galaxy, which is not attempted to model here. The temperature in the external medium is set to T = 3 â 107 K. The cooling time in the shocked cluster gas is approximately 100 Myr. The jet is expected to propagate through the whole volume in 10 Myr. So, cooling by bremsstrahlung marginally influences the state of the gas. This was taken into account (comp. (5) and sect. 2). Thus the calculation is not scalable anymore, formally. But given the smallness of the effect, scaling should be possible, in practice. Because the atmosphere is isothermal, the pressure varies in the same way as the density. In order to keep the system in hydrostatic equilibrium, gravity by an assumed dark matter distribution had to be applied. The gravitational potential, necessary to prevent the King atmosphere from exploding is:
DM

=

3 k T log 1 + (r/a) 2µmH

2

,

(7)

Figure 2. Bow shock shapes for the 3D run at t = 2.04 Myr. On the top panel the shape in one meridional plane is compared to the elliptical R = 7.5 1 - (z + 1)2 /169. The bottom panel compares the bow shock shape for different meridional planes to the parabola R = 2.3(26 - Z )1/3 (thick line). The bow shock is axisymmetric in the middle, where it can be well represented by an ellipse. For |Z | > 10, the bow shock is bumpy and not axisymmetric. That part is called cigar-like.

where µ is the number of particles per proton mass.

3.2. Early evolution & verification of b oundary conditions
A testrun was order to verify ution is shown of a bow shock performed on a ten times smaller grid, in the boundary conditions. This early evolin Fig. 4. The images show the formation and backflows, for both jets.


4

M. Krause: Very Light Jets II

Figure 4. The very early evolution of the 3D (left) and Mach number (right) are shown at of the Mach number. The bottom row shows Mach number, which includes the radial velo

simulation. The meridional slices at = 0, , joined at the axis, of number density 0.03 Myr (top) and 0.06 Myr (bottom). Values below 0.1 were omitted in the plot that sometimes small numerical artefacts at the axis are present, especially in the city.

In that phase, the jet almost ignores the stochastic nature of the density. The bow shock has a round and regular shape, and the density varies smoothly along its surface. However, the evolution on the two sides is different. The jet on the left-hand side, where the average background density is smaller, evolves faster: It produces a bigger bow shock, and a faster backflow at equal times. At t = 0.06 Myr, the two bow shocks are nearly joined together. As the later evolution shows, a single round bubble forms, soon after. The upper and lower halfs of the pictures fit quite good together, although there is sometimes a suspect spike at the apex of the bow shock on both sides. This numerical artefact reminds us about the imperfect treatment of the axis. The effect is rather small, and even absent for most of the simulation time. This supports the choice of boundary conditions on the axis (compare section 2.1).

3.3. Medium term evolution
A snapshot of several quantities at t = 0.32 Myr is shown in Fig. 6. The jet plasma takes the lowest density, and the highest temperature values. At that time, the evolution continued in the same way as in the early phase: the right-hand side, propagating into the on average denser medium, develops a broader cocoon, and is slower. The once separated bow shocks have united. The shape of the bow shock in that phase is oval, a sign of the blastwave phase. The bow shock shows two extensions in jet dir-

ection. Due to their later appearance, these parts of the bow shock will be called cigar-like. The aspect ratio of the bow shock is nearly 2, thanks to the extensions, which contribute approximately 0.5 to the aspect ratio at that time. The plots of the Mach number show a slightly supersonic backflow. The radial Mach number is positive for motion away from the axis, the toroidal Mach number changes sign at R = 0 for motion out of (into) the plane. The flow in R and Z directions is well ordered, whereas in the azimuthal direction, arbitrary fluctuations are seen. The jet beam shows both rotation (same colour above and below R = 0) and translation away from the axis (colour change at R = 0). The latter supports the chosen approach as being able to represent beam motions away from the axis of symmetry. The pressure plot shows a high pressure at shocks in the beam, especially on the axis. While this is in principal correct, the exact amount of the pressure could be influenced by the choice of the cylindrical coordinate system. The first shock in the beam is stronger on the left-hand side. This follows from the higher pressure there, but also from the higher inclination angle to the axis. Since the beams have identical conditions, why are the shocks not identical? The shocks in the beam are driven by whatever hits it. From the density plot, this seems to be the cocoon on the left-hand side, and the entrained cluster gas on the right-hand side. The backflows collide approximately at the center, forming a region of enhanced pressure.


M. Krause: Very Light Jets II

5

Figure 5. Slices of the tracer at t = 0.32 Myr for different axial positions. The shocked ambient gas is shown dark. Compressed regions are darker, rarefied and diluted regions appear lighter. Cells that contain less then roughly 10% ambient gas, i.e the beam plasma, are shown in white.

This region is asymmetric at the time shown here. Due to the stronger backflow from the right-hand side, the region has moved to the left. The higher pressure there drives a stronger shock into the beam. This explains why the first shocks are different on both sides. A passive tracer variable was advected with the flow, set initially to unity in the shocked ambient gas, to zero in the right beam, and to -1 in the left one. This tracer enables differentiation between the beam and the ambient gas. Slices at constant axial position of the tracer at t = 0.32 Myr are shown in Fig. 5. The right side has an approximately axisymmetric cocoon, while the left one is clearly not (compare slices at Z = -3.5 kpc; Z = -1.5 kpc). The slice at Z = -0.83 kpc demonstrates how the cocoons manage to merge together: The left one is smaller and slips inside the right one. This also brings shocked ambient gas inside the cocoon, which can still be spotted in the Z = 0.5 kpc slice. The stronger shock on the axis on the left-hand side causes the beam to widen, due to high pressure. On the right-hand side, less energy has been converted into heat, the beam stays narrow and delivers more power to the terminal shock, where the pressure is consequently higher. Fingers of shocked cluster gas, reaching down to the beam surface, are present just as they are in the unipolar 2.5D simulations (compare e.g. Krause & Camenzind 2001). They are generated in the following way: Kelvin-Helmholtz instabilities are excited at the boundary between cocoon and shocked external medium. The backflow advects those instabilities, while they are amplified. Hence, they are biggest in the symmetry plane. The material in the symmetry plane could in principal

flow outward, creating bumps in the bow shock, or inward. The simulation clearly shows no sign of outward motion. Instead, the gas is transported far down into the radio cocoon in geometrically thin fingers. Soon, their extension falls below the resolution limit of the simulation, and the gas mixes with cocoon plasma. However, in reality, the two phases may remain separate. The magnetic field of the radio plasma further supports the separation of the two phases. In principle, in that way new fuel could be channeled to the central source. However, the dynamics of the central kpc is beyond the scope of this work. The central high pressure region pushes the entrained shocked cluster gas to the right, where it slips between beam and cocoon, thereby widening the right cocoon. It has already been pointed out above that at the time the simulation is shown in Fig. 6, 0.32 Myr, the left jet converts more kinetic energy into heat than the jet on the right-hand side. At that time the tip of the jet has already advanced approximately 20 % further towards the left than towards the right. The difference of the average density is only 2.5 %, which cannot explain the fast propagation of the left jet, considering the estimate by the one dimensional force balance: vhead vbeam , where is the ratio of beam to head area. One important result is therefore that the nonlinear dynamics amplifies the effect of the different density on both sides by more than a factor of ten. The situation at that time is unstable, and the later timesteps show that the right jet catches up, and outruns the left one not later than at t = 0.95 Myr. The situation stays that way until the end of the simulation. On average, the right jet is approximately 10% faster than the left one, at late times. This is in conflict with naive intuition, but is


6

M. Krause: Very Light Jets II

Figure 6. Meridional slices for = 0/ at 0.32 Myr. The quantity plotted is indicated on top of the individual figures, units can be found next to the bar. The radial Mach number is defined to be positive away from the Z-axis. A positive toroidal Mach number is intended to mean motion into the plane for r > 0, and out of the plane for r < 0. The beam shows rotation in some places and translation in some others. Jet beam material can be decerned by its high entropy index.

readily explained, from the right jet to the left, where the beam slowing

considering that shifts the central stronger oblique the left jet down,

the stronger backflow pressure enhancement shocks are driven into as pointed out above.

3.4. Long term evolution
The final snapshot at t = 2.04 Myr is shown in Fig. 7. The plots show the usual picture of a hydrodynamic jet simulation, largely consistent with FR II radio galaxies: The cocoon is now nicely placed around the jet beam. The aspect ratio of the bow shock is 3.6. The Kelvin-Helmholtz instabilities show up prominently. They still grow towards the center, and develop into long fingers at the innermost

positions. The pressure shows a regular spacing of shock compression and rarefaction zones in the beam. High pressure regions are small and show up only at the end of the beams, where the Mach disk is located. The oblique shocks in the beam are now weaker than at t = 0.32 Myr. This follows from the smaller angle with the jet axis. The central region, with a diameter of roughly 10 kpc, is now dynamically calm. No large Mach numbers are observed there, and the pressure is approximately constant. The density takes intermediate values. This is now a relaxed region where jet plasma and shocked cluster gas are mixed (mixing may not happen in nature, see above).


M. Krause: Very Light Jets II

7

Figure 7. The same as Fig. 6, but for t = 2.04. (Only some of the variables are shown.)

3.4.1. The shap e of the b ow sho ck
Figure 2 shows the bow shock shape for the final snapshot in detail. It has an axisymmetric part in the middle, where it can be well represented by an ellipse. The center of this ellipse is not the origin of the grid, but is shifted to the left by one kpc. Such a shift of the center is also found in the larger 2.5D simulation. The elliptical shape ends at |Z | 10 kpc where two cigar-like extensions join the bow shock. These extensions are 3D in nature with several bumps. They can be represented, on average, by a parabola of rank three.

3.4.2. The law of motion of the b ow sho ck
In all the plots in the long term evolution, the two bow shock phases are easily discernable. The inner bow shock structure is a remnant from the blastwave phase. In its early phase, the radius of the bow shock should have obeyed equation (2), i.e. r t0.6 . (This part of the bow shock is far inside the core radius wherefore the density profile is roughly constant.) In the following, this law is checked for early and later evolution using the bow shock position in the Z = 0 plane. For every hundredth timestep, the positive radial bow shock position for = 0 and Z = 0


8

M. Krause: Very Light Jets II

was saved, up to the final timestep 220200, corresponding to t = 2.04 Myr. These values are plotted against time in Fig. 3. The positions were fitted with a function of the type: a + btc . The constant a was allowed for in order to take into account an artificial offset due to the initial conditions imposed by a jet with a final radius and nondeveloped structure at t = 0. The fit was carried out globally, and separately for the early and the late evolution. The result is given in Table 1. The global fit gives an exTable 1. Fit parameters for the bow shock position. The star denotes a fit with fixed a = 0. time range [Myr] global [0 : 0.5] [1 : 2] [1 : 2] a 1.02 0.41 1.25 0 b 4.44 4.52 4.21 5.45 c 0.58 0.38 0.60 0.49

This is approximately 40% of the true jet power (it would be 50% taking b from the global fit, and 130% for the late fit with a = 0). One can also use the bow shock position and velocity at e.g. t = 2.04 Myr in order to compute its power. From (2) one gets: L = 100 0 (r
b ow

)2 (v

b ow

/3)3 .

(11)

This yields a jet power of L 1046 erg/s, accurate to a factor of a few because of the slow bow shock velocity and the low spatial resolution. The true jet power is 1.04 â 1046 erg/s. It follows that the information on the jet power is conserved in the law of motion of the central part of the bow shock at least with an accuracy of a factor of a few.

3.4.4. Emission maps
The emission due to optically thin bremsstrahlung (e.g. Shu 1992a) was computed, mapped onto a Cartesian grid, and integrated for different viewing angles (Figs. 8 and 9). The general X-ray emission properties of shocked ambient gas in jet simulations have been discussed by Clarke et al. (1997), which has been updated recently by Zanni et al. (2003). The idea is that the gas is pushed aside by the jet cocoon. Depending on its compression, it will or will not form X-ray deficits at the location of the cocoon, and bright shells at the edges. Averaging the gas distribution and neglecting flows in the shocked ambient gas, the critical parameter is the relative shell thickness , defined as the width of the shocked ambient gas region divided by the local bow shock radius. Then, the ratio of observed flux to the flux of the initial condition is given by: f /f = F
-1

ponent c of 0.58, quite close to the exponent analytically derived for the blastwave phase (0.6). However, the behaviour is quite different at late times compared to the early evolution. At late times, the exponent is in the range 0.5 to 0.6 (compare below), whereas for the early evolution, c is close to 0.4, which would be the value expected for a supernova blastwave (initial energy input with no additional supply, compare (1)): r= 15E0 t 4 0
2 1/5

(8)

This is due to the initial conditions of the simulation: in order to get a propagating jet, one has to start with a short propagating jet. The energy of that initial jet is quite high. It can be estimated by calculating the energy stored in the 2 2 initial injection box: E0 Rj hj vj , where h is the height of the the cylindrical region. This amounts to approximately 1057 erg. The energy is also given by the parameter b, according to equation (8), and since the density is known: E0 = 2.35 â 1080 0 b5 10
58

(2 - )

-2

,

(12)

erg.

(9)

The difference of a factor of ten could be attributed to increased heat production in the very early evolution, where the unphysical initial condition, which is always present in such simulations, relaxes into a regular jet structure. At late times, the initial conditions may be regarded as relaxed. In that case, the fit parameter a could be set to zero. This would change the best fitting exponent from 0.6 to 0.5. An exponent lower than 0.6 would be expected due to the increase of the aspect ratio. However, the effect is quite small.

3.4.3. Implications on the jet p ower
From the b parameter of the late evolution, one can calculate the jet power, according to equation (2): L
kin

= 2.23 â 1067 0 b5 3.75 â 10

45

erg.

(10)

where F depends on the pre-shock and post-shock temperatures, and is close to unity. Therefore, X-ray deficits could be observed for a shell thickness above 38%, for sources located at 90 inclination. Here, the shell thickness is comparatively low for the whole simulation. Therefore, at high inclinations the X-ray surface brightness never falls below that of the undisturbed King atmosphere. This picture implies that the deficit is pronounced for low inclinations, since most of the gas is shifted aside. Indeed, we find these deficit regions in Fig. 9 for angles of 10 and lower. The deficit is much more evident at the later time. This is due to the different morphology: At t = 0.32 Myr (Fig. 8) the bow shock shape is essentially elliptical with only a small cigar-like extension. In that phase, the gas is compressed more isotropically (compare Fig. 6). The thickness of the shocked ambient gas shell is roughly the same in all directions. At the later time, the prominent cigar has displaced the gas sideways, producing the deficit when viewed pole on. Another reason is the decreasing density. The region in front of the jet head is furthermost from the center and the gas is therefore thinnest. The two phases of the bow shock, cigar and elliptical (comp. sect. 3.4.1), show up prominently in the emission


M. Krause: Very Light Jets II
0.9 0.8 0.7 erg/s/cm2 erg/s/cm2 0.6 0.5 0.4 0.3 0.2 0.1 -4 0.9 0.8 0.7 erg/s/cm2 erg/s/cm2 0.6 0.5 0.4 0.3 0.2 0.1 -4 0.8 0.7 0.6 erg/s/cm2 0.5 0.4 0.3 0.2 0.1 -4 0.8 0.7 0.6 erg/s/cm2 0.5 0.4 0.3 0.2 0.1 -4 0.9 0.8 0.7 erg/s/cm2 erg/s/cm2 0.6 0.5 0.4 0.3 0.2 0.1 -4 -3 -2 -1 0 X [kpc] 1 2 3 4 -3 -2 -1 0 X [kpc] 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 -6 -4 -2 0 Z [kpc] 2 4 6 1 2 3 4 erg/s/cm2 -3 -2 -1 0 X [kpc] 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 -6 -4 -2 0 Z [kpc] 2 4 6 1 2 3 4 erg/s/cm2 -3 -2 -1 0 X [kpc] 1.6 1.4 1.2 1 0.8 0.6 0.4 0.2 0 -4 -3 -2 -1 0 Z [kpc] 1 2 3 4 1 2 3 4 0.2 0 -4 -3 -2 -1 0 Z [kpc] 1 2 3 4 0.8 0.6 0.4 -3 -2 -1 0 X [kpc] 1.2 1 1 2 3 4 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 -4 -3 -2 -1 0 Z [kpc] 1 2 3 4

9

Figure 8. Bremsstrahlung emission maps for the 3D run at t = 0.32 Myr. From top to bottom, the viewing angle is 0 , 10 , 30 , 60 , and 90 . The left column shows the emission map, the middle one a vertical slice through the center, and the right one a horizontal slice through the center of the emission map. For comparison, the undisturbed cluster emission is given. The bow shock morphology is reflected in the ring like structures of differing sizes and surface brightness profiles.

maps. They form circular and elliptical rings, respectively, depending on the viewing angle. Where the rings partially overlap, they are brighter, producing the impression of ring segments (e.g. Fig. 8, 10 and 30 ., Fig. 9, 10 ). The structures can also intersect on the line of sight, producing bright spots (Fig. 9, 300 ). The pole-on figures show at least

two rings: one from the cigar phase and one from the inner elliptical bow shock part. Initially, the edges are brightened for any viewing angle. This changes at late times for high inclination. Here, the edges are prominently brightened only for the direction perpendicular to the jet axis. The reason is the declin-


10

M. Krause: Very Light Jets II
0.8 0.7 0.6 erg/s/cm2 erg/s/cm2 -8 0.6 0.55 0.5 0.45 erg/s/cm2 0.4 0.35 0.3 0.25 0.2 0.15 0.1 -8 0.45 0.4 0.35 0.3 0.25 0.2 0.15 -8 0.45 0.4 0.35 0.3 0.25 0.2 0.15 -8 0.45 0.4 0.35 0.3 0.25 0.2 0.15 -8 -6 -4 -2 0 X [kpc] 2 4 6 8 -6 -4 -2 0 X [kpc] 0.45 0.4 0.35 erg/s/cm2 erg/s/cm2 0.3 0.25 0.2 0.15 0.1 0.05 -20 -10 0 Z [kpc] 10 20 2 4 6 8 -6 -4 -2 0 X [kpc] 0.45 0.4 0.35 erg/s/cm2 erg/s/cm2 0.3 0.25 0.2 0.15 0.1 0.05 -25 -20 -15 -10 -5 0 Z [kpc] 5 10 15 20 25 2 4 6 8 -6 -4 -2 0 2 X [kpc] 4 6 8 erg/s/cm2 -6 -4 -2 0 X [kpc] 1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 -10 0.4 0.35 0.3 0.25 0.2 0.15 0.1 -15 -5 0 Z [kpc] 5 10 2 4 6 8 0.5 0.4 0.3 0.2 0.1 0 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 -8 -6 -4 -2 0 Z [kpc] 2 4 6 8

erg/s/cm2

erg/s/cm2

-10

-5

0 Z [kpc]

5

10

15

Figure 9. The same as Fig. 8 for t = 2.04 Myr. Some artefacts of the rotation procedure are visible at large Z values.

ing density profile that is superimposed on the jet features. Edges and rings are typically enhanced by a factor of a few in surface brightness. They should therefore be detectable for some sources in this phase.

4. 2.5D simulation 4.1. Simulation Setup
In order to study the jet evolution on larger scale, a 2.5D simulation was performed with similar initial conditions to the 3D simulation in the previous section. The simulation was run for 20 Myrs of simulation time; during that time the jet reached an extent of 110 kpc which corresponds

to 220 jet radii. In the radial direction, the jet reached 31 kpc. The jet radius was set to 0.5 kpc and the resolution was set to 20 cells per jet radius. This corresponds to [4400 â 1240] cells in total. With that resolution global parameters like the bow shock velocity on the Z-axis or energy and momentum conservation are accurate to 10% (Krause & Camenzind 2001). e,0 , a and were chosen slightly different to the previous model, because the sideways expansion of the bow shock should be able to escape the constant density part during the simulation time. The values are: e,0 = mp /cm-3 , a = 10kpc and = 0.75. The temperature was again set to T = 3 â 107 K. The jet is injected with a density of jet = 10-4 â e,0 , the sound


M. Krause: Very Light Jets II

11

4.2.1. The shap e of the b ow sho ck
The bow shock surface was extracted and fitted by elliptical functions as in the previous section. The result is shown in Fig.10. As in the 3D case, the bow shock can be represented by an ellipse in the center and by a parabola near the tip. This parabola now has rank two while it was a rank three parabola in the 3D case. The difference is due to the different jet nature. The beam in the 3D simulation stays intact all the way to the tip of the bow shock. Therefore it can deliver its momentum to a smaller area. In the 2.5D simulation the beam becomes turbulent, entraining cocoon and shocked ambient gas, delivering the momentum to a larger area. Also here, the center of the ellipse is shifted to the right, by 1.3 jet radii. This corresponds to the armlength asymmetry, i.e. the right jet is longer for t > 3 Myr, by a few percent. This confirms the findings of the 3D simulation.
Figure 10. Shape of the bow shock at 1 Myr (top) and at 20 Myr (bottom) for the 2.5D simulation. The fit functions are: R = 6.61 â 1 - (Z + 0.23)2 /53.28 (1 Myr), and R = 3.66 â Z + 55.6, R = 3.91 â 55.5 - Z , and R = 30.7 â 1 - (Z + 0.74)2 ./1277 (20 Myr). This means that the bow shock evolves from almost elliptical to elliptical + parabola extensions.

4.2.2. Beam stability
The beam in the 3D simulation was stable, whereas the beam in the 2.5D simulation becomes turbulent. The critical factor for this behaviour is not so much the dimensionality (that should work in the opposite direction, although the stability is probably somewhat enhanced by the choice of the cylindrical coordinates) but the density contrast. Lighter jets automatically moderate their Mach number (Krause 2003). This can be understood from the pressure equilibrium within the whole jet. Because in a very light jet, the sound speed is high and the expansion speed is low, the pressure is almost equal anywhere in the jet (compare Fig. 12). The jet beam adjusts to that pressure by varying the strength of its oblique shocks. At the actual location of the shocks, the pressure is higher. This sums up to a decline according to a power law with exponent -8 in the histogram and can be seen in the pressure distribution (Fig. 12). The sound speed in the beam is therefore above 10% the speed of light, which makes it hard for the jet to be highly supersonic in a non-relativistic simulation (the typical Mach number in the beam is between one and two). The problem is even more severe at earlier times because the average pressure decreases with time (see below). It is well-known that jet beams at low Mach number are disrupted quickly (e.g. Bodo et al. 1994, 1995, 1998).

speed in the jet was set to 20%c, and the jet's Mach number to M = 3. The high internal sound speed and the low Mach number was chosen because a parameter study (Krause 2003) has shown that very light jets adjust their pressure quickly to the cocoon pressure by oblique shocks (no expansion or contraction) and that high Mach numbers cannot be sustained in a non-relativistic very light jet. No cooling was taken into account, and hence the simulation is scalable as outlined in Zanni et al. (2003). With these parameters, the simulation took about a month on an SX-5 supercomputer. It was therefore not possible to do this in 3D, for now.

4.2. Results
We present logarithmic density plots of the simulation results for four different simulation times (5, 10, 15, 20 Myr) in Fig. 11. The morphology that appears in these figures is a continuation of previous simulations that could not reach the size shown here. The 5 Myr figure shows the state that was reached by Krause (2003), and extensively discussed therein. In this early phase, the bow shock is spherical, its radius following an expansion law given by the force balance equation (1). At later times the cocoon transforms via a conical phase towards a cylindrical one. This is a remarkable result because classical double radio galaxies also often have cylindrical cocoons.

4.2.3. Pressure evolution
The average jet pressure is the driving force of the inner elliptically shaped part of the bow shock. Figure 12 shows that the pressure in the jet system monotonically decreases with radius. Close to the axis, the pressure is higher because of shocks in the beam region. In the shocked ambient gas region, a new equilibrium of gravity and pressure appears. The smallest pressure values are located at the bow shock, roughly 20% below the average. In the previous section, it was shown that the inner part of the bow shock roughly follows the spherical


12

M. Krause: Very Light Jets II

Figure 11. Four snapshots of the 2.5D simulation. The logarithm of the number density is shown. The times for the snapshots are indicated on top of the individual figures. The same part of the grid is shown in each case. The jet forms first a spherical bow shock, associated with a spherical cocoon. The cocoon then transforms via a conical state to a cylindrical one. The bow shock's aspect ratio (length/width) also grows, up to 1.8.

Figure 12. Left: Pressure distribution for the 2.5D simulation after 20 Myr. Top right: Axial pressure average over the radius. Bottom right: Pressure histogram within the jet (i.e. beam, cocoon, and shocked ambient gas). It can be represented by a broken power law with humps around the most frequent value. The power law indices are 5 and -8, respectively. 90% of the volume has a pressure in the range [1 - 2] â 10-9 dyn/cm2 . The noise represents the statistical error.

propagation law (1) inside the core radius. The accuracy of this law will be checked in the following also for the larger 2.5D simulation. In this case, the bow shock has propagated more than three core radii in the sideways direction. Notice that the formulae derived for blastwaves with strong shocks are still applicable here, as shown in the appendix. In the spherical approximation, the average pressure inside the whole jet is given by 1 : pj = ( - 1)
1

Lt - Mv 2 /2 , Vj

(13)

Here, = 5/3 is the adiabatic index, and Vj is the jet volume (including beam, cocoon and shocked ambient gas). The power L includes all sources of energy, i.e. the flux of kinetic and internal energy through the jet channel, the flux of internal energy entering through the surface of the bow shock, and the energy lost by work against the gravitational field. The simulation takes these effects properly into account. In order to test the pressure evolution against deviations from the simple spherical law these power sources are calculated as follows: Via the jet beam, a constant power is injected, given 2 3 by: Lj = Lj,kin + Lj,int = Rj vj (1 + ( ( - 1)M 2 /2)-1 ) = 45 9 â 10 erg/s.

The formula neglects the increasing but small part of kinetic energy stored in the beam.


M. Krause: Very Light Jets II
1.5

13

1 power [1046 erg/s]

0.5

0

-0.5

-1 0 2 4 6 8 10 12 time [Myr]

Lb,int Lgrav 0-5 Myr fit 10-20 Myr fit 14 16 18 20

Figure 13. Approximate internal energy flux through the bow shock (Lb,int ) and energy loss rate due to gravitational uplifting (Lgrav ) over simulation time. Lint,e = 4 r2 pe v /( - 1) was calculated using the bow shock's radius and sideways velocity from the simulation. The fitted functions are: 0.57 t0.37 (0-5 Myr), and 1.51 t-0.08 (10-20 Myr). The gravitational energy loss rate levels off at 30% of the beam energy flux.

10 dyn/cm , %

10

ity (p, r ), and the exponent of a local power law approximation to the sideways bow shock propagation (r t , = 3/( + 5)). But for a King profile, it takes some time until adjusts to the value it should have due to the local power law approximation. The latency is caused by the memory of the flow of its history due to the fixed amount of energy and mass at a given time. Lb,int is shown in Fig 13 together with fits at early and late times. It should first rise proportional to t4/5 . The fit gives an exponent of 0.37 which again reflects the effect of the initial conditions. The fitted exponent for Lb,int of -0.08 after 15 Myr is in agreement with a local power law exponent for the pressure of = -1.84, where = 0.8 has been adopted from the measured bow shock propagation for that time span. The local in the fitted region is in the range -1.7 to -2.0. The rate of change of gravitational energy, Lgrav , is computed from the simulation data, averaging over one million years and is plotted in Fig 13. Lgrav rises from one to thirty percent of the jet beam's energy flux, where it seems to converge towards the end of the simulation. The gravitational energy loss rate and the power entering as internal energy through the bow shock make up a similar contribution to the energy within the jet system as the energy flux through the beam. Using L = Lj + Lb,int - Lgrav , where Lb,int and Lgrav now denote the time-averaged value at a given time, (13) can be evaluated, where v and t are given by (1): v= Lt2 Mr 3 t= L (14)
r 0 1/3

-9

2

1

0-5 Myr fit 15-20 Myr fit relative deviation average pressure (spherical approximation) average pressure (simulation) 1 time [mio yrs] 10

M(r )dr

,

(15)

Figure 14. Average jet pressure over simulation time. The stars show the values measured in the simulation, crosses mark the pressure according to a spherical approximation, plusses show the relative difference between the former (in %). Corresponding symbols are connected with solid lines. The lines show fits to the pressure measured from the simulation. The fits are: 32.84 t-0.69 (0-5 Myr, solid line), and 227.95 t-1.66 (15-20 Myr, dashed line). The best fit for the spherical approximation in the range 15-20 Myr is: 126.05 t-1.50 (not shown).

and r denotes the sideways extent of the bow shock. Figure 14 shows this analytical estimate together with the data from the simulation. Here, r was related to time via measurement from the simulation. The agreement is quite good, in general. The analytical formula follows the slope of the simulation data, but underestimates it by up to 20%. Towards the end of the simulation, the pressure evolution can be approximated by a power law with exponent -1.66. For the same time range, the spherical approximation is well approximated by a power law with exponent -1.50. Note that these are only local approximations to curved functions that have not yet reached the asymptotic power law regime.

Through the jet's bow shock, the power Lb,int = 2 4 rb vb pKing (rb )/( - 1) enters the jet system in the spherical approximation. rb and vb are the bow shock's position and velocity, which is measured at the position of the greatest cylindrical radius. pKing denotes the initial pressure distribution of the isothermal King profile. In the spherical approximation, Lb,int should be proportional to t(+3)-1 , where would be the exponent of an isothermal power law distribution in pressure and dens-

4.2.4. Sideways motion of the inner b ow sho ck part
From the pressure evolution (previous section) one should expect that the sideways motion of the inner bow shock part roughly follows the spherical approximation for the whole simulation time. The density profile used here has the asymptotic power law approximations: limr0 () = 0 , and limr () r-9/4 . Therefore, the bow shock should expand with r t0.6 at the beginning, steepening


14
0 - 5 Myr fit 15 - 20 Myr fit 15 - 20 Myr fit (spherical approximation) spherical approximation measured radial bow shock position

M. Krause: Very Light Jets II
measured axial bow shock extention 0 - 3 Myr fit 17 - 20 Myr fit

100

10

Zmax-Zmin [kpc] 10 1 t [Myr] 10

r [kpc]

1 t [Myr]

10

Figure 15. Bow shock radius at Z=0 versus time (squares) with fits and compared to spherical approximation, including all power sources (see text, plus-signs). The three fits are: 4.66 t0.53 (0-5 Myr), 2.99 t0.78 (15-20 Myr), 2.36 t0.81 (1520 Myr, spherical approximation).

Figure 16. Bow shock extention on the axis versus time. The two fits are: 10.45 t0.49 (0-3 Myr), and 3.49 t1.15 (17-20 Myr).

towards r t1.09 , at least as long as it remains spherical. The radial bow shock position was determined every 0.4 Myr (Fig. 15). The resulting function has a curved nature. Following the arguments of sect. 3.4.2, the bow shock propagation was locally fitted by a function of type a + b tc . The resulting parameters for the different regions are shown in Table 2. Usually, a = 0, since only the very late evolution of the jet is studied. Only for the time span up to 5 Myr a fit with a = 0 has been included because here it is possible that effects from the initial condition still dominate the propagation. For comparison, also fits
Table 2. Fit parameters for the bow shock position in the 2.5D simulation. The star denotes a fit with fixed a = 0. Two stars denote fits to the spherical approximation with fixed a = 0. time range [Myr] [0 : 5] [0 : 5] [0 : 5] [15 : 20] [15 : 20] a 2.00 0 0 0 0 b 2.67 4.66 3.22 2.99 2.36 c 0.76 0.35 0.67 0.78 0.81

From the increasing aspect ratio, an exponent lower than the one of the spherical approximation should be expected. The simulation shows that the effect is small.

4.2.5. Axial b ow sho ck propagation
According to self-similar hydrodynamic jet models (compare e.g. Alexander 2002, and references therein), the velocity of the jet head should be proportional to vbow,ax = t(+2)/(+5) , where is the exponent of a local power law for the external density ( r ). For a jet head with constant area, the velocity should be given by the one dimensional estimate: vbow,ax = vbeam t2/(+2) . The last proportionality holds for > -2, only. The bow shock propagation in axial direction is shown in Fig 16. Its evolution can be represented by two power laws that are much closer to the expectation from the self-similar models (compare Table 3) than to the onedimensional estimate. For the time range up to three milTable 3. axial bow pected in denotes Comparison between the fitted exponent e for the shock velocity (vbow,ax te ) and the exponent exthe self-similar and one-dimensional (1D) estimate. denotes velocity growth faster than any power law. measured -0.51 0.15 self-similar -0.4 0.09 1D 0

to the detailed spherical approximation are given, computed by application of (15). According to that, the exponent for the first five million years should be 0.67. Using the pure power law, an exponent of 0.35 is achieved in the simulation data. Allowing for the radial offset gives a best fit exponent of 0.76. Since the exponent of 0.35 is much below any expectation, it follows that the initial condition is still important in that phase, and the fit with offset is more appropriate. The concurrence of the curves increases with time and for the last five million years, the exponent for the power law fit of the simulation data (0.78) differs from that of the spherical approximation by only 0.03.

time range [Myr] [0:3] [17:20]

lion years, the bow shock is still spherical. The exponent should therfore be exactly -0.4. But the initial conditions set up the jet with a finite amount of energy from the beginning which would cause an exponent of -0.6, analogous to the supernova case. The measured exponent of -0.51 shows that the jet is just in transition between these two phases. Towards the end of the simulation the jet head velocity grows slightly faster than self-similar. The good


M. Krause: Very Light Jets II

15

agreement with the self-similar models is probably due to the unstable beam, at least in the later phases, where the momentum is distributed over a large and increasing area. Jet models with intact beams (at higher jet density) have been shown to have head advance speeds more consistent with the one-dimensional estimate than with the self-similar models (Carvalho & O'Dea 2002a).

4.2.6. Asp ect ratio
Defining the aspect ratio to be the length divided by the width of the bow shock, it grows proportionally to t0.34 towards the end of the simulation. This is clearly a non-selfsimilar feature, and may be partially due to the squeezing of the cocoon by the gravity of the swept up gas, as discussed recently by Alexander (2002).

4.2.7. Evolution of the contact discontinuity.
The temporal evolution of the contact discontinuity is shown in Fig. 17. Both the maximum and the average radius are shown. The position of the contact discontinuity was determined by a passive tracer variable that is advected with the flow. Due to mixing, the cocoon plasma does not keep its initial value of zero. The contact surface is therfore defined to be at the largest radius with the tracer having less than 10% of the ambient medium value. The average and maximum values of the contact discontinuity evolve in a dissimilar way, which is due to the geometry changes discussed above. Other than at the earliest times, the contact discontinuity is always decelerating. For some time (up to t = 1 Myr) this is sufficient to stabilise against the RayleighTaylor instability. To be Rayleigh-Taylor stable, the deceleration of the contact discontinuity must overcome the gravity of the assumed dark matter halo, which is given by gD = 3 k T r/a = 1.94 10 µmH a 1 + r2 /a2
-6

value, which is the case for a fraction of the first Myr only. Another complication is that the contact surface is not smooth but Kelvin-Helmholtz fingers always penetrate through. The conclusion is that, at least in the central region, swept up gas is constantly entrained into the cocoon. The entrainment rate is shown in Fig. 19. It is slightly rising ( t0.32 ). The entrained mass is compared to the mass of the gas that filled the cocoon volume before the jet activity. This mass fraction is linearly rising, and would reach unity after 1.6 Gyr, if it would continue in the same way. The width of the cocoon relative to that of the bow shock radius () is also given in Fig. 17. For self-similar evolution this value should be constant, and detailed selfsimilar models (Heinz et al. 1998; Kaiser & Alexander 1999) place it in the range of 0.8 to 0.9. Such a high value is reached at no time. The relative width decreases with time to 0.4 (0.26) in the case of the maximum (average) value. This result is in agreement with simulations by Zanni et al. (2003). For the late phase, the average value for levels off. This may indicate a phase of nearly self-similar behaviour, only disturbed by the gravity of the swept up gas, as discussed by Alexander (2002). The slowly growing aspect ratio also supports this view. The small value of may be due to the weak bow shock, as discussed in more detail by Zanni et al. (2003).

5. Comparison to observations 5.1. Cygnus A
This section refers to published data on Cygnus A (Carilli & Barthel 1996; Smith et al. 2002). For illustration, an overlay of the smoothed X-ray and radio data is shown in Fig. 20 below the synthetic X-ray image of the 2.5D simulation. The X-ray-radio comparison shows the jet-IGM interaction, impressively. Coincident with the radio cocoon, there is a deficit of X-ray emission (disregarding for the moment the bright filaments in the center and the emission from the beam and the core). As discussed above, and detailed in Clarke et al. (1997) and Zanni et al. (2003), the parameter , i.e. the thickness of the shocked ambient gas region over the bow shock radius, should exceed 0.38. Although the 2.5D simulation presented here differs from Cygnus A because of the instable beam, it should be appropriate for the overall morphology. In the simulation, the shocked ambient gas has elliptically shaped X-ray isophotes. Smith et al. (2002) report elliptically shaped isophotes from the cocoon boundary out to 66/(h/0.7) kpc (i.e. the outer boundary of their annulus four). Further out, the spherical fit is the better one. Identifying the bow shock's radius with that position, together with a cocoon width of 32/(h/0.7) kpc, results in = 0.76, in agreement with the discussed limit. Taking the average value for the position of the contact discontinuity, the 2.5D simulation produces = 0.74 at t = 20 Myr, showing that

r/a cm/s 1 + r2 /a2

2

(16)

The central region is represented by the maximum value, which shows the strongest deceleration at the beginning. We determine a characteristic deceleration from a fit to the maximum value of the contact discontinuity up to 1 Myr. Note that this is also an upper limit for the whole simulation. The result is compared to gravity in Fig. 18. This shows that the contact surface should be Rayleigh-Taylor unstable for t > 2 Myr, which is consistent with the density plots. The linear growth time l/g = 2.5 Myr (l/25 pc)(10-6 cm s-2 /g ), where ( l is the wavelength of the instability and g the total acceleration) is typically below the simulation time for wavelengths greater than the resolution limit. The details of the acceleration of the contact discontinuity are more complex than in the global discussion of the previous paragraph. Figure 18 shows that it is constantly shaken with values of the order 10-5 cm s-2 . For stability, the global deceleration has to exceed this


16

M. Krause: Very Light Jets II

10 r [kpc]

relative radius

average maximum 0-1 Myr fit

0.8 0.7 0.6 0.5 0.4 0.3 0.2

average maximum

1 0.1 1 t [Myr] 10

0

2

4

6

8 10 12 14 16 18 20 t [Myr]

Figure 17. Position of the contact discontinuity over time. The contact discontinuity is determined by a passive tracer. The maximum and the average value is shown. The left figure shows absolute values, the right one shows the value relative to the maximum sideways bow shock radius. The fit in the left figure is a power law with exponent 0.28.

0.0001 1e-05

deceleration (fit) gravity
2

2 1.5 1 10 cm/s 0.5 0 -0.5 -1 -1.5 -2 average maximum 0 2 4 6 8 10 12 14 16 18 20 t [Myr]

cm/s

2

1e-07 1e-08 0.1 1 t [Myr] 10

Figure 18. Left: Comparison of the deceleration (determined by a global fit for the time span up to one Myr) of the contact surface against the local gravity (at the average position of the contact surface) of the dark matter halo. Right: Detailed acceleration averaged over 0.2 Myr (because of the spatial resolution) for both, the maximum and the average position of the contact surface.

such a high value is indeed possible for sources of that size. Using the knowledge of the sideways bow shock position, its aspect ratio (length to width) may be derived, which turns out to be 1.1 (using the hot spot separation of 147/(h/0.7) kpc from Carilli & Barthel (1996), it is 1.1 for the eastern and 1.2 for the western jet). From the low aspect ratio, Krause (2003) has concluded a density ratio of < 10-3 . However, also the = 10-4 2.5D simulation presented here reaches an aspect ratio of 1.8 at the end of the simulation. Since aspect ratios can only grow (Krause 2003), and the simulated jets are smaller than the observed jet in Cygnus A, it follows that the real jet still encounters more mass than the simulated ones. This is certainly due to a shallower density profile in the atmosphere of Cygnus A of = 0.51 (Smith et al. 2002) compared to = 0.75 employed for the simulations (the latter was the old value from the ROSAT data). At the largest extent of the bow shock in the 2.5D simulation

-5

1e-06

this would increase the ambient density by a factor of three, which already might result in the desired reduction of the aspect ratio. A still lower central value for the density contrast cannot be precluded. However, a density contrast of 10-4 would also be supported by the large cocoon width (Rosen et al. 1999). Another constraint for the density contrast was presented by Krause (2003), who showed that > 3 â 10-6 , because the jet is no longer in the spherical bow shock phase. It has been shown above that the sideways, pressuredriven part of the bow shock satisfies (1) with good accuracy even at an aspect ratio of 1.8. Given the low aspect ratio of Cygnus A, this equation should be very accurate here, and can be used to derive the jet's power and age. In order to do that, one needs information on the gas distribution before the jet activity. One requirement is that the distribution should join the density distribution far from the radio lobes, smoothly. Assuming it was a King


M. Krause: Very Light Jets II

17

35 30 25
sun

entrained mass 2-20 Myr fit

4.5 4 3.5 % 3 2.5 2 1.5

mass fraction 2-20 Myr fit

10 M

20 15 10 5 0 0 2 4 6 8 10 12 14 16 18 20 t [Myr]

8

0

2

4

6

8 10 12 14 16 18 20 t [Myr]

Figure 19. Entrained gas mass in the cocoon over time. Left: Absolute values, the fit function is: 0.677(t/Myr) 1.32 . Right: Percentage with respect to the mass of the initial condition in the same volume, the fit function is: 1.5 + 0.06(t/Myr).

distribution (6), as found for other clusters of galaxies, the following equation has to be satisfied:
e,0

weak shock. The shock conditions for a weak shock yield (e.g. Shu 1992b): v
b ow

= 1.64 (a/10 kpc)

-1.53

mp /cm

3

(17)

=

The core radius a (for the density profile before the jet activity) cannot have been much greater than 20/(h/0.7) kpc. This follows from the enhanced X-ray emission next to the radio cocoon. Such enhancements have been shown analytically to appear in atmospheres steeper than r -1 (Alexander 2002, and references therein). In the 2.5D simulation, the line-of-sight integrated X-ray emission is strongest next to the radio lobes for positions of the contact discontinuity in excess of 60% of the core radius. There is no obvious lower limit to the core radius. However, it will be shown in the following that the exact value of a is not needed for the present discussion. Approximating the gas density e,0 with a constant, one infers the following conditions for the jet's power and age: L= t= 100 27
2 3 e,0 rb ow vb ow



2

The measurement indicates a temperature jump T1 /T0 of 1.7, and errors limit it to smaller than 4. Also, the bow shock velocity should exceed the speed of sound in the preshock gas: cs = 1313 T0 /5.4 keV km/s. This results in a sideways bow shock velocity of: v
b ow

2 kB T0 T0 = 1313 mp 5.4 keV 2 8 T1 T1 7 7 = + -+ - 5 T0 8 T0 8

km s 15 . 64

= 2217 ±

1948 904

km/s .

Inserting this into (18) and (19) and assuming an average pre-jet density of e,0 = 0.05mp /cm3 yields L = 4.4 ±255 3. â1047 erg/s and t = 17 8 Myr. 12 The parameters can also be estimated assuming a King distribution in the pre-jet era. It follows from (1): L= I (rbow /a)3 9a J (rbow /a)2 3a2 J (rbow /a) t= rbow vbow I (rbow /a) 4
3 3 e,0 rb ow vb ow y

(18) (19)

(20) (21)

3rbow . 5vbow

An interesting feature of these equations is that the source age is independent of the external density, which remains true for the normalisation of arbitrary density distributions. Applied to the simulation data, this method yields quite accurate values, as demonstrated in sect. 3.4.3. In order to apply it to Cygnus A, the bow shock's sideways radius is taken to be 66 kpc (see above), and the velocity from the temperature jump measured from the X-ray data at that position. According to Smith et al. (2002), the temperature of the preshock gas is at about 7-8 keV at large radii, falls to 5.4 keV immediately before and then rises significantly to 9.2 keV at the position where the bow shock was proposed to be situated above. It is therefore a

I (y ) =
0 z

x2 (1 + x2 ) I (y )y dy ,

0.765

dx

J (z ) =
0

where e,0 is given by (17). The result is shown graphically in Fig. 21. It turns out that all the parameters are only weakly dependent on the core radius. The average value for each curve is given in the following summary:
. L = 7.9 ±4433 â1047 erg/s 6. t = 24 11 Myr 16

E = 5.9 ±

14.9 3.8

1062 erg ,


18

M. Krause: Very Light Jets II

Figure 20. Top: Line-of-sight integrated X-ray emissivity for the 2.5D simulation. Bottom: smoothed X-ray map of Cygnus A (Chandra archive, credits: NASA / UMD / A.Wilson et al.) with 6cm radio contours (VLA, credits: NRAO/AUI, Chris Carilli and Rick Perley, (compare Carilli & Barthel 1996)) overlayed. h denotes the Hubble constant in units of 100 km/s/Mpc.

where E is the total energy released by the jets. As discussed above, in order to obtain the power of the beam, the flux of internal energy through the bow shock has to be subtracted and the gravitational energy increase has to be added. A reasonable estimate of the internal energy contained within the bow shock region before the jet activity can be obtained by assuming an isothermal King atmo-

sphere. Adopting a core radius of 10 kpc, the X-ray data far from the center is consistent with a central pressure of p0 = 10-9 erg/cm3 . This leads to a total of 1060 erg, which is negligible here. The gravitational energy increase is not so easy to estimate. But for the present discussion it is sufficient to say that the 2.5D simulation indicates that up to 30% of the jet power is used to lift up gas. Adopting


M. Krause: Very Light Jets II

19
B ,HS

60 55 50 45 40 35 30 25 20 15 10 0 80 70

vbow= 1313 km/s vbow= 2217 km/s vbow= 4165 km/s

= 5 â 1044 (Rj /0.57 kpc)2 2 (uB /u

).

2

4

6

8 10 12 14 16 18 20 a [kpc]

Here, is the bulk Lorentz factor = vj /c, and uB ,HS = 9â10-10 erg/cm3 is the magnetic energy density measured in the hot spots by Wilson et al. (2000). It is again an upper limit, since the magnetic field in the jet will be lower than in the hot spot, where it is shock compressed. It is clear that the derived power can also not be raised by the magnetic energy flux, as long as we assume a nonrelativistic jet. Therefore, the jets in Cygnus A have to be relativistic. The kinetic power for a relativistic jet may be written as follows. L
kin,R

age [Myr]

erg/s]

60 50 40 30 20 10 0 0

vbow= 1313 km/s vbow= 2217 km/s vbow= 4165 km/s



= 8.8 â 1046 (Rj /0.55 kpc)2 (104 R,0 ) (e,0 /mp cm-3 )(1 - (h)-1 ) erg/s
R

2 = 2 Rj R e (1 - (h)

-1

) c

3

= j h2 /e ,

2

4

6

8 10 12 14 16 18 20 a [kpc]

Figure 21. Detailed results for age (top), and power (bottom) for Cygnus A's jets, assuming h = 0.7 and the sideways bow shock velocity indicated on the figures.

where h = 1 + e/c2 is the specific enthalpy. R denotes the relativistic generalisation of the density contrast. This number should be similar to the value for the nonrelativistic (compare Rosen et al. 1999). The central cluster density e,0 has a weak dependency on the core radius which is subsumed in a variation of R,0 , which is regarded to have a value below 10-4 . This can be combined with constraints on the relativistic Mach number M = /s s , where the index s denotes the values for the sound speed (e.g. Marti et al. 1997). The sound speed is given by c2 = p/h = ( - 1)(1 - 1/h)c2 , which yields s for the specific enthalpy: h= M 2 + 2 2 -2 , g= . 2 + g 2 2 M -1 (22)

power [10

47

a value of 20%, the power through the jet channel is determined to L > 1.9 â 1047 erg/s. Neglecting this rather uncertain contribution, the lower limit to Cygnus A's jet power is 1.6 â 1047 erg/s. This has to be compared to the power in the jet channel. For a non-relativistic jet: L
kin 2 = R j j v 3 j 45

= 5 â 10 (Rj /0.57 kpc)2 (104 0 )( (vj /(c/2))3 erg/s .

e,0

/mp cm

-3

)

The jet's radius was adopted from Carilli & Barthel (1996) using the latest value for the Hubble constant, h = 0.72 (Spergel et al. 2003), and the subscript zero denotes values in the center. Hence, even when using most optimistic numbers, the kinetic jet power falls short of the lower limit derived above, by a factor of ten or more. This means that the internal energy cannot be responsible for the rest, since this would yield a subsonic jet. Therefore the energy should be in the (toroidal) magnetic field: L
B 2 = 2 Rj uB 2 v j

Since Cygnus A shows a stable and laminar jet, the Mach number has to exceed unity. Most likely the Mach number is greater than that. Carilli & Barthel (1996) infer a Mach number of M 8 (from the oblique shocks) and M < 13 (from the opening angle). The total power derived in this way is plotted against the Lorentz factor in Fig. 22 for different values of the central density contrast and Mach number. The lower power limit derived above results in a Lorentz factor around 10, the central value gives 37 to 39. An important result is that the jet's Lorentz factor should exceed ten for reasonable assumptions. According to (22) M should then exceed 14, in order to get h > 1, which might point to a significant Alfv´ speed. en If one would regard a Lorentz factor above twenty as unreasonable, the jet power could be further constrained to less than 3 â 1047 erg/s. The derived jet power exceeds the total emitted radio power by at least a factor of 100 (Carilli & Barthel 1996), and might be compatible with estimates by Kaiser & Alexander (1999) (2â1046 erg/s) and Zanni et al. (2003) (3 â 1046 erg/s), who do not state their errors. The simulations have shown that the KelvinHelmholtz instability produces large fingers in the central regions of the cocoon. These may be identified with the


20
5 1.6 M= 1, M= 8, M=13, M= 1, M= 8, M=13,

M. Krause: Very Light Jets II

47

erg/s]

4

kin,R

[10

3

R.0 R.0 R.0 R,0 R,0 R,0=10

=10-4 -4 =10-4 =10-5 =10-5 =10 -5

2

1

0 2 4 6 8 10 12 14 16 18

Figure 22. Power through the jet channel of Cygnus A over the Lorentz factor compared to the lower limit derived in the text. Curves for differing Mach numbers are almost identical.

belts of cooler (4 keV) gas inside the cocoon of Cygnus A. Additionally, the contact discontinuity might be RayleighTaylor unstable (Alexander 2002). The 2.5D-simulation shows that the deceleration is stronger at the beginning, and falls below the critical value for most of the simulation time. Also, local shaking dominates over the global deceleration. Therefore the contact discontinuity is expected to be Rayleigh-Taylor unstable.

5.2. Hercules A
Radio and ROSAT data on Hercules A are given by Gizani & Leahy (2003, 2004). This source is an intermediate one, showing huge radio lobes like FR II sources, but no edge brightening. It seems to be in a similar state to the 2.5D simulation towards the end, presented above. The jet is unstable, probably because it cannot reach the required Mach number in order to stabilise it up to the tip of the cocoon. The whole cocoon seems to be turbulent, and entraining external gas across the contact discontinuity. This might be the reason for the lack of an X-ray deficit as observed in Cygnus A. However, the ROSAT data may lack the necessary resolution. What it does show is that the X-ray gas is elongated in the direction of the radio source. This can be interpreted as a wide region of shocked ambient gas, as in Cygnus A.

suggests that the situation may be similar to the ten degree plot of Fig. 9: A low density jet may have blown an approximately spherical bubble. After reaching the critical radius (see above), a jet with a narrow beam, i.e. high thrust bored a cigar-like extension into the shell. Now, we view the cigar-like extension from such an angle that it partly overlaps with the big shell, producing enhanced emission at those parts of the X-ray rings. An interesting point in this interpretation is that when modeling the X-ray enhanced regions by two elliptical, overlapping rings, these are elongated in the same direction, and depro jecting to a circular ring results in an inclination of roughly 37 for both. This inclination is compatible with the estimate of 45 from the radio data by Blanton et al. (2003). A question concerning this interpretation arises from the temperature map of the X-ray shells. The smaller shell seems to be colder than the larger one. While the uniform temperature of each of the elliptical shells is consistent with their proposed nature, the cigar part, corresponding to the smaller elliptical ring, is always hotter in the simulations than the inner bow shock part. This problem might be alleviated by the fact that the jet in this source is presently not active. Therefore, the radio plasma that was present in the cigar may have left it in the backflow, leaving behind an underpressured region. The gas in the shell would then expand somewhat into the cigar cavity, thereby lowering its pressure. For that reason, the expansion velocity of the leading bow shock would slow down. Assuming a cylindrical shell for the cigar part, the adiabatic expansion (T V -8/3 ) would require a change of the inner radius of the shell from e.g. 80% to 70% of the outer one, in order to explain a temperature change by a factor of two. Following the smaller X-ray shell that has been interpreted as cigar part above, there are emission line filaments (Blanton et al. 2001). These might be due to weak radiative shock waves inside of the shell that might have been excited by the re-expansion described above.

LB +L

5.4. Radio galaxies at higher redshift
Radio galaxies with redshifts greater than z 0.6 show associated large scale emission line regions that are aligned with the radio jets (McCarthy 1993). The aligned emission line regions are often located where one would expect the radio cocoon, and are stronger towards the center (e.g. Best et al. 1996). It has been shown in this paper that the central cocoon regions can entrain several percent of the gas mass that was situated there before the jet event. The X-ray belts in Cygnus A have been interpreted above as evidence for this entrainment process. These belts are the coldest X-ray gas found in the Cygnus A system. It is therefore possible that such entrained gas reaches thermal instability, if the cluster's gas density is higher, or the source has some more time to cool. The cooling time for these belts is some 100 Myr. One can therefore imagine a

5.3. Ab ell 2052
X-ray, radio, emission line, and other data on this central cluster galaxy can be found in Blanton et al. (2001, 2003). It is an old radio galaxy, the radio emission is diffuse, and no sign of hot spots or beams has been found. The diffuse radio emission is devoid of X-rays, and is surrounded by bright X-ray shells. Blanton et al. (2001, 2003) argue that the shells consist of two oppositely situated bubbles, as found in other X-ray clusters. However, the study of the dependency of the X-ray emission on the viewing angle


M. Krause: Very Light Jets II

21

situation where the entrained gas cools down to emission line temperatures whereas the other shocked ambient gas does not cool. As suspected above for Abell 2052, the cooling could also proceed within X-ray filaments that contain weak, radiative shocks. This would explain why the line ratios for some sources (preferentially the smaller) are consistent with shock excitation (e.g. Inskip et al. 2002), without relying on the compression of pre-existing clouds by the bow shock. Wether this picture can cope with other observational data has yet to be explored.

6. Summary and conclusions
Bipolar jet simulations have been presented in 2.5D and 3D. For the 3D simulation, a cylindrical grid was employed, and the boundary conditions were given. It has been shown that the influence of the axial boundary remains acceptably small. Because of small disturbances in the ambient medium, the backflows are located at different distances from the jet axis, and permeate each other in the center. There, a turbulent mixing region emerges that entrains shocked ambient gas across the contact discontinuity via Kelvin-Helmholtz and Rayleigh-Taylor instabilities. After a certain time, both simulations show two bow shock phases. The shape of the outer one is well fitted by a parabola. The inner part has an elliptical shape and stays axi-symmetric, even in the 3D-simulation. It follows the blastwave equation of motion with good accuracy (a few percent difference in the exponent of a local power law fit), even when the aspect ratio is high. The aspect ratio keeps growing for the whole simulation time, i.e. a self-similar regime is not reached. The viewing angle-dependent emission maps show that the different parts of the bow shock may form circular and elliptical rings, and produce elliptical isophotes when viewed at high inclination. Varying the inclination can also produce X-ray deficits. At late times of the 2.5D simulation, the beam is unstable, barely reaching the tip of the bow shock. This could be expected from stability analysis, because a non-relativistic very light jet cannot keep a high Mach number, which is necessary for stability. Regarding the propagation of the tip of the bow shock, an armlength asymmetry of a few percent was measured in both simulations. For late times, the jets are faster on the side with the on average higher density. The X-ray structure in Abell 2052 was shown to be explainable by two elliptical rings, one from the outer bow shock part of a former jet and one from the inner one. From the ellipticity, an inclination of 37 is derived consistent with other estimates in the literature. Based on the simulation results, the low inclination is one reason for the pronounced X-ray deficit inside of the rings. The X-ray gas in Hercules A and Cygnus A shows elliptical isophotes elongated in the direction of the radio jets. The ellipticity decreases at a certain distance from the center. The location where the ellipticity drops was identified as being the bow shock position. For Cygnus A,

the quality of the data is sufficient to determine the sideways bow shock position to 66 kpc. At the same position, a temperature jump can be infered from the X-ray data, implying a sideways bow shock velocity between 1300 km/s and 4200 km/s, i.e. a Mach number between one and three. The width of the radio cocoon is only a quarter of the bow shock width, which is consistent with the observation of an X-ray deficit, and shown to be possible in the 2.5 D simulation presented here. However, in self-similar models, this fraction is usually above 80%, which is a shortcoming of these models. Applying the simulation results and using the width and velocity of the sideways bow shock, a jet power of > 1047 erg/s and an age of 24 Myr is derived. This result was found to be inconsistent with the jet being non-relativistic, and a lower limit for the Lorentz factor of ten was infered. Explaining the belts of low temperature X-ray gas within the radio cocoon of Cygnus A by entrained gas over the contact discontinuity, it was suggested that such gas could cool further to form emission line regions, if the density of the ambient gas would be somewhat higher. This scenario could be realised at higher redshift and explain the origin of the gas producing the alignment effect.
Acknow ledgements. This work was supported by the Deutsche Forschungsgemeinschaft (Sonderforschungsbereich 439). The computations have been carried out on the NEC-SX5 of the HLRS in Stuttgart (Germany). I thank P. Strub for help with the X-ray data, and C. Carilli for providing the radio data.

References
Alexander, P. 2002, MNRAS, 335, 610 Basson, J. F. & Alexander, P. 2003, MNRAS, 339, 353 Best, P. N., Longair, M. S., & Rottgering, H. J. A. 1996, MNRAS, 280, L9 Bicknell, G. V., Saxton, C. J., & Sutherland, R. S. 2003, PASA, 20, 102 Blanton, E. L., Sarazin, C. L., & McNamara, B. R. 2003, ApJ, 585, 227 Blanton, E. L., Sarazin, C. L., McNamara, B. R., & Wise, M. W. 2001, ApJ, 558, L15 Bodo, G., Massaglia, S., Ferrari, A., & Trussoni, E. 1994, A&A, 283, 655 Bodo, G., Massaglia, S., Rossi et al. 1995, A&A, 303, 281+ Bodo, G., Rossi, P., Massaglia, S. et al. 1998, A&A, 333, 1117 Carilli, C. L. & Barthel, P. D. 1996, A&AReview, 7, 1 Carvalho, J. C. & O'Dea, C. P. 2002a, ApJS, 141, 337 --. 2002b, ApJS, 141, 371 Clarke, D. A., Harris, D. E., & Carilli, C. L. 1997, MNRAS, 284, 981 Gizani, N. A. B. & Leahy, J. P. 2003, MNRAS, 342, 399 --. 2004, ArXiv Astrophysics e-prints Heinz, S., Reynolds, C. S., & Begelman, M. C. 1998, ApJ, 501, 126+ Inskip, K. J., Best, P. N., Rawlings, S. et al. 2002, MNRAS, 337, 1381


22

M. Krause: Very Light Jets II
10 d=0, =-1 d=1, =-1

Kaiser, C. R. & Alexander, P. 1999, MNRAS, 305, 707 Krause, M. 2002, A&A, 386, L1 --. 2003, A&A, 398, 113 Krause, M. & Camenzind, M. 2001, A&A, 380, 789 Krause, M. & Camenzind, M. 2002, in High Performance Computing in Science and Engeneering '01, eds.: Krause, E. and J¨ ager, W., Springer, 329+ --. 2003, New Astronomy Review, 47, 573 Marti, J. M. A., Mueller, E., Font, J. A. et al. 1997, ApJ, 479, 151 McCarthy, P. J. 1993, A&AReview, 31, 639 Reynolds, C. S., Heinz, S., & Begelman, M. C. 2001, ApJ, 549, L179 Rosen, A., Hughes, P. A., Duncan, G. C., & Hardee, P. E. 1999, ApJ, 516, 729 Saxton, C. J., Sutherland, R. S., Bicknell, G. V. et al. 2002, A&A, 393, 765 Shu, F. H. 1992a, Physics of Astrophysics, Vol. I (University Science Books) --. 1992b, Physics of Astrophysics, Vol. II (University Science Books) Smith, D. A., Wilson, A. S., Arnaud, K. A. et al. 2002, ApJ, 565 Spergel, D. N., Verde, L., Peiris, H. V. et al. 2003, ApJS, 148, 175 Wilson, A. S., Young, A. J., & Shopbell, P. L. 2000, ApJ, 544, L27 Zanni, C., Bodo, G., Rossi, P. et al. 2003, A&A, 402, 949 Ziegler, U. & Yorke, H. W. 1997, Computer Physics Communications, 101, 54

1 r/r0 0.1 0.01 0.01

0.1

1 t/t0

10

100

Figure for late energy radius.

A.1. Numerical solutions of the force balance equation spherical blast waves. The blast wave with constant injection changes the power law slope at the critical The one with finite energy falls back.

static equilibrium, this equation can be nondimensionalised and rearranged in the following way: 2 2 (t/t ) 0 r r0
+5

= (t/t0 )d -

r r0

+3

.

(A.2)

The scales r0 and t0 are given by: r0 = t
(d-3)/5 0

9( + 5)/5cs t 9 5
1/2

0 1/5

=

2 +3

( + 5)

3/10

0 L

1/5

cs .

App endix A: Late phase of spherical blastwaves
For which radii is it possible to use the blastwave approximation, i.e. the solutions of (1)? This equation was derived under the assumption of vanishing external pressure, which is equal to a strongly supersonic shock. In the simulations we get weak bow shock Mach numbers very soon. Then, the external pressure (pext ) and possibly also the gravitation (remember that the atmosphere is held in equilibrium of pressure and gravity) has to be taken into account. This will be discussed in the following for the case of isothermal power law atmospheres, applicable to the presented work. The force ballance equation for a shell, driven by energy input E (t) = Ltd into an environment with gas mass r profile M(r) = 0 4 r2 (r)dr, where the density distribution is given by (r) = 0 (r/r0 ) , can be written as follows: (Mv ) = S (p t
int

For typical values of galaxy clusters, say 0 = 10-26 g/cm3 , constant sound speed of cs = 1000 km/s, and L = 1047 erg/s for d = 1 the typical scale is Mpc, reducing by a factor of a few when the jet turns off and the energy stays constant at some 1062 erg (d = 0). The force balance equation was numerically integrated, an example is shown in Fig. A.1. The general result is that blast waves with finite energy fall back after reaching the critical radius, the formal solution being oscillatory. The ones with constant energy injection break and change the slope of their power law at that radius. Up to that radius, the propagation follows the same power law that can be derived using the strong bow shock hypothesis. The conclusion is that for the simulated jets, and the typical real sources, the critical radius has not yet been reached, and the blastwave equation of motion can be expected to be a good approximation as long as non-spherical effects remain small.

-p

ext

)-

GMDM M . r2

(A.1)

Here, S = 4 r2 , G is the gravitational constant, MDM denotes the gravitating Mass profile of the dark matter halo, and the internal pressure is given by pint = 2(E (t) - Mv 2 /2)/3V , where V = 4 r 3 /3. Using hydro-