TOWARD A COMPREHENSIVE MODEL FOR FEEDBACK BY ACTIVE GALACTIC NUCLEI: NEW INSIGHTS FROM M87 OBSERVATIONS BY LOFAR, FERMI, AND H.E.S.S.

The characteristic frequency ν_s at which a CR electron (CRe) emits synchrotron radiation scales with the magnetic field strength B and the CRes' cooling time τ as ν_s∝B⁻³τ⁻² (assuming the strong-field limit, B > 3.2 (1 + z)² μG, for which inverse Compton (IC) cooling is negligible). Hence, observations at lower radio frequencies should enable us to probe older fossil CRe populations. Those are expected to accumulate at lower energies, which would imply a low-frequency spectral steepening. Equivalently, such observations should unveil the time-integrated non-thermal action of AGN feedback that is expected to be spread out over larger spatial scales. However, recent LOFAR observations of Virgo A (de Gasperin et al. 2012; see also Figure 1), the radio source of M87, provided two unexpected findings. (1) The observations do not reveal a previously hidden spectral steepening toward low frequencies and show the extended radio halo to be well confined within boundaries that are identical at all frequencies (de Gasperin et al. 2012). The possibility that this may point to a special time, which is characterized by a recent AGN outburst after long silence, is unlikely because statistical studies of the radio-mode feedback in complete samples of clusters imply that the duty cycle of AGN outbursts with the potential to heat the gas significantly in CC clusters is ⩾60% and could approach unity after accounting for selection effects (Bîrzan et al. 2012). (2) While the core region displays a radio spectral index² of α_ν = 0.6 ± 0.02 (consistent with an injection CR spectral index of α = 2α_ν + 1 = 2.2 ± 0.04), the low-frequency halo spectrum outside the central region (of diameter ∼2 arcmin) is considerably steeper with α_ν = 0.85. As suggested by de Gasperin et al. (2012), the most probable explanation of this steepening at low frequencies is strong adiabatic expansion in combination with a superposition of differently synchrotron-aged CRe populations. Adiabatic expansion across a volume expansion factor f = V₂/V₁ = ρ₁/ρ₂ > 1 from the compact source to the more dilute halo regions shifts all CR momenta by the factor p₂ = f^−1/3p₁. CRes that are leaving the lobes have spent different time intervals in regions of strong magnetic fields, which leads to a distribution of spectra with different break momenta. Superposition of those spectra steepens the spectrum down to low frequencies if the cooling time of the most strongly cooled CRes has been comparable to the transport time from the core to the current emitting region.

Zoom In Zoom Out Reset image size

In Figure 2, we show the cooling time of CRes (see Pfrommer et al. 2008 for definitions). CRes of energy ∼200 (400) MeV that are responsible for ∼25 (100) MHz radio emission in a 20 μG field have energy loss timescales of ∼150 (75) Myr, which drops further by the factor f^1/3 after accounting for adiabatic expansion. CRes can cool through IC interactions with various seed photon fields, including the cosmic microwave background, starlight, and dust-reprocessed infrared emission. Within the central radio halo region (r < 34 kpc), the star-formation-induced photon fields at least contribute as much as the cosmic microwave background (Pinzke et al. 2011), which motivates our choice of a seed photon energy density that is twice that of the cosmic microwave background. Synchrotron cooling in equipartition magnetic fields of 10–20 μG (de Gasperin et al. 2012) and adiabatic expansion substantially decreases the cooling times further during the advective transport of CRes in the radio lobes. Those rise in the cluster atmosphere (either buoyantly or driven by the jet) and uplift dense gas from the central regions (Werner et al. 2010). Those rising motions of light relativistic fluid cause strong downdrafts at the sides of the lobes, which induce Rayleigh–Taylor and Kelvin–Helmholtz instabilities that disrupt the lobes and excite turbulence, which mixes the CRes with the dense uplifted ambient plasma. This picture appears to be supported by radio maps of M87 (Owen et al. 2000; de Gasperin et al. 2012) that show a diffusely glowing radio halo surrounding the lobes, which themselves show the characteristic mushroom shape of a Rayleigh–Taylor instability (see Figure 1). This radio map demonstrates that the radio halo emission fills the entire CC region (in projection), indicating that CRs are distributed over 4π steradian by means of turbulent advection and CR streaming. The comparably sharp halo boundary, which exactly coincides with the emission region at 1.4 GHz (Owen et al. 2000), suggests a magnetic confinement of CRs to the core region. Coulomb interactions of the released CRes with the dense ambient ICM cause them to be thermalized after a maximum cooling time of τ ≃ 40 Myr, which corresponds to the radio halo age (see Figure 2). Thus, this naturally explains the absence of a fossil CRe signature in the low-frequency radio spectra, as well as the spatially similar morphology of Virgo A at all frequencies.

Zoom In Zoom Out Reset image size

Diffusive shock acceleration operates identically on relativistic particles of the same rigidity ($R=\tilde{p}c/z e$, where $\tilde{p}$ is the physical momentum and z is the ion charge), and unlike electrons, protons do not suffer an injection problem into the acceleration process and have negligible radiative losses. Hence, we expect CR protons (CRps) to dominate the pressure of the lobes³ and to exhibit a similar injection spectrum as the CRes in the central regions (α ≈ 2.2). Similarly to the CRes, CRps should also be released from the lobes and mixed into the ambient ICM. Since adiabatic expansion leaves their spectral shape invariant and CRps have negligible radiative losses, we do not expect a spectral steepening of the CRp population in the halo. A necessary consequence of this picture are hadronic interactions of CRps with the ambient gas. This should produce a steady gamma-ray signal from M87 (Pfrommer & Enßlin 2003), with a spectral index above GeV energies of α_γ = α ≈ 2.2.

To test this hypothesis, we turn to the gamma-ray regime. A joined unbinned likelihood fit with a power law (${\it dN}/{\it dE}\propto E^{-\alpha _\gamma }$) to the Fermi and H.E.S.S. data of M87's low state yields F(> 1 GeV) = 1.35(± 0.35) × 10⁻⁹ photons cm⁻² s⁻¹ and α_γ = 2.26 ± 0.13 (Abdo et al. 2009), which is consistent with our expectation if the gamma rays are produced in hadronic CRp–p interactions. Moreover, a comparison of the Fermi light curve for the last 14 months with the initial 10 month data revealed no evidence of variability in the flux above 100 MeV (Abramowski et al. 2012; Abdo et al. 2009).⁴ Gamma-ray upper limits on M87 set by EGRET (Reimer et al. 2003) are consistent with the steady flux seen by Fermi. This demonstrates that the EGRET and Fermi observations of M87, which together span more than two decades, show no apparent increase in flux. Additionally, the low flux state at very high energies E > 100 GeV does not provide any hint of variability (Aharonian et al. 2006), although the low signal statistics renders such a test very difficult. Encouraged by the spectral similarity of radio and gamma-ray observations and the constant gamma-ray flux, we will now assume that the gamma rays are tracing hadronic CR interactions.

To estimate the CR pressure that is implied by the observed gamma-ray emission, we choose a constant CR-to-thermal pressure ratio, X_CR = P_CR/P_th (which we require for self-consistency, as we will see in Section 3). We assume an isotropic one-dimensional CR distribution function in momentum space,⁵

$$\begin{equation} \frac{d^4N}{dp\,d^3x}=C p^{-\alpha }\theta (p-q), \end{equation} \tag{ 1 }$$

where θ denotes the Heaviside step function, C denotes the normalization, $p={\tilde{p}}/mc$ and q denote the normalized momentum and low-momentum cutoff, and α = 2.26. We decided in favor of the spectral index of the gamma-ray emission, which likely probes the prevailing CRp distribution in the radio halo region. In contrast, the radio emission of the central region may signal current acceleration conditions that, albeit generally very similar, may be slightly different from the ones that have accelerated the CRs in the past, which have now been released.

We use a fit to the spherically symmetric thermal pressure profile as inferred from X-ray observations (Matsushita et al. 2002; Böhringer et al. 1994). The thermal pressure is P_th = nkT = n_ekT/μ_e, where μ_e = n_e/n = μ (X_H + 1)/2 ≃ 0.517 for a fully ionized gas of primordial element composition with a hydrogen mass fraction X_H = 0.76 and mean molecular weight μ = 0.588. The electron density and temperature profiles⁶ of the central region of the Virgo cluster (r ≲ 100 kpc corresponding to 20') are given by fits to the X-ray profiles and have the form

$$\begin{eqnarray} n_e (r) &=& \left[ \sum _{i=1}^2 n_i^2 \, \left(1+\frac{r^2}{r^2_i} \right)^{-3 \beta _i} \right]^{1/2}, \end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray} T_e (r) &=& T_0 + (T_1 - T_0)\, \left[ 1+\left(\frac{r}{r_t} \right)^{-1}\right]^{-1}, \end{eqnarray} \tag{ 3 }$$

where (n₁, n₂) = (0.15, 0.013) cm⁻³, (r₁, r₂) = (1.6, 20) kpc, (β₁, β₂) = (0.42, 0.47), (kT₀, kT₁) = (1, 3) keV, and r_t = 13.5 kpc.

The sharp drop of radio morphology of Virgo A's halo suggests that the CRps are also confined to a central region of radius r_max = 34 kpc. Hence, we integrate the gamma-ray flux due to hadronic CRp–p interactions from within a spherical region bounded by r_max. We use a semi-analytic formalism (Pfrommer & Enßlin 2004) that has been modified to constrain the CR-to-thermal pressure ratio (Pfrommer et al. 2008). We obtain X_CR = 0.38 for a vanishing low-momentum cutoff, q = 0. Accounting for a non-zero low-momentum cutoff of the CRp distribution as a result of Coulomb cooling, we obtain X_CR = 0.31 (0.34) for q = 0.58 (0.25), which corresponds to a Coulomb cooling time of CRps of τ = 40 Myr, assuming an ambient electron density of n_e = 10⁻¹ (10⁻²) cm⁻³. Since there is considerable uncertainty about the shape of the low-momentum regime of the CR distribution function, we decided to adopt the most conservative value of X_CR = 0.31 for the remainder of this work. We note that a fresh supply of CRps could lower q and therefore increase the CR pressure.

