Buoyancy Instabilities in Galaxy Clusters: Convection Due to Adiabatic Cosmic Rays and Anisotropic Thermal Conduction
Abstract
Using a linear stability analysis and two and threedimensional nonlinear simulations, we study the physics of buoyancy instabilities in a combined thermal and relativistic (cosmic ray) plasma, motivated by the application to clusters of galaxies. We argue that cosmic ray diffusion is likely to be slow compared to the buoyancy time on large length scales, so that cosmic rays are effectively adiabatic. If the cosmic ray pressure is of the thermal pressure, and the cosmic ray entropy (; is the thermal plasma density) decreases outwards, cosmic rays drive an adiabatic convective instability analogous to Schwarzschild convection in stars. Global simulations of galaxy cluster cores show that this instability saturates by reducing the cosmic ray entropy gradient and driving efficient convection and turbulent mixing. At larger radii in cluster cores, the thermal plasma is unstable to the heat fluxdriven buoyancy instability (HBI), a convective instability generated by anisotropic thermal conduction and a background conductive heat flux. The HBI saturates by rearranging the magnetic field lines to become largely perpendicular to the local gravitational field; the resulting turbulence also primarily mixes plasma in the perpendicular plane. Cosmicray driven convection and the HBI may contribute to redistributing metals produced by Type 1a supernovae in clusters. Our calculations demonstrate that adiabatic simulations of galaxy clusters can artificially suppress the mixing of thermal and relativistic plasma; anisotropic thermal conduction allows more efficient mixing, which may contribute to cosmic rays being distributed throughout the cluster volume.
Subject headings:
convection — cooling flows — galaxies: active — galaxies: clusters: general — magnetic fields1. Introduction
Microscopic transport of heat and momentum in dilute plasmas, like those in clusters of galaxies, is primarily along magnetic field lines (Braginskii, 1965). This anisotropic transport dramatically affects the convective stability of the plasma; convective stability is no longer determined by the entropy gradient (Schwarzschild, 1958). Instead, a plasma is unstable to buoyant motions irrespective of the background entropy and temperature gradients (Balbus, 2000; Quataert, 2008). Cosmic rays diffusing along magnetic field lines also affect the convective stability of the plasma (Chandran & Dennis, 2006; Dennis & Chandran, 2008). Nonlinear simulations show that these instabilities driven by anisotropic thermal and cosmic ray transport can change the magnetic field configuration, and the background temperature and density profiles in the plasma, but they do not drive efficient convection (e.g., Parrish & Stone, 2007; Parrish & Quataert, 2008; Parrish, Stone, & Lemaster, 2008; Sharma, Quataert, & Stone, 2008). The instabilities saturate largely by rearranging the magnetic field configuration, thereby slowing down the instability and reaching a state of marginal stability to linear perturbations. By contrast, hydrodynamic convection in stellar interiors redistributes energy efficiently to make the plasma nearly adiabatic and thus marginally stable to convection.
One of the key astrophysical motivations for studying the transport properties of dilute plasmas in the presence of cosmic rays is to understand the dynamical and thermal structure of clusters of galaxies. The radiative cooling time ( Gyr) is much less than the Hubble time ( Gyr) in cluster cores. Thus, it was expected that the intracluster medium (ICM) would cool rapidly, resulting in large rates ( yr) of mass cooling to form cold gas and stars (e.g., Fabian, 1994). However, Xray observations have failed to detect copious emission from the expected cold plasma component in cluster cores (e.g., Peterson et al., 2003). The lack of cooling flows implies that cooling is balanced by some source of heating, e.g., heating by thermal conduction from large radii (Bertschinger & Meiksin, 1986; Zakamska & Narayan, 2003), heating by jets and bubbles blown by a central AGN (Binney & Tabor, 1995; Ciotti & Ostriker, 2001), or heating by cosmic rays (e.g., Rosner & Tucker, 1989; Loewenstein, Zweibel, & Begelman, 1991; Chandran & Rasera, 2007). Although thermal conduction may operate at large radii, it appears that the plasma at small radii must be heated by a feedback process which efficiently selfregulates. This is required to avoid the fine tuning of thermal conductivity required in models that include only conduction (e.g., Guo & Oh, 2008; Conroy & Ostriker, 2008).
Cosmic rays from a central AGN have been invoked to prevent catastrophic cooling of the plasma in cluster cores, either directly via Alfvén waves driven by cosmic rays heating the plasma (e.g., Loewenstein, Zweibel, & Begelman, 1991; Guo & Oh, 2008), or indirectly via convection driven by cosmic rays, the dissipation of which heats the plasma (Chandran & Rasera, 2007; Rebusco et al., 2006). Although there is ample evidence for the presence of cosmic rays in radio emitting bubbles in clusters (e.g., Bîrzan et al., 2004), it is unclear how/whether cosmic rays can be spread throughout the cluster volume at a sufficient level for these heating mechanisms to work. Simple hydrodynamic jets do not couple their energy to most of the ICM and instead simply drill through it, without heating and without transporting cosmic rays (and metals) throughout the ICM (Vernaleo & Reynolds, 2006).
In this paper we show that cosmic rays, which are likely to be centrally concentrated in clusters, can drive efficient convection and mixing if the cosmic ray pressure is not negligible compared to the plasma pressure (); this is likely the case in and around radio bubbles. We argue that on the large scales that likely dominate the turbulent dynamics in the ICM, cosmic rays are effectively adiabatic rather than diffusive (i.e., the cosmic ray diffusion time is longer than the buoyancy time). As a result, the cosmic rays can drive a Schwarzschildlike adiabatic convective instability. We present a linear analysis demonstrating that, while magnetic reorientation can shut off diffusive (“isobaric”) cosmic ray instabilities, it cannot shut off the adiabatic buoyancy instability driven by a negative cosmic ray entropy gradient. We then present twofluid (plasma and cosmic rays) numerical simulations with thermal conduction and cosmic ray diffusion along magnetic field lines.
We do not include plasma cooling in this paper, nor do we include the heating of the thermal plasma that arises from the excitation of short wavelength Alfvén waves by streaming cosmic rays. These are both significant omissions and preclude our results from being an accurate representation of the plasma in cluster cores. However, the main focus of this paper is not to solve the cooling flow problem per se, but rather to isolate and understand the transport and turbulence properties of cluster plasmas with realistic physics (e.g., anisotropic conduction, convection, and cosmic rays). By neglecting cooling, our calculations implicitly assume that some unspecified source of heating is preventing the rapid cooling of the ICM. A study of cluster cores with cooling and anisotropic conduction will be presented in a separate paper.
The remainder of this paper is organized as follows. In §2 we present the basic equations used in our analysis and derive the dispersion relation for linear buoyancy waves in the presence of thermal plasma and cosmic rays. In §3 we describe our numerical simulations and show that a steep cosmicray entropy gradient drives convective motions at small radii in cluster cores, and that anisotropic thermal conduction drives convection at intermediate radii that rearranges the magnetic field structure in clusters. Readers interested in just the numerical results can skip the linear stability calculation in §2. In §4 we summarize and discuss the implications of our work.
2. Basic Equations
To describe cosmic rays and thermal plasma in the ICM, we use the twofluid model of Drury & Völk (1981), modified to include anisotropic transport, gravitational acceleration (, where is the gravitational potential), and a cosmic ray source term to build up cosmic ray pressure. The equations of this model are
(1) 
(2) 
(3) 
(4) 
and
(5) 
where is the Lagrangian time derivative,
(6) 
is the heat flux, is the plasma temperature,
(7) 
is the diffusive flux of cosmicray energy (multiplied by ), is the cosmic ray energy source term, is the mass density, is the common bulkflow velocity of the thermal plasma and cosmic rays, is the magnetic field, , and are the thermalplasma and cosmicray pressures, is the parallel thermal conductivity, is the diffusion coefficient for cosmicray transport along the magnetic field, and =5/3 and =4/3 are the adiabatic indices of the thermal plasma and cosmic rays, respectively.
As mentioned in §1, we do not include radiative cooling in the energy equation (eq. [4]). In addition, we do not include the effects of cosmic ray streaming relative to the thermal plasma: in particular we neglect Alfvén wave heating of the thermal plasma and Alfvén wave streaming in the cosmic ray energy equation (e.g., Loewenstein, Zweibel, & Begelman, 1991). This physics will be included in the future together with plasma cooling. For the present paper we focus on the physics of convective instabilities in clusters using a simplified but physically reasonable model.
2.1. Linear Stability Analysis
We take all quantities to be the sum of an equilibrium value plus a smallamplitude fluctuation: , etc. We take the equilibrium velocity to vanish, set ( for the cluster, and is chosen such that local magnetic field lies in the plane), and take , , and to be functions of alone. We employ a local analysis in which all fluctuating quantities vary as with , where is the scale on which the equilibrium quantities vary. We consider the limit and work in the Boussinesq approximation, ( is the sound speed). We also do not include the perturbed cosmic ray source term (eq. [5]) in our linear analysis since its form is uncertain.
In terms of the plasma displacement, , the perturbed magnetic field () can be combined with equation (2), to give
(8) 
where is the totalpressure perturbation. Dotting equation (8) with and using the nearincompressibility condition (), we find to leading order in that
(9) 
and
(10) 
where , and is the square of the Alfvén speed.
In terms of the perturbation to the magneticfield unit vector, , the perturbed heat flux can be written as , where we have dropped the term involving the perturbed thermal conductivity since it is smaller than the other terms by a factor of . Taking the dot product of the heat flux with we find that
(11) 
where and is the component of . Using equations (4) and (11), and using , we find that
(12)  
where is the rate at which thermal conductivity smoothes out temperature fluctuations along the magnetic field, and
is the square of the BruntVäisälä frequency.
In the same way that we obtained equations (11) and (12), for the cosmic rays we find that
(13)  
where is the rate at which diffusion smoothes out variations in along the magnetic field, and
is the square of the BruntVäisälä frequency associated with the cosmic ray pressure.
Adding equations (12) and (13), making use of equation (9), and noting that the righthand side of equation (9) is much smaller than the individual terms on the righthand sides of equations (12) and (13), we find that
(14) 
where
,
and
Upon substituting equation (14) into equation (10), we obtain an equation for the plasma displacement alone,
(15)  
We set , where and . After taking the dot product of equation (15) with and , we obtain two equations which can be written in matrix form as
(16) 
Setting the determinant of the matrix on the lefthand
side of equation (16) equal to zero, we obtain the dispersion
relation
,
where
is a factor that depends on the magnetic field geometry and wavenumber, , and is the component of .
While one solution to equation (16), , describes waves unaffected by buoyancy, corresponds to the modes modified by buoyancy,
(17) 
In the limit that , equation (17) reduces to equation (13) of Quataert (2008). For nonzero , we consider two limiting cases of equation (17). The first is a highly diffusive limit for cosmic rays, in which . In this case, equation (17) reduces to
(18) 
This is a generalization of the HBI (heatflux driven buoyancy instability; Quataert 2008) and MTI (magnetothermal instability; Balbus 2000) in the presence of diffusive cosmic rays (e.g., Chandran & Dennis, 2006; Dennis & Chandran, 2008).
The second limit we consider is that of adiabatic cosmic rays, with . In this case, the cosmic rays diffuse slowly compared to the buoyancy time and behave nearly adiabatically; equation (17) then reduces to
(19) 
This implies that if the cosmic ray pressure is significant and if cosmic ray entropy () is sufficiently peaked and decreasing outward, the plasma will become unstable to a Schwarzschild type buoyancy instability. The entropy of cluster plasma is stably stratified according to the Schwarzschild criterion. The thermal plasma response at small scales is nonetheless governed by the temperature gradient and not the entropy gradient. For typical cluster parameters, although the global conduction timescale may be longer than the buoyancy timescale (see Fig. 1), anisotropic conduction determines the local buoyant response (i.e., for ) via HBI/MTI depending on the sign of temperature gradient. The adiabatic cosmicray instability (ACRI) described by equation (19) is different from the MTI and HBI in that it cannot saturate by magnetic field reorientation (the cosmicray driving in eq. [19] does not depend on the field configuration term that shows up in the thermal driving). It must lead to vigorous convection that changes the cosmic ray entropy profile to become nearly adiabatic.
2.2. The Cosmic Ray Diffusion Coefficient
A cosmic ray particle streaming through a magnetized plasma is efficiently scattered in pitch angle by magnetic fluctuations with wavelengths comparable to its Larmor radius. An unavoidable source of magnetic fluctuations is the selfexcited streaming instability (e.g., Kulsrud & Pearce, 1969). The effect of pitchangle scattering due to Alfvén waves is that the bulk speed of cosmic rays relative to the thermal plasma is close to the Alfvén speed, i.e., , where is the cosmic ray drift velocity relative to the thermal plasma, is pitchangle scattering rate, and is cosmic ray gradient scale (Kulsrud, 2005). The cosmic rays stream along the magnetic field direction, and down the cosmic ray pressure gradient. In addition to streaming with Alfvén wave packets, cosmic rays also undergo momentumspace diffusion, which leads to spatial diffusion along the field lines with . Thus, if cosmic ray scattering is efficient and , the diffusion timescale over scales comparable to is much longer than the Alfvén crossing time.
With the above model for cosmic ray scattering one can selfconsistently calculate the cosmic ray diffusion coefficient in terms of the plasma parameters (e.g., Loewenstein, Zweibel, & Begelman, 1991). In the present calculations we do not explicitly include the effects of cosmic rays streaming with respect to the plasma. Instead, the diffusion coefficient in equation (7) should be interpreted as an effective diffusion coefficient, taking into account both microscopic diffusion and streaming along turbulent magnetic field lines. The cosmic ray “diffusion” time due to cosmic rays streaming along random magnetic field lines can be crudely bounded by the Alfvén crossing time (). Together with the fact that the true microscopic diffusion time is much longer than the Alfvén crossing time (as argued above), this motivates our choice of for the cosmicray diffusion coefficient in most of our numerical simulations, where is a factor of order unity. Somewhat arbitrarily, we take , but our results are insensitive to so long as .
This estimate of () is consistent with the measured Galactic cosmic ray diffusion coefficient for GeV particles, cms (Berezinskii et al., 1990), using typical values for magnetic field strength and cosmic ray scale height in the Galaxy. However, at higher energies the diffusion coefficient increases as , where is the cosmic ray energy (e.g., see Engelmann et al. 1990). For a cosmic ray energy distribution function steeper than (in the Milky Way it scales as from 1 to GeV), the cosmic ray pressure in equation (2) will be dominated by the lowest energy cosmic rays, and a diffusion coefficient seems appropriate for the fluid description of cosmic rays considered here. With this choice, the ratio of the buoyancy timescale () to the cosmic ray diffusion timescale () is . In clusters (Govoni & Feretti, 2004) so that we expect cosmic rays to be adiabatic on relatively large length scales and thus to be susceptible to cosmic ray driven convection when their entropy gradient is sufficiently large. Our linear stability analysis differs from that of Chandran & Dennis (2006) and Dennis & Chandran (2008), in that we allow for a background heat flux. In addition, while they obtained the dispersion relation allowing the thermal plasma and cosmic rays to be either both in the diffusive limit or both in the adiabatic limit, we have argued that the more relevant case in clusters is likely that of diffusive thermal plasma and adiabatic cosmic rays. Because of uncertainties in , we have carried out numerical simulations for different values of . We find that for = 10 cm s, the cosmic rays are effectively adiabatic at all scales, and even for as high as cm s, the cosmic rays are adiabatic on large scales, kpc (see sec. 3.2.4).
3. Numerical Simulations
We have extended the ZEUSMP MHD code (Hayes et al., 2006; Stone & Norman, 1992a, b) to include thermal conduction along magnetic field lines (Sharma & Hammett, 2007), and have added cosmic rays as an additional fluid diffusing along magnetic field lines. We numerically solve equations (1)(7). As mentioned earlier, we do not include plasma cooling. The cosmic ray energy equation, thermal conduction in equation (4), and the cosmic ray pressure gradient in the equation of motion are implemented in an operator split fashion with appropriate source and transport terms. Both thermal conduction and cosmic ray diffusion are subcycled.
We have tested the code extensively. The anisotropic cosmic ray diffusion equation is analogous to anisotropic thermal conduction. We use the method of Sharma & Hammett (2007) which preserves positivity of . Appendix shows a 1D shock tube test, adapted from Pfrommer et al. (2006), which shows that adiabatic evolution of cosmic rays is accurate. The thermal conductivity is chosen to be the Spitzer value,
(20) 
Based on the discussion in §2.2 we choose for most of our calculations, where is the local Alfvén speed; to test the dependence on , we also carry out simulations with a constant value of cm s.
How cosmic rays are produced and distributed in the ICM is still poorly understood (e.g., Eilek, 2003). Thus, we use a simple phenomenological source term to drive the cosmic ray pressure in the inner parts of the cluster. The cosmic ray energy source term in equation (5) is based on Guo & Oh (2008),
(21) 
We take kpc and , which leads to a centrally
peaked cosmic ray entropy, and take the cosmic ray energy injection
rate to be erg
s.
The cosmicray pressure fraction and its dependence on radius are not that wellconstrained observationally. For a few clusters (e.g., Virgo, Perseus, and Fornax) observational constraints indicate that averaged over the cluster core (Pfrommer & Enßlin, 2004; Churazov et al., 2008). It is possible, however, that larger cosmicray pressure fractions arise in clusters in which AGN feedback is expected to play a particularly strong role, such as Hydra A or Sersic 15903 (see, e.g., Zakamska & Narayan 2003; Chandran & Rasera 2007); cosmicrays may also be particularly important at small radii, and quite likely contribute significantly to the pressure in and around radio bubbles (which are associated with a deficit of Xray emission from the thermal plasma; e.g., Bîrzan et al. 2004). In our simulations, we choose parameters in equation (21) such that the cosmic ray pressure is smaller than the plasma pressure even at late times. For a larger , we find that the cosmic ray pressure builds up faster than the rate at which cosmic ray pressure can be transported outwards via convection (or diffusion, but the latter is slow), and the cosmic ray pressure can become larger than the plasma pressure. In calculations with cooling, it is likely that a larger cosmic ray injection rate could remain consistent with : if the gas is allowed to cool then a large part of the cosmic ray energy may be channeled into plasma heating (and thus cooling) without building up a larger cosmic ray pressure.
To assess how effectively turbulence generated by the HBI or ACRI mixes the plasma, we solve for the advection of a passive scalar density (e.g., a proxy for metallicity) , using
(22) 
Appendix A shows the behavior of the passive scalar density in a shock tube test.
3.1. Simulation Parameters
The simulations are carried out in spherical (,,) geometry with the inner boundary at kpc and the outer boundary at kpc. Strict outflow boundary conditions are applied to the radial velocity at the inner and outer radial boundaries. The plasma pressure and density are held fixed at the outer boundary to prevent spurious oscillations; the plasma cooling time at is longer than the Hubble time. All other plasma and cosmic ray variables are copied on ghost zones at both the inner and outer radial boundaries. Reflective boundary conditions are applied at the boundaries (). Periodic boundary conditions are applied in the direction. A logarithmic grid is chosen in the radial direction, while the grid is uniform in and . Our fiducial run uses a grid, with and . With these choices, , , and . We also carried out a number of 2D (axisymmetric) simulations, and one higher resolution, , 2D simulation for convergence studies.
We typically initialize a weak ( everywhere) splitmonopole magnetic field with , although in one calculation (CRM), we use a monopole field to compare simulations with and without net magnetic flux. In runs with cosmic rays, the initial cosmic ray pressure is small (0.005 times the plasma pressure at ) and varies as ; the cosmic ray pressure builds up in time via the source term in eq. [5]). The initial pressure is chosen such that the plasma is in dynamical equilibrium (). Initial thermal equilibrium is not imposed.
As in Guo & Oh (2008), we use cluster parameters relevant for Abell 2199. The gravitational potential () is the sum of the dark matter potential (),
(23) 
where is the characteristic dark matter mass, kpc is the scale radius (Navarro, Frenk, & White, 1997), and the potential due to the central cD galaxy (),
(24) 
where kpc, g cm (see Kelson et al. 2002 for cD galaxy NGC 6166).
The ideal gas law is used with , where () is the mean molecular weight per thermal particle (electron) and () is the total (electron) number density. We assume a fully ionized plasma with hydrogen mass fraction , helium mass fraction , such that and . Since we do not include cooling, the metallicity appears only in the conversion of pressure into plasma temperature used to calculate the thermal conductivity in equation (20).
To roughly match the observations (Johnstone et al., 2002), the initial temperature increases linearly from 1.5 keV at to 4.6 keV at (see Fig. 2); the electron number density is fixed to 0.0015 cm and the temperature is fixed to keV at at all times. The initial density, obtained from imposing hydrostatic equilibrium, is quite a bit larger than the observed density at . This is because the density obtained from hydrostatic equilibrium is extremely sensitive to the form of the temperature profile, for which we use a simple linear fit. However, since we do not include cooling, the large inner density does not significantly affect our results.
3.2. Results
Label  Dim.  initial B  CR angle  

