STOPPING COOLING FLOWS WITH COSMIC-RAY FEEDBACK

X-ray spectra of hot virialized gas in galaxy groups and clusters have firmly established that the rate at which gas cools to low temperatures is much lower than that expected from uninhibited radiative cooling in subsonically inflowing cooling flows (e.g., Peterson et al. 2001). As a result of this realization, many theoretical and computational studies have sought to prevent cooling inflows by introducing a gas heating mechanism usually exploiting the phenomenal accretion energy released by massive black holes in the cores of group and cluster-centered elliptical galaxies (e.g., McNamara & Nulsen 2007). In spite of this almost unlimited source of accretion energy, no specific mechanism or mode of gas heating has been accepted that is successful for long cosmic times and for cooling flows of all scales: galactic, group, and cluster.

In their recent review of the cooling flow puzzle, McNamara & Nulsen (2007; hereafter MN07) emphasize two types of active galactic nucleus (AGN)-related heating mechanisms: PdV, work done when cavities in the hot gas are inflated by cosmic rays (CRs) and energy-dissipating shock (or other) waves that deliver energy from the central energy source throughout the cluster gas. CRs, often observed inside cluster cavities from radio synchrotron emission, are thought to be deposited or produced by nonthermal jets that project out from the central AGN. However, our recent study of the energetics of cluster cavities inflated by CRs (Mathews & Brighenti 2008b, hereafter MB08) shows that the (shock) heating produced as young cavities expand due to locally enhanced CRs is more than compensated by a global cooling as the entire cluster gas expands more or less permanently to accommodate the partial pressure of the CRs as they diffuse into the cluster gas. Most of the PdV work expended as cavities form is stored in potential energy of displaced gas that moves out in the cluster potential. The dissipative heating beneath buoyant cavities discussed in MN07 appears to be absent in the cavities computed in MB08, and this is consistent with our earlier estimate (Mathews et al. 2006). The total (bulk and turbulent) kinetic energy is always very small and cannot be an important source of dissipative heating (MB08).

The second currently favored heating mechanism, large-scale wave dissipation, is also problematic. Since the hot gas density profiles in groups and clusters is generally flatter than ρ ∝ r⁻², successive shocks preferentially heat gas near the cluster center until it exceeds observed temperatures (Mathews et al. 2006). Overheating the central gas is a persistent difficulty with previous theoretical attempts to resolve the cooling flow problem by using central AGN heating scenarios (Brighenti & Mathews 2002b, 2003) since the cluster gas closest to the source of AGN energy typically has the lowest temperatures observed in the hot cluster gas. Finely tuned heating scenarios in which the gas inflow is assumed to be essentially stopped are inconsistent with secular changes in the gas that evolve without limit. For example metal enrichment by Type Ia supernovae in the central galaxy would eventually increase the local hot gas metal abundances far above those observed.

As an alternative solution to the cooling flow problem, we proposed that gas receiving AGN energy is buoyantly transported outward. In these mass-circulating flows gas flows in both radial directions simultaneously. The combined flows maintain a cooling flow appearance similar to observed X-ray images (Mathews et al. 2003, 2004; Brighenti & Mathews 2006).

In this Letter we describe the evolution of virialized hot gas in the Virgo cluster in which successive cavities are formed with CRs that buoyantly transfer gas from near the cluster core to large radii where its radiative cooling time is very much longer. The average net gas mass that flows past every cluster radius is very small although significant inflows and outflows occur locally in the gas. In buoyant outflows generated by CRs, the cluster gas need not be directly heated by the CRs, although some modest heating results from the dissipation of weak shocks that propagate away from newly formed cavities (Mathews & Brighenti 2007; MB08). Unlike most AGN heating scenarios, buoyant mass-circulating flows generated by CRs can preserve the observed cluster temperature and density profiles for at least several Gyr. This resolution of the cooling flow problem, using mass transfer by CR buoyancy rather than heating, is attractive because almost all group and cluster-centered galaxies are observed to produce CRs in jets, radio lobes, or in central nonthermal radio sources.