We checked that the radio emission from secondary CRes, which are inevitably injected by hadronic CRp–p reactions, has a subdominant contribution to the observed radio emission. While the total radio emission due to secondaries is well below the observed flux density of the radio halo of M87, future high-frequency observations at ≳ 10 GHz of the external halo parts may be able to probe the hadronically induced radio emission (see Section 4.2 for more details). At these frequencies, the hadronically induced emission component cannot any more be neglected and complicates inferences on the CRe pitch-angle distribution in different spectral aging models of radio-emitting CRes.

In general, cluster merger and accretion shocks are also expected to accelerate CRs, which add to the CR population injected by the central AGN. Non-observations of gamma-ray emission at GeV and TeV energies limit their contribution to X_CR < 0.017 in the Coma and Perseus clusters (Arlen et al. 2012; Aleksić et al. 2012) and typically to less than a few percent for the next best targets (Ackermann et al. 2010; Pinzke et al. 2011). A joint likelihood analysis searching for spatially extended gamma-ray emission at the locations of 50 galaxy clusters in 4 yr of Fermi-LAT data does not reveal any CR signal and limits the CR-to-thermal pressure ratio to be below 0.012–0.014 depending on the morphological classification (Ackermann et al. 2013). We hence neglect their contribution in the center, where the AGN-injected population apparently dominates.

Our findings for X_CR are in agreement with constraints on the non-thermal pressure contribution in M87 that have been derived by comparing the gravitational potentials inferred from stellar and globular cluster kinematics and from assuming hydrostatic equilibrium of the X-ray emitting gas (Churazov et al. 2010). Depending on the adopted estimator of the optical velocity dispersion, these authors find a non-thermal pressure bias of X_nt = P_nt/P_th ≈ 0.21–0.29 (which is even higher in M87 for more sophisticated estimates that are less sensitive to the unknown orbit anisotropy). There are several effects that could (moderately) increase X_nt. (1) A more extended distribution of P_nt in comparison to P_th could hide a larger non-thermal pressure contribution (since only the pressure gradient enters into the equation for hydrostatic equilibrium), (2) the presence of an X-ray "invisible" cool gas phase that is coupled to the hot phase could bias X_nt low (Churazov et al. 2008), and (3) the tension in the orbit anisotropy estimates of the stars and globular clusters may call for triaxial potential models (Romanowsky & Kochanek 2001).

So far, we implicitly assumed that the pressure of radio lobes is provided by a relativistic component. While it is unknown what physical component provides the dominant energetic contribution to radio lobes, here we will review the reasons for our assumption. Equipartition arguments for the radio-emitting lobes show that the sum of CRes and magnetic fields can only account for a pressure fraction of ≃ 10% in comparison to the external pressure, with which the lobes are apparently in approximate hydrostatic equilibrium (Blanton et al. 2003; de Gasperin et al. 2012). However, these estimates rely on power-law extrapolations from the high-energy part of the CRe distribution (with Lorentz factor γ ≳ 10³, which was visible so far at radio frequencies ν ≳ 300 MHz) to lower energies that dominate the CRe energy budget, leaving the possibility of an aged fossil CRe population that would dominate the pressure. LOFAR observations have now unambiguously demonstrated the absence of such a hypothetical fossil CRe population in M87 (de Gasperin et al. 2012). Since CRps should be accelerated at least as efficiently as the visible CRes and since CRps have negligible radiative cooling rates (in comparison to CRes), we make the very plausible assumption in this work that CRps provide the dominating pressure component of lobes (and work out a few smoking-gun tests for this scenario in Section 4). Such a CRp population could originate from entrained material as a result of the interaction of a Poynting-flux-dominated jet with the surrounding plasma. This entrained material can get accelerated by processes such as internal shocks or Kelvin–Helmholtz turbulence to give rise to a spine-layer jet geometry at larger scales that is needed to explain the multi-frequency data of some AGNs (e.g., in NGC 1275; Aleksić et al. 2010).

The alternative scenario of a hot but subrelativistic (Maxwellian) proton population is unlikely because of three reasons. (1) CRp Coulomb cooling timescales are too short (see, e.g., Enßlin et al. 2007), implying a burst-like heating once the content has been set free as a result of disrupting the radio lobes by the action of Kelvin–Helmholtz and Rayleigh–Taylor instabilities. However, there is no signature in the X-ray observations of such a localized heating around the disrupted lobes (Million et al. 2010). (2) Moreover, LOFAR observations show no sign of such a ∼100 MeV electron distribution, requiring a helicity fine tuning for the acceleration process. (3) Adiabatic expansion of a composite of relativistic and non-relativistic (thermal) populations favors the CR over the non-relativistic population since P_CR/P_th∝f^1/3, due to the softer equation of state of CRs. Here we assumed an ultra-relativistic CR population with an adiabatic exponent of γ = 4/3 and f = V₂/V₁ = ρ₁/ρ₂ is the volume expansion factor starting from the densities at which the two populations got energized.

Advection of CRs by the turbulent gas motions produces centrally enhanced profiles. In contrast, CR streaming along magnetic field lines implies a net outward flux of CRs that causes the profiles to flatten (Enßlin et al. 2011). A necessary requirement for CR streaming to work is a centrally peaked profile of the CR number density. Do we have any evidence for this in M87? Radio maps of Virgo A (Owen et al. 2000; de Gasperin et al. 2012) suggest that turbulent advection, driven by the rising lobes, is strong enough to maintain a sufficiently strong inward transport of CRs. Alternatively, the smoothly glowing radio halo that fills the central region for r < r_max may also be caused by escaping CRs from the lobes already at the early stages of their rise. Following observational and theoretical arguments, the current distribution of the CR pressure appears to follow a centrally peaked profile.

The CR Alfvén wave heating term is given by (Wentzel 1971; Guo & Oh 2008)

$$\begin{equation} \mathcal {H}_\mathrm{CR}= -{\boldsymbol{\upsilon}} _A {{{\cdot }}} \boldsymbol{\nabla } P_\mathrm{CR}= -{\boldsymbol{\upsilon}} _A {{{\cdot }}} \left(\boldsymbol { {e}}_r\nabla _r \langle P_\mathrm{CR}\rangle _\Omega + \boldsymbol { {e}}_l\frac{\delta P_\mathrm{CR}}{\delta l}\right), \end{equation} \tag{ 4 }$$

where ${\boldsymbol{\upsilon}} _A = \boldsymbol { {B}}/\sqrt{4 \pi \rho }$ is the Alfvén velocity, ρ is the gas mass density, and we decompose the local CR pressure gradient into the radial gradient of the azimuthally averaged CR pressure profile, $\langle P_\mathrm{CR}\rangle _\Omega$, and a fluctuating term, $\boldsymbol { {e}}_l \delta P_\mathrm{CR}/\delta l \equiv \boldsymbol{\nabla } P_\mathrm{CR}- \boldsymbol { {e}}_r \nabla _r \langle P_\mathrm{CR}\rangle _\Omega$, which in general also points in a non-radial direction, defined by the unit vector $\boldsymbol { {e}}_l$. In steady state, CR streaming transfers an amount of energy to the gas per unit volume of

$$\begin{equation} \Delta \varepsilon _\mathrm{th} = -\tau _A {\boldsymbol{\upsilon}} _A {{{\cdot }}} \boldsymbol{\nabla } P_\mathrm{CR}\approx P_\mathrm{CR}= X_\mathrm{CR}P_\mathrm{th}, \end{equation} \tag{ 5 }$$

where we used the Alfvén crossing time $\tau _A = \delta l/\upsilon _A$ over the CR pressure gradient length, δl. Comparing the first and last terms of this equation suggests that a constant CR-to-thermal pressure ratio is a necessary condition if CR streaming is the dominant heating process, i.e., the thermal pressure profile adjusts to that of the streaming CRs in equilibrium so that we will take the radial profile of P_th as a proxy for P_CR.

The observed radial profile of P_th averages over the local pressure inhomogeneities within a given radial shell, especially toward the center where the azimuthally averaged radial profile plateaus out. However, such a constant P_th at the center is in strong disagreement with the fluctuating thermal pressure maps (see Figure 7 of Million et al. 2010). Those pressure fluctuations could have their origin in shocks, turbulence, perturbations of the gravitational potential, or entropy variations of advected gas that has not had time to come into pressure equilibrium with the surroundings. To illustrate one of these possibilities, we consider pressure fluctuations due to weak shocks and expand the Rankine–Hugoniot jump conditions in the limit of small Mach number, $\mathcal {M}_1 - 1 = \varepsilon \ll 1$:

$$\begin{equation} \frac{P_2}{P_1} = \frac{2 \gamma _\mathrm{th} \mathcal {M}_1^2 - (\gamma _\mathrm{th}-1)}{\gamma _\mathrm{th}+1} \approx 1 + \frac{4 \gamma _\mathrm{th}}{\gamma _\mathrm{th}+1}\,\varepsilon = 1 + \frac{5}{2}\,\varepsilon , \end{equation} \tag{ 6 }$$

where subscripts 1 and 2 denote upstream and downstream quantities in the shock rest frame and we assume γ_th = 5/3 in the third step. Adopting a constant CR-to-thermal pressure ratio, we obtain an estimate for the CR heating due to pressure fluctuations as a result of a sequence of weak shocks, i.e., the second term on the right-hand side of Equation (4),

$$\begin{equation} \frac{\delta P_\mathrm{CR}}{\delta l} \approx -\frac{5}{2}\, \frac{\varepsilon \, X_\mathrm{CR}P_\mathrm{th}}{r}. \end{equation} \tag{ 7 }$$

Here we assume that the characteristic size scale of the pressure fluctuations is of order the cluster-centric radius, which can be realized, e.g., for spherically symmetric shocks whose curvature radius increases as r. We note that some of the other mechanisms that may be responsible for the observed pressure fluctuations can in principle have a different scaling with radius and hence have a modified radial dependence of the fluctuating term in $\mathcal {H}_\mathrm{CR}$, in particular at the center. It will be subject to future work to quantify the contribution of those processes to the CR heating rate.