CR  3D  splitM  90  
NCR  3D  splitM  0  0 
CRM  3D  monopole  90  
CR2D  2D  splitM  90  
CR2Ddbl  2D  splitM  90  
CR28  2D  splitM  90  
CR29  2D  splitM  90  
CR30  3D  splitM  30  
CR30ad  3D  splitM  30 
Table 1 summarizes the properties of our simulations. Although we list a number of calculations in Table 1, we focus on two 3D simulations: CR (the fiducial run) and NCR. Cosmic rays are not included in NCR; NCR thus serves as a control run that allows us to isolate the effects of cosmic rays. The aim of the rest of the simulations is to understand certain aspects of the physics in more detail. As described earlier, for all runs except CRM we initialize a splitmonopole magnetic field. For convergence studies, we carried out a two dimensional (axisymmetric) version of CR, CR2D, and compared it with a run with double the resolution, CR2Ddbl. All runs, except CR28 and CR29, use ; a fixed value for is chosen for CR28 and CR29 to test the influence of on our results. In runs CR30 and CR30ad, the cosmic ray source term is only applied within 30 of the pole to study angular diffusion of cosmic rays with (CR30) and without (CR30ad) thermal conduction along field lines; in the rest of the calculations, cosmic ray injection is spherically symmetric.
Our initial profiles are in dynamical equilibrium, but not in thermal equilibrium. The plasma will remain static if thermal conduction and cosmic rays are not included, since the magnetic field is very weak and the plasma is stably stratified according to the Schwarzschild criterion. However, in the presence of thermal conduction, the background temperature and density change in time; in addition, the magnetic field lines are reoriented by the HBI. When the cosmic ray source term is applied, the cosmic ray entropy gradient can drive convection and mixing due to the ACRI. In this section we discuss the influence of these effects on the structure of the ICM.
Figure 1 shows different timescales in the initial state: the isothermal sound crossing time (solid line) is the shortest timescale; the plasma cooling time (dotted line) is included for comparison with the other timescales, although cooling is not included in the simulations; the growthtime of the HBI (longdashed line) varies slowly as a function of radius; the cosmic ray injection timescale (shortdashed line), i.e., the timescale for the cosmicray source term to increase the cosmic ray pressure by an amount comparable to the plasma pressure, increases rapidly with radius since the cosmic ray source term is centrally concentrated; the thermal conduction timescale (short dotdashed line) has a maximum at intermediate radii and is comparable to the HBI timescale at both and ; finally, the cosmic ray diffusion timescale (long dotdashed line; shown for cms) is shorter than the buoyancy timescale only within kpc. Since the diffusion time is , this implies that cosmic rays are effectively adiabatic for smaller , i.e., in all runs except CR29 (Table 1).
The Fiducial Run (CR)
We discuss the fiducial simulation (CR) in detail and compare it with the simulation without cosmic rays (NCR). Both simulations CR and NCR show similar properties for radii kpc, the radius outside of which the cosmic ray pressure is always small. Since there is no cooling to balance heating by thermal conduction in the initial state, the initial thermal properties will be modified on the conduction timescale.
Figure 2 shows the angle averaged temperature and density profiles as a function of radius for different times for runs CR and NCR. Figure 2 shows that for run NCR, the temperature becomes isothermal near the inner and outer radii due to thermal conduction along initially radial magnetic field lines. Although Figure 1 shows that the HBI timescale is always shorter than the conduction time, the two are of the same order at both and where thermal conduction makes the plasma isothermal before the HBI can grow significantly. The HBI, which is active within 30 kpc 100 kpc, reorients the magnetic field lines to be primarily perpendicular to the radial direction (see Fig. 3). This creates a thermal barrier and a large temperature gradient at these radii. The formation of such a thermal barrier is not forced by the boundary condition since the temperature at is floating; if the conduction timescale were shorter than the HBI timescale at all radii, the plasma would become isothermal (4.6 keV) at all radii. Run CR shows similar behavior at larger radii but differs substantially at small radii. Since the cosmic ray pressure provides a substantial fraction of the pressure support to balance gravity within 30 kpc, the plasma density and pressure at small radii in run CR is substantially smaller than in NCR. The plasma temperature is similar in magnitude for CR and NCR but shows a minimum at 20 kpc for CR at late times.
Figure 4 shows the cosmic ray entropy profile (top) and the ratio of cosmic ray pressure to plasma pressure (, bottom) as a function of radius for different times; we define cosmic ray entropy as since this is the quantity whose gradient determines the buoyant response of adiabatic cosmic rays (see, e.g., eq. [19]). Since the cosmic ray source term is chosen to be a strong function of radius ( for kpc), it drives convection due to the ACRI when the cosmic ray pressure builds up and becomes comparable to the plasma pressure. At later times, cosmic ray injection does not continuously increase the inner cosmic ray pressure with time; instead, a cosmic ray driven convection front spreads radially outwards. Convection drives the cosmic rays to be adiabatic ( constant) in regions where the cosmic ray pressure is not negligible compared to the plasma pressure.
Figure 5 shows 2D contour plots of turbulent velocity (absolute value of velocity) in the plane at 9 Gyr for CR and NCR. Velocities are similar for kpc where cosmic ray pressure is negligible and fluid motions are driven by the HBI. The maximum turbulent velocity in the HBIdominated regions is km s. Run CR shows much larger turbulent velocities ( 100 km s) for kpc; note that the turbulent velocities induced by the ACRI are consistent with mixing length theory, with , where is the power carried by convection and is the resulting convective velocity. The large turbulent velocities in the presence of the ACRI are sufficient to prevent catastrophic cooling according to the models of Chandran & Rasera (2007). By contrast, the turbulent velocities are extremely small at kpc for NCR because the plasma temperature gradient is wiped out by conduction before the HBI can drive any turbulence.
While the turbulent velocity vectors are roughly isotropic at kpc for CR, they are aligned primarily perpendicular to the radial direction at large radii where the HBI dominates. This is consistent with the result that while the HBI saturates by reorienting magnetic field lines perpendicular to gravity (e.g., Parrish & Quataert, 2008), the ACRI drives roughly isotropic convection irrespective of the magnetic field geometry.
Figure 3 shows the timeaveraged (from to ) angle between the magnetic field vector and the radial direction as a function of radius for runs CR and NCR; the angle is defined with respect to the radial direction, such that it is for the initial splitmonopole field. The average angle is similar for kpc, where the cosmic ray pressure is negligible, for both simulations. For run CR, radii kpc are convectively stirred by the ACRI and the average angle between the magnetic field and the radial direction is , close to the value expected for a uniform, random magnetic field unit vector (). For both CR and NCR at 30 kpc kpc, the magnetic field is nearly perpendicular to the radial direction because of the HBI. The average angle between the magnetic field unit vector and the radial direction at these radii is . For kpc and for kpc in run NCR, the HBI is weak since thermal conduction wipes out the temperature gradient before the HBI can grow significantly; as a result, the field lines are not perpendicular to the radial direction.
Figure 6 shows angle and time averaged (from to ) kinetic and magnetic energy profiles as a function of radius for runs CR, CRM, and NCR. Both the magnetic and kinetic energies are generally amplified, although to varying degrees. The kinetic energy is very effectively amplified at small radii by the action of the ACRI in run CR; the turbulent Mach number is (see also Fig. 5). The magnetic energy as a function of radius in CR is quite striking in that there is no magnetic energy enhancement at small radii: there is a slight bump in magnetic energy at 20 kpc but the final magnetic energy is smaller than the initial magnetic energy for kpc. In turbulent dynamos one often finds the turbulent magnetic energy to be of the order of turbulent kinetic energy (e.g., Cho & Vishniac, 2000); this is clearly not the case for the run CR. Magnetic field amplification can, however, be subtle; e.g., in local shearing box simulations of the magnetorotational instability with no net magnetic flux, magnetic field amplification (and associated stress and turbulence) occurs only for Prandtl numbers exceeding unity (e.g., Fromang et al., 2007). We do not have explicit viscosity and resistivity and it is possible that the effective Prandtl number in our simulation is small, so that dissipation of the initially split monopolar field at the grid scale dominates over magnetic field enhancement by convective turbulence. To better understand this, we have done a simulation with cosmic rays with a net initial magnetic flux (run CRM), but with all other properties of the simulation the same. Figure 6 shows that in this case, the magnetic energy density is much larger in the inner regions as compared to CR, although the magnetic energy is still times smaller than the kinetic energy. In this case, the magnetic field strength in the center of the cluster is amplified to 0.11 G by convective motions. The kinetic energy density profiles are almost identical for CR and CRM, indicating that the properties of the convection are not very different in the two cases, although the efficiency of magnetic field amplification differs dramatically.
At larger radii (30 kpc 100 kpc), the HBI causes amplification of the magnetic and kinetic energies in both CR and NCR. The magnetic energy is enhanced by a factor at kpc. This level of field amplification – a factor of , primarily of the and components – is what is required to reorient initially radial magnetic fields into fields that are perpendicular to the radial direction (as was seen in previous HBI and MTI simulations; Parrish & Quataert 2008; Parrish & Stone 2007; Sharma, Quataert, & Stone 2008). The kinetic energy is also enhanced at these radii because of turbulence driven by the HBI.
In order to study the numerical convergence of our fiducial run, we
carried out an axisymmetric 2D run analogous to CR – CR2D – and an
axisymmetric run with double the resolution – CR2Ddbl (see Table
1); studying convergence directly with the 3D run CR
would have been computationally prohibitive, requiring
times more cpu time. The results from runs CR2D and CR2Ddbl are
nearly identical to each other; in particular, angle averaged plots
such as those shown in Figures 2, 4,
3, and 6 are quite similar in the two
cases. We have not included these 2D results in every figure, but to
illustrate the basic convergence result, Figure 3 shows
that the average angle between the magnetic field unit vector and the
radial direction is quite similar for CR2D and CR2Ddbl (which are
both similar to the 3D run CR). It is interesting to note that the
2D runs differ from the 3D run CR in one important respect: the
amplification of the magnetic field at small radii is significantly
larger in 2D than in 3D; the magnetic field energy density at small
radii in the 2D runs is comparable to the 3D run with a net magnetic
flux (CRM) in Figure 6.
Diffusion of a Passive Scalar
Figure 7 shows , the passive scalar density, in the inner 40 kpc for CR and NCR at 9 Gyr, as well as the projection of the magnetic field unit vectors. The passive scalar density is initialized to be a large number () for kpc (corresponding to four radial zones) and is negligible () for kpc. The goal of initializing a passive scalar is to study mixing due to turbulence. Observations of clusters reveal a metallicity distribution that is more spatially extended than the light distribution of the central galaxy. This may indicate turbulent transport of metals in clusters (e.g., Rebusco et al., 2006; Rasera & Chandran, 2008). As expected, run CR with large turbulent velocities at kpc also results in efficient mixing. Mixing is negligible for NCR because the inner radii ( kpc) are isothermal and are thus not stirred by the HBI. For a direct comparison with observations, one must include a time dependent, spatially distributed source term in the passive scalar equation which represents metal enrichment due to Type Ia supernovae (e.g., Rebusco et al., 2005); this is beyond the scope of the present paper. Nonetheless, our results indicate that cosmicray driven convection is an efficient mechanism for mixing plasma in clusters of galaxies.
Figure 8 shows the angle averaged passive scalar density as a function of radius for run CR at 1/3, 1, 3, 9 Gyr. Also shown is a Gaussian fit (at 9 Gyr) with a diffusion coefficient of cms. The passive scalar density at 9 Gyr is flatter than the Gaussian fit at 20 kpc, implying that diffusion due to convection driven by cosmic rays corresponds to an effective diffusion coefficient somewhat larger than cms. Beyond kpc, the cosmic ray pressure is unimportant, and turbulence is driven by the HBI. This change in the source of the turbulence accounts for the rapid decline in the passive scalar density at large radii in Figure 8.
To isolate just the mixing induced by turbulence driven by the HBI, we have taken run NCR at Gyr) and initialized a passive scalar density peaked at kpc, the radius where the temperature gradient is large and the HBI is active (see Fig. 2). The passive scalar is initialized to be at two grid points at kpc, and , and negligible () elsewhere. Figure 9 shows snapshots of the passive scalar density at later times. The passive scalar diffuses more rapidly in the direction. This is because the turbulent velocities due to the HBI are larger in the direction perpendicular to gravity than they are in the radial direction (just as the magnetic field components perpendicular to gravity are preferentially amplified). To estimate the diffusion coefficient in the and directions, we compare how much the passive scalar has spread in the two directions; the full width at half maximum (FWHM) for along and at 9 Gyr is , 40 kpc, respectively. For comparison, the FWHM at 9 Gyr for for run CR shown in Fig. 8 is kpc. Thus the diffusion coefficient due to the HBI alone is (perpendicular to gravity) and (parallel to gravity) times smaller than the diffusion coefficient due to the ACRI in run CR. Although these precise numerical values likely depend on the detailed parameters of our simulations, the fact that the HBI primarily induces turbulence and mixing in the plane perpendicular to gravity is a generic result.
Heat Flux modified by the HBI and ACRI
Many 1D models of clusters parameterize thermal conduction by its ratio to the Spitzer value (Zakamska & Narayan, 2003; Chandran & Rasera, 2007; Guo & Oh, 2008). However, our simulations show that, because of plasma instabilities that operate in clusters (e.g., HBI and ACRI), a reduction of the conductivity by a fixed factor is not applicable (see also Parrish & Stone 2007; Parrish & Quataert 2008; Sharma, Quataert, & Stone 2008). At large radii in cluster cores, where the cosmic ray pressure is negligible, the HBI can orient field lines perpendicular to the radial direction, but at small radii where cosmic rays can be significant, the magnetic field may be significantly more radial. For example, for run CR the average angle of the magnetic field relative to the radial direction for kpc is (see Fig. 3), corresponding (roughly) to a reduction factor of . For 30 kpc kpc, however, the turbulence is dominated by the HBI, and the average angle between the magnetic field and the radial direction is , corresponding to a reduction factor of .
An even more subtle result is that the HBI can change the background temperature gradient by forming thermal barriers (see Fig. 2), which will be absent with isotropic conduction. For example, the temperature gradient at 30 kpc 100 kpc for CR and NCR at late times is times larger than the initial temperature gradient. Thus, the HBI not only reduces the conductive heating by a factor of , it also increases it by making the temperature gradient larger by a factor . Approximating the conductivity of a magnetized plasma by a constant factor with respect to the Spitzer value misses all of this interesting dynamics. Whether this is important in real clusters with radiative cooling and various sources of heating remains to be seen.
Runs with larger (Cr28 & Cr29)
We have also carried out 2D simulations with larger cosmic ray diffusion coefficients (), since the value of the cosmic ray diffusion coefficient is poorly constrained (see §2.2). Run CR28 uses cm s, the cosmic ray diffusion coefficient estimated for GeV cosmic rays in the Galaxy. Run CR29 uses cm s. All other parameters and initial conditions are same as the fiducial run. Since we are comparing these 2D simulations with CR, which is a 3D simulation, we have verified that the run CR2D gives results similar to the 3D results presented here; in particular, the profiles for and are identical to the profiles for CR shown in Figure 4.
Figure 10 shows the angle averaged cosmic ray entropy () and the ratio of the cosmic ray pressure to the plasma pressure () as a function of radius for CR28 and CR29. The profiles for CR28 and CR (see Fig. 4) look similar; cosmic ray entropy and pressure are slightly more radially spread out for CR28. The profiles for CR29 are quite different. The entropy and pressure ratio profiles have a peak at kpc; this is the radius beyond which the cosmic ray diffusion time is longer than the buoyancy timescale (see Fig. 1). The ratio is smaller in CR29 as cosmic rays are spread out over a larger volume by diffusion. The cosmic ray entropy increases outwards for kpc since cosmic ray diffusion, and not convection driven by cosmic rays, dominates the outward cosmic ray transport. This is different from CR and CR28 where cosmic ray diffusion is subdominant. Even in CR29, cosmic rays are effectively adiabatic for kpc; cosmicray driven convection is absent at these radii, however, because the cosmic ray pressure is not large enough to drive the ACRI (see Fig. 10). At smaller radii, the cosmic rays are nearly isobaric because of rapid diffusion, and the system is formally unstable to the CR mediated version of the MTI (eq. [18]). However, because the field lines are nearly radial at these radii, the growthrate of the CRMTI is quite slow and we do not see any indications that it develops in our simulations.
To explicitly study the possibility of the ACRI setting in at larger radii in the cluster core, we carried out a simulation with the cosmic ray source term (eq. [21]) three times larger than in run CR29. This larger source term increases the CR pressure at large radii, and at late times is large enough to drive the ACRI. The turbulent velocities are km s at kpc, where the cosmic rays are effectively adiabatic in spite of the large . Thus, even in the presence of rapid cosmic ray diffusion, the ACRI can set in at large radii where the cosmic rays are adiabatic, provided the cosmic ray pressure is sufficiently large; this may naturally occur in the vicinity of cosmicray filled buoyant bubbles. More generally, our results demonstrate that so long as , cm s, cosmic rays will behave effectively adiabatically throughout the cluster core and bulk transport by convection and other mechanisms will dominate the diffusive transport.
Runs with Cosmic ray Sources at the Poles (CR30 & CR30ad)
It is very unlikely that cosmic rays in clusters are injected spherically symmetrically. Instead, the injection likely occurs preferentially in the polar direction. To study the resulting physics in this case, we carried out 3D simulations in which the cosmic ray source term is applied only within of the pole: CR30 and CR30ad. Except for this difference all parameters for run CR30 are the same as run CR. Run CR30ad differs from CR30 in that the plasma is adiabatic, i.e., thermal conduction is not included. One of the aims of these simulations is to show the dramatic differences that result from including anisotropic thermal conduction (relative to a more typical adiabatic simulation). Cluster plasmas are observed to be stable to adiabatic convection because the entropy increases outwards (e.g., Piffaretti et al., 2005). However, convection in an anisotropically conducting plasma depends on the temperature gradient, and not the entropy gradient, and the system is unstable independent of the sign of the temperature gradient. This makes it much easier to mix a thermally conducting plasma than an adiabatic plasma. In clusters, this implies that turbulence produced by external means, e.g., the ACRI, wakes due to galaxy clusters, etc., may be an effective way of mixing the thermal plasma.
Figure 11 shows contour plots ( snapshot) of the ratio of cosmic ray pressure to plasma pressure () at 1/3, 1, 3, 6 Gyr, for CR30 and CR30ad. For CR30, the cosmic rays become unstable to the ACRI in the polar region. The resulting turbulence is able to convectively mix the plasma, not only in the unstable radial direction, but also in the marginally stable direction. Instead of cosmic rays being confined only to the cone, convection effectively mixes plasma in both the radial and angular directions. In addition, at radii beyond 20 kpc where the temperature gradient is appreciable (see Fig. 2), the HBI mixes material primarily in the direction, as seen by the oriented fingers at late times in Figure 11. For CR30, at 3 Gyr in the equatorial region for kpc. Although cosmic rays are still dominant near the pole, convection brings a nonnegligible amount cosmic rays into the equatorial region. In comparison, there is no sign of convective overshoot in CR30ad because convective motions in the thermal plasma near the equator are strongly stabilized by a large positive plasma entropy gradient. The polar cosmic ray dominated plasma does expand somewhat in both and as it becomes overpressured. However, the value of at 3 Gyr is everywhere in the equatorial plane for CR30ad.
4. Summary & Astrophysical Implications
The Xray emitting plasma in clusters of galaxies is hot (
keV) and dilute ( cm), so that the transport of
heat and momentum along magnetic field lines can be energetically and
dynamically important. In addition, jets launched by a central AGN
produce relativistic plasma (cosmic rays), which are observed in part
as bubbles of radio emission associated with deficits of thermal Xray
emission (Xray cavities; e.g., Bîrzan et al. 2004). In this paper we have
studied the transport properties of an ICM composed of cosmic rays and
thermal plasma. We have argued that cosmic ray diffusion is likely to
be slow because of scattering by selfgenerated Alfvén waves (§2.2); as a result, the cosmic rays are adiabatic on
moderately large length scales 110 kpc.
More concretely, cosmic rays are adiabatic on the scale of cluster
cores, so long as their parallel diffusion coefficient satisfies
cms.
It is now well established that anisotropic conduction and anisotropic cosmic ray diffusion can dramatically modify buoyancy instabilities in low collisionality systems, producing qualitatively new instabilities such as the MTI, HBI, and their cosmic ray counterparts (e.g., Balbus, 2000; Chandran & Dennis, 2006; Parrish & Stone, 2007; Quataert, 2008; Parrish & Quataert, 2008). Nonlinear studies of these instabilities (including those in this paper) have demonstrated that they saturate by approaching a state of marginal stability to linear perturbations, just as in hydrodynamic convection. However, in a magnetized plasma, there is an additional degree of freedom that is not present in hydrodynamic convection, namely the local direction of the magnetic field. The primary mechanism by which these diffusive buoyancy instabilities saturate is by rearranging the magnetic field lines, so that the linear growth rate becomes extremely small (see Fig. 3). This is different from entropy gradient driven convection in adiabatic fluids, which saturates by producing convection that wipes out strong entropy gradients. Even nonlinearly, most of the energy flux in systems unstable to the MTI and HBI is transported by thermal conduction, rather than convection. Moreover, the saturation of these instabilities is quasilinear in the sense that the saturated magnetic energy is proportional to the initial magnetic field energy (Sharma, Quataert, & Stone, 2008).
In this paper we have shown analytically and through numerical simulations that when cosmic rays have appreciable pressure, , and an outwardly decreasing entropy (), they can drive strong convection and mixing in clusters of galaxies. This adiabatic cosmic ray instability (ACRI) in the central regions of clusters of galaxies is a cosmic ray analogue of hydrodynamic convection familiar in the context of stars and planets. In particular, the nonlinear saturation of adiabatic cosmic ray convection is similar to that of hydrodynamic convection, and thus quite different from the saturation of the MTI and HBI. Our simulations of cluster cores also provide insight into the global saturation of the HBI at radii in clusters where the cosmic ray pressure is negligible ( kpc in our models), and thus where the only convective instability is that driven by the background conductive heat flux. More specifically, the primary results of this paper include:

