Anisotropic Thermal Conduction and the Cooling Flow Problem in Galaxy Clusters
Abstract
We examine the longstanding cooling flow problem in galaxy clusters with 3D MHD simulations of isolated clusters including radiative cooling and anisotropic thermal conduction along magnetic field lines. The central regions of the intracluster medium (ICM) can have cooling timescales of Myr or shorter—in order to prevent a cooling catastrophe the ICM must be heated by some mechanism such as AGN feedback or thermal conduction from the thermal reservoir at large radii. The cores of galaxy clusters are linearly unstable to the heatfluxdriven buoyancy instability (HBI), which significantly changes the thermodynamics of the cluster core. The HBI is a convective, buoyancydriven instability that rearranges the magnetic field to be preferentially perpendicular to the temperature gradient. For a wide range of parameters, our simulations demonstrate that in the presence of the HBI, the effective radial thermal conductivity is reduced to % of the full Spitzer conductivity. With this suppression of conductive heating, the cooling catastrophe occurs on a timescale comparable to the central cooling time of the cluster. Thermal conduction alone is thus unlikely to stabilize clusters with low central entropies and short central cooling timescales. High central entropy clusters have sufficiently long cooling times that conduction can help stave off the cooling catastrophe for cosmologically interesting timescales.
Subject headings:
convection—galaxies: clusters: general—instabilities—MHD—plasmas—Xrays: galaxies: clusters1. Introduction
Xray observations of the intracluster medium (ICM) of relaxed galaxy clusters show a centrally peaked surface brightness distribution. The observed temperatures and densities are high enough that the plasma can have a cooling time much less than 500 Myr (Sarazin, 1986). The standard isobaric cooling flow model predicts mass dropping out of the Xray emitting ICM at rates in excess of 100–. However, Xray spectroscopic observations with Chandra and XMMNewton have ruled out classical cooling flows of material below 1 keV as it would copiously emit lines such as Fe XVII that are not observed (Peterson & Fabian, 2006). Therefore, a mechanism is required to heat the ICM to avert this cooling catastrophe.
The plasma in the ICM has temperatures of 1–15 keV and number densities of – cm. The magnetic field in the ICM is estimated to be in the range of 0.1–10 G depending on where the measurement is made (Carilli & Taylor, 2002). Under these conditions, the Coulomb mean free path is many orders of magnitude larger than the gyroradius; e.g., at keV, , and G, the mean free path is kpc, while the electron gyroradius is cm. The mean free path is, however, smaller than the temperature gradient scale length. As a result, a fluid description of the ICM plasma, e.g. MHD, is appropriate, but the effects of anisotropic heat and momentum transport must also be included. The BraginskiiMHD equations (Braginskii, 1965) are a modification of the standard MHD equations to include anisotropic transport due to the magnetic field.
AGN feedback and conduction from the thermal bath at large radii are two of the most often discussed mechanisms for heating cool cluster cores. This work shall only briefly consider the former. The latter, thermal conduction, has been studied by many authors (e.g., Binney & Cowie, 1981; Narayan & Medvedev, 2001; Zakamska & Narayan, 2003; Guo et al., 2008). Because of uncertainties associated with the magnetic geometry, the heat flux is often parameterized as an effective thermal conductivity given as a fraction, , of the ideal (fieldfree) Spitzer heat flux. Previous work, however, has not considered the dynamic consequences of anisotropic conduction. The plasma in galaxy clusters is unstable to two different convective instabilities driven by anisotropic thermal conduction along magnetic field lines. The first, the magnetothermal instability, or MTI (Balbus, 2000; Parrish et al., 2008), occurs when the temperature gradient and gravity are in the same direction, as is true at large radii in galaxy clusters. Well inside the cooling core, the heatfluxdrivenbuoyancy instability, or HBI, occurs where the temperature gradient is in the opposite direction of gravity (Quataert, 2008; Parrish & Quataert, 2008). The HBI occurs in the cooling core of the ICM where the temperature increases outward. Nonlinear simulations of the MTI and HBI have shown that they significantly modify the thermal conductivity, because they saturate by rearranging the magnetic field geometry. Thus, determining the effective conductivity in the ICM requires considering the backreaction of the anisotropic heat flux on the magnetic field geometry.
In this work, we examine the stability of the cooling cores of galaxy clusters using threedimensional MHD simulations including anisotropic thermal conduction and cooling. In particular we assess the interplay between cooling and the HBI. This paper is organized as follows. In §2 we summarize the physics of the HBI and the thermal instability and their nonlinear saturation in local calculations. In §3 we then describe the equations of BraginskiiMHD and the numerical tools we utilize to solve them. §4 presents our fiducial cluster model (based on Abell 2199) in detail and examines its evolution over cosmic time. In §5, we examine a variety of variations of our fiducial model including cluster halo parameters, magnetic field strength and geometry, and the central cluster entropy. We also show a few experiments with a simple AGN heating model (§5.7). Finally, in §6 we discuss the implications of this work for the cooling flow problem and highlight some of our plans for future work.
In parallel to our work described here, Bogdanovic et al. (2009) has conducted a similar study using similar numerical methods. Their simulations start with different initial conditions and cover a different part of parameter space but reach broadly similar conclusions.
2. Physics of the HBI and Cooling
The linear physics of the HBI and its nonlinear evolution are outlined in Quataert (2008) and Parrish & Quataert (2008), respectively. We briefly review them here for clarity. The heatfluxdriven buoyancy instability is a convective instability driven by a background heat flux with the temperature gradient as the source of free energy. In contrast, the entropy gradient drives the more familiar Schwarzschild convection. For an arbitrarily oriented magnetic field, the dispersion relation of the HBI is (Quataert, 2008):
(1)  
where is the BruntVäisälä frequency,
(2) 
and
(3) 
where is the Alfvén speed, and
(4) 
is the frequency for conduction to act on a given scale, with the unit vector directed along the magnetic field and the thermal diffusivity
For a weak magnetic field, equation (1) has unstable solutions for either sign of the temperature gradient. The case of corresponds to the HBI. For a weak, initially vertical magnetic field (, ), the growth rate of the HBI is given by
(5) 
where is with respect to gravity. Qualitatively, one can picture the HBI as being driven by regions of converging and diverging perturbed magnetic field lines. Regions of converging magnetic field are conductively heated and become buoyant. In local simulations, the HBI generates MHD turbulence and a modest magnetic dynamo that amplifies the field. The most prominent method by which the HBI saturates is via a significant reorientation of the magnetic field geometry. The HBI takes a largely vertical field and reorients it to become largely horizontal. This fact is crucial for cluster cores since this reorientation of the magnetic field can significantly reduce the heat transport across an HBI unstable region.
In addition to driving the HBI, thermal conduction also has a significant impact on thermal instability as originally shown in Field & Saslaw (1965). The Field criterion states that wavelengths longer than
(6) 
are thermally unstable, where is the cooling function discussed later. For modes with wavelengths shorter than , the conduction time is shorter than the cooling time, and local perturbations are stabilized. In the ICM, the plasma is often locally stable to thermal instability, but unstable to global modes even in the presence of conduction (Kim & Narayan, 2003). Equation (6) was derived under the assumption of isotropic conduction. The results are similar for anisotropic conduction, except that conduction can only stabilize perturbations with short wavelengths along the magnetic field.
An interesting way to examine the physics of cooling in clusters comes from the recent work of Voit et al. (2008), who examined the role of the central entropy of the ICM as an indicator of feedback and star formation in galaxy clusters. The entropy is defined as . In Voit et al. (2008), the entropy profile for a cluster is fit using
(7) 
where is approximately the central entropy, is the power law normalization and is the power law exponent. They find that low entropy clusters, those with keV cm, have stronger H emission, star formation indicators, and AGN activity than higher entropy clusters, keV cm (Cavagnolo et al., 2008). As a matter of terminology, clusters with an inwardly decreasing temperature are referred to as coolcore or relaxed clusters. Clusters with an isothermal or inwardly increasing temperature profile are referred to as non coolcore clusters.
These observational results can be qualitatively understood in light of the Field criterion (eq. [6]). When cooling is pure Bremsstrahlung, the Field length becomes a function of entropy, scaling as . Thus, cooling can take place on small length scales when the entropy is low. The HBI decreases , making smaller wavelengths unstable to cooling. Fully nonlinear simulations, such as those presented here, are needed to understand the combined dynamics of cooling and the HBI.
3. Method
3.1. Equations of MHD with Anisotropic Heat Conduction
We solve the usual equations of ideal MHD with the addition of a vector heat flux, , and a gravitational acceleration :
(8) 
(9) 
(10)  
(11) 
where the symbols have their usual meaning. The total energy is given by
(12) 
and the internal energy, . We assume throughout.
The anisotropic heat flux is given by
(13) 
where the thermal diffusivity is given by the Spitzer value (Spitzer, 1962) and is a unit vector in the direction of the magnetic field. The Spitzer thermal diffusivity is given by
(14) 
Note that is a thermal diffusivity, and it depends inversely on the density. The Spitzer conductivity is which has only the wellknown dependence and no density dependence.
The energy equation also includes heating () and cooling () terms. The cooling function we adopt is from Tozzi & Norman (2001) with the functional form
(15) 
with units of erg cm s. The temperature dependence is a fit to cooling dominated by Bremsstrahlung above 1 keV and metal lines below 1 keV with
(16) 
where , , and , for a metallicity of , with units of . The majority of our runs do not include additional heating terms. For these runs in equation (10).
In runs with heating, we adopt a heating profile of the form
(17) 
where is the scale radius of heating and the normalization is chosen as
(18) 
where is the total thermal heating input. This model is motivated by a simple description of AGN feedback. We discuss simulations with heating in more detail in §5.7.
3.2. Timescales
We now discuss several of the key timescales in galaxy clusters in order to provide some intuition for the important physical processes (See Table 1). We examine these timescales for volumeaveraged quantities in our fiducial atmosphere of A2199. The generation of the fiducial atmosphere is discussed in §4.1. The volumeaveraged sound speed is approximately cm s, corresponding to a sound crossing time over 50 kpc of Myr. A weak magnetic field of 1 nG corresponds to an Alfvén crossing time over 50 kpc of years.
Timescale  Symbol  Time 