We choose conservative values for the magnetic field strengths that are about half as much as the minimum energy estimates presented in de Gasperin et al. (2012). We adopt $B = 20\,\mu \mathrm{G}\,\sqrt{n_e / 0.1\ \mathrm{cm}^{-3}}$, which yields B = 24 (5) μG for the central value (at the halo boundary, r_max = 34 kpc). The density scaling yields a constant Alfvén velocity $\upsilon _A \approx 130\ \mathrm{km\ s}^{-1}$. The plasma beta factor β = P_th/P_B increases from 20 (at the center) to 47 (at r_max), which corresponds to a magnetic pressure contribution of X_B = 0.05 (center) and 0.02 (r_max).

Radiative cooling of cluster gas is dominated by thermal bremsstrahlung and collisional excitation of heavy ions to metastable states and the subsequent line emission following de-excitation (the latter of which dominates for particle energies ≲ 3 keV). These processes are characterized by the cooling rate

$$\begin{equation} \mathcal {C}_\mathrm{rad} = n_e n_t \Lambda _\mathrm{cool}(T, Z), \end{equation} \tag{ 8 }$$

where the target density n_t is the sum of ion densities and we defined the cooling function Λ_cool(T, Z) that depends on temperature T and metal abundance Z. The cooling function is fairly flat between the range of temperatures of kT₀ = 1–2 keV we encounter in the central regions, bounded by the radio halo. The X-ray inferred metallicity map shows substantial spatial variations on small scales ranging between Z₀ = 0.7 and 1.3 Z_☉, which are of similar amplitude as the overall radial trend in Z (Million et al. 2010). Hence, we decided to bracket the cooling function for these ranges of temperature and metallicity, which yields Λ_cool(T₀, Z₀) ≈ (1.8–2.8) × 10⁻²³ erg cm³ s⁻¹ while assuming collisional ionization equilibrium (see Figure 8 of Sutherland & Dopita 1993).

Figure 3 shows the radial profile of the radiative cooling rate $\mathcal {C}_\mathrm{rad}$ and that of the CR heating rate $\mathcal {H}_\mathrm{CR}$ assuming $\upsilon _A=\mathrm{const}$. For radii r ≳ 2 kpc, CR heating balances cooling already for the "smoothed" gradient pressure term ($\nabla _r \langle P_\mathrm{CR}\rangle _\Omega$). The central pressure plateau causes the CR pressure gradient and hence the CR heating to drop. To model the large visible fluctuations in P_th (and hence in P_CR according to Equation (5)), we assume that they are caused by weak shocks of typical Mach number $\mathcal {M}=1.1$ with a size scale that is of the order of the radius (i.e., we adopt ε = 0.1 in Equation (6)). This fluctuating pressure term appears to be critical in balancing the cooling rate with CR heating in the central region for r < 2 kpc. We conclude that CR Alfvén wave heating can globally balance cooling, hence justifying our assumption of a constant CR-to-thermal pressure. In Appendix B, we show that a varying $\upsilon _A$ (within physical bounds) does not change our main conclusion about the ability of CR heating to balance radiative cooling globally although there exists the possibility that heating dominates cooling at the external halo regions, which would cause these to adiabatically expand.

Zoom In Zoom Out Reset image size

Particularly notable is the non-trivial result that the radial shapes of $\mathcal {C}_\mathrm{rad}\propto n^2$ and $\mathcal {H}_\mathrm{CR}\propto \nabla _r P_\mathrm{th}$ closely track each other. In fact, if CR heating dominates, P_th should be determined by the pressure profile of the streaming CRs, and hence this resemblance of the radial heating and cooling rate profiles is an important signature of self-consistency. The sharp boundaries of the radio halo of Virgo A may suggest that CRs are confined to within its boundaries, perhaps as a result of azimuthally wrapped magnetic fields. Hence, the CR heating rates become very uncertain outside the maximum halo radius.

We will argue that this similarity of the CR heating and radiative cooling profiles is not a coincidence in M87, but rather generically expected in models of CR Alfvén wave heating due to the presence of a physical attractor solution for non-equilibrium initial conditions (albeit constrained to the time average). Global equilibrium requires matching heating and cooling rates at any radius, i.e., similar radial profiles as well as a matching overall normalization. The latter is a consequence of the self-regulating property of AGN feedback. If $\mathcal {H}_\mathrm{CR}<\mathcal {C}_\mathrm{rad}$ at some small radius, this triggers runaway cooling and eventually accretion onto the central black hole. Phenomenologically, this is connected to the launch of relativistic jets (McNamara & Nulsen 2007) that carry a fraction of the accretion energy to kiloparsec scales and inflate radio lobes with CRs (see Section 2.3) that boost the CR heating rate to become larger or equal to the radiative cooling rate (Guo & Oh 2008). If $\mathcal {H}_\mathrm{CR}> \mathcal {C}_\mathrm{rad}$, the CR heating overpressurizes this region on the Alfvénic crossing time. This heated patch cannot radiate away the thermal energy and expands quasi-adiabatically in the cluster atmosphere until $\mathcal {H}_\mathrm{CR}\simeq \mathcal {C}_\mathrm{rad}$ (since we assume here that the Alfvén velocity is smaller than the sound speed; see Section 3.1). This implies a transfer of gas mass to larger radii and a flattening of the density profile.

To demonstrate that the radial heating and cooling profiles generically match in this model, we use a constant gas-to-dark-matter ratio inside the dark-matter scale radius, i.e., n∝r⁻¹, which is characteristic for the central density profile in CCs (and which applies to the density profile of the Virgo cluster for 2 kpc ≲ r ≲ 100 kpc; see Equation (2)). Outside the temperature floor, we adopt a temperature profile with a small positive slope, $T\propto r^{\alpha _t}$ (with α_t ≲ 0.3 typically), so that the pressure profile scales as $P_\mathrm{th}=n kT\propto r^{\alpha _t-1}$. As a result, we obtain a scaling of the cooling and CR heating rates with radius of $\mathcal {C}_\mathrm{rad} \propto r^{-2}$ and $\mathcal {H}_\mathrm{CR}\propto d (r^{\alpha _t-1})/dr \propto r^{\alpha _t-2}$. Those are identical for a constant temperature (α_t = 0). For a smoothly rising temperature profile (α_t ≲ 0.3) the heating rate profile is tilted and slightly favored over the cooling rates at larger radii, provided they match up at some inner radius (see Figure 3). A corollary of this consideration is that the onset of cooling is smoothly modulated from the outside in, as we will see in the following.

CR heating appears to be able to keep the gas in global thermal equilibrium in M87. However, this does not necessarily imply local thermal stability. We consider small isobaric perturbations to this equilibrium configuration, which maintain pressure equilibrium at each radius. A necessary condition for maintaining hydrostatic equilibrium is a short dynamical time τ_dyn in comparison to the CR heating time τ_A. We obtain the following condition on the dynamical timescale

$$\begin{equation} \tau _\mathrm{dyn} \simeq \sqrt{\frac{R}{g}} \approx \frac{R}{c_s} \stackrel{!}{<} \frac{R}{\upsilon _A} \approx \tau _A, \end{equation} \tag{ 9 }$$

where c_s is the sound speed and g = GM(< R)/R² is the gravitational acceleration. This inequality can be rewritten to yield $c_s / \upsilon _A \simeq \sqrt{\beta \gamma _\mathrm{th}/2} > 1$, where β is the plasma beta factor and γ_th = 5/3. For magnetic field estimates adopted here, we obtain $c_s / \upsilon _A \simeq 4$ (6.3) at the center (at the halo boundary, r = 34 kpc) so that this inequality is always fulfilled.

In Appendix A we derive a general instability criterion for a fluid that is cooling by thermal bremsstrahlung and line emission and heated by a generic mechanism with a power-law dependence on ρ and T while allowing for thermal and CR pressure support. We begin the discussion with a transparent order-of-magnitude calculation that provides insight into the physics and complements the formal treatment in Appendix A. The generalized Field criterion (Field 1965; Balbus 1986) for thermal instability in the presence of isobaric perturbations is given by

$$\begin{equation} \left.\frac{\partial (\mathcal {L}/ T)}{\partial T}\right|_P < 0, \end{equation} \tag{ 10 }$$

where $\mathcal {L}$ is defined as the net loss function per unit mass such that $\rho \mathcal {L}=\mathcal {C}-\mathcal {H}$, where $\mathcal {C}$ and $\mathcal {H}$ are the cooling and heating rates per unit volume.

Here, we assume the existence of small-scale tangled magnetic fields that result, e.g., from magnetohydrodynamic turbulence. A cooling gas parcel necessarily pulls in the flux-frozen magnetic fields alongside the contracting gas. The CRs are bound to stay on a given flux tube and will experience adiabatic compression during the onset of collapse. The streaming motion of CRs into the increasingly denser gas parcel is driven by the free energy stored in the large-scale CR pressure gradient. On the Alfvénic crossing timescale that is of relevance for our stability analysis, the CR population does not suffer non-adiabatic losses such as hadronic interactions (e.g., Enßlin et al. 2007). In particular, CR streaming is an adiabatic process for the CR population (Kulsrud 2005), implying locally P_CR∝ρ^γ, where γ is the adiabatic exponent of the CR population. We emphasize that this local adiabatic behavior of P_CR does not hold for the equilibrium distribution of CRs that is subject to non-adiabatic gain and loss processes, which modify the "CR entropy" (i.e., its probability density of the microstates).⁷ We obtain the scaling of the CR heating rate with temperature,

$$\begin{equation} \mathcal {H}_\mathrm{CR}= - {\boldsymbol{\upsilon}} _A {{{\cdot }}} \boldsymbol{\nabla } P_\mathrm{CR}\propto \rho ^{\gamma +1/3} \propto T^{-\gamma -1/3}. \end{equation} \tag{ 11 }$$