If the cosmic ray entropy decreases outwards and if , convection driven by the ACRI sets in. In the saturated state, the cosmic ray entropy profile becomes nearly constant in the region with significant cosmic ray pressure (Fig. 4). The resulting turbulent velocities are consistent with mixing length theory, with km s where is the total power supplied to cosmic rays and is the pressure height scale of cosmic rays. The ACRI generates turbulent motions more effectively in cluster cores than the HBI alone.

The ACRI drives roughly isotropic convection with the average angle between the field lines and the radial direction ; by contrast, the HBI generates magnetic field lines that are primarily in the and directions (Fig. 3), shutting off the radial conduction of heat. The effective radial conductivity of a cluster plasma thus depends sensitively on which of these instabilities operates at a given location, and may not be adequately approximated as a fixed fraction of the Spitzer value throughout the cluster.

We have quantified the mixing of a passive scalar by the ACRI and HBI: the ACRI produces roughly isotropic mixing with a turbulent diffusion coefficient cms (Fig. 7); mixing length theory predicts that . At larger radii, only the HBI operates and the mixing is primarily in the and directions, rather than in the radial direction (Fig. 9). Both the ACRI and the HBI may contribute to mixing metals in clusters by redistributing, in both radius and angle, metals produced by Type 1a supernovae. Some observations of metallicity gradients in clusters have inferred mixing at levels comparable to those found here (e.g., Rebusco et al., 2006).