Sound Crossing  45  
Alfvén Crossing (1G)  
Conduction  20  
HBI Growth  120  
Cooling 
For a more typical magnetic field of 1 G, the Alfvén speed is cm s, leading to a crossing time across 50 kpc of 1.1 Gyr. We will see shortly that the Alfvén timescale is the longest timescale in the problem. For this cluster the central magnetic beta, the ratio of thermal pressure to magnetic pressure is for G. Even for G, the central magnetic beta is significantly greater than unity. Further out in the core, where the thermal pressure has dropped, the beta parameter is typically of order several hundred.
Also of interest is the heat conduction timescale. Our model cluster has a volumeaveraged thermal diffusivity of , which yields the scaledependent conduction time
(19) 
The HBI growth time in the fast conduction limit is
(20) 
As we will see later, the HBI has an opportunity to grow significantly compared to the average time between major mergers, roughly 5 Gyr, for a typical cluster. The final timescale of interest is the cooling time at the center of the cluster
(21) 
This cooling time estimate is in general too long (by almost a factor of 2), as it does not account for the increase in the cooling rate as temperature decreases and the density and line emission increase. Nonetheless, the cooling time is much longer than the HBI growth rate for the fastest growing, shortwavelength modes. Thus, the HBI growth is likely to play a significant role in the thermal evolution of the cluster core.
3.3. Numerical Tools
For our simulations we use the Athena MHD code (Gardiner & Stone, 2008; Stone et al., 2008) combined with the anisotropic conduction methods of Parrish & Stone (2005) and Sharma & Hammett (2007). In particular, we use harmonic averaging of the conductivity and the monotonized central difference limiter on transverse heat fluxes to ensure stability. The heating, cooling, and thermal conduction are operator split and subcycled with respect to the MHD timestep. The cooling simulations are implemented with a temperature floor of keV, below which UV lines become important, and the cooling curve fit is no longer accurate. This temperature floor prevents the cooling catastrophe from going to completion.
Most of the simulations in this work are done on uniform Cartesian grids of (128). One highresolution run is done at (256). We use modified reflecting boundary conditions for all MHD variables, in which the pressure and density are extrapolated in the ghost zones, but the other variables are reflected. This prevents the gravitational source term from introducing a spurious acceleration at the boundary. For the temperature boundary condition, we introduce a thermal bath at with a fixed temperature . For almost all of our runs kpc. This thermal bath physically represents the massive thermal energy available in the ICM outside the cluster core. Since we are simulating the entire cluster in a Cartesian geometry, there is no inner boundary condition at the cluster center.
To seed multiple modes of the HBI and break symmetry, we add Gaussian white noise perturbations to the initial velocity field such that the applied perturbation is everywhere a fixed fraction of the sound speed, typically %. In the absence of the HBI, cooling, or an applied perturbation, we find that we can hold hydrostatic equilibrium to better than one part in .
All of our runs, unless otherwise noted, are run to a time of 9.5 Gyr, a large fraction of the age of the universe. In a small number of runs, the simulations do not go to completion as a result of very severe cooling flows concentrating large quantities of mass into the central few zones of the grid. These exceptions are noted.
4. Fiducial Simulation
4.1. Initial Conditions
To introduce the phenomenology of the HBI in cluster cores, we begin by discussing a simple fiducial calculation based on observations of the galaxy cluster Abell 2199 as discussed in Zakamska & Narayan (2003) and Johnstone et al. (2002). This fiducial run is identified as T1 in the table of Runs (Table 2). Our initial conditions for the cluster core are obtained by constructing a spherically symmetric atmosphere in both hydrostatic and thermal equilibrium. We have found that it is quite advantageous to start from thermal equilibrium. Runs without thermal equilibrium experience thermal fronts propagating from the boundaries. Starting in thermal equilibrium produces a much more physical initial condition. The equations for this equilibrium are
(22) 
(23) 
where is the initial effective thermal conductivity and the heating function is neglected for our fiducial case. Our potential is chosen to be a softened NFW potential with a dark matter density distribution given by
(24) 
where is the scale mass, is the scale radius, and is the softening radius (Navarro et al., 1997). The potential softening is important for numerically maintaining a very accurate hydrostatic equilibrium. The potential is given by
(25)  
The plasma is modeled as a fully ionized ideal gas with and .
Run  ()  (kpc) 