Here we identified the characteristic size of our unstable clump with an exponential pressure scale height h so that the local gradient length scales as δl ∼ h∝ρ^−1/3, and we assumed $\upsilon _A=\mathrm{const}$. We suppress the dependence on pressure, which is constant for isobaric perturbations. The radiative cooling rate scales with temperature as

$$\begin{equation} \mathcal {C}_\mathrm{rad} = \rho ^2 \left[\Lambda _0 Z^2 T^{1/2} + \Lambda _l(T,Z)\right] \propto T^{-3/2} + \eta T^{-2+\chi }, \end{equation} \tag{ 12 }$$

where the first and second term in the brackets denote the bremsstrahlung and line-emission cooling. Here we introduce the logarithmic temperature slope of the line cooling function χ and the ratio of line-to-bremsstrahlung cooling rate η. For an ultra-relativistic CR population with γ = 4/3, we find a stable fixed point in the bremsstrahlung regime since the slope of the heating rate (=−5/3) is steeper than that of the cooling rate (=−3/2). Increasing (decreasing) the temperature of a perturbed parcel of gas implies a stronger cooling (heating) that brings the parcel back to the background temperature. However, in the line-cooling regime (at (1–3) × 10⁷ K), which is characterized by η ≳ 1 and χ ≲ 0, the slope of the cooling rate (≲ − 2) is steeper than that of the CR heating rate (=−5/3), implying instability.

This simple order-of-magnitude calculation is substantiated in Figure 4, which shows the instability criterion derived in Equation (A11) and restricted to CR Alfvén wave heating as the sole source of heat in the central region of the Virgo cluster. A gas in thermal equilibrium ($\rho \mathcal {L}=0$) is thermally unstable if

$$\begin{equation} \mathcal {D} \equiv \frac{1}{2} + \eta \chi - \frac{(5/3 - \gamma) (1 + \eta)}{1+X_\mathrm{CR}\gamma } < 0. \end{equation} \tag{ 13 }$$

We compute the logarithmic temperature slope of the line cooling function χ(T) and the ratio of line-to-bremsstrahlung cooling rate η(T) for a gas at solar metallicity and assume collisional ionization equilibrium (Sutherland & Dopita 1993). In Figure 4, we show $\mathcal {D}$ as a function of temperature and adopt values appropriate for M87 (γ = 1.37 and X_CR = 0.31; see Appendix A.2). For most of the temperature range, we find that $\mathcal {D}<0$ and the gas is subject to thermal instability. However, there are a few "islands of stability" where $\mathcal {D}>0$.

Zoom In Zoom Out Reset image size

The location of such an "island of stability" corresponds to a sufficiently positive slope of the cooling function Λ_cool as a function of T (i.e., it is located at the low-temperature side of every line-cooling complex). Increasing the metallicity Z increases the cooling rate and sharpens the line complexes. This implies more stable "islands of stability," i.e., they protrude more above the "ocean of instability," but they also extend over a slightly smaller temperature range. However, increasing Z does not significantly change their localization in temperature (as long as the CR heating can compensate the increasing cooling rate that accompanies the elevated Z).

Apart from parameters associated with radiative cooling (e.g., metallicity), $\mathcal {D}$ also depends on cluster-immanent parameters associated with the CR distribution and magnetic fields, i.e., the CR adiabatic exponent γ, the CR-to-thermal pressure ratio X_CR, and $\upsilon _A$. Of those, only X_CR and the possible density dependence of $\upsilon _A$ can in principle vary substantially among clusters. To address the robustness of the shape of the instability criterion $\mathcal {D}$ with respect to variations in X_CR, in Figure 4 we additionally show $\mathcal {D}$ for a 10 times lower CR pressure (which would, however, not suffice to maintain global thermal equilibrium in M87). Most importantly, the shape of $\mathcal {D}$ versus temperature barely changes, except that the second "island of stability" vanishes for low values of X_CR. In Appendix B, we address how a varying $\upsilon _A$ impacts the local stability analysis. In particular, for the relation $\upsilon _A\propto \rho ^{1/2}$, which represents the extreme case of collapse solely perpendicular to $\boldsymbol { {B}}$, we find more prominent "islands of stability" while still obeying global thermal equilibrium (within the range of uncertainty).

Before we can discuss the fate of the cooling gas, we will be considering the typical length scale of a region that becomes subject to thermal instability. We expect cooling gas parcels only to appear in systems whose dimension is greater than a critical length scale λ_crit, below which CR Alfvén wave heating smoothes out temperature inhomogeneities. We will estimate λ_crit by considering thermal balance for a cool parcel of radius r embedded in a medium with an equilibrium CR pressure P_CR = X_CRP_th. CR streaming transfers an amount of energy to a given parcel at a rate ${\sim} r^3 f_s \upsilon _A (P_\mathrm{CR}/r)$, where f_s is a suppression factor that depends on the magnetic topology connecting the parcel with its surroundings. Line emission and bremsstrahlung cooling can radiate energy from the cloud at a rate ∼r³n_en_tΛ_cool(Z, T). Cooling and CR heating are thus approximately balancing each other for systems with a radius

$$\begin{equation} \lambda _\mathrm{crit} = \frac{f_s \upsilon _A P_\mathrm{CR}}{n_e n_t \Lambda _\mathrm{cool}} = \frac{f_s X_\mathrm{CR}}{\mu _t}\,\frac{\upsilon _A kT}{n_e \Lambda _\mathrm{cool}}, \end{equation} \tag{ 14 }$$

where μ_t = n_t/n = (3X_H + 1)/(5X_H + 3) ≃ 0.48 assuming primordial element composition.

In Figure 5 we show λ_crit as a function of radius for the parameters of M87 used throughout, i.e., $\upsilon _A = 130\ \mathrm{km\ s}^{-1}$, X_CR = 0.31, Λ_cool(T₀, Z₀) of Section 3.2, and the observed electron temperature and density profiles of Equations (2) and (3). The CR heating rate without azimuthal fluctuations (f_s = 1) yields a critical wavelength λ_crit > r at all radii, which implies that the unstable wavelengths cannot be supported by the cluster core.

Zoom In Zoom Out Reset image size

Since CRs stream anisotropically along magnetic fields, this can suppress the heating flux into a given volume element depending on the magnetic connectivity to the surroundings. In the case of an isotropically tangled small-scale magnetic field, CR heating along a given direction is suppressed by a factor f_s ≈ 0.3. As a result, the minimum size of an unstable region at 10 kpc is λ_crit = 5–8 kpc, depending on metallicity. However, such a large region is unlikely to be exclusively heated from one direction, which makes it challenging for cooling multi-phase gas to form unless a given region is sufficiently magnetically insulated with a suppression factor f_s ≪ 1. Since we have demonstrated in Figure 3 that CR heating globally balances cooling, a magnetic suppression of heating into one region implies an overheating in another. Because the gas is thermally unstable for temperatures in the range of (1–3) × 10⁷ K, the additional energy would be radiated away on the radiative cooling time.

The previous considerations enable us to construct a comprehensive "cooling flow" model for M87, which may serve as a blueprint for other CC clusters. We distinguish three different cases, each of which may be realized by a different magnetic suppression factor. However, quantifying their relative abundance requires appropriate numerical simulations.

Case A. In the absence of strong magnetic suppression of the CR heat flux, a gas parcel is thermally stabilized by CR heating because the system cannot support the unstable wavelengths, which are larger than the cluster radius.

Case B. For strong magnetic suppression of the CR heat flux, gas at temperatures T = (1–3) × 10⁷ K is thermally unstable, starts to cool radiatively, becomes denser than the ambient gas, and begins to sink toward the center. Since CRs are largely confined to stay on individual field lines, which are flux frozen into the contracting gas, the CR population is adiabatically compressed. This implies an increasing CR heating rate $\mathcal {H}_\mathrm{CR}\propto \rho ^{\gamma +1/3}$ (Equation (11)), which, however, increases at a slower rate than the cooling rate for T = (1–3) × 10⁷ K (see Section 3.3). The CRs transfer energy to the gas on a timescale comparable to the Alfvén crossing time over the CR pressure gradient length, $\tau _A = \delta l/\upsilon _A \sim 75$ Myr for δl = 10 kpc.

If the gas parcel remains magnetically insulated, it cannot be supplied with a fresh influx of CRs, which causes the CR heating to lag the radiative cooling, implying eventually the formation of multi-phase gas and potentially stars, provided that favorable conditions for their formation are met. Such a magnetic suppression could be realized through processes that dynamically form coherently aligned magnetic sheaths, which are wrapped around the gas parcel as a result of either the action of magnetic draping (Lyutikov 2006; Ruszkowski et al. 2007; Dursi & Pfrommer 2008), the heat-flux-driven buoyancy instability (Quataert 2008; Parrish & Quataert 2008), or through radial collapse, which amplifies the azimuthal components of the magnetic field, in combination with coherent motions of the cooling gas parcel so that the azimuthal magnetic field becomes causally unconnected to its surroundings.

The comparably large critical wavelength of the instability, λ_crit, at large radii (even for substantial magnetic suppression) may make it difficult for thermal instability to condense out large gas clouds (which requires a stable magnetic insulation on the correspondingly long timescale τ_A). Instead, it is more likely that cooling (multi-phase) clouds form in the central regions at radii ≲ 3 kpc because of the small λ_crit. For this case, the presence of small pockets of multi-phase gas at larger radii will then necessarily imply that they would have to be uplifted by eddies driven by the rising radio lobes.

Case C. For moderate magnetic suppression of a given gas parcel at T = (1–3) × 10⁷ K, CR heating will be able to approximately match radiative cooling. That means we can consider the gas parcel to be in global thermal equilibrium but subject to local thermal instability. If this global thermal equilibrium can be maintained during the cooling of a necessarily large region of wavelength λ ≲ r, and if the perturbations are in pressure equilibrium with the surroundings, then the results of our stability analysis apply. In particular, the gas collapse is halted at around 10⁷ K by the first "island of stability" situated at T_floor = (0.7–1) × 10⁷ K. Hence, for case C, we predict a temperature "floor" at the cluster center with typical values distributed within T_floor. The insensitivity of the shape of the first "island of stability" on the various parameters entering the instability criterion suggests that we found a universal mechanism that sets a lower temperature floor in CC clusters (as long as the CR pressure is sufficient to provide global thermal stability).