The procedure and assumptions employed here are nearly identical to those described in detail by MB08 to which the reader is referred. The flow equations are

$$\begin{equation} { \partial \rho \over \partial t} + {\bf \nabla }\cdot \rho {\bf u} = 0 \end{equation} \tag{ 1 }$$

$$\begin{equation} \rho \left({ \partial {\bf u} \over \partial t} + ({\bf u \cdot \nabla }){\bf u}\right) = - {\bf \nabla }(P + P_c) - \rho {\bf g} \end{equation} \tag{ 2 }$$

$$\begin{equation} {\partial e \over \partial t} + {\bf \nabla \cdot u}e = - P({\bf \nabla \cdot u}) - (\rho /m_p)^2 \Lambda (T,z) \end{equation} \tag{ 3 }$$

$$\begin{equation} {\partial e_c \over \partial t} + {\bf \nabla \cdot u}e_c = - P_c({\bf \nabla \cdot u}) + {\bf \nabla \cdot }(\kappa {\bf \nabla }e_c) + {\dot{S}}_c, \end{equation} \tag{ 4 }$$

where artificial viscosity terms are suppressed. Gas and CR pressures are related to their corresponding energy densities by P = (γ − 1)e and P_c = (γ_c − 1)e_c respectively where for simplicity we assume γ = 5/3 and γ_c = 4/3 for fully relativistic CRs. The energy spectrum of the CR particles N(E) is not considered in detail so the CRs are characterized only by their integrated energy density e_c ∝ ∫EN(E)dE. Furthermore, the CR particles are unspecified and can consist of electrons or protons in any combination.

Unlike the discussion in MB08, we now include radiative cooling. We adopt the radiative cooling coefficient n_ionsn_eΛ_sd(T, z) erg cm⁻³ s⁻¹ from Sutherland & Dopita (1998) which, as expressed in Equation (3), (ρ/m_p)²Λ erg cm⁻³ s⁻¹, requires that Λ = 1.1[(4 − 3μ)(2 + μ)/25μ²]Λ_sd where μ = 0.61 is the molecular weight and the coefficient is n_ions/n_p = 1.1. We do not compute the metallicity in the gas throughout the flow but simply adopt a uniform metallicity near solar z = 0.75 z_☉ since most of the gas that cools has been metal-enriched in the cluster core. The computed flows are not sensitive to the gas metallicity. Computational time steps are required to be less than the constant-pressure radiative cooling time in all zones. Since some cooling occurs in off-center zones near the z-axis, it is useful to remove cooled gas from the computational grid. If a grid zone cools below T_min = 5 × 10⁵ K, gas is removed at constant pressure by reducing the zone density until the temperature resets to T_max = 5 × 10⁶ K. Gas removal at constant pressure from a zone preserves the pressure gradients relative to nearby zones so the gas flow to or from neighboring zones is unaffected. Moreover, the total thermal energy in the zone is also preserved so thoroughly cooled gas removed in this way does not alter the total thermal energy of uncooled gas. We assume that CRs are not removed with the cooled gas.

Gas is accelerated by pressure gradients in both the thermal gas and CRs. The two fluids are assumed to be coupled by the presence of small magnetic fields in the cluster gas. This coupling can be efficient even when the energy density in the magnetic field is much less than that in the thermal gas, which is consistent with the small fields observed in clusters, 1–10 μG (Govoni & Feretti 2004). Consequently, magnetic pressure and stresses are not explicitly considered.

Equation (4) above describes the advection of CRs by locally flowing gas and the diffusion of CRs through this gas. Little is known about the diffusion coefficient κ, but it is likely that it varies inversely with the gas density since the magnetic field is probably larger in denser gas. As in MB08, we consider a range of CR diffusion coefficients that depend on a single density parameter n_e0:

$$\begin{eqnarray} \kappa = \left\lbrace \begin{array} {@{}ll@{}} 10^{30} \;{\rm cm}^2 \;{\rm s}^{-1}l &:\quad n_e \le n_{e0}\; {\rm cm}^{-3} \nonumber\\ 10^{30}(n_{e0}/n_e)\; {\rm cm}^2\; {\rm s}^{-1}l &:\quad n_e > n_{e0} \;{\rm cm}^{-3}.\nonumber\\ \end{array} \right. \end{eqnarray}$$