It is considerably easier to mix thermal plasma in the presence of anisotropic thermal conduction, since the plasma is formally always buoyantly unstable and thus already prone to mixing! By contrast, treating the plasma as adiabatic (i.e., ignoring thermal conduction) results in an artificially stabilizing entropy gradient in cluster plasmas.
^{12} As a concrete example of these effects, we have demonstrated that cosmic rays initially injected into the polar regions can be partially mixed to the equator by convective overshooting in the ACRI and HBI unstable regions (Fig. 11); this effect is largely absent in simulations that treat the plasma as adiabatic. If wave heating due to cosmic ray streaming or heating due to Coulomb interactions is important in clusters (e.g., Guo & Oh 2008), a mechanism similar to that described here may be crucial in redistributing cosmic rays throughout the cluster volume. More generally, to study the mixing produced by external sources of turbulence such as galactic wakes or cosmicray filled bubbles, we suspect that anisotropic thermal conduction must be accounted for, so that the buoyant response of the thermal plasma is correctly represented.
Having summarized our primary results, we now describe several caveats and directions for future research. First, to generate the ACRI, we have injected cosmic rays using a subsonic source term at small radii. In reality a significant fraction of the cosmic rays produced by an AGN are expected to be produced in a supersonic jet shocking against the ICM. The spatial distribution of cosmic rays produced by jets is poorly understood. The intuition drawn from our simulations should apply as long the source of cosmic rays ultimately produces a centrally concentrated bubble of relativistic plasma that expands subsonically.
Our calculations intentionally do not include plasma cooling. Instead of trying to solve the cooling flow problem, our goal has been to study the basic physics of buoyancy instabilities in the combined relativistic + thermal plasma, implicitly assuming that some heating process is preventing catastrophic cooling of the plasma. Our calculations also do not include anisotropic ion viscosity which is times smaller than electron thermal conductivity. Finally, we do not treat the effects of cosmic ray streaming with respect to the thermal plasma from first principles, although our choice of the cosmic ray diffusion coefficient qualitatively accounts for limits on cosmic ray streaming produced by selfexcited Alfvén waves (see §2.2). In future work, we intend to include all of the above effects, which will provide a more quantitative model of plasma in cluster cores.
Finally, we note that in a full cosmological context, galaxy clusters will contain a large number of galaxies and other dark matter substructure. The motion of such bound objects through the ICM will reorient the magnetic field and generate downstream turbulence. The interplay between this turbulence and that generated by the instabilities studied in this paper is worth investigating in detail in future work. This interaction may create a magnetic dynamo in the ICM that is more effective than that produced by the HBI alone: galaxies moving through the ICM will comb out the magnetic field lines in the radial direction, while the HBI will amplify the field and generate a strong perpendicular magnetic field component from the seed radial field created by galactic wakes.
Appendix A Numerical Tests
a.1. Cosmic ray shock tube
We have tested the adiabatic implementation of cosmic rays with a 1D shock tube problem discussed in Pfrommer et al. (2006) and in Rasera & Chandran (2008). Like Rasera & Chandran (2008), we also use 1024 grid points. Thermal conduction and cosmic ray diffusion are absent for this test problem. The left () and right () states are given by (, , , ) = (1, 0, , ) and (, , , ) = (0.2, 0, , ), respectively. Figure 12 shows profiles at (solid line) and at (shorter than the crossing time; points). The profiles match very well with the analytic result and with Fig. 3 of Rasera & Chandran (2008). In addition we also test advection of a passive scalar density governed by equation (22). The passive scalar density shows a discontinuity at the location of the contact discontinuity; volume integrated is not conserved but volume integrated is.
Footnotes
 affiliation: Chandra Fellow
 affiliation: Chandra Fellow
 We have dropped a term on the righthand side of expression for equal to , which is small for ; in the opposite limit with . Notice that .
 A physically more realistic model would be to include a feedback source term for cosmic rays where, instead of a fixed cosmic ray luminosity, a fixed fraction of the instantaneous mass accretion rate is converted into cosmic ray power; this level of detail is unnecessary for studying the basic physics of cosmicray driven convection but will be included in future calculations with radiative cooling.
 footnotetext: Halfangle around the polar axis over which CR source is applied.
 footnotetext: The fiducial run.
 footnotetext: Although monopolar, everywhere, including the boundaries, since the origin is excluded from the computational domain.
 footnotetext: Resolution for CR2Ddbl is 256 128, double that of CR2D.
 footnotetext: Plasma is adiabatic () for this run.
 By the antidynamo theorem, the amplification in 2D must be transient. The dynamical time in clusters is so long, however, that the “transient” can last a Hubble time!
 A purely thermal electronion plasma can also show an adiabatic, rather than diffusive, response even in the presence of rapid electron thermal conduction; this occurs if electrons and protons are not collisionally coupled on the buoyancy timescale. We find, however, that even in the outer parts of clusters, the electronproton energy exchange time is shorter than the buoyancy timescale and thus the MTI/HBI limits are appropriate.
 Even including isotropic thermal conduction reduces the stabilizing effect of the entropy gradient; it does not, however, capture the MTI/HBI, which are driven by anisotropic thermal conduction along magnetic field lines.