(keV)  (keV)  (G)  (keV cm) 

R1  390  20  2  5  , radial  15.5  
T1 
390  20  2  5  , tangled  22.4  
T1HB  390  20  2  5  , tangled  22.4  
T1256  390  20  2  5  , tangled  22.4  
T2  520  26  4  9.5  , tangled  31.1  
T3  390  20  1  6  , tangled  5.46  
E1  390  20  3  5  , tangled  43.6  
E2  390  20  4  5  , tangled  83.1  
E2NC 
390  20  4  5  , tangled  83.1  
E2HB  390  20  4  5  , tangled  83.1  
E3  390  20  4.5  5  , tangled  122.  
E3NC 
390  20  4.5  5  , tangled  122.  
H1 
390  20  2  5  , tangled  14.0  
I1 
390  20  2  5  , radial  1.5  
Iso1  390  20  4.5  4.5  , tangled  52.5  
Iso2  390  20  4.5  4.5  , tangled  154. 
We solve equations (22) and (23) as a two point boundary value problem with the constraints of matching at both the inner () and outer () boundary. As these are a thirdorder system of equations, we choose a further symmetry constraint, namely, that the heat flux vanishes at the center. This system of equations is slightly stiff, but generally is soluble with a shooting method with a good initial guess. We only solve these ODEs to establish our initial condition. The subsequent evolution of the equilibrium is calculated using the full system of partial differential equations (PDEs) for MHD, equations (8)–(11).
For our fiducial model, we use a physical initial magnetic field geometry, that of tangled magnetic fields. First, in constructing the twopoint boundary value for our initial conditions, we set . Second, the initial fields are tangled with a Kolmogorov spectrum using the method outlined in detail §4.2 of Parrish et al. (2008). We initialize the fields in Fourier space as
(26) 
where is chosen as the wavenumber corresponding to 2–4 times the grid scale. We also choose a low cutoff corresponding to wavelengths of 1/2 to 1/4 of the domain size. We randomize the phase and use the Fast Fourier Transform to calculate the vector potential in real space. We choose , the appropriate space scaling for the 1D power spectrum of Kolmogorov turbulence. Our last step in initializing the magnetic field is to difference the vector potential to calculate a manifestly divergencefree initial field. The normalization of the magnetic field for our fiducial run is such that G.
Our model cluster, Abell 2199 has an inner temperature of roughly 2 keV and a temperature of 5 keV near 200 kpc (Johnstone et al., 2002). The gravitational potential is fit to an NFW profile with a scale radius of kpc, a softening radius of kpc, and a mass of . Figure 1 shows the fiducial atmosphere that results from these parameters. The central density in our thermal equilibrium model is slightly lower than the observed value.
4.2. Evolution of the Fiducial Case
The evolution of our fiducial cluster model prominently illustrates the physics of the HBI and cooling in the galaxy cluster core. This evolution is best understood through a variety of diagnostics. First, in Figure 2 we show the time evolution of the magnetic and kinetic energies in the cluster. There is a very weak magnetic dynamo in the first 2 Gyr or so due to the HBI. The initial drop in magnetic energy is due to reconnection. After Gyr, the core loses central pressure support, inflow begins, and the kinetic energy increases. Correspondingly, the magnetic energy is amplified through fluxfreezing leading to a maximum increase of . We define the magnetic energy amplification at the final simulation time, , as
(27) 
where the angle brackets denote a volume average. The kinetic energy dominates the magnetic energy at all times. Figure 2 shows that the HBI is not a strong source of magnetic field amplification in cluster cores.
The real hallmark of the HBI is the reorientation of the magnetic field geometry. Figure 3 shows the evolution of the volumeaveraged angle of the magnetic field with respect to the radial direction from its initial tangled state () to a final saturated state of 74.6. The angles given here and in Table 3 are volume averages over the entire cluster. The reorientation of the magnetic field is quite dramatic and takes place on just a few Gyr. Concomitantly, Figure 4 shows the evolution of the azimuthallyaveraged radialtemperature profile. We typically bin the temperature into 5 kpc radial bins. As the magnetic geometry evolves to be more azimuthal, the thermal conduction from outside the core begins to shut off and the temperature starts to fall. At Gyr, the central temperature has hit the cooling floor of keV—an effective proxy for the cooling catastrophe. We do not remove gas from the grid after hitting the temperature floor, and thus, we never see a true ”cooling flow.” We define the time of the cooling catastrophe, , to be when the average temperature of the inner bin has reached the temperature floor.
We can explore the magnetic field amplification in slightly more detail.
We bin the magnetic field into radial bins, and then for each bin, , we calculate the local amplification of the magnetic energy as
(28) 
Figure 5 shows the amplification as a function of radius at several times. The HBI amplifies the energy most efficiently in the middle of the core, just beyond 100 kpc. The high fields at late times in the central region are primarily due to fluxfreezing as the density increases in the core.
Finally, Figure 6 shows a post hoc calculation of the heat fluxes. See Parrish et al. (2008) §5.3 for a full discussion of the heat flux diagnostics. We calculate the heat fluxes as a shell average within the radial range of kpc. We begin by defining a fiducial heat flux to be the radial flux through the shell if the conduction were isotropic at the Spitzer value, namely:
(29) 
This value is the same as the heat flux with anisotropic conduction and purely radial magnetic field lines. We calculate the conductive heat flux and normalize to the fiducial value to calculate the Spitzer fraction, or effective conductivity, defined as
(30) 
where is given by equation (13). We also define a flux due to mass advection
(31) 
where angle brackets denote shell averages. The final component of the heat flux is the convective heat flux given by
(32)  
where refers to the local deviation from the mean of a quantity, e.g. . The first terms of equations (31) and (32) are the dominant terms.
The initial Spitzer fraction for tangled magnetic fields is due to the average over the random field geometry. From run to run there is some variation in this initial value as there is mode power on scales larger than our averaging volume. As the HBI grows, the heat flux is reduced significantly, eventually saturating at (Fig. 6). This dramatic reduction in heat flux leads to the cooling catastrophe. As the core cools and loses pressure support, the advective heat flux increases as mass is transported inwards. It is interesting to note that the convective heat flux is very small, especially during the HBI phase of the evolution.
For comparison purposes, we also run our fiducial model with purely isotropic conductivity and radial magnetic fields (run I1 in Table 1). The HBI is not present for purely isotropic conduction. Figure 7 shows the evolution of the azimuthallyaveraged temperature profile—the cluster reaches an almost isothermal temperature profile with no hint of the cooling catastrophe. At fixed pressure, the thermal instability can lead to either runaway heating or runaway cooling. This is easy to see by examining the form of the cooling term and conduction. In the Bremsstrahlung regime, cooling scales as at fixed pressure, while the conduction term scales as . For the case illustrated here, as the temperature is perturbed upwards, conductive heating increases much more rapidly than cooling. Thus, the thermal runaway drives the cluster towards an almost isothermal profile. In this simulation, there is a small amount of noise which randomizes the field a very small amount, hence in Table 3.
To summarize the fiducial case, we begin with a cluster core in hydrodynamic and thermal equilibrium. If thermal conduction were isotropic, the central temperature would latch onto the heating branch of the thermal instability, continuing to rise. Instead, the HBI begins to act on a 100 Myr timescale, reorienting the magnetic field geometry, and reducing the effective thermal conductivity from the outer part of the cluster. As the cluster center becomes denser and cooler, a thermal runaway proceeds, leading to a cooling catastrophe on a timescale comparable to the initial cooling time in the core of the cluster.
5. Variation of Parameters
In order to fully understand the behavior of the HBI in galaxy clusters, we now turn to an exploration of parameter space. Table 2 lists the various runs in which we vary the cluster properties, magnetic field strength and geometry, and the central entropy of the cluster. The following sections describe each of these experiments.
Run  max  min 
(Gyr) 
(Gyr) 


