text:
# HI Signatures of Reionization

## 1. Exploring the Dark Ages

The Dark Ages are ended by the appearance of the first stars and/or quasars, that reheat the diffuse cosmic baryons and then reionize the Universe (see Rees 1999). The reionization is defined as the epoch when the volume–filling factor of HII regions is $`1`$, and the Universe becomes transparent to ionizing radiation (for an extensive review see Loeb & Barkana 2001). The absence of a Gunn–Peterson (GP) effect in the line of sight of distant quasars, put the epoch of reionization ($`z_{reion}`$) at redshifts $`>5`$. To date, there are little additional constraints on the physics of the cosmic baryons at such high redshifts. This is also due to the complexities involved by any theoretical model that must include the nature and the birthrate of the first luminous objects, their spectrum and emissivity, their feedback into the surrounding Intergalactic Medium (IGM), etc.

A significant step forward has been recently obtained with Keck spectroscopy of the most distant SLOAN quasars (Becker et al. 2001; Djorgovski et al. 2001). The presence of a sudden increase in the Ly$`\alpha `$ opacity between redshift 5 and 6 may indicate that the Universe is approaching reionization at $`z6`$. However, Barkana (2001) pointed out that, due to the high Ly$`\alpha `$ opacity, these observations are also consistent with a post–reionization phase, and that a conclusive proof requires the observation of similar GP absorptions along several additional line of sights. Alternatively, the observation of smoking–gun signatures has been proposed. Among them: dectection of scattered Ly$`\alpha `$ emission around sources beyond $`z_{reion}`$ (Loeb & Rybicki 1999); features in galaxy number counts at $`z>z_{reion}`$ (Barkana & Loeb 2000). These observations would be effective probes of the physics of cosmic baryons and of the nature of the first ionizing sources at the same time, but they must wait for NGST and/or SIRTF, which will be operating at the end of this decade. Signatures of the reionization are expected also in the CMB. MAP and Planck can probe the optical thickness at $`z_{reion}`$ in the 30–150 GHz range. In particular, fluctuations on $`0.1`$ degree reflect the topology of reionization. However, such measures are highly degenerate with cosmological parameters.

It is clear that the investigation of the Universe at $`z10`$ is a difficult task with present–day facilities. In this situation, the observational window offered by the redshifted 21cm line can open much of the Universe to a direct study of the reheating and reionization epochs. In the following we will discuss some observations of HI at high $`z`$ that can shed light on many aspects of these processes.

## 2. The History of the Cosmic Baryons

In this section we outline a simple scenario for the thermal history of the diffuse cosmic baryons in the high–$`z`$ Universe, and we show under which condition they emit in the 21cm line.

During the Dark Ages ($`z>10`$), most of the baryons of the Universe are distributed uniformly. The typical density contrast on scales $`>1h^1`$ Mpc is still in the linear regime. The temperature of the diffuse baryons, experiencing only adiabatic cooling since the decoupling from the cosmic background radiation, is about $`T_K=0.026\times (1+z)^2`$ K (Couchman 1985), which, at redshift $`z10`$, is one order of magnitude smaller than the CMB temperature $`T_{CMB}=2.73\times (1+z)`$. The spin temperature $`T_S`$ of the cosmic baryon is defined via the ratio of the population of the triplet and the singlet state of the spin transition originating the 21cm emission. In absence of any coupling, the spin of the IGM is in thermal equilibrium with the background radiation, and $`T_S=T_{CMB}`$ on a very short timescale. To observe the 21cm emission from HI against the CMB, we need some process that decouple $`T_S`$ from $`T_{CMB}`$. Collisions are ineffective, since the density contrast on Mpc scales is too low (Madau, Meiksin & Rees 1997).

There is, however, an efficient mechanism that makes HI visible in the redshifted 21cm line: the Wouthuysen–Field effect (WT). In this process, a Ly$`\alpha `$ photon field mixes the hyperfine levels of HI in its ground state via intermediate transition to the $`2p`$ state (see Appendix in Tozzi et al. 2000). The process effectively couples $`T_S`$ to the color temperature $`T_\alpha `$ of a given Ly$`\alpha `$ radiation field (Field 1958). The color temperature is easily driven toward the kinetic temperature $`T_K`$ of the diffuse IGM due to the large cross section for resonant scattering (Field 1959). Therefore the spin temperature can be written as:

$$T_S=\frac{T_{\mathrm{CMB}}+y_\alpha T_K}{1+y_\alpha },$$
(1)
where $`y_\alpha 3.6\times 10^{13}P_\alpha /T_K`$ is the Ly$`\alpha `$ pumping efficiency, and $`P_\alpha `$ is the total rate at which Ly$`\alpha `$ photons are scattered by an hydrogen atom.

Thus, the diffuse IGM at very high redshift (between recombination and full reionization) can be observed in the redshifted 21cm radiation against the cosmic background in presence of a Ly$`\alpha `$ background (likely produced by the same sources that reionize the Universe). The expected emission is less than 10 mK (Shaver et al. 1999; Tozzi et al. 2000). However if $`T_ST_KT_{CMB}`$, the HI is visible in absorption with a signal significantly larger than the maximum detectable in emission (see Hogan & Rees 1979). Thus, the most favourable condition for detection is given by a cold IGM in a strong Ly$`\alpha `$ field.

At the same time, the Ly$`\alpha `$ photon field also reheats the diffuse gas, driving $`T_K`$ towards larger values. The thermal history of the diffuse IGM results from the competition between adiabatic cooling and reheating due to the photon field, and it can be simply described as:

$$\frac{dT_K}{dz}=\frac{2\mu }{3}\frac{\dot{E}}{k_B}\frac{dt}{dz}+2\frac{T_K}{(1+z)},$$
(2)
where $`\dot{E}`$ is the heating rate due to recoil of scattered Ly$`\alpha `$ photons and, possibly, to additional contribution from an X–ray background (here $`\mu =16/13`$). The heating is quite efficient, since only 0.004 eV per particle are needed to heat the IGM above the CMB at $`z10`$. Therefore, Ly$`\alpha `$ resonant scattering or photoelectric heating by soft X–ray photons can reheat the IGM in 10% of the Hubble time (Madau, Meiksin & Rees 1997), and the WT mechanism leads preferentially to emission in the 21cm line.

The signal from high–$`z`$ HI can be predicted by solving Eq. 1 and Eq. 2, after assuming a value for the Ly$`\alpha `$ radiation field. A reference value is the thermalization rate $`P_{th}=(7.6\times 10^{12}s^1)(1+z)/10`$, which is the critical value that would drive $`T_S`$ towards $`T_K`$, and, at the same time, would heat the baryons at a rate of about 200 K per Gyr (see Madau, Meiksin & Rees 1997).

## 3. 21cm Signatures of Reionization

Between $`5<z_{reion}<20`$ the Universe is expected to experience a phase transition from the neutral to the reionized phase. A first effect of the reionization is the sudden disappearance of the 21cm emission from HI. In fact, Shaver et al. (1999) showed that with a bandwidth of 5 MHz, and an observation of 24 hours, the reionization signal (Figure 1) is detected at 100$`\sigma `$ independent of telescope size. For this measure sensitivity is not an issue, the challenge concerns signal contamination. In fact, the sharpness of this feature, expected in the 70–240 MHz range, makes it recognizable from galactic and extragalactic foregrounds (Di Matteo et al. 2001). In addition, Gnedin & Ostriker (1997) computed the 21cm emission in a specific $`\mathrm{\Lambda }`$CDM cosmology within an N–body simulation, and showed that the fluctuations in the 21cm emission will dominate the fluctuations in the CMB by two orders of magnitude (Figure 2). The global, all–sky signature of the reionization, can be easily detected by present–day radio telescopes also in the case of a non-homogeneous, patchy process, as in the case of isolated, strong ionizing sources.