This prediction of a generic lower temperature floor at T_floor is in agreement with high-resolution X-ray observations of Virgo (Matsushita et al. 2002; Million et al. 2010; Werner et al. 2010), suggesting a viable solution for the puzzling absence of a large amount of lower temperature gas in Virgo that is expected for an unimpeded cooling flow. Despite earlier claims of a universal temperature profile in CC clusters (e.g., Allen et al. 2001), a temperature floor at a fixed fraction of one-third of the virial radius (Peterson et al. 2003) appears not to be supported by the latest data. Instead, using a large complete sample of 64 clusters with available high-quality X-ray data, Hudson et al. (2010) find a broad distribution of the central temperature drop of T₀/T_vir = 0.2–0.7 in CC clusters (see Figure 4 of Hudson et al. 2010). Note that T₀ is the temperature estimate of the central annulus and thus an upper limit to the true central temperature since T₀ depends on angular resolution and sensitivity of the cluster observation. The lower temperature in their modified cooling flow model is in the range (0.5–2) × 10⁷ K and clusters around 10⁷ K. It is generally smaller than T₀, which is expected if there are multi-temperature components present (see Figure 16 of Hudson et al. 2010). These findings are in perfect agreement with our presented model, which implies a lower temperature floor around (0.7–1) × 10⁷ K as an attractor solution but certainly can stray from this equilibrium value and assume higher central temperature values for an intermittent period of overheating. A temperature floor similar to that observed in M87 of around 5 × 10⁶ K (Werner et al. 2010) is found in comparable high-resolution observations of the Perseus cluster, Centaurus cluster, 2A 0335+096, A262, A3581, and HCG 62 (Sanders & Fabian 2007; Sanders et al. 2008, 2009, 2010). Clearly, the lack of a larger amount of gas that is radiatively cooling below 5 × 10⁶ K implies the presence of heating.⁸

However, metastable gas in the high-density tail of density fluctuations can cool below T ≃ 7 × 10⁶ K and hence leave this "island of stability." These dense fluctuations again become subject to thermal instability, sourcing some of the observable multi-phase gas. The collapse could again be halted at the second and third "island of stability" just above 10⁶ K and 10⁵ K. However, since cooling clouds at these temperatures are associated with higher densities and smaller length scales, the amount of CR heating now depends critically on the magnetic connectivity of the cloud with the surroundings since that determines whether there will be a sufficient influx of CRs responsible for sustaining the heating rate. These theoretical considerations compare favorably to the observational result of a multi-wavelength study of the emission-line nebulae located southeast of the nucleus of M87, which reveals multi-phase material spanning a temperature range of at least five orders of magnitude, from ∼100 K to ∼10⁷ K (Werner et al. 2013). This material has a small total gas mass and has most likely been uplifted by the AGN from the center of M87.

We emphasize that the success of our model came about without any tunable parameters and only used input from non-thermal radio and gamma-ray observations. The model implies that gas in the core region of Virgo is in global equilibrium but susceptible to local thermal instability. Hence, this provides a physical foundation of the phenomenological prescription used in numerical simulations that study the interplay of thermal instability and feedback regulation (McCourt et al. 2012; Sharma et al. 2012). Apparently, thermal instability can produce spatially extended multi-phase filaments only when the instantaneous ratio of the thermal instability and free-fall timescales falls below a critical threshold of about 10 (Sharma et al. 2012; Gaspari et al. 2012b). There are also notable differences since these prescriptions did not have the modulating effect of our "islands of stability." We may speculate that our model may then inherit some of their successes including quantitative consistency with data on multi-phase gas and star formation rates but holds the promise to cure the problem of the absence of a (physical) temperature floor in these numerical experiments.