R1  15.0  65  0.13  0.82  4.0 
T1  2.5  74  0.07  1.4  2.7 
T1HB  1.05  74  0.06  1.4  2.7 
T1256 
5.9  77  0.03  1.4  2.4 
T2 
1.5  0.06  1.1  3.2  
T3  2.5  0.08  0.20  1.5  
E1  2.3  75  0.07  3.0  3.5 
E2  2.3  76  0.05  5.9  5.7 
E2NC  2.2  53  0.48  5.9  4.0 
E2HB  0.47  74  0.08  5.9  7.0 
E3  2.2  77  0.05  9.3  
E3NC  4.1  50.6  0.55  9.3  5.9 
H1 
0.30  72  0.09  1.3  
I1  1.3  24  1.0  0.82  none 
Iso1  6.1  67  0.17  2.6  1.7 
Iso2  1.74  76  0.06  12.8 
Table 3 lists the saturation properties of these runs. The magnetic field amplification, , is given as a volumeaverage over the cluster. The maximum of the magnetic field angle, , is the maximum in time of the volumeaveraged magnetic angle. Likewise, , is the minimum in time of the shellaveraged heat flux.
5.1. Radial Magnetic Fields
To assess the importance of the initial magnetic field geometry, we consider a split monopole radial magnetic field such that . This geometry is useful for illustrating the effects of the HBI on the equilibrium state. We choose a mean magnetic field of 1 nG in order to minimize the effects of magnetic tension on the equilibrium state. The resulting initial condition has a slightly higher central pressure and density than our fiducial case. In general, the atmosphere is in the rapid conduction limit on all but the largest scales. The choice of radial magnetic fields gives conduction the best chance of thermally stabilizing the cluster.
The evolution of the HBI in the radial field simulation is quite similar to that of our fiducial case. Initially the atmosphere is in thermal equilibrium with radial fields. Figure 8 shows the evolution of the volumeaveraged magnetic field from the initially radial () geometry. The HBI mode grows rapidly on a timescale of Myr and reorients the magnetic field to have in just over 4 Gyr. Figure 9 shows the resultant evolution of the temperature profile. The atmosphere initially latches onto a heating branch of the thermal instability, peaking at a central temperature of just over 3 keV at 1.1 Gyr. In the absence of the HBI, it would continue to evolve to an almost isothermal state as in run I1. Instead, however, the cluster undergoes a cooling catastrophe after 4 Gyr.
The driver of this cooling catastrophe is easily seen from Figure 10, which plots the heat fluxes as a function of time at 100 kpc. As the HBI rearranges the magnetic geometry, the Spitzer fraction plummets to . Having reduced the contact with the thermal bath at large radii, cooling becomes dominant in the core and the central temperature starts to rapidly decrease. This cooling drives mass inflow to small radii giving rise to the large inward advective flux at late times. At all times, the convective heat flux is small compared to both the saturated conductive heat flux and the mixing length estimate. Figure 11 shows the overall evolution of this cluster in a 2D slice taken at . As the magnetic field lines wrap in the azimuthal direction, the central temperature decrease is easily observed.
The cooling catastrophe’s vigor is driven by two complementary properties of the cooling curve. First, at fixed pressure, the cooling increases as the temperature declines at the center of the cluster. Second, below 1 keV, line emission from metals like iron and oxygen becomes increasingly important, scaling as . Thus, once gas has cooled below 1 keV, the cooling is much harder to reverse. Observations show very few clusters with central temperatures below 1 keV.
5.2. Strong Magnetic Fields
In the previous sections, we have demonstrated the evolution of the HBI in clusters for weak magnetic fields. Under these conditions, magnetic tension is not significant. We now consider a more realistic magnetic field of G. See Carilli & Taylor (2002) for a review of cluster magnetic field measurements. When the tension force becomes comparable to the buoyancy time, the HBI growth is suppressed at small scales. For example, for G in our fiducial cluster, magnetic tension becomes important on scales smaller than 5–20 kpc depending on the location in the cluster.
The cluster evolution of this case (labeled run T1HB) is similar to the weaker magnetic field tangled case. The maximum field strength we can simulate in this constant field model is limited since a modest field at the center can be dynamically important and low several scale heights out from the center. For run T1HB we find that there is a modest amount of numerical reconnection that dissipates some of the initial magnetic energy. A very weak dynamo leads to a maximum magnetic energy increase of only 5%. The maximum magnetic field angle and minimum heat flux are quite similar to the lower field case. In addition, the cooling catastrophe occurs at almost exactly the same time. Thus, a constant 1 G magnetic field provides very little stabilization for our fiducial cluster.
It is interesting to examine the magnetic field geometry in this simulation from a different perspective, namely, how the average magnetic field angle varies versus radius. This is shown in Figure 12.
Initially, the angle is distributed as a random variable about 60. At 1.3 Gyr, the HBI has had somewhat more than 10 growth times to significantly increase the average angle and decrease the conductivity. By 2.2 Gyr, radial infall from inhomogeneously cooling material has straightened the magnetic field out within 30 kpc, enhancing the radial heat flux. This idea was suggested by Balbus & Reynolds (2008) as a mechanism for opposing the HBI and slowing the cooling catastrophe. Unfortunately, by the time the magnetic field has been straightened out, the plasma within 20 kpc has a mean temperature of 0.2 keV and a mean density of . This corresponds to an increase in cooling from the initial equilibrium of . Meanwhile, the average angle in this region is about 50, implying that the effective conductivity has only increased a factor of 1.7 over the initial value—clearly not enough to stabilize the equilibrium. Similar behavior is found for all of our runs that reach a cooling catastrophe.
5.3. Dependence on Cluster Parameters
We now consider the effect of varying the cluster parameters, as exemplified by runs T2 and T3 in Table 2. We will only compare the tangled magnetic field geometries since that is the most physically relevant setup. Run T2 is modeled on Abell 2390, a hot, massive cluster. This NFW halo has a larger mass and scale radius. The core temperature rises from a central temperature of 4 keV to an outer temperature of 9.5 keV at 300 kpc. We choose the softening radius to be approximately of the scale radius as we did for our fiducial case. The dependence of our results on cluster mass and temperature is small. The HBI growth time and central cooling time for our model of A2390 are 130 Myr and 1.1 Gyr, respectively, similar to our fiducial run. This run evolves in a similar way to run T1, reaching the cooling catastrophe in 3.2 Gyr with similar saturated parameters. Due to the vigor of the cooling catastrophe in this more massive case, the run does not reach 9.5 Gyr, as most of the mass has collapsed to the central few zones.
In run T3, we examine the effect of changing the initial thermal profile but not the NFW parameters. This run is the same as the fiducial case but the temperature now initially varies from 1 to 6 keV. The physics of how the magnetic geometry is modified remains very similar; however, the cooling catastrophe occurs at an earlier time of 1.5 Gyr. The primary reason for this is that the higher initial central density and lower temperature make this cluster cool faster to the temperature floor. Thus, we see only a very weak dependence on the cluster parameters, mostly being driven by the initial location on the cooling curve.
5.4. Dependence on Central Entropy
Motivated by the discussion of the role of central entropy in §2, we have undertaken a parameter study in central entropy. The runs T1, E1, E2, and E3 form a monotonic progression of central entropy from 22.4 to 122 keV cm. We have effected the entropy variation by modifying the initial central temperature while maintaining thermal equilibrium. We fit our cluster entropy profiles with a powerlaw of the form of equation (7). The parameter is typically within 10–20% of the central entropy. The primary difference evident in examining Table 3 is that the time of the cooling catastrophe increases as increases–reaching a maximum of just over 10 Gyr for the highest entropy case with nG magnetic fields (run E3). This phenomenon is easy to understand since the central cooling time itself increases as increases.
It is interesting to compare these runs to several highentropy simulations without conduction, runs E2NC and E3NC. These runs have cooling but no thermal conduction. First, it is clear that the central cooling time estimated by equation (21) is an overestimate compared to the actual time of the cooling catastrophe. For example, run E3 has a predicted cooling timescale of 9.3 Gyr but in fact reaches a cooling catastrophe without conduction in 5.9 Gyr. Second, these runs without conduction show that even though the effective conductivity is reduced by the HBI, the time to a cooling catastrophe is significantly longer than in the absence of conduction. In the case of run E3, thermal conduction delays the cooling catastrophe from 5.3 Gyr to Gyr.
We now consider runs T1HB and E2HB which have higher magnetic fields of 1 G. In the former, the low entropy run, magnetic tension does very little to stabilize the HBI resulting in a cooling catastrophe at almost the same time as the 1 nG case. In the latter, the high entropy run, magnetic tension plays a more significant role and increases the time of the cooling catastrophe from 5.7 Gyr (the low field case) to 7.0 Gyr (the higher field case).
Figure 13 shows the evolution of the temperature profile for run E2HB. The cooling is especially slow when the temperature is high. In fact, the time to reach the cooling catastrophe, Gyr, is longer than the typical time between major mergers for galaxy clusters –5 Gyr (Cohn & White, 2005).
How does a cluster reach the high entropy states for which the cooling catastrophe can be avoided? Cluster mergers may provide shocks that heat the cluster and boost its entropy. In addition, strong AGN feedback may increase the cluster entropy enough to slow the cooling instability, even in the presence of the HBI.
5.5. Isothermal Initial Conditions
As a further exploration of the important physics in clusters, we present the results of two simulations, runs Iso1 and Iso2, that are initialized with isothermal temperature profiles at 4.5 keV. Our usual approach of imposing thermal equilibrium (eq. [23]) does not work in this case, and the central density becomes a free parameter for constructing the hydrostatic equilibrium. By definition, the initial state is not in thermal equilibrium. The other initial parameters are the same as our fiducial case.
Run Iso1 has an initial central density of , which corresponds to an entropy of 52.5 keV cm and a cooling time of 2.6 Gyr. By virtue of this short cooling time and the action of the HBI, the core experiences the cooling catastrophe at 2.1 Gyr. By lowering the central density by a factor of 5 to , run Iso2 has a central entropy of keV cm and a cooling time of 12.8 Gyr. After 9.5 Gyr, this run has developed a slightly relaxed core ( keV), but is very far from the cooling catastrophe.
Starting from the nonequilibrium initial condition, we see very similar qualitative behavior to our equilibrium model. Runs with short cooling times develop both the cool core and the HBI quickly, and runs with long cooling times are much more stable. In our isothermal clusters the state is merely transient, as even clusters with long central cooling times and high entropy would eventually evolve into relaxed, coolcore clusters.
5.6. Resolution Dependence
We perform one high resolution simulation of our fiducial case, T1256, which is a tangled field simulation. We terminate this run at 6.3 Gyr (2/3 of the normal run time) in light of the large processor time required. All of the qualitative results discussed thus far hold up at higher resolution. We do see some minor differences resulting from the increased resolution: smaller scales are now available on which the HBI can act. The HBI both amplifies the magnetic energy slightly more and is able to reach an even larger average magnetic field angle. The final volumeaveraged magnetic field angle of corresponds to an incredibly azimuthally wrapped magnetic field with a tiny effective conductivity of . This precipitates the cooling catastrophe on a slightly shorter timescale.
5.7. Experiments with Central Heating
Surveys of galaxy cluster cores consistently find Xray cavities or bubbles filled with radio emitting plasma or cosmic rays (e.g., Bîrzan et al., 2004; Dunn & Fabian, 2006). These structures indicative of feedback are especially prevalent in clusters with short central cooling times, roughly the same population of low entropy clusters mentioned previously. In an effort to understand the interplay between the HBI and heating, we proceed with a preliminary analysis of heating in cluster cores.
A number of groups have proposed cosmic rays from a central AGN as a heating mechanism. In particular, streaming cosmic rays can excite Alfvén waves which nonlinearly Landau damp to heat the plasma (Loewenstein et al., 1991). Guo & Oh (2008) have demonstrated in 1D models that a combination of parameterized cosmic ray feedback and conduction can prevent significant cooling for a Hubble time. For our test problems, our heating luminosity is parameterized as in equation (17) with the initial normalization set by equation (18), motivated by Chandran (2005). These heating functions are generic and do not discriminate among cosmic ray or mechanical energy injection.
The feedback dynamics of the cluster core is qualitatively simple. As the core cools, the accretion rate onto the supermassive () black hole at the center of the cD galaxy increases. As increases, the feedback heating increases, slowing the cooling. If heating becomes too effective, the accretion ceases, and a feedback loop is established. Thus, a simple static heating model is insufficient, and we instead implement a rough timevariable version. Namely, we sum the cooling luminosity within the heating effective radius, , at ,
(33) 
We then calculate the cooling luminosity in a similar way at every timestep and scale the initial heating luminosity to the current cooling luminosity as
(34) 
This methodology is not ideal, but given that we cannot resolve the Bondi radius on our grid, it is preferable to extrapolating the central density and temperature down to the Bondi radius; and it ensures we have an approximate feedback mechanism. It should be noted that this heating model can be numerically unstable when the cooling instability has progressed, and large feedback heating is added within a small region. We will improve this treatment in future work.
As an example calculation, we take our fiducial atmosphere and add a total initial heating luminosity of erg s with a characteristic radius of kpc. We construct an initial thermal equilibrium such that heating, cooling, and thermal conduction all balance. This run, labeled H1, has a very interesting thermal evolution as is shown in Figure 14.
The HBI proceeds very slowly for this simulation since the initial average magnetic field is 3.5G, strong enough to exert significant tension. What is especially interesting about this run is that the centrally concentrated feedback heating drives a minimum in the temperature profile at kpc after 5.4 Gyr of evolution. This type of profile, with a minimum slightly offset from the center, is actually observed in some clusters (Sanderson et al., 2006). The cooling flow region seems to be simply pushed outwards from the center of the cluster, perhaps evidence that an additional volumetric heating component is needed. The heating power at 6 Gyr has increased only modestly to a new value of erg s. Due to the combination of heating and the HBI the cooling catastrophe occurs much later than the initial central cooling time of 1.3 Gyr. Unfortunately, we are not able to follow this run to completion as the sharp temperature discontinuity combined with the rapid cooling that follows leads to numerical instabilities. While not necessarily thermally stable, this run demonstrates that cluster cores with heating and the HBI can remain stable for longer than the time between major mergers, although not necessarily a Hubble time. In general, it appears that the combination of magnetic fields of G (to slow the HBI) and a modest amount of AGN feedback significantly slow the cooling catastrophe, even for lowentropy clusters.
6. Discussion and Conclusions
The plasma in the intracluster medium of galaxy clusters is dilute and magnetized with a mean free path large compared to the gyroradius. Under these conditions, heat transport is anisotropic along magnetic fields. This results in the ICM being unstable to the heatfluxdriven buoyancy instability in regions where the temperature increases outward (Quataert, 2008). The cores of galaxy clusters also often have short cooling times of Myr. In order to understand the thermal evolution of galaxy clusters with cooling and the HBI, we have performed threedimensional timedependent MHD simulations of galaxy clusters cores.
Isolated galaxy clusters evolved with magnetic fields, anisotropic conduction, and cooling share a number of common properties. We begin with a cluster that is in both hydrostatic and thermal equilibrium. After Myr, the HBI begins to rearrange the magnetic field geometry in the cluster core; the magnetic field saturates with an average angle between the magnetic field and the radial direction of . Second, as the magnetic geometry is rearranged to be tangential to the temperature gradient, the magnetic field exerts a thermally insulating effect, reducing the effective radial thermal conductivity to % of the Spitzer value. Finally, having reduced the thermal conduction from the outer parts of the cluster, and lacking another heat source, the core proceeds to a cooling catastrophe on a timescale comparable to the initial central cooling time.
We have studied a number of different parameter variations and find several interesting trends. Motivated by the observational work of Voit et al. (2008), we have explored the effects of different initial entropies. For larger central cluster entropies, the time of the cooling catastrophe is delayed to more than 9.5 Gyr for our highest entropy cluster (122 keV cm). In addition, we find that stronger magnetic fields, G, can suppress the HBI via magnetic tension forces. The onset of the cooling catastrophe can thus be delayed in these higher magnetic field calculations.
We have also carried out initial calculations of the effects of heating on the ICM using a parameterized heating function in which the total heating power is proportional to the total rate of cooling in the central 20 kpc. Despite some difficulty with numerical instabilities inherent in the method, we find that modest heating rates of erg s can substantially delay the cooling catastrophe to Gyr, longer than the time between major cluster mergers. There is, however, some evidence that centrally concentrated heating may simply move the cooling catastrophe further out in the core (Fig. 14).
Observations and simplified theoretical models both suggest that there are two distinct quasistable cluster states (Voit et al., 2008; Guo & Oh, 2008; Cavagnolo et al., 2009). High entropy, fairly isothermal cluster cores have long growth times for both the Field instability and the HBI. The longer central cooling times require less conductive heating to balance cooling. These thermal states are longlived even in the absence of mergers. By contrast, low entropy relaxed clusters (coolcore clusters) with short central cooling times are unstable to both the Field instability and the HBI. As we have shown, this population of clusters cannot be stabilized by conduction alone and must have an additional feedback mechanism, plausibly the central AGN but potentially other sources. Observations of cool core clusters with central entropies keV cm show a number of feedback indicators, including H emission indicative of cool gas at K, radio emission indicative of AGN feedback, and optical color gradients indicative of central star formation in the BCG (Cavagnolo et al., 2008; Voit et al., 2008). High entropy clusters, in which conduction is more important, generally show none of these feedback indicators. Our simulations of these high entropy clusters show that they are thermally stable for cosmologically long timescales, and that conduction provides a significant stabilizing effect, e.g. runs E2 and E3 (see, Tables 2 and 3).
It may be possible for clusters to transition between these two populations. Relaxed clusters may be promoted to high entropy clusters by significant heating, such as a major merger or an especially energetic feedback event. Recent work shows that disruption of cool cores in a merger is possible at cosmologically early times but difficult at late times (Burns et al., 2008). Alternatively, isothermal moderate entropy clusters can eventually become relaxed coolcore clusters over long timescales.
A key task for future work is to better understand the proposed heating mechanisms for low entropy clusters. In particular, bubbles and jets from an AGN are far from geometrically isotropic. Not only must the heating be locally efficient, but the heating must then be distributed by some mechanism azimuthally around the cluster core to prevent a cooling catastrophe. The enhanced azimuthal heat transport from the HBI may play a significant role in redistributing local AGN heating throughout the cluster core.
In future work, we will examine these heating mechanisms in more detail including the relevant physics. For the case of buoyant bubbles, there are many unanswered questions about the disruption time of the bubbles. This shredding is governed by RayleighTaylor and KelvinHelmholtz instabilities. In the full BraginskiiMHD treatment, momentum is transported by ions anisotropically along magnetic field lines. If the bubbles are indeed draped by magnetic fields, then the RT and KH instabilities will be modified by an anisotropic Braginskii viscosity. Cosmic rays may play a role in directly heating the plasma by exciting Alfvén waves (Guo & Oh, 2008). Finally, galaxy wakes in a full cosmological context can also provide heating to the ICM (Conroy & Ostriker, 2008). They may also increase the importance of thermal conduction by competing with the HBI to reorient the magnetic field.
A key lesson of this work is that it is difficult to characterize the ICM plasma as having a single thermal conductivity parameterized by a constant . Buoyancy instabilities such as the HBI and MTI directly modify the magnetic geometry and selfconsistently evolve the system to a new state that may enhance or suppress the effective conductivity. In the cores of galaxy clusters, the HBI suppresses thermal conduction from the large heat reservoir at large radii. In the absence of AGN feedback or very large magnetic fields in cluster cores, it appears that conduction alone cannot solve the cooling flow problem.
Footnotes
 affiliation: Astronomy Department & Theoretical Astrophysics Center, 601 Campbell Hall, University of California, Berkeley, CA 94720; iparrish@astro.berkeley.edu
 affiliation: Chandra/Einstein Fellow
 affiliation: Astronomy Department & Theoretical Astrophysics Center, 601 Campbell Hall, University of California, Berkeley, CA 94720; iparrish@astro.berkeley.edu
 affiliation: Astronomy Department & Theoretical Astrophysics Center, 601 Campbell Hall, University of California, Berkeley, CA 94720; iparrish@astro.berkeley.edu
 affiliation: Chandra/Einstein Fellow
 The literature is not consistent regarding the use of and . We will use to represent a true diffusion coefficient and to represent a conductivity in units erg cm s K.
 Evaluated for kpc as volumeaveraged quantities.
 Softening radius of NFW halo (eq. [24])
 Fiducial Run
 No conduction
 No conduction
 Simulation includes additional heating (see §5.7).
 Isotropic conduction only
 Initial cooling time at the innermost radii ( kpc).
 Time of Cooling Catastrophe, when the inner gridpoint has reached the cooling floor (see §4.2).
 Run time is less than 9.5 Gyr.
 Run time is less than 9.5 Gyr.
 Run time is less than 9.5 Gyr.