Another all–sky, distinctive signature is expected at lower fequencies. A narrow but deep absorption feature would appear in correspondence of the short but significant interval when the Ly$`\alpha `$ field reach the thermalization value while $`T_S<T_K`$. The signal computed in the case of heating only by Ly$`\alpha `$ photons, is shown in Figure 3. Such a strong absorption feature ($`40`$ mK) marks the transition from a cold and dark universe to a universe populated with radiation sources, but not yet reionized. The effect is weakly dependent on the reheating epoch and on the adopted cosmology. On the other hand, the amplitude of the signal will be strongly dependent on the timescale $`\tau `$ on which the Ly$`\alpha `$ field reaches the thermalization value. In any case, the maximum of the absorption is always larger than 10 mK (see Tozzi et al. 1999).

Thus, the reheating and reionization epochs are expected to leave a double signature superimposed on the CMB in the 100–200 MHz frequency range. The detailed appearance of such a feature will depend on the spectrum of the ionizing sources. A physical scenario, including a well defined population of sources, can be built with the code CUBA (F. Haardt and P. Madau, see http://pitto.mib.infn.it/ haardt/cosmology.html), which provide the evolution of the background radiation to be plugged into Eq. 2.

A tomography of the IGM at high $`z`$ can be obtained scanning this frequency range. Apart from the all–sky, global signatures that are potentially detectable with a single–dish telescope, an array with a sensitivity of few $`\mu `$Jy can detect also the fluctuations in the neutral baryons on scales $`>1`$’, corresponding to few comoving Mpc. In the linear regime, the fluctuations induced in the brightness temperature will be directly proportional to $`\mathrm{\Delta }\rho /\rho `$, allowing a straightforward reconstruction of the perturbation field at that epoch. In Figure 4 we show results for two cosmologies, a tilted critical CDM universe (tCDM), and an open $`\mathrm{\Omega }_0=0.4`$ universe (OCDM). In both cases the fluctuations are normalized to reproduce the local abundance of clusters of galaxies. In OCDM the fluctuations are larger by a factor of 3 since for a given local normalization at $`z=0`$, the amplitude of the perturbations at high $`z`$ is correspondingly larger than that in tCDM. In both figures the density field has been evolved with a collisionless N–body simulation using the Hydra code (Couchman, Thomas, & Pearce 1995), where the baryons are assumed to trace the dark matter distribution without any biasing.

If reheating is provided by a sparse distribution of quasars, 21cm emission on Mpc scales will be produced in the quasar neighborhood (outside the HII bubble) as the medium surrounding it is heated to $`T_S=T_K>T_{\mathrm{CMB}}`$ by soft X-rays from the quasar itself. The emission region is followed by an absorption ring, since the Ly$`\alpha `$ photons reach regions where $`T_K<T_{CMB}`$. The size and intensity of the detectable 21cm region will depend on the quasar luminosity and age. The radio map resulting from a quasar ‘sphere of influence’, 10 Myr after it turns on at $`z=8.5`$ (tCDM) is shown in Figure 5 (left). The absorption region is limited to a small edge. Moreover, in Figure 5 (right), we show a quasar with the same Ly$`\alpha `$ luminosity but with a spectrum absorbed at energies larger than the Lyman limit, as in the case with an intrinsic $`N_H>10^{22}`$ cm<sup>-2</sup>. Consequently the HII region is reduced, and the X–ray warming front is well behind the light radius. This occurrence leads to a larger absorption ring where the signal reaches $`40`$ mK. In addition, imaging the gas surrounding a quasar in 21cm emission could provide a direct means of measuring intrinsic properties of the source, like the emitted spectrum and the opening angle of quasar emission.

The proposed observations are an obvious target for the Square Kilometer Array (SKA, see http://www.ras.ucalgary.ca/SKA/science/science.html). If $`2/3`$ of SKA were reserved for a compact configuration ($`7`$ km in diameter) the IGM fluctuations would be detectable at 150 MHz with 1’ resolution and a bandwidth of 1 MHz, as also quasar emission/absorption shells. In Figure 2 of Tozzi et al. (2000), we showed that a 5$`\sigma `$ detection of fluctuations on scales larger than 2’ can be reached with less than 100 hours of integration.

A potential problem of this kind of observation is the extremely low frequencies involved. To explore the range $`10<z<20`$, where the reheating phase is expected, we need good sensitivity and resolution down to 70 MHz and lower, where technical problems due to the structure of the ionosphere arise and a dedicated instrument is needed. This bandwidth will be observable with the Low Frequency Array (LOFAR). Such observations will have greater probability of detecting the reheating epoch at high $`z`$, taking advantage of the strong signal expected during the short epoch when the cold HI reaches its maximum optical depth against the CMB.

As we said above, we do not rely only on future radio facilities. All–sky features or $`30`$’ scale fluctuations (corresponding to $`1`$ mK) would be detectable with a small but dedicated telescope, with diameter of about 250 m. A near–future radio facility like the Canadian Large Adaptive Reflector (CLAR, see S. Coté, these Proceedings) would provide a 5$`\sigma `$ detection in 300 hours of integration for a typical density fluctuation.

## 4. Conclusions

The reionization is a crucial epoch in the evolution of cosmic structures. The same ingredients that characterize it (for example: nature of the first sources, the role of black holes in galaxy formation, the role of feedback, the chemical enrichment of the IGM) will deeply affect the evolution of the Universe in the following epochs. Getting direct information about these crucial early stages of the luminous Universe is an obvious goal for the next–decade cosmology. To date, several strategies have been proposed, but most of them must wait for the next generation of astronomical facilities.

The observations of high–$`z`$ HI in the redshifted 21cm, on the other hand, offer the possibility to investigate the reionization epoch using future but also present–day radio facilities. All–sky signals from HI at redshift $`z10`$ are recognizable as discontinuities superimposed on the CMB in the wavelength range $`100200`$ MHz. In particular, the epoch of reheating can be seen as a deep ($`40`$ mK) absorption feature against the CMB, at the corresponding redshifted 21cm line. The density perturbation field at a redshift $`z5÷20`$ can be reconstructed looking for mK fluctuations at $`1^{}`$$`5^{}`$ resolution in the radio sky. Finally, luminous quasars can be seen by identifying peculiar, ring–shaped signals whose morphology depends on the source’s age, luminosity and geometry. Such observations represent an exciting scientific and technological challenge, and will constitute a unique investigation of the Dark Ages, for many aspects complementary to future observations at other wavelengths.

## References

Barkana, R., & Loeb, A. 2000, ApJ, 539, 20

Barkana, R. 2001, New Astr. in press, astro-ph/0108431

Becker, R.H., et al. 2001, AJ in press, astro-ph/0108097

Couchman, H. M. P. 1985, MNRAS, 214, 137

Couchman, H. M. P., Thomas, P. A., & Pearce, F. R. 1995, ApJ, 452, 797

Di Matteo, T., et a. 2001, ApJ in press, astro-ph/0109241

Djorgovski et al. 2001, ApJL in press, astro-ph/0108069

Field, G. B. 1958, Proc. I.R.E., 46, 240

Field, G. B. 1959, ApJ, 129, 551

Gnedin, N. Y., & Ostriker, J. P. 1997, ApJ, 486, 581

Hogan, C.J., & Rees, M.J. 1979, MNRAS, 188, 791

Loeb, A., & Rybicky, J. B. 1999, ApJL, 524, 527

Loeb A., & Barkana, R. 2001, ARA&A, 39, 19

Madau, P., Meiksin, A., & Rees, M. J. 1997, ApJ, 475, 429 (MMR)

Rees, M.J. 1999, Physics Reports, astro-ph/9912345

Shaver, P. A., et al. 1999, A&A, 345, 380

Tozzi, P., et al. 1999, Proc. of the 19th Texas Symposium, astro-ph/9905199

Tozzi, P., Madau, P., Meiksin, A., & Rees, M. J. 2000, ApJ, 528, 597


format: text

text:
# The hydrodynamics of dead radio galaxies

## 1 Introduction

There is now overwhelming observational and theoretical evidence that powerful radio galaxies possess highly collimated and relativistic twin jets of matter that emerge from a central active galactic nucleus (AGN). In the powerful Faranoff-Riley type II radio galaxies (FR-II; Fanaroff & Riley 1974), the jets are thought to remain relativistic along most of their length before terminating in a series of shocks resulting from interaction with the surrounding material. A substantial amount of theoretical and numerical (hydrodynamic) work suggests that, after passing through the terminal shock, the spent jet material inflates a broad cocoon which encases the jets (e.g., see Scheuer 1974; Norman et al. 1982; Lind et al. 1989; Begelman & Cioffi 1989; Cioffi & Blondin 1992; Hooda, Mangalam & Wiita 1994).

In the early life of a powerful source, the cocoon is probably overpressured (Begelman & Cioffi 1989) with respect to the surrounding ambient material (either the interstellar medium \[ISM\] of the host galaxy, or the intracluster medium \[ICM\] of the host cluster). The cocoon then undergoes supersonic pressure-driven expansion into the ambient medium, sweeping the ambient medium into a shocked shell. A contact discontinuity separates the relativistic cocoon material from the shocked ambient material. From the onset of activity, Kelvin-Helmholtz (KH) instabilities act to shred the contact discontinuity. If the ambient medium has a density profile characterized by $`\rho r^\alpha `$ where $`\alpha >2`$, the contact discontinuity will accelerate and hence will be Rayleigh-Taylor (RT) unstable. However, in the more physical case of $`\alpha <2`$ the contact discontinuity will possess an initial deceleration, which is sufficient to stabilize the contact discontinuity against RT modes. Only at later times, once the cocoon comes into approximate pressure balance with its surroundings, will the deceleration of the contact discontinuity no longer exceed the local gravitational acceleration, thereby rendering it RT unstable. Eventually, AGN activity will cease and the jets will turn off. Although this stage of a radio galaxy’s life has been little studied, it seems likely that the combined action of the KH and RT instabilities will transform the remnant of the cocoon into ‘plumes’ that rise in the potential of the galaxy/cluster under the action of buoyancy forces (Gull & Northover 1973; Churazov et al. 2000; Brüggen & Kaiser 2000).

Observationally, radio galaxies are often seen to be associated with galaxy or cluster cooling flows. This raises an obvious and interesting set of questions. To what extent is radio galaxy activity a natural response to the cooling flow phenomenon? Do radio galaxies act back on their environment to partially offset the cooling flow, i.e., are radio galaxies a dramatic manifestation of le Chatelier’s principle? Is radio galaxy activity a crucial component for our understanding of galaxy clusters? Since it takes only a small fraction of the cooling flow mass to fuel a powerful AGN, it is easy to see how radio galaxy activity can result from a cooling flow. However, the impact of a radio galaxy on its environment is much less well understood.

In this paper, we use hydrodynamic simulations to study the evolution of a radio galaxy situated at the center of a galaxy/cluster. In particular, we follow the evolution of the source at times after the jets have terminated. We examine the interaction of the radio galaxy with the ISM/ICM and assess the environmental impact of the AGN activity. Section 2 presents our set of hydrodynamic simulations and focuses on some numerical issues. In Section 3, we briefly discuss the ‘active’ phase of the source. Section 4 presents our results for dead sources (i.e., the passive phase after the jets have turned off), and Section 5 discusses some astrophysical implications as well as the limitations of our calculations. Section 6 draws our conclusions together.

## 2 The hydrodynamic simulations

### 2.1 Basics and simulation setup

We model a source with back-to-back supersonic jets situated at the center of a spherical gaseous atmosphere. To permit parameter studies of various high resolution simulations, we make the assumption of axisymmetry and neglect magnetic fields. We discuss these assumptions in Section 5. The simulations were performed using the ZEUS-3D code (Stone & Norman 1992a,1992b; Clarke, Norman & Fiedler 1994) in its 2D (axisymmetric) pure hydrodynamic mode. ZEUS-3D is an explicit Eulerian finite difference code, formally of second order, which uses an artificial viscosity to reproduce shocks. Major advantages of this code are its flexibility in the choice of the computational grid, and the fact that it has been extensively tested. The simulations were performed on a Sun Microsystems Ultra 60 workstation.

Unlike many previous simulations, we model both jets and make no assumption regarding reflection symmetry in the plane normal to the jets. This allows us to study phenomena associated with the collision of the two backflows. Our simulations were performed in spherical polar coordinates $`(r,\theta ,\varphi )`$, and the computational domain was the region $`r(0.1,10)`$, except for one run (Run 4 below) which had the domain $`r(0.1,30)`$. All gas within the calculation was assumed to have an adiabatic index of $`\gamma =5/3`$. The ambient galaxy/cluster gas was assumed to be initially in an isothermal configuration with a (adiabatic) sound speed $`c_\mathrm{s}=1`$ and a density profile given by a $`\beta `$-model with $`\beta =0.5`$:

$$\rho (r)=\frac{1}{[1+(r/r_0)^2]^{3/4}},$$
(1)
where the core radius was set to $`r_0=2`$. The gravitational potential, $`\mathrm{\Phi }`$, was set so as to keep this ambient material in hydrostatic equilibrium

$$\frac{\mathrm{d}\mathrm{\Phi }}{\mathrm{d}r}=\frac{1}{\rho }\frac{\mathrm{d}p}{\mathrm{d}r},$$
(2)
which gives (to within an additive constant)

$$\mathrm{\Phi }=\frac{c_\mathrm{s}^2}{\gamma }\mathrm{ln}(\rho ).$$
(3)
It is assumed that the gravitational potential is dominated by a stationary distribution of background dark matter, i.e., the self-gravity of the gas is negligible. This reasonable assumption prevents us from having to solve Poisson’s equation as part of the hydrodynamic simulations.

We inject back-to-back jets into this background medium via the use of inflow boundary conditions at $`r=0.1`$. The jets are initially conical with a $`15^{}`$ half opening angle, and in pressure equilibrium with the ambient medium. They are given an initial density of $`\rho _{\mathrm{jet}}=0.01`$ and Mach number (with respect to the sound speed of the injected jet material) of $`=10`$. As discussed in the Appendix, these parameters are chosen so that the Kelvin-Helmzholtz growth rate at the cocoon-ICM contact discontinuity in this non-relativistic simulation is approximately the same as the corresponding growth rate in the real relativistic situation. At time $`t=1`$, the jet activity is stopped and the inner boundary is made reflecting. The source is then allowed to evolve passively.

We can relate quantities within the code to physical quantities by fixing the parameters of the background medium. Suppose we set $`r_0=100\mathrm{kpc}`$, $`c_{\mathrm{ISM}}=1000\mathrm{km}\mathrm{s}^1`$ and a central number density of $`n_0=0.01\mathrm{cm}^3`$ — values representative of a rich galaxy cluster. Then, one code unit of time corresponds to 50 Myr. The total kinetic luminosity of the source is then $`9.3\times 10^{45}\mathrm{erg}\mathrm{s}^1`$. These parameters are relevant for powerful sources in rich galaxy clusters such as Cygnus A (Arnaud et al. 1984; Reynolds & Fabian 1996) and 3C 295.

Unless otherwise stated, this paper shall use the above scalings when converting our simulation results into quantities that can be compared directly with real systems. However, in the present universe at least, many radio galaxies are found in smaller galaxy clusters or poor groups. Taking objects like Hydra A and Virgo A as prototypes, the typical relevant ISM/ICM parameters are $`r_0=10\mathrm{kpc}`$, $`c_{\mathrm{ISM}}=500\mathrm{km}\mathrm{s}^1`$ and $`n_0=0.1\mathrm{cm}^3`$ (e.g., see the detailed analysis of the Chandra X-ray Observatory data for Hydra A by David et al. 2000). In this case, one time unit corresponds to 10 Myr and the total kinetic luminosity of the source is $`1.16\times 10^{44}\mathrm{erg}\mathrm{s}^1`$. Of course, it must be appreciated that these $`\beta `$-models are very crude parameterizations of the real ISM/ICM density structure.

### 2.2 The need for adequate numerical resolution

During the active phase, the large scale evolution and properties of the cocoon are strongly influenced by much smaller scale shock structures that appear in the jet. As will be discussed in more detail below, a conical shock appears in the jet which sprays the jet thrust over a large working surface (which expands almost self-similarly). Failure to resolve this shock structure in a simulation leads to the jet thrust being deposited over a very small area which leads to the jet rapidly ‘drilling’ a narrow cavity through the ICM.

Initially, we performed simulations with $`r_{\mathrm{out}}=10`$ at three different resolutions: $`n_\mathrm{r}\times n_\theta =300\times 300`$ (Run 1), $`600\times 600`$ (Run 2), and $`1200\times 1200`$ (Run 3). The resolution was slightly enhanced at small $`r`$ and along the jet axes. We find that Run 1 does indeed fail to resolve the jet-shock structure. On the other hand, Run 2 and Run 3 both well resolve the internal shocks within the jet and produce cocoons that expand in almost a self-similar manner, at least during the early active phase. The similarity of these two simulations gives us some confidence that we have achieved a degree of numerical convergence.

We note one additional numerical technicality. Our simulations possess a relatively high dynamical range in $`r`$. However, to maintain adequate resolution on large scales, we choose not to use a logarithmic (“scaled”) grid for the $`r`$ variable. Instead, we space our grid using a more general geometric progression (i.e. “ratioed” scaling). While allowing us to maintain adequate resolution on the scales that matter to our problem, it has the undesirable effect of forcing highly elongated grid elements to exist in the very central regions of the simulation. One might then worry about the effects of an effective anisotropic artificial viscosity in this inner region. However, this is acceptable since we are not interested primarily in the detailed hydrodynamic properties of these inner regions. One can view the innermost region of our grid as a buffer zone, which connects smoothly to the interesting region of the grid, while protecting it from the hard inner boundary. Experiments with high resolution simulations using a logarithmic $`r`$-grid give us confidence that this is not a problem.

For the rest of this paper, we discuss results from Run 2 (i.e., $`r_{\mathrm{out}}=10`$ with a grid size of $`n_\mathrm{r}\times n_\theta =600\times 600`$), and one other simulation. This final simulation (Run 4) has a larger outer radius ($`r_{\mathrm{out}}=30`$) and a larger computational grid ($`n_\mathrm{r}\times n_\theta =1200\times 600`$). In the spatial region covered by both Run 2 and Run 4, the resolutions are approximately (but not exactly) the same.

## 3 The active phase

There is already a substantial body of numerical work on the hydrodynamics, and magnetohydrodynamics, of active radio galaxy jets. Thus, we shall discuss only briefly our simulation results for the active phase of our source. In a companion paper (Reynolds, Heinz & Begelman 2000) we focus on the X-ray appearance of such active sources.

Our results during the active phase are very similar to those obtained by Lind et al. (1989), who performed axisymmetric hydrodynamic simulations of a single jet propagating through a uniform medium. The initially conical jet is focused by the external pressure and suffers a series of oblique shocks. Near the working surface with the background medium, a strong annular shock deflects most of the flow into a wide fan. Thus, the jet thrust is deposited over a large area, thereby reducing the advance speed of the jet below the value that would be obtained if the jet just ‘drilled’ into the background medium. This is the 2-dimensional equivalent of the so-called “dentist drill effect” (Scheuer 1974). After passing through the annular shock, the spent jet material enters a mildly supersonic backflow which gradually decelerates through a series of weak shocks. Powerful Kelvin-Helmholtz (KH) instabilities exist along the contact discontinuity between the backflowing jet material and the shocked/compressed ambient material. These KH instabilities, which become more powerful as the backflow decelerates toward the region where the two backflows collide, act to mix ambient material into the cocoon. These instabilities can be seen clearly in Fig. 1, which shows specific entropy maps ($`s=p/\rho ^{5/3}`$; note that we neglect logarithms and constants in our definition of the specific entropy) for Run 2 at times $`t=0.5`$, $`t=1.0`$, $`t=3.0`$ and $`t=5.0`$. During the active phase, the structure expands approximately self-similarly, at least as long as the cocoon is still sufficiently overpressured as to undergo supersonic lateral expansion into the ambient material.

It is worth noting explicitly that the backflow in the cocoon has a speed comparable to the sound speed of the cocoon material. In code units, the advance speed of the working surface of the jet is $`3`$ (i.e., a Mach number of 3 with respect to the sound speed of the undisturbed ambient material). This should be compared with a backflow speed of $`2030`$ which exists in a large portion of the cocoon. Translated into terms applicable to real radio galaxies, the backflow speed would be a significant fraction of the speed of light (since the spent jet material is a relativistic plasma with a sound speed of $`c/\sqrt{3}`$).

## 4 The passive phase

As discussed in Section 2, we turn off the jets after $`t=1`$ and allow the existing structure to evolve. Both Run 2 and Run 4 display very similar behavior during the passive phase although, of course, Run 4 can be followed to much later times due to the larger spatial domain of the calculation (see Fig. 2). We will describe our findings based on Run 4. It should be noted that Brüggen & Kaiser (2000) and Churazov et al. (2000) have also simulated the late evolution of radio galaxies. Their simulations make the assumption that the source evolution before the jets turn off is governed by simple supersonic expansion of the lobes and consequently, their initial conditions for the passive phase are based on relatively simple geometries; namely, spherical and ellipsoidal bubbles placed in isothermal atmospheres. However, turbulence and powerful circulatory motion seems to be present at all times inside the active cocoon, which leads to significant excitation of the Kelvin-Helmholtz instability before the source turns off. Whether, and by how much, the cocoon is overpressured at the moment that the jet activity ceases can be estimated by the simple analytic models of Begelman & Cioffi (1989).

Thus, we extend the work of Churazov et al. (2000) and Brüggen & Kaiser (2000) by simulating the active phase of the source and then shutting the jets off, rather than modelling the late stages of evolution by letting a static bubble evolve under the action of buoyancy. We also simulate a substantially larger portion of the cluster atmosphere while maintaining the required high resolution, thereby allowing us to track the source evolution at later times. We consider these aspects to be important extensions of these previous works.

### 4.1 The stages of an inactive source’s life

As would be expected from simple qualitative arguments, turning off the jets results in the almost immediate disappearance of internal shocks in the jet channels, which subsequently collapse. The very rapid backflow motions in the cocoon slow down accordingly, due to the lack of high-inertia material and the disappearance of the terminal shock, which previously re-directed the flow. While the radial lobe expansion during the active phase was governed mostly by the jet thrust, the absence of ram pressure from the jets in the inactive phase results in a rapid slowing of the cocoon head to subsonic speeds. If the cocoon is still significantly overpressured when the jet activity ceases, this subsonic stage will be preceded by a Sedov expansion phase.

During the early inactive phase of our simulations, the source is still slightly overpressured with respect to the ambient medium and we do indeed find a Sedov phase, during which the cocoon expands in order to achieve pressure balance. When pressure balance is reached (at $`t1.5`$), the driving force on the shocked/compressed ICM shell ceases and this shell rapidly decelerates to subsonic speeds. Once the lobes reach rough pressure equilibrium with their surroundings, the dynamics of the lobes closely resemble a rising bubble. In the initial stages of this late inactive, or buoyant, stage, the denser IGM begins to settle back into the center, squirting the light cocoon material along the major axis of the cocoon. Eventually, the ICM core will reform and the cocoon will be pinched off into two separate, buoyantly rising plumes.

Material that had been entrained previously is now carried along by each plume, partly due to its own inertia and partly due to the significant circulation inside the plume. This vorticity continues to drag ambient material into the bubble at the plume’s back end, where strong eddies produce a jelly-fish-like structure (Fig. 3). At the head of the plume, a dense cap of swept up material forms, which is pushed ahead.

Continuing with the bubble analogy, we can write the buoyant velocity $`v_\mathrm{c}`$ of the plume as

$$v_\mathrm{c}𝒞\sqrt{\frac{4r_{\mathrm{bub}}g(\rho _\mathrm{a}\rho _\mathrm{c})}{\rho _\mathrm{a}}}$$
(4)
where $`r_{\mathrm{bub}}`$ the radius of the bubble, $`g`$ is the local gravitational acceleration, and $`\rho _\mathrm{a}`$ and $`\rho _\mathrm{c}`$ are the densities of the ambient and bubble materials respectively. The top left panel of Fig. 4 shows the volume averaged plume height (i.e., its distance from the center). The plumes have approximately constant velocity of $`v_\mathrm{c}0.8`$. Together with the cocoon size of $`r_{\mathrm{bub}}(t=4)1`$ and the gravitational acceleration at that position, this indicates that $`𝒞`$ is slightly less than unity.

### 4.2 Energetics and Thermodynamics

One goal of these simulations is to determine a set of physical parameters that are only weakly dependent on the details of the simulations (most importantly, they should be weakly dependent on the numerical resolution and the dimensionality). In particular, we are interested in the overall matter distribution, energetics, and thermodynamics of the system. For this discussion, it is interesting to separate the system into different structural components, i.e., the cocoon/plumes, and the ambient medium (including the shocked/compressed shell that bounds the radio source).

Defining these regions is not trivial. Since ZEUS-3D is an Eulerian code we cannot, a priori, separate cocoon plasma from ambient plasma by tracking fluid elements. Due to the turbulent cocoon surface, a simple geometric identification of the cocoon is also not an option. However, since the jet material is of very high entropy compared to the background gas (especially once it has passed through the thermal shock), we use the specific entropy $`s=p/\rho ^{5/3}`$ as a discriminant.

#### 4.2.1 The entropy threshold method and numerical mixing

For the simulation parameters we chose, the background material has specific entropy (in code units)

$$s=\frac{p}{\rho ^{5/3}}=\frac{3}{5}[1+(r/r_0)^2]^{1/2}.$$
(5)
The injected material, on the other hand, has entropy $`p/\rho ^{5/3}1300`$ and cocoon material, which has gone through the terminal shock, can have entropy as high as $`s10^5`$. We chose to define the cocoon as the volume inside of which the entropy satisfies $`s_\mathrm{c}>10`$. At all but the very earliest times in our simulations, the shock driven into the ambient material is too weak to raise its entropy above this threshold. Conversely, an examination of entropy maps from the simulations shows that, until about $`t8`$, the contact discontinuity between the cocoon and ambient material is sharp and well tracked by the $`s_\mathrm{c}=10`$ contour. In fact, until $`t5`$, the location of the contact discontinuity is insensitive to the chosen entropy threshold, provided that threshold is in the range $`10<s_\mathrm{c}<100`$. For times $`t>8`$ numerical mixing reduces substantial parts of the cocoon/plumes below the $`s_\mathrm{c}=10`$ threshold. At these late times, identification of a well-defined contact discontinuity separating the plume from the ambient material becomes problematic.

To better illustrate the effects of numerical mixing at later times, we calculate the entropy distribution of Run 4. Figure 5 shows this function at different times. Due to numerical mixing the second peak at high entropy (i.e., the cocoon material) diffuses to lower entropies. It is clear that after a certain amount of time an entropy cut is not a reasonable estimator for cocoon material anymore, namely, when numerical diffusion has lowered the entropy inside the cocoon to levels comparable with the ambient entropy.

Numerical mixing is an irreversible process and hence creates entropy. On the other hand, ‘mixing’ across the contact discontinuity in real radio galaxies is not necessarily irreversible as can be seen from the following argument. Suppose that the ambient ISM/ICM material is not mixed on the microscopic level. Instead, assume that intact dense ISM/ICM filaments and clouds are simply swept up into the light cocoon and kept distinct due to being on different magnetic structures. The fact that this is a thermodynamically reversible process can be seen by noting that, given a sufficiently long time, buoyancy effects could separate the ISM/ICM gas from the cocoon plasma. Thus, the entropy structure of the simulated cocoon becomes increasingly less physical as time proceeds. Indeed, once numerical mixing has significantly affected a substantial part of the simulated cocoon (which occurs at $`t8`$), the entropy structure of the cocoon may be a poor reflection of reality.

#### 4.2.2 Definitions of energies

To explore the energetics of our simulation in a quantitative manner, we compute the following quantities:

1. Total mass of cocoon,

$$M_{\mathrm{cocoon}}(t)=_C\rho dV,$$
(6)
2. Total internal energy of cocoon (ambient) material,

$$E_{\mathrm{int},\mathrm{coc}(\mathrm{amb})}(t)=\frac{1}{\gamma 1}_{C(A)}pdV,$$
(7)
3. Total kinetic energy of cocoon (ambient) material

$$E_{\mathrm{kin},\mathrm{coc}(\mathrm{amb})}(t)=\frac{1}{2}_{C(A)}\rho v^2dV,$$
(8)
4. Total gravitational potential energy of cocoon (ambient) material

$$E_{\mathrm{pot},\mathrm{coc}(\mathrm{amb})}(t)=_{C(A)}\rho \mathrm{\Phi }dV,$$
(9)

where $`C`$ is the region in which $`s>10`$, and $`A`$ is the region in which $`s10`$. We then reference all energies to their values at the initial time:

$$\mathrm{\Delta }E=E(t)E(0).$$
(10)
#### 4.2.3 The active and early-inactive phases

The bottom panels of Fig. 4 shows how these energies vary with time (upto $`t=10`$) for Run 4. During the active phase, the energy injected by the jets is transformed into internal and kinetic energy of both the cocoon and ambient material, as well as the potential energy of the ambient material. All of these energy forms increase linearly with time as the injected energy is shared roughly according to equipartition. After the source turns off, the kinetic energy of the cocoon rapidly decreases as the rapid motions associated with the jet cease. During the short Sedov phase that follows (between $`t=1`$ and $`t1.5`$), the internal energy of the cocoon also undergoes a rapid decrease (adding mainly to the internal energy of the ambient medium), as the overpressured cocoon expands adiabatically to achieve pressure equilibrium with the ambient material. After the Sedov phase, the cocoon source then enters the buoyant phase.

At all times during the buoyant phase in our simulations, the rising buoyant plumes entrain and lift material (mainly in the trailing regions). This contributes to an increase in the gravitational potential energy of both the plume and ambient material. Concurrently, an approximately spherical pulse propagates radially outwards in the ISM/ICM atmosphere. This pulse is the remnant of the strong shock that bounded the supersonic cocoon at early times. As will be discussed below, this pulse mediates a general expansion (inflation) of the ISM/ICM atmosphere which results due to a heating of the core regions by the jet activity.

#### 4.2.4 The entropy structure at late times

The radio galaxy influences the entropy of the ISM/ICM atmosphere by driving shocks into it. At very long times after the jet activity has ceased, the ICM/ISM will be in a relaxed hydrostatic equilibrium in which the entropy is increasing outwards (or else there will be convective turnover until such a state is achieved). A parcel of gas which passed through a very strong radio-galaxy induced shock (which may have initially been in the centralmost regions of the galaxy/cluster) will have a high final specific entropy and achieve its new equilibrium state some large distance from the cluster core. If the same parcel of gas suffers a weaker shock, the entropy increase will be less and it will settle at a lower point in the final ISM/ICM atmosphere.

The density, temperature, and entropy plots in Fig. 6 show that, after the episode of radio galaxy activity, the ISM/ICM core is left on a higher adiabat than it was prior to the activity. This is easy to understand. Since the initial radio galaxy expansion was supersonic, the material which initially had the lowest specific entropy (i.e. at the center of the initial ISM/ICM distribution) has all been shock heated. The final ISM/ICM core contains a mixture of (mildly) shocked gas, and unshocked ISM/ICM which has adiabatically flowed inwards from more distant regions (and thus possesses a higher specific entropy) to displace highly shocked material. The highly shocked material is expelled to the outer regions of the galaxy/cluster. For the specific case represented by our simulation, the excess entropy over the undisturbed cluster gas peaks at 20%. This translates into an excess entropy of about 14 $`\mathrm{keV}\mathrm{cm}^2`$.

Although this is a small change in the entropy of the cluster core, one might think that repeated radio galaxy activity would lead to a significant accumulation of excess entropy in the cluster core. This depends upon the assumed entropy profile for the ISM/ICM. Within the $`\beta `$-model parameterization used here, the entropy profile of the ISM/ICM is rather shallow. Thus, even if one could imagine completely shocking a large part of the central ISM/ICM, the entropy increase of the core seen in the final equilibrium would still be rather modest since material just outside of the strongly shocked region will simply flow inwards and take its place. If the initial entropy profile is much steeper (as it is in some strong cooling flow models) then the same radio galaxy could impact the entropy of the core rather more significantly. Clearly, further work is needed on the contribution of radio galaxies to the ISM/ICM entropy in which the entropy injection is considered in relation to particular detailed cooling flow models.

Using ROSAT data, Ponman, Cannon & Navarro (1999) and Lloyd-Davies, Ponman & Cannon (2000) have found evidence for an “entropy floor” in groups and clusters. This suggests from an observational point of view that these systems possess excess entropy (as compared to the entropy profile expected from simple scaling relations and detailed simulations; however, see alternative viewpoint of Bryan 2000). Lloyd-Davies, Ponman & Cannon (2000) claim an excess entropy of 70–140 $`\mathrm{keV}\mathrm{cm}^2`$ in a sample of 20 clusters and groups. Interestingly, these authors also conclude that this entropy injection occurred at low densities ($`13\times 10^4\mathrm{cm}^3`$) and hence they suggest that it occurred prior to the collapse/formation of the cluster. In principle, microscopic mixing of the ICM and radio plasma inside of the cocoon/plumes might also lead to ICM entropy injection at low density. However, our simulations suggest that buoyancy effects tend to transport any such mixed material completely out of the cluster center. Thus, if the determination of the density at which the entropy injection occurred is robust, it does indeed seem likely that the bulk of the entropy injection occurred prior to cluster formation and is not due to more recent radio galaxy activity.

#### 4.2.5 Energetics at late times

The energetics of the very late stage of Run 4 are reported in Fig. 6. The large scale pulse mentioned in §4.2.3 manifests itself in each of the quantities plotted in Fig. 6 at $`r20`$. The top-left panel shows the difference in enclosed energies as a function of radius between late times and the initial state $`t=0`$. Comparing the values of this function on each side of the pulse may suggest to the casual reader that the wave contains most of the injected energy. This is a curious result since we have already noted appreciable entropy and internal energy increases associated with the irreversible heating of the ISM/ICM atmosphere.

The bottom panels of Fig. 6 elucidate the true nature of this pulse. In the bottom-left panel, we note that the major contributor to the increase of enclosed energy comes from gravitational potential energy (rather than kinetic and internal energy as would be expected if this were a sound wave carrying away the injected energy). Secondly, displaying the enclosed energy difference in mass coordinates (Fig. 6 bottom-right panel) shows that the mass contained within the inner edge of the wave-packet has acquired an appreciable excess of energy during the evolution corresponding to at least half of the injected energy.

Thus, it can be seen that the ISM/ICM atmosphere has been inflated (or puffed up) due to the heat deposited in the core of the atmosphere by the jet activity. The spherical pulse (travelling at the sound speed) drives mass to a higher location in the potential thereby mediating this inflation. That this inflation has occurred is obvious from the reduced core density in the final state (Fig. 6 top right). Examining the enclosed energy in mass coordinates lets one see that this inflation actually accounts for a substantial fraction of the injected energy, implying that the injected energy has been heated the ISM/ICM with a high efficiency ($`\eta 0.5`$).

## 5 Discussion

### 5.1 Cooling flows and radio galaxy heating

As discussed above, we find that a large fraction ($`\eta 0.5`$) of the energy injected by the jet activity is thermalized in the core regions of the ISM/ICM. As mentioned in the Introduction, this raises the possibility that radio galaxy activity can have significant effects on the energetics of the ISM/ICM and, therefore, act as negative feedback for cooling flow activity.

In particular, let us address the case of cluster cooling flows. Peres et al. (1998) compared the radio and X-ray luminosity of 58 clusters in a flux limited ROSAT sample. They found strong radio sources in about $`10\%`$ of the clusters within the sample. However, using the results of Peres et al. (1998) and summing over all members of the sample (not just those with radio sources), we find that the radio luminosity of the sample (defined as the sum over all sources of $`\nu L_\nu `$ measured at 1.4 GHz) is about 1% of the sample bolometric X-ray luminosity from within the “cooling radius”. To discuss the relevance of radio galaxy heating on the cluster, we must relate the radio luminosity to the kinetic luminosity of the sources.

Simple dynamical models of expanding radio lobes suggest that the maximum synchrotron luminosity of homogeneous, space-filling, lobes is only a few percent of the kinetic luminosity (e.g., see Bicknell, Dopita & O’Dea 1997). Departures from equipartition always act to decrease the synchrotron efficiency, and the efficiency further decreases as the source expands (Begelman 1996, 1999). The only way to increase this efficiency is to concentrate the radio emission into a small volume of the radio lobe (e.g., strong shocks within the lobe). In the absence of such inhomogenieties, it seems likely that the radio luminosity is at most a few percent of the kinetic luminosity. Therefore, using the Peres et al. (1998) results quoted above, the sample-averaged kinetic luminosity is comparable to the X-ray luminosity from within the cooling radius. The results of our work suggest that a large fraction of this jet power can be thermalized in the ICM.

Thus, from a purely energetic point of view, radio galaxies can balance the radiative cooling associated with cooling flows. However, it is clear that the ICM of real galaxy clusters does not experience a steady-state situation in which jet-induced heating balances radiative cooling. As has already been noted, most of the integrated radio luminosity in the Peres et al. (1998) sample is contributed by the powerful radio galaxies that reside in only 10% of the clusters. The arguments of the previous paragraph would suggest that, in these clusters, the radio galaxy heating would vastly exceed the radiative cooling. Conversely, radiative cooling would certainly seem to dominate in the remaining 90% of the systems (unless they possess powerful, but invisible, nuclear outflows; Binney 1999). Thus, it is clear that even if radio galaxy heating is relevant, the heating/cooling balance can be vastly different from cluster to cluster, or at different points of a given cluster’s lifetime.

It is beyond the scope of this paper to confront, in any direct manner, the plethora of new data that the Chandra X-ray Observatory and XMM-Newton are providing on the properties of the ICM and cooling flows. Here, we simply note some of the challenges that future models of cooling-flows and embedded radio-galaxies must meet.

Chandra finds spectroscopic evidence for cooling flows (albeit with a mass deposition rate smaller than previously estimated with ROSAT; Peres et al. 1998) even in clusters that contain energetic radio galaxies. For example, Fabian et al. (2000) find cool gas in the Perseus cluster, even though the central radio galaxy is clearly interacting with the ICM is a rather vigorous manner. Curiously, the cool gas seems to form a shell around the radio lobes and occupies precisely the location where we would expect shock/compressionally heated gas to reside. Furthermore, Chandra has observed sharp temperature discontinuities in clusters that do not appear to host powerful radio galaxies. Firstly, a few clusters (e.g., Abell 2142, Markevich et al. 2000) have large scale temperature discontinuities which have become known as “cooling fronts”. These may be discontinuities in the cooling flow associated with a group/cluster merger event. Secondly, Mazzotta et al. (2001) has found a bubble of hot gas in the cluster MKW3S. They associated this with some previous outburst of nuclear activity from the central galaxy. However, if such a bubble had buoyantly risen from the centre of the cluster, we would expect KH instabilities to disrupt it to a larger extent than observed. One would need to invoke some mechanism, such as large scale magnetic fields, to stabilize the large scale KH instabilities and prevent disruption of the bubble.

XMM-Newton has also produced a further mystery. By performing detailed emission line spectroscopy of several cooling flow clusters (e.g., Abell 1835, Peterson et al. 2001; Abell 1795, Tamura et al. 2001) with the reflection grating spectrometer (RGS) on XMM-Newton, one can probe the detailed temperature distribution of the ICM. Multiple temperatures were observed in concordance with cooling flow models, with components ranging from the virial temperature down to $`1`$ keV. However, no emission lines were observed from gas below $`1`$ keV whereas cooling models would predict large amounts of emission from such gas. Fabian et al. (2001) discusses various solutions for this puzzle and favours a model in which there are large abundance inhomogenieties affecting the predicting cooling rates. Another solution is to postulate a heating mechanism which targets cool ($`1`$ keV) gas. Examples of such a mechanism are the “percolation” model of Begelman (2001), or the “reconnection” model of Norman & Meiksin (1996). Heating models based upon simple hydrodynamic shocks (such as presented in this work or that of David et al. 2001) do not possess such temperature dependent targeting. In order to match the observations, the putative heating mechanism must meet the further challenge of heating the cool gas rapidly up to the virial temperature.

Finally, we note that neither Chandra nor XMM-Newton have yet found clear evidence for strong shocks associated with radio galaxies. Using the characterization of Paper I, all cluster radio sources observed so far appear to be in the weak-shock or the no-shock regime. Although more work is required before making a rigorous statement, there appears to be a disagreement between the observed frequency of strong shocks and current hydrodynamic models of radio galaxy evolution. This might suggest that radio jets are less powerful than we believe (i.e. the synchrotron efficiency is higher than a few percent, implying that the radio emitting regions fill a small portion of the radio lobes). This problem was first pointed out by David et al. (2001).

### 5.2 The limitations of this work

The simulations that we present in this paper are just the first step in addressing the long-term evolution of radio galaxies. In order to make the problem tractable, we have made several simplifying assumptions. Here, we briefly discuss how these assumptions might affect our results.

Our biggest simplification is that we have reduced the problem to two spatial dimensions by assuming strict axisymmetry about the jet axis. Given the inherent axisymmetry of the system, and the fact that observed radio cocoons are approximately axisymmetric, this is not an outrageous assumption. However, there are several undesirable consequences of assuming axisymmetry. Firstly, we prevent the formation of non-axisymmetric KH and RT modes. All of the available free energy is then channeled into the axisymmetric and RT modes. The axisymmetric KH modes are most efficient at mixing (since they present the most effective barrier to the flow), and hence the axisymmetry assumption maximizes the amount of mixing between the ambient and cocoon material. Secondly, the assumed axisymmetry prevents us from including the effects of any realistic random ISM/ICM motions (i.e., “weather”). These ICM motions will buffet the buoyantly rising plumes and will become crucial in determining the morphology of very old sources.

We have also chosen to ignore the presence of magnetic fields (i.e., we perform a hydrodynamic, rather than a magnetohydrodynamic, calculation). Magnetic fields in the radio cocoon might well be in equipartition with the particle energies and hence could be dynamically significant. Dynamically important fields that are tangled in an isotropic manner on small scales behave like a relativistic gas and so can be handled within the context of a pure hydrodynamic simulation. However, anisotropic tangled fields, or large scale magnetic fields, will introduce new features to the dynamics which require a full magnetohydrodynamic treatment. Similarly, the importance of neglecting the cocoon’s relativistic equation of state can only be assessed once calculations are performed using relativistic hydrodynamic codes.

## 6 Conclusions

Radio galaxies are known to be dynamical systems that evolve on timescales of $`10^710^8\mathrm{yr}`$. Hence, we expect many galaxies and clusters of galaxies to be host to dead, or relic, radio sources. It is possible, and maybe even probable, that this past radio activity has influenced the energetics and thermodynamics of the hot, space-filling medium in these hosts. Thus, a study of the environmental impact of dead radio galaxies is important and, given the high quality X-ray data now coming from the Chandra X-ray Observatory and XMM–Newton, very timely.

We have used high resolution hydrodynamic simulations to investigate one particular scenario relevant to dying radio galaxies — i.e., a rather powerful radio-loud AGN situated at the center of a galaxy or cluster of galaxies whose activity ceases abruptly. In particular, we have simulated back-to-back jets propagating in a $`\beta `$-model galaxy/cluster atmosphere. We then shut down the jet activity and let the resulting structure evolve. To make the problem more tractable, we assume axisymmetry about the jet axis (thereby reducing the problem to two spatial dimensions) and neglect the action of magnetic fields. This is a significant extension of the work of Churazov et al. (2000) and Brüggen & Kaiser (2000) since we simulate the active phase of the source and then shut off the jets, rather than modelling the late stages of evolution by simply letting a static bubble evolve under the action of buoyancy.

In the early lifetime of the simulated source (while the radio jets are still active), the shocked jet material forms a cocoon which is bounded by a shell of shocked ICM. Even during the active phase, hydrodynamic instabilities start to shred this cocoon. However, only after the jet activity ceases do KH and RT instabilities destroy the integrity of the cocoon. Thereafter, the old cocoon material forms two plumes which rise in the cluster potential through the action of buoyancy forces. These plumes entrain a significant amount of cooler material from the ICM core and lift this material high up into the cluster atmosphere. At very late times, the galaxy/cluster atmosphere settles back into hydrostatic equilibrium but with a core specific entropy that has been enhanced by $`20\%`$ over its initial value. This entropy enhancement is due to the shocking of the lowest entropy material by the strong shock which bounds the radio cocoon during the early active phase. We find that a large fraction ($`0.5`$) of the injected energy is thermalized in the ISM/ICM. Comparing late times to the initial state, we find that the ISM/ICM atmosphere has become inflated in order to maintain hydrostatic equilibrium after the thermalization of the jets energy. Such a large thermalization efficiency raises the possibility that radio galaxies are important in the overall energy budget of the ISM/ICM.

## Acknowledgements

CSR appreciates support from Hubble Fellowship grant HF-01113.01-98A. This grant was awarded by the Space Telescope Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA under contract NAS 5-26555. This work has also benefited from support by NASA grant NAG5-7329. Finally, we appreciate support from NASA under LTSA grant NAG5-6337, and the National Science Foundation under grants AST-9529170 and AST-9876887.

## Appendix

In this appendix, we present a brief derivation of the relativistic dispersion relation for the Kelvin-Helmholtz instability in cylinderical geometry. We then use this dispersion relation to connect our simulations to the behavior of real sources. Similar relativistic KH analyses have been performed by many authors (Blandford & Pringle 1976; Ferrari, Trussoni & Zaninetti 1978; Ray 1981; Birkinshaw 1984; Zaninetti 1985). We repeat the derivation here for completeness.

The continuity and momentum equations of relativistic fluid dynamics for a $`\gamma `$-law gas are

$`{\displaystyle \frac{}{x^i}}(\rho u^i)`$ $`=`$ $`0`$ (11)
$`hu^k{\displaystyle \frac{u_i}{x^k}}{\displaystyle \frac{p}{x^i}}+u_iu^k{\displaystyle \frac{p}{x^k}}`$ $`=`$ $`0`$ (12)
$`u^k{\displaystyle \frac{}{x^k}}(p\rho ^\gamma )`$ $`=`$ $`0,`$ (13)
where $`\rho `$ and $`p`$ are the comoving density and pressure, respectively, $`u^i`$ are the components of the four velocity, and $`\gamma `$ is the polytropic index. Here, $`h=\rho +\gamma p`$ is the relativistic enthalpy. We wish to use these equations to study the relativistic Kelvin-Helmholtz instability.

Consider first-order perturbations of these equations:

$`p`$ $`=`$ $`p_0+\mathrm{\Delta }p`$ (14)
$`\rho `$ $`=`$ $`\rho _0+\mathrm{\Delta }\rho `$ (15)
$`u^i`$ $`=`$ $`u_0^i+\mathrm{\Delta }u^i.`$ (16)
For simplicity, we shall drop the subscript ‘0’ in the subsequent equations. Assume that the unperturbed quantities are time-independent and spatially uniform. The linearized fluid equations can then be cast in the form

$$\left(\frac{h}{\gamma p}1\right)u^iu^k\frac{^2\mathrm{\Delta }p}{x^kx^i}=\frac{^2\mathrm{\Delta }p}{x^ix_i}$$
(17)
For definiteness, we choose cylinderical polar coordinates and say that the unperturbed velocity field is along the $`z`$-direction with velocity $`V`$. This last equation can then be written as

$$\mathrm{\Gamma }^2\left(\frac{h}{\gamma p}1\right)\left(\frac{^2}{t^2}+2V\frac{^2}{tz}+V^2\frac{^2}{z^2}\right)\mathrm{\Delta }p=\left(^2\frac{^2}{t^2}\right)\mathrm{\Delta }p.$$
(18)
Now consider the following geometry which is appropriate for the analysis of the Kelvin-Helmholtz instability in radio galaxy cocoons. Consider a contact discontinuity situated at the surface $`r=R`$. For $`r<R`$, we have a fluid (denoted by subscript 1) which has an unperturbed velocity $`V_1`$ in the $`z`$-direction. For $`r>R`$, we have a fluid (denoted by subscript 2) which has an unperturbed velocity of $`V_2`$. We search for perturbed modes of the form $`\mathrm{\Delta }p=f(r)\mathrm{exp}[i(kz+m\varphi \omega t)]`$, imposing a regularity condition at $`r=0`$ and only permitting outgoing waves in the region $`r>R`$. Matching the pressure perturbations at $`r=R`$ gives the following dispersion relation:

$$\frac{J_m(K_1R)}{K_1J_m^{}(K_1R)}\mathrm{\Gamma }_1h_1(\omega kV_1)^2=\frac{H_m^{(1)}(K_2R)}{K_2H_m^{(1)}(K_2R)}\mathrm{\Gamma }_2h_2(\omega kV_2)^2,$$
(19)
where

$$K_1^2=\mathrm{\Gamma }_1^2\left(\frac{h_1}{\gamma _1p}1\right)(\omega kV_1)^2k^2+\omega ^2,$$
(20)
with a similar expression for $`K_2`$. Without loss of generality, we can suppose that the fluid in the region $`r>R`$ is at rest ($`V_2=0`$, $`V_1=V`$). It can be seen from this last expression with a little algebra that, in limit of small wavelengths and only mildly relativistic wave speeds, the general solution has the form

$$\omega =kf(h_1/h_2,,\gamma _1,\gamma _2)$$
(21)
where $`=V/c_{\mathrm{s},2}`$, where $`c_{\mathrm{s},1}`$ is the sound speed in the region $`r<R`$.

We can now use this general solution to comment on the connection between our simulations and the hydrodynamics of real sources. Since most of the interesting hydrodynamics is associated with the development of the KH instability, we should endeavour to approximate the expected real-life KH growth rates within our simulations.

The cocoons in real radio galaxies will, at least initially, be dominated by plasma with a relativistic equation of state (i.e. $`\gamma _1=4/3`$). On the other hand, our simulations are inherently non-relativistic (i.e., $`\gamma _1=5/3`$). While fully relativistic simulations are needed to address, in full, the effect of a relativistic equation of state, we can use the above dispersion relation to judge the effect on the KH growth rates of neglecting the relativistic nature of the cocoon plasma. Numerical solutions of eqn (19) show that the KH growth rates and wavelengths are little affected by changing from $`\gamma _1=4/3`$ to $`\gamma _1=5/3`$. Hence, neglecting the relativistic equation of state is not a serious flaw of our simulations.

Noting that the Mach number of the cocoon backflow is determined by the hydrodynamics of the system (and is always approximately unity), the remaining parameter relevant to the KH growth rate is the ratio $`h_1/h_2`$. In real systems, the ambient thermal material will be non-relativistic (with sound speed $`c_{\mathrm{s},2}`$) while the cocoon (especially near the jet head) will be dominated by relativistic plasma ($`h_1p`$). Since $`cc_{s,2}`$, the ratio of relativistic enthalpies will be

$`{\displaystyle \frac{h_1}{h_2}}`$ $`=`$ $`{\displaystyle \frac{p}{\rho _2c^2}}`$ (22)
$`=`$ $`\left({\displaystyle \frac{c_{s,1}}{c}}\right)^2.`$ (23)
Plausible parameters give $`c_{s,1}30004000\mathrm{km}\mathrm{s}^1`$, giving $`h_1/h_212\times 10^4`$. Any pollution of the relativistic plasma in the cocoon (due to mixing) will reduce this ratio appreciably. For comparison, in our simulations we have

$`{\displaystyle \frac{h_1}{h_2}}`$ $`=`$ $`{\displaystyle \frac{\rho _1}{\rho _2}}`$ (24)
$`=`$ $`\left({\displaystyle \frac{c_{\mathrm{s},2}}{c_{\mathrm{s},1}}}\right)^2,`$ (25)
and by measuring the sound speeds, we determine that we achieve $`h_1/h_2510\times 10^4`$. We conclude that while it would be desirable to push this ratio lower by a factor of a few, our simulations are exploring the appropriate parts of parameter space.