{'fraction_non_alphanumeric': 0.0295862898875472, 'fraction_numerical': 0.01192995559981742, 'mean_word_length': 4.671805130618969, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The hydrodynamics of dead radio galaxies

We present a numerical investigation of dead, or relic, radio galaxies and the environmental impact that radio galaxy activity has on the host galaxy or galaxy cluster. We perform axisymmetric hydrodynamical calculations of light, supersonic, 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 passively. The dead source undergoes an initial phase of pressure driven expansion until it achieves pressure equilibrium with its surroundings. Thereafter, buoyancy forces drive the evolution and lead to the formation of two oppositely directed plumes that float high into the galaxy/cluster atmosphere. These plumes entrain a significant amount of low entropy material from the galaxy/cluster core and lift it high into the atmosphere. An important result is that a large fraction (at least half) of the energy injected by the jet activity is thermalized in the ISM/ICM core. The whole ISM/ICM atmosphere inflates in order to regain hydrostatic equilibrium. This inflation is mediated by an approximately spherical disturbance which propagates into the atmosphere at the sound speed. The fact that such a large fraction of the injected energy is thermalized suggests that radio galaxies may have an important role in the overall energy budget of rich ISM/ICM atmospheres. In particular, they may act as a strong and highly time-dependent source of negative feedback for galaxy/cluster cooling flows.


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 ⋆ Hubble Fellow while at University of Colorado 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 ρ ∝ r −α where α > 2, the contact discontinuity will accelerate and hence will be Rayleigh-Taylor (RT) unstable. However, in the more physical case of α < 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;.
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.

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,1992bClarke, 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, θ, φ), 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 γ = 5/3. The ambient galaxy/cluster gas was assumed to be initially in an isother-mal configuration with a (adiabatic) sound speed cs = 1 and a density profile given by a β-model with β = 0.5: where the core radius was set to r0 = 2. The gravitational potential, Φ, was set so as to keep this ambient material in hydrostatic equilibrium which gives (to within an additive constant) 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 ρjet = 0.01 and Mach number (with respect to the sound speed of the injected jet material) of M = 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 r0 = 100 kpc, cISM = 1000 km s −1 and a central number density of n0 = 0.01 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 × 10 45 erg 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 r0 = 10 kpc, cISM = 500 km s −1 and n0 = 0.1 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 × 10 44 erg s −1 . Of course, it must be appreciated that these β-models are very crude parameterizations of the real ISM/ICM density structure.

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 rout = 10 at three different resolutions: nr × n θ = 300 × 300 (Run 1), 600 × 600 (Run 2), and 1200 × 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., rout = 10 with a grid size of nr × n θ = 600 × 600), and one other simulation. This final simulation (Run 4) has a larger outer radius (rout = 30) and a larger computational grid (nr × n θ = 1200 × 600). In the spatial region covered by both Run 2 and Run 4, the resolutions are approximately (but not exactly) the same.

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 2dimensional 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/ρ 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 ∼ 20 − 30 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/ √ 3).

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  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  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 Figure 1. Specific entropy maps for t = 0.50, t = 1.0, t = 3.0 and t = 5.0 for Run 2. As the same greyscale levels are used for all four panels, one can clearly see the decrease in the specific entropy of the cocoon it rises and entrains ambient thermal plasma. Note the change of scale in the first three of these plots.
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.

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 t ∼ 1.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 vc of the plume as where r bub the radius of the bubble, g is the local gravitational acceleration, and ρa and ρ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 vc ∼ 0.8. Together with the cocoon size of r bub (t = 4) ∼ 1 and the gravitational acceleration at that position, this indicates that C is slightly less than unity.

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 sepa-rate 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/ρ 5/3 as a discriminant.

The entropy threshold method and numerical mixing
For the simulation parameters we chose, the background material has specific entropy (in code units) The injected material, on the other hand, has entropy p/ρ 5/3 ∼ 1300 and cocoon material, which has gone through the terminal shock, can have entropy as high as s ∼ 10 5 . We chose to define the cocoon as the volume inside of which the entropy satisfies sc > 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 t ∼ 8, the contact discontinuity between the cocoon and ambient material is sharp and well tracked by the sc = 10 contour. In fact, until t ∼ 5, the location of the contact discontinuity is insensitive to the chosen entropy threshold, provided that threshold is in the range 10 < sc < 100. For times t > 8 numerical mixing reduces substantial parts of the cocoon/plumes below the sc = 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 t ∼ 8), the entropy structure of the cocoon may be a poor reflection of reality.

Definitions of energies
To explore the energetics of our simulation in a quantitative manner, we compute the following quantities: . Normalized volume distribution of material with given specific entropies (computed over the whole computational domain). The solid line corresponds to very early times t = 0.02, the dotted line corresponds to t = 2 and the dash-dotted line corresponds to t = 4. The second broad peak in the distribution, which corresponds to the cocoon material, gradually moves to lower entropies as it is polluted by low entropy ISM/ICM gas.
(ii) Total internal energy of cocoon (ambient) material, (iii) Total kinetic energy of cocoon (ambient) material (iv) Total gravitational potential energy of cocoon (ambient) material where C is the region in which s > 10, and A is the region in which s ≤ 10. We then reference all energies to their values at the initial time:

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 t ≈ 1.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 equilib-rium 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.

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 keV 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 β-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 keV cm 2 in a sample of 20 clusters and groups. Interestingly, these authors also conclude that this entropy injection occurred at low densities (1 − 3 × 10 −4 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.

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 r ≈ 20. 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 (η ∼ 0.5).

Cooling flows and radio galaxy heating
As discussed above, we find that a large fraction (η ∼ 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 νLν 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(Begelman , 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 coolingflows 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).

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.

CONCLUSIONS
Radio galaxies are known to be dynamical systems that evolve on timescales of 10 7 − 10 8 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 Xray 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 backto-back jets propagating in a β-model galaxy/cluster atmo-sphere. 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  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.

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 γ-law gas are where ρ and p are the comoving density and pressure, respectively, u i are the components of the four velocity, and γ is the polytropic index. Here, h = ρ + γ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: 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 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 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 V1 in the z-direction. For r > R, we have a fluid (denoted by subscript 2) which has an unperturbed velocity of V2. We search for perturbed modes of the form ∆p = f (r) exp[i(kz +mφ−ω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: where with a similar expression for K2. Without loss of generality, we can suppose that the fluid in the region r > R is at rest (V2 = 0, V1 = 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 ω = kf (h1/h2, M, γ1, γ2) where M = V /cs,2, where cs,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. γ1 = 4/3). On the other hand, our simulations are inherently non-relativistic (i.e., γ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 γ1 = 4/3 to γ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 h1/h2. In real systems, the ambient thermal material will be non-relativistic (with sound speed cs,2) while the cocoon (especially near the jet head) will be dominated by relativistic plasma (h1 ≈ p). Since c ≫ cs,2, the ratio of relativistic enthalpies will be h1 h2 = p ρ2 c 2 (22) = cs,1 c 2 .
Plausible parameters give cs,1 ∼ 3000 − 4000 km s −1 , giving h1/h2 ∼ 1 − 2 × 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 = cs,2 cs,1 2 , and by measuring the sound speeds, we determine that we achieve h1/h2 ∼ 5 − 10 × 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.