References
 Allen, S. W., Dunn, R. J. H., Fabian, A. C., Taylor, G. B., & Reynolds, C. S. 2006, MNRAS, 372, 21
 Balbus, S. A. 2000, ApJ, 534, 420
 Berezinskii, V. S., Bulanov, S. V., Dogiel, V. A., & Ptuskin, V. S., Astrophysics of Cosmic Rays, Amsterdam:NorthHolland, 1990, Ed. Ginzburg, V. L.
 Bertschinger, E. & Meiksin, A. 1986, ApJ, 306, L1
 Binney, J. & Tabor, G. 1995, MNRAS, 276, 663
 Bîrzan, L., Rafferty, D. A., McNamara, B. R., Wise, M. W., & Nulsen, P. E. J. 2004, ApJ, 607, 800
 Braginskii, S. I., Reviews of Plasma Physics, Vol. 1, Consultants Bureau, 1965
 Chandran, B. D. G. & Dennis, T. J. 2006, ApJ, 642, 140
 Chandran, B. D. G. & Rasera, Y. 2007, ApJ, 671, 1413
 Cho, J. & Vishniac, E. T. 2000, ApJ, 538, 217
 Churazov, E., Forman, W., Vikhlinin, A., Tremaine, S., Gerhard, O., & Jones, C. 2008, MNRAS, 388, 1062
 Ciotti, L. & Ostriker, J. P. 2001, ApJ, 2001, 551, 131
 Conroy, C. & Ostriker, J. P. 2008, ApJ, 681, 151
 Dennis, T. & Chandran, B. D. G. 2009, ApJ, 690, 566
 Drury, L. & Völk, H. 1981, ApJ, 248, 344
 Eilek, J. A. 2003, Proc. of The Riddle of Cooling Flows in Galaxies and Clusters of Galaxies, E13, Ed. T. H. Reipeich, J. C. Kempner, & N. Soker
 Engelmann, J. J. et al. 1990, A&A, 233, 96
 Fabian, A. C. 1994, ARA&A, 32, 277
 Fromang, S., Papaloizou, J., Lesur, G., & Heinemann, T. 2007, A&A, 476, 1123
 Govoni, F. & Luigina, F. 2004, Int. J. Mod. Phys. D, 13, 1549
 Guo, F. & Oh, S. P. 2008, MNRAS, 384, 251
 Hayes, J. C. et al. 2006, ApJS, 165
 Johnstone, R. M., Allen, S. W., Fabian, A. C., & Sanders, J. S. 2002, MNRAS, 336, 299
 Kelson, D. D. et al. 2002, ApJ, 576, 720
 Kulsrud, R. M. & Pearce, W. P. 1969, ApJ, 156, 445
 Kulsrud, R. M., Plasma Physics for Astrophysics 2005, Princeton Univ. Press
 Loewenstein, M., Zweibel, E. G., & Begelman, M. C. 1991, ApJ, 377, 392
 Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 490, 493
 Parrish, I. J. & Stone, J. M. 2007, ApJ, 664, 135
 Parrish, I. J. & Quataert, E. 2008, ApJ, 677, L9
 Parrish, I. J., Stone, J. M., & Lemaster, M. N. 2008, ApJ, 688, 905
 Peterson, J. R. et al. 2003, ApJ, 590, 207
 Pfrommer, C. & Enßlin, T. A. 2004, A&A, 413, 17
 Pfrommer, C. et al. 2006, MNRAS, 367, 113
 Piffaretti, R. et al. 2005, A&A, 433, 101
 Quataert, E. 2008, ApJ, 673, 758
 Rasera, Y. & Chandran, B. 2008, ApJ, 685, 105
 Rasera, Y., Srivastava, K., & Chandran, B. 2008, ApJ, 689, 825
 Rebusco, P. Churazov, E., Böhringer, H., & Forman, W. 2005, MNRAS, 359, 1041
 Rebusco, P., Churazov, E., Böhringer, H., & Forman, W. 2006, MNRAS, 372, 1840
 Reynolds, C. S., McKernan, B., Fabian, A. C., Stone, J. M., & Vernaleo, J. C. 2005, MNRAS, 357, 242
 Rosner, R. & Tucker, W. H. 1989, ApJ, 338, 761
 Schwarzschild, M., Structure and Evolution of Stars, Princeton University Press, 1958
 Sharma, P. & Hammett, G. W. 2007, J. Comp. Phys., 227, 123
 Sharma, P., Quataert, E., & Stone, J. M. 2008, MNRAS, 389, 1815
 Stone, J. M. & Norman, M. L. 1992, ApJS, 80, 753
 Stone, J. M. & Norman, M. L. 1992, ApJS, 80, 791
 Vernaleo, J. C. & Reynolds, C. S. 2006, ApJ, 645, 83
 Zakamska, N. & Narayan, R. 2003, ApJ, 582, 162