Here we consider n_e0 = 6 × 10⁻³ and 6 × 10⁻⁶ cm⁻³. The cluster gas dynamics are not very sensitive to n_e0 or to other κ(n_e) variations that have been explored (MB08).

The flow equations are solved in two-dimensional cylindrical coordinates (R, z) using a ZEUS-like code (Stone & Norman 1992) but with the additional CR equation. As in MB08 the computational grid consists of 100 equally spaced zones out to 50 kpc plus an additional 100 zones in both coordinates that increase logarithmically out to ∼1 Mpc, filling a large hemisphere. Also following MB08, a spherical gravitational potential for the Virgo cluster is found by assuming that the gas is in hydrostatic equilibrium with the temperature and modified density profiles T_G(r) and ρ_mG(r) observed by Ghizzardi et al. (2004). The density profile was found by integrating the equation of hydrostatic equilibrium using T_G(r) and the total gravitational acceleration from a 3.2 × 10⁹ M_☉ black hole, the M87 stellar mass ρ_*(r) from Brighenti & Mathews (2002a), and an NFW profile with M_vir = 1.30 × 10¹⁴M_☉ and concentration 8.60. The resulting modified observed density profile ρ_mG(r) agrees very well with that of Ghizzardi et al. ρ_G(r) for ≳1 kpc but avoids an unphysical maximum in g(r) near ∼4 kpc if g(r) is found from T_G(r) and ρ_G(r). All flow calculations begin with T_G(r) and ρ_mG(r).

CRs are introduced at intervals of t_cyc = 2 × 10⁸ yr in bursts lasting t_cav = 2 × 10⁷ yr. These timescales are roughly consistent with the lifetimes of visible X-ray cavities and radio synchrotron electrons. If the CR luminosity during injection times is ${\dot{E}}_c$, the time-averaged CR luminosity is

$$\begin{equation} \langle L_{\rm cr} \rangle = 2(t_{\rm cav} / t_{\rm cycle}) {\dot{E}}_c, \end{equation} \tag{ 5 }$$

assuming that CRs are injected into both cluster hemispheres by symmetric double jets. For comparison with the cavity calculations in MB08, all cavities are assumed to form in a source region 10 kpc along the z−axis described by a normalized Gaussian profile. The CR source term in Equation (4) is therefore

$$\begin{equation} {\dot{S}}_c = {\dot{E}}_c \tau (t) {e^{-(({\bf r - r}_{\rm cav})/r_s)^2} \over \pi ^{3/2} {r_s}^3} {\rm erg}\; {\rm cm}^{-3}\;{\rm s}^{-1}, \end{equation} \tag{ 6 }$$

where r is the distance from the origin, r_s = 2 kpc and τ(t) is unity during injection times and zero otherwise. All calculations described here continue for 3 Gyr. For each assumed κ(n_e0) we vary ${\dot{E}}_c$, seeking that value for which the gas temperature and density profiles continue to agree best with profiles observed in Virgo. Except for radiative cooling and a variable ${\dot{E}}_c$, the parameters chosen here are identical to those discussed for single cavities by MB08.

The left column of Figure 1 illustrates the evolution of Virgo cluster gas as a pure radiatively cooling flow without CRs. The long-dashed lines show the radial gas density, temperature and pressure profiles after evolving for 3 Gyr from the observed profiles. The evolved profiles are determined by ordering ρ, T and P in each computational zone by increasing distance from the origin r = (R² + z²)^1/2, then averaging in appropriate bins Δr. After 3 Gyr, ρ(r), and T(r) have evolved away from the initial (observed) profiles (solid lines). The ratio of gas entropy to its initial value (subscript 0), log(S/S₀) = log(T/T₀) − (2/3)log(ρ/ρ₀), decreases monotonically toward the cluster center where the gas cools in the central 1.5 kpc. The discrepancy in integrated mass of (uncooled) gas between the initial and evolved density profiles in the uppermost panel, ΔM = M(r, 3Gyr) − M(r, 0), peaks at ∼15 kpc.¹ Between 2 and 3 Gyr the quasi-steady cooling rate in the core is $\langle {\dot{M}} \rangle = 86$ M_☉ yr⁻¹ for the entire Virgo cluster (both hemispheres)²