References
 Balbus, S. A. 2000, ApJ, 534, 420
 Balbus, S. A., & Reynolds, C. S. 2008, ApJ, 681, L65
 Binney, J., & Cowie, L. L. 1981, ApJ, 247, 464
 Bîrzan, L., Rafferty, D. A., McNamara, B. R., Wise, M. W., & Nulsen, P. E. J. 2004, ApJ, 607, 800
 Bogdanovic, T., Reynolds, C. S., Balbus, S. A., & Parrish, I. J. 2009, ApJ Submitted
 Braginskii, S. I. 1965, Reviews of Plasma Physics, ed. M. A. Leontovich, Vol. 1 (New York: Consultants Bureau), 205
 Burns, J. O., Hallman, E. J., Gantner, B., Motl, P. M., & Norman, M. L. 2008, ApJ, 675, 1125
 Carilli, C. L., & Taylor, G. B. 2002, ARA&A, 40, 319
 Cavagnolo, K. W., Donahue, M., Voit, G. M., & Sun, M. 2008, ApJ, 683, L107
 —. 2009, ArXiv:0902.1802
 Chandran, B. D. G. 2005, ApJ, 632, 809
 Cohn, J. D., & White, M. 2005, Astroparticle Physics, 24, 316
 Conroy, C., & Ostriker, J. P. 2008, ApJ, 681, 151
 Dunn, R. J. H., & Fabian, A. C. 2006, MNRAS, 373, 959
 Field, G. B., & Saslaw, W. C. 1965, ApJ, 142, 568
 Gardiner, T. A., & Stone, J. M. 2008, Journal of Computational Physics, 227, 4123
 Guo, F., & Oh, S. P. 2008, MNRAS, 384, 251
 Guo, F., Oh, S. P., & Ruszkowski, M. 2008, ApJ, 688, 859
 Johnstone, R. M., Allen, S. W., Fabian, A. C., & Sanders, J. S. 2002, MNRAS, 336, 299
 Kim, W.T., & Narayan, R. 2003, ApJ, 596, 889
 Loewenstein, M., Zweibel, E. G., & Begelman, M. C. 1991, ApJ, 377, 392
 Narayan, R., & Medvedev, M. V. 2001, ApJ, 562, L129
 Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 490, 493
 Parrish, I. J., & Quataert, E. 2008, ApJ, 677, L9
 Parrish, I. J., & Stone, J. M. 2005, ApJ, 633, 334
 Parrish, I. J., Stone, J. M., & Lemaster, N. 2008, ApJ, 688, 905
 Peterson, J. R., & Fabian, A. C. 2006, Phys. Rep., 427, 1
 Quataert, E. 2008, ApJ, 673, 758
 Sanderson, A. J. R., Ponman, T. J., & O’Sullivan, E. 2006, MNRAS, 372, 1496
 Sarazin, C. L. 1986, Reviews of Modern Physics, 58, 1
 Sharma, P., & Hammett, G. W. 2007, Journal of Computational Physics, 227, 123
 Spitzer, L. 1962, Physics of Fully Ionized Gases (Physics of Fully Ionized Gases, New York: Interscience (2nd edition), 1962)
 Stone, J. M., Gardiner, T. A., Teuben, P., Hawley, J. F., & Simon, J. B. 2008, ApJS, 178, 137
 Tozzi, P., & Norman, C. 2001, ApJ, 546, 63
 Voit, G. M., Cavagnolo, K. W., Donahue, M., Rafferty, D. A., McNamara, B. R., & Nulsen, P. E. J. 2008, ApJ, 681, L5
 Zakamska, N. L., & Narayan, R. 2003, ApJ, 582, 162