If the gamma-ray emission of the low state of M87 results from hadronic CRp–p interactions, then (1) the emission should be extended with a characteristic angular profile worked out by Pfrommer & Enßlin (2003) for their isobaric CR model, (2) the gamma-ray emission should not show any sign of time variability, and (3) the emission should be characterized by the distinct spectral signature of decaying pions around 70 MeV (in the differential spectrum). The triad of these observations will provide strong support for the existence of CRps mixed into the ICM and will allow us to measure their pressure contribution. Spatially resolving the gamma-ray emission helps in estimating the large-scale CR pressure gradient and its contribution to the CR-heating rate. However, this does not shed light on the contribution from small-scale pressure fluctuations. The comparably coarse angular resolution of Fermi (68% containment angle of 12' above 10 GeV) precludes such a measurement in M87. In contrast, the superior resolution of current and future Cerenkov telescopes (which is about a few arcminutes with H.E.S.S. II above TeV energies and expected to be five times better with CTA, approaching theoretical limits of 03 at 1 TeV; see Wagner et al. 2009; Hofmann 2006) is a promising route toward measuring the pressure gradient of the high-energy CRps, which, however, carry only a small amount of the total pressure that is dominated by the GeV-energy regime.

We note that a critical prediction of this model is the similarity of the spectral indices of the gamma-ray emission above GeV energies (which corresponds to that of the CR population) and the injection index of the CRe population, as probed by low-frequency radio observations of the core regions of Virgo A. Similar observations of nearby AGNs in CC clusters are necessary to show the universality of the presented mechanism.

The hadronic CRp–p reaction also implies the injection of secondary CRes that contribute to the radio synchrotron emission in M87. Assuming steady state of the hadronic injection of secondary CRes and synchrotron cooling, we obtain a spectral index for the secondary radio emission of α_ν = α/2 = 1.13 (adopting the CRp index α = 2.26 inferred in Section 2.2). Using the formalism of Pfrommer et al. (2008), we calculate the hadronically induced radio halo emission within a spherical region of radius 34 kpc (adopting the same parameters as in the rest of the paper, namely, X_CR = 0.31, α = 2.26, q = 0.58, and the model for the magnetic field of Section 3.1). At 1.4 GHz (10.55 GHz), this emission component is expected to have a total flux density of 21.1 Jy (2.15 Jy), which falls short of the observed halo emission by a factor of 6.6 (10).

However, the spatial profile of the expected hadronically induced emission is very different from the observed radio halo emission, which is centrally more concentrated. If we exclude a cylinder centered on the central radio cocoon of radius 2.5 kpc (corresponding to 051), we obtain a flux density at 1.4 GHz (10.55 GHz) for the external hadronic halo of 19.8 Jy (2.0 Jy), which still falls short of the observed emission of 69.8 Jy (3.4 Jy) that subtends the same solid angle (see Figure 6) but makes up a larger flux fraction in comparison to that of the entire halo region. In particular, the hadronically induced emission is expected to make up ∼60% of the observed (external) halo emission at 10.55 GHz.

Zoom In Zoom Out Reset image size

In Figure 6 we also show the continuous injection model (green), which assumes an uninterrupted supply of fresh CRes from the central source (Pacholczyk 1970). The CRes are subject to IC and synchrotron cooling, which introduces a steepening of the power-law spectrum by Δα_ν = 0.5 above a break frequency. The pitch angles of CRes are conjectured to be continuously isotropized on a timescale shorter than the radiative timescale (Jaffe & Perola 1973). Additionally shown is a similar model, in which the source is switched off after a certain time (orange; Komissarov & Gubanov 1994). This causes the spectrum to drop exponentially above a second break frequency, which corresponds to the cooling time since the switch-off (see Appendix D for details). Depending on the model of primary electrons, our model predicts that the hadronically induced radio emission (blue) starts either to dominate the halo emission or to make up a large fraction of the total emission at frequencies ν_trans ≳ 20 GHz. This would reveal itself as a radio spectral flattening at this transition frequency ν_trans. Clearly, the possibility of such an emission component should be accounted for when drawing conclusions about the CRe pitch-angle distribution in spectral aging models of synchrotron and IC emitting CRes, especially if the high-frequency data are sparse so that individual data points impart high statistical weights in the fits.

If CRes and protons have a very different distribution within the radio halo (which could be achieved through enhanced radiative losses of CRes in regions of strong magnetic fields), there may be regions in which hadronically induced emission prevails over the primary emission already at frequencies ν < 20 GHz. Those regions should then also show a radio spectral flattening in comparison to the adjacent primary emission-dominated regions.

Faraday rotation measurements of background sources or polarized regions of Virgo A (central region, halo, or rising lobes) will provide a reliable estimate not only of the spatial distribution of field strengths but also about the magnetic morphology. This is important for predicting the heating rates (that scale with the Alfvén velocity) and the microphysics of CR Alfvén wave heating, which depends on in situ alignments of magnetic fields with the CR pressure gradients. In particular, CR release timescales from the radio lobes may depend on the pre-existing morphology of the magnetic fields that dynamically drape around a super-Alfvénically rising lobe (Lyutikov 2006; Ruszkowski et al. 2007). This magnetic draping effect stabilizes the interfaces against the shear that creates it in the first place (Dursi 2007). However, magnetic draping can only occur if the ratio of the integral length scale of the magnetic power spectrum λ_B and the curvature radius R of the radio lobe at the stagnation line obeys the relation (see Supplementary Information of Pfrommer & Dursi 2010)

$$\begin{eqnarray} &&\frac{\lambda _B}{R} \gtrsim \mathcal {M}_A^{-1} = \left(\mathcal {M}\sqrt{\frac{\gamma _\mathrm{th}\beta }{2}}\right)^{-1} \simeq 0.25\, \mathcal {M}^{-1}\,\left(\frac{\beta }{20}\right)^{-1/2}, \nonumber\\ \end{eqnarray} \tag{ 15 }$$

where $\mathcal {M}$ and $\mathcal {M}_A$ are the sonic and Alfvénic Mach numbers of the rising lobes, γ_th = 5/3, and β is the plasma beta factor. The violation of this relation could cause Kelvin–Helmholtz and Rayleigh–Taylor instabilities to disrupt the lobes and to release CRs into the ICM. Higher order Faraday rotation measure statistics are sensitive to the statistics of the projected magnetic field distribution. Complementary to this are magnetic draping studies at spiral galaxies within this central region (Pfrommer & Dursi 2010) that provide a means of measuring the projected in situ geometry of the magnetic field at the galaxies' position. Additionally, those could be of paramount importance to provide indications for a spatially preferred direction of magnetic fields that could hint at the presence of a magnetic buoyancy instability.

Cosmological cluster simulations of CR feedback by AGNs suggest a substantial modification of the X-ray emissivity and the thermal SZ effect, with modifications of the integrated Compton-y parameter within the virial radius up to 10% (Sijacki et al. 2008; Vazza et al. 2013). High-resolution SZ observations can even reveal the composition of the plasma within radio lobes (Pfrommer et al. 2005). Different dynamically dominating components within radio lobes, which can generally range from relativistic to hot thermal compositions, produce different SZ contrasts at the position of these lobes. A relativistic filling implies an SZ cavity at the lobe positions. This is because CRps do not Compton upscatter cosmic microwave background photons while CR leptons that are in pressure equilibrium with the ambient ICM scatter only a tiny fraction of background photons to X-ray energies in comparison to a situation where the lobes would be filled with a slightly hotter than average thermal population. In fact, the latter case would make the lobes almost indistinguishable from the ambient cluster SZ effect. Both signatures are also different from the hypothetical alternative of a hot but sub-relativistic particle population. This would cause an SZ effect that is spectrally distinguishable from the ordinary thermal SZ effect since it crosses through a null above the canonical 217 GHz (but within the microwave band). High-resolution "SZ spectroscopy" would be ideally suited to detect the difference between these cases (Pfrommer et al. 2005). However, the large dynamic range required for such a detection in Virgo A will make such a measurement very challenging.

Thermal conduction of heat from the hot outer cluster regions to the cooling cores has been proposed to solve the problem (Narayan & Medvedev 2001; Zakamska & Narayan 2003). However, there are a number of drawbacks to this proposal. (1) The strong temperature dependence of thermal conductivity (κ∝T^5/2) makes this a very inefficient process for groups and small clusters. (2) Thermal conduction ceases to be important if gas cools to low temperatures and, as such, requires in some clusters a conductivity above the classical Spitzer value, as well as a high degree of fine tuning (Voigt & Fabian 2004; Conroy & Ostriker 2008). (3) The heat-flux-driven buoyancy instability (Quataert 2008) causes the magnetic field lines to align with the gravitational equipotential surfaces of the cluster in the nonlinear stage of the instability (Parrish & Quataert 2008). At first sight, this should suppress the inward transport of heat by more than an order of magnitude and should reinforce the cooling flow problem. The instability, however, is easily quenched by the presence of external turbulence at the level of ∼1% of the thermal pressure (Parrish et al. 2010; Ruszkowski & Oh 2010), leaving the question of how much thermal conduction is suppressed by the presence of magnetic fields in realistic cosmological situations. Overall, thermal conduction could certainly aid in mitigating the fast onset of cooling at intermediate cluster radii (Voit 2011) but is unlikely to play a major role in the cluster centers. In particular, we show in Appendix C that CR heating dominates over thermal conduction for the entire CC region of the Virgo cluster.

Over the past years, AGN feedback has become the lead actor because of its promising overall energetics and the observational confirmation that the nuclear activity appears to be tightly coupled to a self-regulated feedback loop (Bîrzan et al. 2004). Guided by high-resolution Chandra X-ray observations, which emphasize sharp density features in high-pass filtered maps as a result of the gas being pushed around by AGN jets and lobes, many papers have focused on the dissipation of PdV work done by the inflating lobe on the surroundings. The detailed heating mechanisms include sound waves, weak shocks, and turbulence (see McNamara & Nulsen 2012 and references therein). Since sound waves and weak shocks are excited in the center by the AGN, their outward-directed momentum flux in combination with the small dissipation efficiency implies that only a small fraction of this energy will be used for heating the central parts (e.g., Gilkis & Soker 2012). Those parts, however, have the lowest cooling time. Additionally, the implied heating rate is very intermittent in space and time with an unknown duty cycle. Nevertheless, these processes will certainly happen and (if the outburst energy is large enough) may locally even dominate during some times (Nulsen et al. 2005). As we argue in this paper, these processes are, however, unlikely to be the most important heating source because they do not provide a generic solution to the lower temperature floor observed in CCs⁹ and they are not favored on energetic grounds due to their small heating efficiency.

Assuming a relativistic filling of the lobes (for which we provide arguments in Section 2.3), three-quarters of the total lobe enthalpy is stored in internal energy after the end of the momentum-dominated jet phase. These light relativistic lobes rise buoyantly and get subsequently disrupted by Rayleigh–Taylor or Kelvin–Helmholtz instabilities. Once released, those CRps stream at about the Alfvén velocity with respect to the plasma rest frame and excite Alfvén waves through the "streaming instability" (Kulsrud & Pearce 1969). Damping of these waves transfers CR energy and momentum to the thermal gas at a rate that scales with the product of the Alfvén velocity and the CR pressure gradient (Guo & Oh 2008). While very plausible, there was the open question of whether the CRp population has been mixed into the central ICM and whether it is abundant enough to provide sufficiently strong heating to offset the cooling catastrophe. In this paper, we answer both questions affirmatively for the case of M87 on the basis of a detailed analysis of non-thermal observations. Moreover, in addressing the issue of local thermal stability in this model, we put forward an attractive solution to the problem of an apparent temperature floor in CC clusters.

In particular, we observe in M87 that the lobes get disrupted at about one scale height so that the density drops from the radius where the lobes start to rise buoyantly to the current position at the radio halo boundary by a factor of about f = ρ₁/ρ₂ ≃ 10. Hence, the internal (relativistic) energy is adiabatically cooled from U₁ = 3P_thV to U₂ = U₁f^{1 − γ} ≃ 1.5P_thV at the time of release. A fraction of CRs gets caught in the advective downdrafts that are driven by the buoyantly rising relativistic lobes. This re-energizes the CRs adiabatically at the expense of mechanical (e.g., turbulent) energy, thereby partially replenishing the energy initially lost to the surroundings and making it available for CR heating at the dense ICM regions with short cooling times. We end up with a fraction of (1.5–3) P_thV in CR and (1–2.5) P_thV in mechanical energy. Since streaming is an adiabatic process for the CR population (Kulsrud 2005), a fraction of the CR population can pass through several cycles of outward streaming and downward advection, suggesting the physical realization of a larger CR energy fraction closer to 3 P_thV and a correspondingly smaller mechanical energy fraction. The central cooling region delineated by the radio halo is characterized by a sharp radio boundary, which suggests magnetic confinement of CRs that can be realized by an azimuthal field configuration. Hence, most of the CR energy will get dissipated within the halo boundary, whereas it is not clear which fraction of the mechanical energy gets dissipated in the central region. Among other uncertainties, this depends on the partitioning of mechanical energy into turbulence and shocks.

The puzzling absence of fossil CRes implied by the LOFAR observations of M87 can be naturally explained by the observed disruption of radio lobes during their rise in the gravitational potential, e.g., through the action of Kelvin–Helmholtz and Rayleigh–Taylor instabilities. This releases and mixes CRes and protons into the ambient dense ICM. We argue that mixing is essential in this picture so that the CRes have direct contact with the dense ambient plasma, which causes CRes to efficiently cool through Coulomb interactions. This implies a maximum CRe cooling timescale of ≃ 40 Myr, which is similar to the age of the Virgo A radio halo. A direct consequence of this interpretation is that CRps are similarly mixed with the dense gas and must collide inelastically with gas protons. This produces pions that decay into a steady gamma-ray signature, which resembles that of the flux in the "low state" of M87 observed by Fermi and H.E.S.S. In particular, the CRp population (witnessed through the gamma-ray emission) shares the same spectral index as the injected CRe population in the core region of Virgo A, which we consider as strong support for a common origin of the radio and gamma-ray emission.¹⁰ Hence, the hadronically induced gamma-ray emission enables us to estimate the CR pressure, which has a CR-to-thermal pressure contribution of X_CR = 0.31 (allowing for a conservative Coulomb cooling timescale of the order of the radio halo age that modifies the low-energy part of the CRp distribution).

We show that streaming CRs heat the surrounding thermal plasma at a rate that balances that of radiative cooling on average at each radius. Hence, the thermal pressure profile is a consequence of the interplay of CR heating and radiative cooling and should match that of the streaming CRs (hence justifying in retrospect our assumption of a constant CR-to-thermal pressure in our analysis). We argue that this similar radial scaling of the heating and cooling rate profiles is not a coincidence in M87 but rather generically expected in models of CR Alfvén wave heating.

However, the global thermal equilibrium is locally unstable if the gas temperature drops below kT ≃ 3 keV (for solar metallicity). The fate of a gas parcel depends on its magnetic connectivity to the surroundings, which may obey a yet unknown distribution function. Nevertheless, we can distinguish three different cases.

• 1.  
  For magnetically well-connected patches without appreciable fluctuations in the CR heating rate, temperature inhomogeneities are smoothed out on scales smaller than the cluster-centric radius. As a result, the system is stabilized because it cannot support the unstable wavelengths.
• 2.  
  Contrarily, for strong magnetic suppression of the CR heat flux, the thermally unstable gas starts to cool radiatively, becomes denser than the ambient gas, and sinks toward the center. After an initial phase when the magnetically trapped CRs within the parcel deposit their energy to the gas, CR heating lags cooling. This eventually causes the formation of multi-phase gas and potentially stars. At large radii, the critical wavelength λ_crit above which thermal instability can operate is comparably large, thereby making it more likely for multi-phase gas to form further inward where λ_crit is smaller. Thus, small clumps of multi-phase gas at larger radii would have to be uplifted by eddies driven by the rising radio lobes.
• 3.  
  For the perhaps most probable case of moderate magnetic suppression, CR heating will be able to approximately match radiative cooling. If this global thermal equilibrium can be maintained during the Lagrangian evolution of the cooling gas parcel that is in pressure equilibrium with its surroundings, then the gas collapse is halted at around 1 keV by an "island of stability." There, the CR heating rate has a steeper dependence on temperature than the radiative cooling rate. Increasing (decreasing) the temperature of a perturbed parcel of gas implies a stronger cooling (heating) that brings the parcel back to the background temperature. This is consistent with the observed temperature floor of 0.7–1 keV seen in the X-ray observations of the center of Virgo. Since this result is insensitive to significant variations of the CR pressure support, this may suggest a viable solution to the "cooling flow problem" in other CC clusters (as long as the CR heating rate is sufficient to provide global thermal stability). However, the high-density tail of density fluctuations that is associated with temperatures kT ≲ 0.7 keV again becomes subject to thermal instability, thereby providing a plausible mechanism for the formation of some of the observed cooling multi-phase medium in and around M87.

We note that the success of the CR heating model came about without any tunable parameters and only used input from radio and gamma-ray observations. We present several testable predictions of this model that include gamma-ray, radio, and SZ observations.

We emphasize the importance of non-thermal observations to establish this physical picture of AGN feedback and to pinpoint the details of this heating mechanism in M87. While tempting, we caution against directly extrapolating our result to other CC clusters. To faithfully do so, we would need similar non-thermal observations to establish that CRs are injected and mixed into the ambient ICM and to estimate the CR pressure contribution that determines the CR heating rate. Since the residency time of the CRps in the central core region is longer than the age of the radio halo, CR heating can take place even without (high-)frequency radio halo emission.

In particular, CR release timescales from the radio lobes are crucial and may depend on the pre-existing morphology of the magnetic fields that dynamically drape around the rising lobes (Ruszkowski et al. 2007). Those draped fields provide stability in the direction of the magnetic field against the shear that creates it if the magnetic coherence scale is large enough for draping to be established in the first place. Future Faraday rotation measure studies and galaxy draping observations (Pfrommer & Dursi 2010) may help to clarify this point.

We hope that this work motivates fruitful observational and theoretical efforts toward consolidating the presented picture and to establish non-thermal observations of the underlying high-energy processes as a driver for the understanding of the evolution of cluster cores and perhaps the modulation of star formation in elliptical galaxies.

Consider a fluid subject to external heating and cooling. The generalized Field criterion (Field 1965; Balbus 1986) for thermal instability in the presence of isobaric perturbations is given by

$$\begin{equation} \left.\frac{\partial (\mathcal {L}/ T)}{\partial T}\right|_P < 0. \end{equation} \tag{ A1 }$$

Here $\mathcal {L}$ is the net loss function in units of erg s⁻¹ g⁻¹ defined such that $\rho \mathcal {L}=\mathcal {C}-\mathcal {H}$, where $\mathcal {C}$ and $\mathcal {H}$ are the cooling and heating rates per unit volume. $\mathcal {L}$ is generally a function of gas mass density ρ, temperature T, and possibly space and time. Earlier local stability analyses neglected CR pressure support (Loewenstein et al. 1991) or were restricted to a cooling function of the interstellar medium at temperatures around 10⁴ K (Wiener et al. 2013b). We study the instability in the presence of a relativistic fluid and generalize the derivation by Conroy & Ostriker (2008) to also include cooling by line emission in addition to thermal bremsstrahlung, which will turn out to be important for putting the complete picture in place (in particular for cooling flows in clusters).

A gas in thermal equilibrium ($\rho \mathcal {L}=0$) obeys the following identity:

$$\begin{equation} \left.\frac{\partial (\mathcal {L}/ T)}{\partial T}\right|_P = \frac{1}{\rho T}\left[ \left.\frac{\partial (\rho \mathcal {L})}{\partial T}\right|_\rho + \left.\frac{\partial \rho }{\partial T}\right|_P \left.\frac{\partial (\rho \mathcal {L})}{\partial \rho }\right|_T \right]. \end{equation} \tag{ A2 }$$

The first and third derivatives on the right-hand side are governed by the heating and cooling processes under consideration, while the second depends on the various sources of pressure support. We consider pressure support provided by thermal gas and a relativistic fluid (dominated by CRps):

$$\begin{equation} P = P_\mathrm{th} + P_\mathrm{CR}= K_1 \rho T + P_\mathrm{CR}(\rho), \end{equation} \tag{ A3 }$$

where K₁ is a constant and P_CR is only a function of density. We use the following identity

$$\begin{equation} \left.\frac{\partial \rho }{\partial T}\right|_P = -\left.\frac{\partial P}{\partial T}\right|_\rho \left(\left.\frac{\partial P}{\partial \rho }\right|_T \right)^{-1} \end{equation} \tag{ A4 }$$

to obtain

$$\begin{eqnarray} \left.\frac{\partial \rho }{\partial T}\right|_P &=& {-}\frac{K_1\rho }{K_1 T + \left.\left(\partial P_\mathrm{CR}/\partial \rho \right)\right|_T} \nonumber\\ &=& {-}\frac{\rho }{T}\frac{P_\mathrm{th}}{P_\mathrm{th} + \rho \left.\left(\partial P_\mathrm{CR}/\partial \rho \right)\right|_T} \nonumber\\ &=& {-} \frac{\rho }{T}\frac{1}{1 + X_\mathrm{CR}\left.\left(\partial \ln P_\mathrm{CR}/\partial \ln \rho \right)\right|_T}, \end{eqnarray} \tag{ A5 }$$

where we defined X_CR ≡ P_CR/P_th.

We now assume that the heating has a power-law dependence on ρ and T and account for cooling by thermal bremsstrahlung and line emission to get

$$\begin{eqnarray} &&\rho \mathcal {L}= \mathcal {C}_\mathrm{rad} - \mathcal {H} = \rho ^2 \left[\Lambda _b (T,Z) + \Lambda _l(T,Z)\right] - \Gamma _0 \rho ^\beta T^{-1/2+\delta }, \nonumber\\ \end{eqnarray} \tag{ A6 }$$

where Λ_b(T, Z) = Λ₀Z²T^1/2 and Λ_l(T, Z) are the cooling functions due to bremsstrahlung and line emission, respectively. Both depend on metallicity and temperature. (For brevity, we omit the arguments in the following.) Λ₀ and Γ₀ are constants, and we define the identical power-law indices δ and β for the heating as in Conroy & Ostriker (2008) to ease comparison with their results. Using the relation for thermal equilibrium ($\rho \mathcal {L}=0$), we have

$$\begin{eqnarray} &&\left.\frac{\partial (\rho \mathcal {L})}{\partial T}\right|_\rho = \Lambda _b\rho ^2T^{-1} \left[ 1-\delta + \frac{\Lambda _l}{\Lambda _b} \left(\left.\frac{\partial \ln \Lambda _l}{\partial \ln T}\right|_\rho + \frac{1}{2} - \delta \right) \right] \nonumber\\ \end{eqnarray} \tag{ A7 }$$

and

$$\begin{equation} \left.\frac{\partial (\rho \mathcal {L})}{\partial \rho }\right|_T = \Lambda _b\rho \, (2-\beta) \left(1 + \frac{\Lambda _l}{\Lambda _b} \right). \end{equation} \tag{ A8 }$$

Combining those with Equations (A2) and (A5), we finally obtain

$$\begin{eqnarray} \left.\frac{\partial (\rho \mathcal {L})}{\partial T}\right|_P &=& \Lambda _b\rho T^{-2} \left[ 1-\delta + \frac{\Lambda _l}{\Lambda _b} \left(\left.\frac{\partial \ln \Lambda _l}{\partial \ln T}\right|_\rho + \frac{1}{2} - \delta \right) \right. \nonumber\\ && -\left. \frac{2-\beta }{1+X_\mathrm{CR}\gamma }\, \left(1 + \frac{\Lambda _l}{\Lambda _b} \right) \right]. \end{eqnarray} \tag{ A9 }$$

To simplify the instability criterion, we define the logarithmic density slope of the CR pressure γ, the logarithmic temperature slope of the line cooling function χ, and the ratio of line-to-bremsstrahlung cooling, η:

$$\begin{eqnarray} \gamma & \equiv & \left.\frac{\partial \ln P_\mathrm{CR}}{\partial \ln \rho }\right|_T,\quad \chi = \chi (T,Z) \equiv \left.\frac{\partial \ln \Lambda _l}{\partial \ln T}\right|_\rho ,\quad \,{\rm and} \nonumber\\ \eta &=& \eta (T,Z) \equiv \frac{\Lambda _l}{\Lambda _b}. \end{eqnarray} \tag{ A10 }$$

Thus, the gas is thermally unstable if the discriminant $\mathcal {D}$ is negative:

$$\begin{equation} \mathcal {D} \equiv 1-\delta + \eta \left(\chi + \frac{1}{2} - \delta \right) - \frac{(2-\beta) (1 + \eta)}{1+X_\mathrm{CR}\gamma } < 0. \end{equation} \tag{ A11 }$$

For η = 0 we recover the result of Conroy & Ostriker (2008) if we identify X_CR ≡ 1/α' in their notation. Clearly, the presence of line cooling (η > 0) can substantially modify the instability criterion, depending on the temperature slope χ.

Now, we specify this instability criterion to our case at hand of CR Alfvén wave heating. We assume the existence of small-scale tangled magnetic fields that result, e.g., from magnetohydrodynamic turbulence. Then, magnetic flux freezing during the contraction of unstable cooling gas ensures adiabatic compression of CRs. On the Alfvénic crossing time that is of relevance for our stability analysis, the CR population does not suffer non-adiabatic losses such as hadronic interactions (e.g., Enßlin et al. 2007). In particular, CR streaming is an adiabatic process for the CR population (Kulsrud 2005), implying locally P_CR∝ρ^γ and yielding

$$\begin{equation} \mathcal {H}_\mathrm{CR}= - {\boldsymbol{\upsilon}} _A {{{\cdot }}} \boldsymbol{\nabla } P_\mathrm{CR}= \Gamma _0 \rho ^\beta T^{-1/2+\delta } \propto \rho ^{\gamma +1/3}. \end{equation} \tag{ A12 }$$

Here we identified the characteristic size of our unstable clump with an exponential pressure scale height h so that the local gradient length scales as δl ∼ h∝ρ^−1/3 and we assume $\upsilon _A=\mathrm{const.}$ In practice, $\upsilon _A$ may inherit a weak density scaling depending on the topology of the magnetic field that is compressed or expanded. We explore the impact of a density-dependent $\upsilon _A$ on the global and local thermal equilibrium in Appendix B but defer dynamical simulations of this interacting CR-magnetohydrodynamic system to future work. Equation (A12) defines the power-law indices β = γ + 1/3 and δ = 1/2 for the heating rate. We emphasize that this local adiabatic behavior of P_CR does not hold for the equilibrium distribution that is subject to gain and loss processes, which modify the "CR entropy" (i.e., its probability density of the microstates). In fact, we argued in the main part of the paper that CR Alfvén wave heating establishes an equilibrium distribution that is characterized by X_CR = P_CR/P_th ≃ const. if the dominant heating source is provided by CRs.

To estimate the CR adiabatic exponent, we remind ourselves of the definition of the CR pressure of a CR power-law population defined in Equation (1), which is given by

$$\begin{eqnarray} P_\mathrm{CR} &=& \frac{m_p c^2}{3}\,\int _0^\infty d p\,\frac{d^4N}{dp\,d^3x} \,\frac{\upsilon }{c}\, p\nonumber\\ &=& \frac{C\,m_p c^2}{6} \, \mathcal {B}_{\frac{1}{1+q^2}} \left(\frac{\alpha -2}{2},\frac{3-\alpha }{2} \right), \end{eqnarray} \tag{ A13 }$$

where $\upsilon /c = p/\sqrt{1+p^2}$ is the dimensionless velocity of the CR particle and $\mathcal {B}_x(a,b)$ denotes the incomplete beta function. The local adiabatic exponent of the CR population is obtained by taking the logarithmic derivative of the CR pressure at constant CR entropy S_CR

$$\begin{eqnarray} \gamma &\equiv& \left.\frac{\partial \log P_\mathrm{CR}}{\partial \log \rho }\right|_{S_\mathrm{CR}}= \frac{\rho }{P_\mathrm{CR}} \left(\frac{\partial P_\mathrm{CR}}{\partial C}\frac{\partial C}{\partial \rho } + \frac{\partial P_\mathrm{CR}}{\partial q}\frac{\partial q}{\partial \rho } \right) \nonumber\\ & = & \frac{\alpha + 2}{3} - \frac{2}{3}\, q^{2-\alpha }\, \frac{\upsilon (q)}{c}\, \nonumber\\ && \times \left[ \mathcal {B}_{\frac{1}{1+q^2}} \left(\frac{\alpha -2}{2},\frac{3-\alpha }{2}\right) \right]^{-1} \simeq 1.37, \end{eqnarray} \tag{ A14 }$$

where we used typical values for M87 of α = 2.26 and q = 0.58 as derived in Section 2.2.

The ratio of line-to-bremsstrahlung cooling at 1 keV for a gas with solar abundance is approximately η ≃ 3.6 (see Figure 8 of Sutherland & Dopita 1993), assuming collisional ionization equilibrium. We can now rewrite Equation (A11) to get a constraint for the temperature slope of the cooling function χ for which the gas can become thermally unstable,

$$\begin{equation} \chi < -\frac{1}{2\eta } + \frac{\left(2-\beta \right)(1+\eta)}{\left(1+X_\mathrm{CR}\gamma \right)\eta }\simeq 0.13, {\rm for 1 keV and }Z=Z_\odot. \end{equation} \tag{ A15 }$$

Hence, if the cooling function is sufficiently flat or even inverted (as it can be the case in the line cooling regime), we have thermal instability. In the bremsstrahlung regime, on the contrary, the gas is in thermal equilibrium, since the criterion for instability

$$\begin{equation} \delta > \frac{X_\mathrm{CR}\gamma -1+\beta }{1+X_\mathrm{CR}\gamma }\simeq 0.8 \end{equation} \tag{ A16 }$$

cannot be fulfilled for CR Alfvén wave heating that does not depend on temperature (δ = 1/2).

So far, we assumed that the Alfvén speed is independent of radius (or density). Whether this is true depends in practice on the magnetic topology of a collapsing (or expanding) fluid parcel. Parameterizing $B\propto \rho ^{\alpha _B}$, we have α_B = 0 for collapse along $\boldsymbol { {B}}$ (implying $\upsilon _{A,\parallel } \propto \rho ^{-1/2}$) and α_B = 1 for collapse perpendicular to $\boldsymbol { {B}}$ (implying $\upsilon _{A,\perp } \propto \rho ^{1/2}$). By construction, we have at the center $c_s/\upsilon _A\simeq 4$. Hence, at the halo boundary of Virgo A, we obtain for the ratio of sound-to-Alfvén speed $c_s/\upsilon _{A,\parallel } =\sqrt{\beta \gamma _\mathrm{th}/2}\simeq 1.5$ and $c_s/\upsilon _{A,\perp }\simeq 26$ (where β is the plasma beta factor and γ_th = 5/3). This neglects a small-scale dynamo that could modify the B − ρ scaling furthermore. Hence, our discussed case of $\upsilon _A=\mathrm{const.}$ is the geometric mean of these extreme cases and implies that field components perpendicular and parallel to the direction of collapse are present. We normalize the CR heating at a radius of r₀ ≃ 1.5 kpc that corresponds to the radial extent of the central radio cocoon, which provides rather homogeneous conditions for field estimates based on the minimum energy conditions. Note that for simplicity, we hold the pressure gradient fixed while varying the Alfvén speed with radius, implicitly assuming a sufficiently large magnetic connectivity in radial direction. Hence, the CR heating rate at r_max = 35 kpc could be at most overestimated by a factor of $\upsilon _A/\upsilon _{A,0} = (\rho /\rho _0)^{-1/2}\simeq 3.8$, which is about the enhancement factor of the CR heating rate with $\upsilon _A=\mathrm{const.}$ over the cooling rate at the halo boundary for Z = 0.7 Z_☉ (see Figure 7). We can interpret this result twofold. Either we have $\mathcal {H}_\mathrm{CR}> \mathcal {C}_\mathrm{rad}$ and the halo boundary is currently adiabatically expanding in the cluster atmosphere, or we have global thermal equilibrium, i.e., $\mathcal {H}_\mathrm{CR}= \mathcal {C}_\mathrm{rad}$.

Zoom In Zoom Out Reset image size

How does a varying Alfvén speed impact on the local stability analysis? We reconsider Equation (A12) and adopt the density dependencies of $\upsilon _A$ discussed in the previous section, i.e., α_B = {0, 0.5, 1}, to obtain

$$\begin{eqnarray} \mathcal {H}_\mathrm{CR} &=& - {\boldsymbol{\upsilon}} _A {{{\cdot }}} \boldsymbol{\nabla } P_\mathrm{CR} \nonumber\\ &=& \Gamma _0 \rho ^\beta T^{-1/2+\delta } \propto \rho ^{\gamma +1/3\pm 1/2}, \mathrm{i.e.}, \nonumber\\ \beta -\gamma &=& \left\lbrace -\frac{1}{6}, \frac{2}{6}, \frac{5}{6}\right\rbrace , \end{eqnarray} \tag{ B1 }$$

where the first and last value of β corresponds to collapse parallel and perpendicular to $\boldsymbol { {B}}$, respectively. Following the discussion surrounding Equation (12), we can understand the character of local thermal stability by comparing the (modified) slope of the CR heating rate to that of the radiative cooling rate. We find that the slope of the parallel CR heating rate is always steeper than that of the radiative cooling rate for T > 2 × 10⁵ K, which implies local instability throughout this temperature regime. Note that this case is characterized by $\mathcal {H}_\mathrm{CR}> \mathcal {C}_\mathrm{rad}$ (Figure 7) so that the instability analysis formally does not apply and would require lowering X_CR until $\mathcal {H}_\mathrm{CR}= \mathcal {C}_\mathrm{rad}$.

In contrast, the "islands of stability" get more prominent for perpendicular collapse in comparison to the standard case of $\upsilon _A=\mathrm{const.}$, and even a new "island of stability" emerges out of the "ocean of instability" at around 5 × 10⁵ K (left panel of Figure 8). It is interesting to note that the latter case either represents the quiescently saturating state of the heat-flux-driven buoyancy instability (Quataert 2008; Parrish & Quataert 2008) or could be obtained through radial collapse, which amplifies the azimuthal components of the magnetic field. Since these processes are likely to generate or amplify a sufficiently large azimuthal field component, the realization of the pure radial instability case has a small probability.

Zoom In Zoom Out Reset image size

In the right panel of Figure 8, we address how a varying Alfvén speed changes the critical instability length scale (λ_crit) that is obtained by balancing the CR heating rate with the radiative cooling rate. On scales λ < λ_crit, CR heating dominates over cooling and wipes out temperature inhomogeneities. Gas parcels can become thermally unstable on scales λ > λ_crit provided that the unstable wavelengths can be supported by the system. For a gas parcel of fixed metallicity and magnetic connectivity outside the center, λ_crit decreases for collapse perpendicular to $\boldsymbol { {B}}$ ($\upsilon _A \propto \rho ^{1/2}$), i.e., smaller regions can become unstable and start to cool until they are stabilized at T = 10⁷ K (see left panel of Figure 8). In the opposite case of collapsing gas along $\boldsymbol { {B}}$ ($\upsilon _A \propto \rho ^{-1/2}$) we have formally local instability for T > 2 × 10⁵ K. However, for gas of metallicity Z = 1.3 Z_☉ and a magnetic suppression factor f_s = 0.3 of the heating rate, cooling is nevertheless stabilized for r ≳ 10 kpc since the unstable wavelengths cannot be supported by the system.