Zoom In Zoom Out Reset image size

Panels in the central column of Figure 1 show the corresponding evolution with CRs having a low CR diffusion coefficient κ with n_e0 = 6 × 10⁻⁶ cm⁻³. In this solution after 3 Gyr the gas density and temperature remain close to the initial, observed profiles. The entropy in 25 ≲ r ≲ 100 kpc is lowered by the accumulated buoyant arrival of gas from the cluster core (cf. MB08) which has lost entropy by radiation. The modest entropy maximum near 10 kpc is probably due to heating by successive weak cavity shocks at this radius. The total mass of hot gas in each spherical annulus changes very little during the evolution, as expected from the small central mass cooling rate between 2 and 3 Gyr, $\langle {\dot{M}} \rangle = 0.13$ M_☉ yr⁻¹ (for both hemispheres). The cooling inflow has been almost completely balanced by the outward circulation of gas from the cluster core to large radii driven by CR buoyancy. The bottom panel shows that the CR pressure P_c after 3 Gyr (short dashed line) is about an order of magnitude less than the gas pressure within 100 kpc, but P_c/P decreases sharply beyond 100 kpc. In this solution the total (spherical) time-averaged CR luminosity is 〈L_cr〉 = 2.70 × 10⁴³ erg s⁻¹. The total CR energy delivered to each cavity is about E_cav = 8.2 × 10⁵⁸ erg.

The right column of panels in Figure 1 shows the identical evolution but with a much larger CR diffusion coefficient (n_e0 = 6 × 10⁻³ cm⁻³). As before, the gas density and temperature profiles remain essentially unchanged after 3 Gyr with a total time-averaged CR luminosity 〈L_cr〉 = 2.60 × 10⁴³ erg s⁻¹. The total gas cooling rate between 2 and 3 Gyr, $\langle {\dot{M}} \rangle = 1.15$ M_☉ yr⁻¹, is only ∼0.01 of the cooling flow rate without CRs.

CR buoyancy efficiently arrests cooling inflows without requiring ad hoc heating or infusions of hot, nonrelativistic gas in the cavities. In the two solutions with CRs gas cools at the center and along the z-axis out to about 3–5 kpc where the free fall time to the center is ≲10⁷ yr. If all this cooled gas were accreted by the central black hole, the efficiency required for accretion energy to produce enough CRs to maintain the observed gas profiles is $\eta _{\rm cr} = \langle L_{\rm cr} \rangle /\langle {\dot{M}} \rangle c^2$. We find η_cr = 3.6 × 10⁻³ and 3.9 × 10⁻⁴ respectively for the low and high κ solutions. But η_cr is probably underestimated since some of the cooled, dusty gas in group-centered elliptical galaxies is heated and buoyantly relocated far out into the hot gas (Temi et al. 2007). The total CR luminosity required to (nearly) stop the Virgo cooling flow, 〈L_cr〉 = 2.60 × 10⁴³ erg s⁻¹, is equal to the (computed) bolometric X-ray luminosity within only ∼56 kpc.

The Fermi Gamma Ray Telescope may confirm that cooling inflows in Virgo and other nearby clusters are stopped with CR buoyancy. Gamma rays are created by neutral pion decays following inelastic p–p interactions between CRs and the thermal plasma. The integrated gamma ray flux (cm⁻² s⁻¹) above energy E_γ,min is

$$\begin{equation} F_{\gamma }(E_{\gamma } > E_{\gamma,{\rm min}}) = {1 \over 4 \pi d^2} \int _{\rm Virgo} dV \int _{E_{\gamma,{\rm min}}}^{\infty } dE_{\gamma } q_{\gamma } (E_{\gamma },r), \end{equation} \tag{ 7 }$$

where d = 16 Mpc is the distance to Virgo and the specific photon emissivity, q_γ(E_γ, r) cm⁻³ s⁻¹ GeV⁻¹, is described by Pfrommer & Ensslin (2004) and Ando & Nagai (2008). Assume for simplicity that all CRs are relativistic protons having a power-law momentum distribution f(p_p, r) ∝ p^−α_p with α = 2.3. Then, after 3 Gyr,the integrated gamma ray photon flux from Virgo above E_γ,min = 2 Gev is F_γ = 2 × 10⁻⁹ or 3 × 10⁻⁹ cm⁻² s⁻¹ respectively for the low and high κ solutions in Figure 1. These fluxes are near the detection threshold of the Fermi telescope, ∼10⁻⁹ cm⁻² s⁻¹ estimated by Ando & Nagai (2008) for extended sources. If half of the CRs are electrons, evidence of the AGN CRs that we require in Virgo may still be marginally visible with the Fermi Telescope. Additional CRs in Virgo may have been produced in shocks and turbulence accompanying mergers or in the accretion shock when baryons first fell into the Virgo potential, but these are not considered here.

For both low and high κ solutions in Figure 1 the CR partial pressure P_c/(P + P_c) increases toward the cluster center. This may seem surprising because all CRs are injected at 10 kpc and diffusion into the dense, central gas is slow because κ decreases. This feature can be understood from the reaction of X-ray cavities on the cluster gas. As a young cavity buoyantly moves away from the cluster center, relatively cooler nearby gas concentrates underneath the cavity, then flows rapidly outward through the cavity along the z-axis (MB08). These jets or filaments of denser thermal gas move out beyond 10 kpc then later fall back under gravity into the cluster core (e.g., Mathews & Brighenti 2008a). It is likely that the central concentration of CRs occurs because they are advected in with each infalling filament, then become trapped by the relatively small κ in the dense central gas. CRs are a likely source of nonthermal pressure that is occasionally observed in the cores of hot galactic and cluster gas.

Figure 2 shows that the distribution of CRs in Virgo after 3 Gyr is concentrated in the core and along the z-axis. In a recent study of warm and cold gaseous filaments in the Perseus Cluster, Ferland et al. (2008, 2009) found that the observed spectrum can be understood if some of the emission lines are excited by collisions with CRs. The enhanced CR distribution in Figure 2 along the z-axis provides a natural source of CRs just where the cool radial filaments form. Although we do not compute the passively evolving magnetic field, it is expected to be stronger in denser thermal filaments and longitudinally oriented along the cool filaments due to the radial expansion of the filament flow along the z-axis.

Zoom In Zoom Out Reset image size

The two-dimensional cylindrical nature of our solution requires that all cavities and CR activity occur along the z-axis. The relatively small number of double–double radio sources (Schoenmakers et al. 2000) indicates that successive CR jets can move along the same direction, as assumed here. Remarkably, the X-ray thermal filament in Virgo (Forman et al. 2007) that was formed by a cavity ∼10⁸ years ago (Mathews & Brighenti 2008a) and its associated radio lobe lies in a different direction than the younger nonthermal jet in the central galaxy M87. Another indication that CR energy flows out in many different directions is the azimuthal variation of successive generations of X-ray cavities in Perseus and other clusters.

It is expected that the mass outflow due to CR buoyancy in Figure 1 also transports iron produced by Type Ia supernova inside the central galaxy into a quasi-spherical region 50–100 kpc in radius as observed in many groups and clusters (De Grandi et al. 2004). The outflows described here can also carry short-lived dust out to 5–10 kpc into the hot gas where its far-infrared emission is observed in many nearby galaxy groups (Temi et al. 2007). Radial entropy variations qualitatively similar to those in Figure 1 have been observed in galaxy groups (Gastaldello et al. 2007; Sun et al. 2009). Finally, the outward mass transport can create an inflection in the density slope as in Figure 1, resembling the double-structured character of many density profiles observed.

Studies of the evolution of hot cluster gas at UC Santa Cruz are supported by NASA and NSF grants for which we are very grateful.