Galaxy Motions, Turbulence and Conduction in Clusters of Galaxies
Abstract
Unopposed radiative cooling in clusters of galaxies results in excessive mass deposition rates onto the central brightest cluster galaxy. However, the cool cores of galaxy clusters are continuously heated by thermal conduction and turbulent heat diffusion due to minor mergers or the galaxies orbiting the cluster center. These process can either reduce the energy requirements for AGN heating of cool cores, or they can prevent overcooling altogether. We perform 3D MHD simulations including fieldaligned thermal conduction and selfgravitating particles to model this in detail. Turbulence is not confined to the wakes of galaxies but is instead volumefilling, due to the excitation of largescale gmodes. We systematically probe the parameter space of galaxy masses and numbers to assess when the cooling catastrophe is prevented. For a wide range of observationally motivated galaxy parameters, we find that the magnetic field is randomized by stirring motions, restoring the conductive heat flow to the core. The cooling catastrophe either does not occur or it is sufficiently delayed to allow the cluster to experience a major merger that could reset the conditions in the intracluster medium. Whilst dissipation of turbulent motions (and hence dynamical friction heating) is negligible as a heat source, turbulent heat diffusion is extremely important; it predominates in the cluster center. However, thermal conduction becomes important at larger radii, and simulations without thermal conduction suffer a cooling catastrophe. Conduction is important both as a heat source and to reduce stabilizing buoyancy forces, enabling more efficient diffusion. Turbulence enables conduction, and conduction enables turbulence. In these simulations, the gas vorticity—which is a good indicator of trapped gmodes–increases with time. The vorticity growth is approximately mirrored by the growth of the magnetic field, which is amplified by turbulence.
keywords:
conduction – cooling flows – galaxies: clusters: general – galaxies: active – instabilities – Xrays: galaxies: clusters1 Introduction
The intracluster medium (ICM) in many galaxy clusters has central cooling times shorter then the Hubble time.
Radiative cooling should lead to large accumulation of cold material
in their centers; however, there is no observational evidence for such gas.
This can be understood if some source of heating balances cooling in the ICM.
The heating mechanisms invoked to explain this overcooling problem
involve AGN “radio mode” heating (e.g., Binney &
Tabor (1995); Churazov et al. (2002); Fabian
et al. (2003); Ruszkowski et al. (2004a, b); Scannapieco
& Brüggen (2008)),
preheating by AGN (McCarthy et al., 2008), cosmic rays from AGN (Guo & Oh, 2008; Sharma et al., 2009a),
supernovae, turbulent mixing (Kim &
Narayan, 2003a; Voigt &
Fabian, 2004; Dennis &
Chandran, 2005),
thermal conduction (Zakamska &
Narayan, 2003; Kim &
Narayan, 2003b), a combination of thermal conduction and AGN (Ruszkowski &
Begelman, 2002)
and dynamical friction (ElZant
et al., 2004; Kim
et al., 2005; Kim, 2007); see Conroy &
Ostriker (2008) and
McNamara &
Nulsen (2007) and references therein for reviews of the above mechanisms).
Conduction alone is unlikely to offer the complete solution to the overcooling problem for the full range of
cluster masses, as its strong temperature dependence implies that it is less effective in lower mass clusters. Furthermore, thermal conduction is well known to be an unstable heating mechanism, either failing to avert a cooling catastrophe, or leading to an isothermal temperature profile (Bregman &
David, 1988; Guo & Oh, 2008; Conroy &
Ostriker, 2008).
Nevertheless, thermal conduction may entirely suppress cooling in non coolcore (NCC) clusters and reduce
the constraints on the required energy injection by AGN in coolcore (CC) clusters (Guo
et al., 2008); indeed, it is different to stabilize massive clusters with AGN feedback alone (Conroy &
Ostriker, 2008), and a second heat source (such as conduction) is generally required. Besides offsetting radiative losses and stemming a cooling catastrophe, conduction can have important implications for establishing the observed bimodality in cluster core entropy (Guo
et al., 2008), and the star formation threshold in brightest cluster galaxies (Voit et al., 2008). Indeed, a sudden increase in conduction (due to say, turbulence from an AGN outburst or a merger) could mediate a cool core to non coolcore transition (Guo & Oh, 2009).
Thermal conduction can be strongly suppressed by magnetic fields that are known to be present in the ICM (e.g., Enßlin et al. (2003), Vogt et al. 2003). However, interest in thermal conduction as a potential heating mechanism was revived by Narayan & Medvedev (2001) who suggested that even in the presence of tangled Bfields, the level of conduction can be an appreciable fraction of the SpitzerBraginskii value. Computing the exact magnitude and distribution of the effective conductivity of the ICM is further complicated by buoyancy instabilities which reorient the magnetic field. When temperature increases in the direction of gravity, as in the cluster outskirts, the magnetothermal instability (MTI; Balbus (2000); Parrish & Stone (2005)) tends to make the Bfields radial, thereby increasing the effective conduction. On the other hand, in cool cores where temperature decreases in the direction of gravity, the heat flux buoyancy instability (HBI; Quataert (2008); Parrish et al. (2009); Bogdanović et al. (2009)) tends to reorient the fields in the direction perpendicular to that of gravity, effectively shutting down thermal conduction.
However, these instabilities do not operate in a static atmosphere. Chandra and XMM observations show that the cluster gas is rarely in perfect hydrostatic equilibrium. Sloshing motions due to minor mergers, AGN, or galaxy motions can continuously and significantly perturb the gas, as has been repeatedly seen in many disparate numerical simulations (Evrard, 1990; Norman & Bryan, 1999; Nagai et al., 2003; Vazza et al., 2009, 2010; Ascasibar & Markevitch, 2006; ZuHone et al., 2010). Current observational evidence for turbulence ranges from the analysis of pressure maps (Schuecker et al., 2004), its effect on resonantline scattering (Churazov et al., 2004), and Faraday rotation maps (Vogt & Enßlin, 2005; Enßlin & Vogt, 2006), as well as constraints on turbulent line widths (Sanders et al., 2009). Low levels of turbulence in the ICM can randomize the field configuration set up by the HBI and restore the heat flow to the core (Ruszkowski & Oh, 2010; Parrish et al., 2010). Both of these works modeled the ICM turbulence via a simple driving mechanism to determine the level of turbulence required to effectively restore thermal conduction. This approach did not allow us to link the level of turbulence to the physical properties of the cluster (such as mechanical luminosity of the central AGN or the properties of cluster galaxies). Also, the driving mechanism led, by construction, to volumefilling turbulence which was very effective in randomizing the magnetic field. While the low amplitude of the required subsonic turbulence is eminently feasible (), the realism of volumefilling turbulence is less clear. For instance, both analytic calculations (Subramanian et al., 2006) and numerical simulations (Iapichino & Niemeyer, 2008) predict that turbulence due to galaxy wakes should not be volumefilling (), as turbulence is largely confined to ’streaks’ behind orbiting galaxies.
In this paper, we extend our previous work and perform threedimensional MHD simulations of the effect of turbulence driven by galaxy motions on the properties of the anisotropic thermal conduction. We show how the trapping of gravity modes excited by the orbiting galaxies can lead to volumefilling turbulence of the right magnitude to restore conductive heat flow. We demonstrate how these subsonic motions generate vorticity and lead to the growth of magnetic field via kinematic dynamo action. We also show that turbulent heat diffusion is an important part of the energy budget.
The paper is organized as follows. In §2 we review basic theoretical expectations for the interaction between turbulence and magnetic fields. In §3 we describe the numerical methods and the setup of the inital conditions. In §4,
we describe our results, including the level and volumefilling nature of turbulence, evolution of the gas temperature, generation
of vorticity and magnetic fields, and nature of heating mechanisms. Conclusions are summarized in §5.
2 Theoretical Expectations: Turbulence and Trapping of gmodes
It is useful to begin by reviewing some basic theoretical expectations for the behavior of turbulence excited in galactic wakes in clusters. In principle, orbiting galaxies can excite galactic wakes by two means: hydrodynamically (as the ICM collides with the ISM of the galaxy), and gravitationally (similar to dynamical friction for collisionless particles). In practice, we shall conservatively assume that ram pressure stripping is efficient in removing the ISM of galaxies and thus that galaxies only exert gravitational influence.
Volumefilling Turbulence The volume filling factor of galaxy wakes is small (for a simple analytic estimate, see Subramanian et al. (2006)). This might seem to imply that the impact of turbulence excited by galactic wakes is confined to a small fraction of the cluster. However, orbiting galaxies can also resonantly excite gmodes, which from a formal WKBJ analysis have the dispersion relation:
(1) 
. The BruntVäisälä frequency for buoyant oscillations is:
(2) 
where is the fluid entropy , and applies if thermal conduction is negligible, while applies if the thermal conduction time is sufficiently short that a displaced blob’s temperature is determined by conductive rather than adiabatic cooling (Sharma et al., 2009b). Note that depend on the entropy and temperature gradient respectively; typically, is about a factor of 2 smaller than .
The above dispersion relation immediately implies that to obtain modes with real , the driving frequency ; otherwise the modes have imaginary radial wavenumber and are evanescent. Physically, one can always achieve a low frequency response by making the mode progressively more tangential, but it is impossible to drive the system at frequencies higher than the maximum response frequency of , corresponding to completely vertical oscillations. This thus implies that waves driven at frequencies can be resonantly excited, and must propagate inward toward the cluster center (as can be seen from their group velocity; Balbus &
Soker (1990)), where they will be trapped, reflected and focused inside the resonance radius where . A linear analysis by Balbus &
Soker (1990) showed that most of the power in modes is in the longest wavelengths
Isotropic Turbulence Turbulence in the fluid has to compete with buoyancy forces arising from stable stratification. One can show that the ratio of tangential and radial velocities is given by (e.g., see discussion in §2 of Ruszkowski & Oh (2010)):
(3) 
where is the eddy turnover frequency at a given scale, and Fr is the Froude number, which compares inertial and gravitational forces ( is the Richardson number). If , then turbulence is fundamentally 2D, and for instance it is difficult to rearrange magnetic fields in the radial direction. However, the level of turbulence required to overcome stable stratification is weak; for typical cluster conditions the critical turbulent velocity is (Sharma et al., 2009b):
(4)  
where is the gravitational acceleration in units of , is a characteristic scale height in units of kpc, and is the critical Richardson number; is typical for hydrodynamic flow.
At first blush, the requirement for might seem to be at odds with the requirement that for modes to be excited. However, note that for homogeneous Kolmogorov turbulence, ; it is therefore conceivable that low frequency gmodes can be excited on large scales, while highfrequency smallscale modes can overcome stabilizing buoyancy forces. Since our background state is not homogeneous, we have to resort to 3D simulations to verify if this expectation is indeed satisfied. This is a major goal of this paper.
Vorticity and Bfield growth Gmodes excite vorticity. An easy way to see this is to examine the vorticity evolution form of the momentum equation for waves (i.e., assuming ) (Lufkin et al., 1995):
(5)  
where is the vorticity, and to note that is in general nonradial, so that is nonzero (indeed, we see in Fig. 2 that since rises toward the center, that gmodes become progressively more tangentially biased there). This implies that vorticity is a good tracer of gmodes, a fact which we shall exploit. It also means that gmodes could conceivably drive an efficient dynamo. There is a wellknown analogy between the vorticity equation:
(6) 
where is the viscosity, and the relation for the magnetic field in the fluxfreezing limit,
(7) 
where is the electricity resistivity. This, together with the fact that the divergence of and both vanish, leads to the expectation that their growth might be related
(8) 
where is the rms turbulent velocity on large scales. The above estimate is consistent with observed G fields (Carilli &
Taylor, 2002), though there are considerable uncertainties
Magnetic tension Magnetic tension can inhibit the HBI (Quataert, 2008). For perturbation scales comparable to the radius (i.e., ) we obtain a critical value:
(9)  
for suppression, where the fiducial values are measured from our simulated cluster at a radius of kpc. Due to the similarity between the field values in equation (8) and (9), it has been suggested that magnetic fields amplified by turbulence can prevent the onset of the HBI (e.g., see discussion in Kunz et al. (2010). Note that their version of equation (9) yields somewhat lower Bfields than ours, for identical parameters. In any case, equation (9) is only an approximation as the derivation assumes the WKB approximation, while the nonlinear saturation of the HBI occurs on global scales). In this paper, we will deliberately ignore this possibility. Observationally, the strength of the magnetic field in the ICM is G and has a large scatter of about an order of magnitude within the ICM and between clusters (Carilli & Taylor, 2002); moreover, there are considerable observational uncertainties in these values, as mentioned above. Numerical simulations show that the HBI still develops for G fields (I. Parrish, priv. comm.), although it can be delayed for increased field strengths. Given the large uncertainty in whether observed field strengths are capable of stabilizing the HBI, past studies of HBI (e.g., Parrish et al. (2009); Bogdanović et al. (2009); Ruszkowski & Oh (2010); Parrish et al. (2010)) focused on the regime where the magnetic tension is unimportant. We also adopt the same approach here, and study if volumefilling turbulence alone can stabilize the HBI. More specifically, we consider plasma and note that, as long as the field is not dynamically important, its exact value does not play a role. In this case, the magnetic field strength scales out of the problem and only serves as a medium to redirect the heat flow via anisotropic thermal conduction. We can therefore study the effects of turbulence alone without the possibly confounding effects of magnetic tension.
Turbulent heating and heat diffusion Turbulence impacts the thermodynamics of the fluid through its effect on thermal conduction, both randomizing and amplifying the magnetic field. Both of these suppress the HBI, and allow thermal conduction at the Spitzer value. However, turbulence can also directly affect the thermal state of the plasma through dissipation of turbulent motions (direct heating), or allowing heat transport via turbulent diffusion (Kim & Narayan (2003b); Dennis & Chandran (2005), and references therein). The heating rate from dissipation of kinetic and magnetic energy is:
(10) 
where is a dimensionless constant of order unity and , the dominant velocity lengthscale, is unknown but almost certainly a function of radius; a reasonable ansatz might be (Dennis & Chandran, 2005), where is some adjustable constant of order unity, and is some minimal lengthscale. On the other hand, the heating rate from turbulent heat diffusion is:
(11) 
where is the specific entropy, and the turbulent diffusivity is:
(12) 
where the second factor of takes into account the damping of radial heat transport by buoyancy forces (Dennis & Chandran, 2005). The fact that is of order the hydrodynamic value even for a magnetized plasma was found in MHD simulations by Cho et al. (2003). Nonetheless, equation (11) should be understood to be only approximate, since it assumes that fluid elements are transported adiabatically, which need not be the case when anisotropic conduction is operating. In reality, both the thermal conduction diffusion coefficient and the turbulent heat diffusion coefficient can be comparable, and either could dominate in a specific situation.
Thermal conduction may indirectly assist with turbulent heat diffusion, as it reduces the impact of buoyancy forces (and thus reduces ). Indeed, simulations by Sharma et al. (2009) show that metal mixing in a stratified plasma is much more efficient once conduction is at play, allowing much broader metallicity profiles, for this very reason. Naively, if we think of gas entropy as a scalar to be advected by turbulent motions, similar conclusions should hold, although of course the interaction between heat transport and dynamics requires detailed simulations. We shall investigate the relative role of all these heating processes in our simulations.
3 Methods
3.1 Initial conditions for the gas
The details of the numerical setup are described in Ruszkowski & Oh (2010; hereafter RO10). Here we summarize key differences.
The cluster parameters used here are similar to those corresponding to coolcore cluster A2199. In addition to the NFW potential
of the cluster halo, we also include the contribution from the central brightest cluster galaxy (BCG), which was not included in RO10. The gravitational potential is described by the sum
of the term due to an NFW profile with a softened core
(13)  
where is the smoothing core radius ( kpc), kpc is the usual NFW scale radius, and the BCG contribution which has a King profile:
(14)  
where , kpc is the core radius for the BCG and km s is its lineofsight velocity dispersion. The parameter in equation (13) determines the cluster mass and is of the order of the total cluster mass, . We then solve the equation of hydrostatic equilibrium assuming the entropy distribution as parametrized by Cavagnolo et al. (2009); see equations (16) & (17) of RO10. Note that we do not include the gravitational contribution from other galaxies (§3.2) in our initial conditions, so the system is not initially in full hydrostatic equilibrium. However, after an initial transient, it rapidly relaxes to a new equilibrium configuration.
The addition of the BCG has two effects. Firstly, due to the increased gravitational acceleration, it results in higher gas densities compared to the models we considered in Ruszkowski & Oh (2010). This allows for a more conservative analysis of the effect of cooling. In fact, the central density here is a factor of times higher, which, combined with a slightly lower assumed central temperature, results in a central cooling time which is nearly 5 times shorter. The higher adopted central density in this paper is in line with that observed in A2199 (Johnstone et al., 2002). Given this more stringent setup, some of the stable models in Ruszkowski & Oh (2010) would actually undergo a cooling catastrophe. Secondly, the change in the gravitational potential has consequences for the BruntVäisälä frequency and the trapping of gmodes, as we discuss below.
The initial distribution of density and temperature is shown in Figure 1. The frequency of circular orbits and the BruntVäisälä frequencies for a hydrodynamic and conducting fluid with this initial density and temperature profile are shown in Figure 2. For a mode with a given value of , modes can be resonantly excited if . This therefore defines an outer trapping radius for such a mode. Note that both and are both strong functions of the gravitational potential. We have directly verified the resonance condition by running simulations both with and without the central cD galaxy; in the latter case, orbiting galaxies fail to excite volumefilling turbulence, which is to be expected since falls inward in this case, and the resonance condition is never satisfied (see also Lufkin et al. (1995); Kim (2007)). Note that fineturning of the resonance condition is not necessary: the resonance is not very sharp (Balbus & Soker, 1990), and in practice galaxies with noncircular orbits excite modes with a variety of harmonics, some of which can potentially fall below .
The magnetic field setup was identical to that in Ruszkowski & Oh (2010): we generate statistically isotropic randomphase complex fields in Fourier space, with 3D Fourier amplitudes given by:
(15) 
as appropriate for Kolmogorov turbulence, where and kpc is a smoothing wavelength. We then apply a divergence cleaning operator in space, and then inverse Fourier transform the field back to real space.
3.2 Initial conditions for the galaxies
The simulations must be initialized with a galaxy population, which has the appropriate spatial distribution, masses and velocities. Rather than relying upon cosmological simulations, we use an empirically grounded approach, which also has the advantage of speed and flexibility. How are the galaxies spatially distributed? From a sample of Kband selected galaxies within 93 clusters and groups, Lin et al. (2004) find that the galaxy number density profile in clusters is well described by the NFW profile (Navarro et al., 1997) with a concentration parameter , with no evidence for cluster mass dependence of the concentration. The theoretical justification for galaxies tracing the NFW profile is somewhat equivocal. If one attempts to use cosmological simulations to set up initial conditions, the radial distribution of subhalos in simulations is wellknown to be less concentrated than the darkmatter, or ’antibiased’ (Nagai & Kravtsov (2005), and references therein). This is due to tidal stripping of subhalos in the central regions. On the other hand, simulations that include galaxy formation allow subhalos to be selected by stellar mass. This generally shows closer agreement with observed profiles, as the stellar mass (which is tightly bound) remains conserved while the dark matter is stripped from outer regions (Nagai & Kravtsov, 2005). Other investigators find that the fraction of such stellardominated halos is small, but caution that numerical resolution effects may preclude robust conclusions at this point (Dolag et al., 2009). Overall, we therefore simply employ the observational result that galaxies trace the NFW profile.
As for the galaxy masses, instead of using a Schechter function, we simply assume (as did, for instance, Subramanian et al. (2006); Kim (2007)) that all galaxies have the same mass. This is for two reasons. Firstly, this allows us to rapidly explore the effect of varying galaxy masses (due, for instance, to different efficiencies of tidal stripping). The assumption of a characteristic mass is reasonable: since dynamical friction scales as , most turbulent motions are induced by galaxies of mass , where most of the mass resides, rather than the more abundant lower mass galaxies. Indeed, we shall find that the induced gas motions are mostly sensitive to the mass of galaxies, and less sensitive to their number (§4.1). Previous hydrodynamic simulations found unchanged results with galaxies drawn from a Schechter distribution, if the characteristic break mass (Kim, 2007). Secondly, it allows us to directly calibrate against lensing estimates for subhalo mass fraction. Unlike Kband surveys, lensing is directly sensitive to total mass, but is generally only sensitive to subhalos with . Natarajan et al. (2009) find from the massive lensing cluster Cl 0024+16 that of the cluster mass can be attributed to substructure with , with typical masses (with a weak radial trend such that galaxies in the outer regions are more massive; see their Fig 6). Their results, including the mass function as a function of radius, is broadly consistent with the results of the Millenium simulation run, except that the typical masses of galaxies is lower in simulations by a factor of . This is subject to the uncertainties of extra binding due to a compact stellar halo mentioned above; note that masses of is also consistent with other observations from lensing (Shin et al., 2008) and galaxy wakes (Sakelliou et al., 2005). Below we explore a grid of models with varying galaxy mass, but never allowing the total substructure mass fraction to rise above . For simplicity in the code, the galaxies are modeled as point masses. Since we are primarily concerned with the excitation of gmodes on scales much larger than kpc galactic scales, we do not expect this simplification to significantly impact our results.
Given these assumptions, the most rigorous way to initialize galaxy velocities is to directly construct the distribution function from the density profile, using Eddington’s formula (Binney & Tremaine, 2008; Kazantzidis et al., 2004). However, velocity anisotropy is only easily incorporated in such models if it has certain parametric forms, as for instance in OsipkovMerritt models. Instead, we construct a selfconsistent velocity model via the local Maxwellian approximation: approximating the velocity dispersion tensor by a multivariate Gaussian at each point, with dispersions given by the solution of the Jeans equation (Hernquist, 1993). This has the virtue of simplicity and flexibility. Note that such models may not be in strict equilibrium, and can demonstrate evolution (Kazantzidis et al., 2004). However, Springel et al. (2005) find the actual amount of relaxation to be small; furthermore, Faltenbacher et al. (2005), who directly simulate the motion of galaxies in clusters, find their velocity distribution is indeed closely Maxwellian, with good agreement between simulation results and equilibrium Jeans equation solutions. We therefore solve the Jeans equation assuming no rotational support or bulk streaming ():
(16)  
where is the velocity anisotropy
parameter
(17)  
is the galaxy number density and is the combined cluster + cD galaxy
gravitational potential. Note that we have not built selfconsistent models and ignore the contribution of galaxies to the gravitational potential; for a large subhalo mass fraction, the system is not in full equilibrium. In practice, this is a small effect, and the galaxy distribution does not evolve significantly over the course of our simulation.
What are appropriate assumptions for
? It may be estimated from observations via Jeans equation
modelling, given knowledge of galaxies positions, the cluster
potential, and line of sight velocities. A detailed study of 10
clusters using a spectroscopic sample of galaxies from SDSS and 2dF
found galaxy orbits to be isotropic within the errors for most
clusters (Hwang &
Lee, 2008). An earlier paper, using ENACS data, found
that the brightest ellipticals do not yield an equilibrium solution,
while other ellipticals, SOs and early spirals have isotropic orbits, and
late spirals prefer radial to isotropic orbits
(Biviano &
Katgert, 2004). Overall, we assume isotropy , and
regard this as our default model. In passing, we note that one could
easily incorporate the effects of the velocity anisotropy by, for example,
considering fits to to measurements in simulations (Hoeft et al. (2004)).
Given that the evidence for orbital anisotropy in observations is marginal to date, we defer the study of the effect of such orbital distributions to future work.
We also note that preferentially radial orbits would enhance the restoration of conduction even further and strengthen our conclusions.
We solve equation (3.2) as an initial value problem, where . Having solved for and , we randomly sample from the
multivariate Gaussian at each position to create a realization of the
velocity field. The above procedure allows us to initialize the
simulations with a realization of galaxies with both masses and sixdimensional
phase space coordinates (position and velocity).
3.3 Simulation
The simulations were performed using the FLASH code (version 3.2). FLASH is a modular, parallel adaptive mesh refinement magnetohydrodynamic
code. Magnetic field evolution was solved by means of a directionally
unsplit staggered mesh algorithm (USM; Lee
et al. (2009)).
The USM module is based on a finitevolume, highorder Godunov scheme combined with constrained
transport method (CT). This approach guarantees divergencefree magnetic field distribution. We implemented the anisotropic conduction unit
following the approach of Sharma &
Hammett (2007). More specifically, we applied monotonized central (MC) limiter to the conductive fluxes.
This method ensures that anisotropic conduction does not lead to negative temperatures in the presence of steep temperature
gradients.
The threedimensional computational domain was approximately 1 Mpc on each side,
enclosing a large fraction of the cluster.
The central regions of the cluster had an enhanced refinement level. The maximum spatial resolution for 6 levels of refinement was
kpc.
The simulations were performed on a 384processor cluster located at the Michigan Academic Computing
Center at the University of Michigan in Ann Arbor and on the Columbia supercomputer at NASA Ames.
4 Results
We performed a total of 16 runs including radiative cooling, anisotropic thermal conduction and selfgravitating particles to emulate the gas “stirring” by galaxies. We considered a uniform grid of parameters: 50, 100, 150, and 200 galaxies characterized by masses of (0.3, 0.6, 0.9, 1.2) M. With our cluster mass of , these parameters corresponds to a mass fraction in galaxies ranging from to a maximum of . For instance, for galaxies with mass , our grid corresponds to . We also performed two control runs: one without the galaxies (and hence without stirring) to isolate the effect of heat buoyancy instability, and one without conduction to isolate the effect of dynamical friction heating by galaxies.
4.1 Gas velocities and volume filling of turbulence
Figure 3 (left panel) shows the evolution of the velocity dispersion measured within 100 kpc from the cluster center. Thin blue (red) lines are for 100 (200) galaxies respectively and for equally spaced masses ranging from to . The mass increases gradually from the lightest to the darkest color. The black dashed line is for the pure HBI case. The HBI case and lightercolored curves are evolved for shorter times. These runs suffer from overcooling and the central temperatures reaches the low temperature threshold at which point the simulation is stopped. The right panel in Figure 3 shows the median velocity within 100 kpc; the color coding corresponds to that in the left panel. It is clear from these figures that there is a clear trend for the velocity dispersion or the median velocity to increase with the typical galaxy mass. A similar, albeit weaker, trend is seen for the galaxy number. This is consistent with the findings of Kim (2007) who found in pure hydrodynamic simulations, that the gas velocity dispersion scales as , where and are the number and mass of galaxies, respectively. Note that a scaling is consistent with dynamical friction in the linear regime, since for dynamical friction.
As and increase, the cooling catastrophe is delayed, and is completely staved off at the upper envelope of these parameters. In this respect our MHD simulations differ markedly from those of Kim (2007), who found that a cooling catastrophe was inevitable in purely hydrodynamic simulations, for all portions of parameter space. We explore these differences further in §4.5.
Note that in our case the velocity dispersion seems to increase more slowly than . Besides the inclusion of MHD in our simulations, differing results could be due to a variety of factors, including the different assumed distribution of galaxies. Note that the stated number of galaxies are distributed over the entire cluster; the number of galaxies in the inner regions which actually result in trapped gmodes is actually considerably smaller, and subject to Poisson fluctuations. Also, the introduction of more galaxies and/or increasing their mass does not cause the velocity dispersion to increase without limit; instead, the growth in velocity dispersion appears to saturate. Kim (2007) also observed this in his hydrodynamic simulations, and attributed it to loss of resonant excitation once density fluctuations become large and the background is nonlinear. We see the same saturation on the same yr timescale, in simulations with driven volumefilling turbulence, where resonant excitation of modes is not an issue (see §3.4 of RO10).
The asymptotic velocities of , while generally insufficient for turbulent heating to be important, is enough to restore thermal conduction and enable turbulent heat diffusion.
The comparison between
the velocity dispersion and median velocity reveals that that both of these quantities are comparable, as might be expected for volume filling turbulence.
4.2 Bias in magnetic field orientation
The evolution of the anisotropy in the orientation of magnetic fields is shown in Fig. 4. The definition of this parameter is similar to that for galaxy velocity defined in Eq. (16). The only difference is that the velocity dispersions are replaced by magnetic field dispersions. Thus, corresponds to isotropic magnetic fields, whilst corresponds to progressively more tangential (radial) fields respectively. Thin blue lines are for 100 galaxies and for equally spaced masses ranging from M to M. The galaxy mass increases gradually from the lightest to the darkest color. Thick red curves are the corresponding lines for 200 galaxies. The black dashed line is for the pure HBI case. As expected, when stirring is weak, the HBI prevails, leading to a systematic tangential bias in the orientation of magnetic fields. This insulates the core against thermal conduction, leading to a cooling catastrophe. The fields become more tangential with time and the cluster eventually suffers from overcooling. On the other hand, for increasingly vigorous stirring (i.e., increasing the individual masses or number of galaxies), the field becomes increasingly isotropic, and a cooling catastrophe is averted.
These results are consistent with the driven turbulence simulations in RO10
4.3 Evolution of Gas Temperature and Entropy
In Figure 5 we show the evolution of temperature profiles. This figure is for the models where heating is more efficient. Specifically, it corresponds to the following pairs of parameters (150,1.2), (200,0.9), (200,1.2), where the first number in the parenthesis is the number of galaxies and the second is the galaxy mass in M. Progressively older profiles correspond to systematically brighter colors. The final time corresponds to 5 Gyr and the curves are plotted every 0.1 Gyr. As is clearly seen in this figure, these models do not lead to the cooling catastrophe. Several features are of interest. The temperature profile does not asymptote toward an isothermal profile, as is generically the case when thermal conduction alone offsets cooling (Bregman & David, 1988; Guo & Oh, 2008; Conroy & Ostriker, 2008). Despite the fact that we have not introduced an additional source of central heating such as an AGN, the cluster is able to remain in a thermally stable CC (i.e., with a central temperature which is lower than at the cooling radius) state via heat transport from the outer heater reservoir alone. Without finetuning, this is impossible to achieve with thermal conduction alone (when the cluster either becomes isothermal or undergoes a cooling catastrophe). Finally, the temperature profile is not always monotonic, but occasionally increases inward—a situation which is thermodynamically impossible if thermal conduction alone is operating. Note that these fluctuations are transient; such reversals are not present in the later stages of the evolution (progressively lighter blue curves correspond to later times). As we shall see in §4.5, all of these features hint that an additional heat transport process is at play: turbulent heat diffusion.
Figure 6 shows the temperature evolution for the parameters where the heating is least efficient. From left to right shown are: (100, 0.3), (50, 0.6), (50, 0.3). Here, the profiles are shown more frequently then in Figure 3 to better capture the evolution of the system just before the imminent cooling catastrophe. For the same reason, we also extend the radial scale to smaller distances from the cluster center to show how the system becomes thermally unstable. It is evident that in all three cases, the cluster quickly evolves toward a cooling catastrophe. In the final stages of the process, the cooling is so fast at the very center that the gas accretion accelerates so much that adiabatic compression in the shells surrounding the center can heat the gas up (e.g., see last profile in the right panel).
Finally, in Fig 7 we show the entropy profiles (where entropy is defined as ) for the strong heating models. The central entropy grows somewhat, consistent with the rise in temperature, and as might be expected if heating by conduction and/or turbulent heat diffusion were taking place. However, these profiles show that turbulent convection/stirring is still a relatively gentle process; we do not see the flat isentropic central profile which might be expected if turbulent convection were extremely efficient. Instead, the fluid always remains stably stratified by entropy, which steadily increases outward at all times.
As discussed above the cluster will develop a cooling catastrophe if the number of galaxies and/or their masses are too small. The time it takes for the
cluster to reach this point (essentially the effective cooling time) is plotted in Figure 8 as a function of galaxy number and mass.
The contours are plotted every Gyr. The models that exhibit the effective cooling time of 6 Gyr (the maximum simulation run time) are thermally stable. We point out that in practice
the models that possess cooling times to 4 Gyr could be considered stable as they are likely to experience cluster mergers that may reset
the conditions in the ICM and further slow down or essentially delay the cooling process. In any case, as can be seen in Figure 5, a substantial
fraction of the models shows appreciably long effective cooling times. As a technical note, we add that the reason for the lack of monotonicity
in some of the contour lines as a function of galaxy number is that a single random seed was used to generate the conditions for a given number
of galaxies and varying galaxy masses.
4.4 Generation of vorticity and magnetic fields
As we have previously seen, modes must be excited for stirring by galaxies to excite volumefilling turbulence. These modes also induce vorticity (equation 2). Vorticity is therefore an excellent tracer of the growth of modes. We compute the evolution of vorticity in the central 100 kpc to assess if modes are indeed generated and trapped. Figure 9 (left panel) shows the evolution of the square of the scaled vorticity for the same set of parameters as in Figure 3 that shows velocity dispersion and median velocity. The scaled vorticity is defined as , where kpc and km s are the reference lengthscale and velocity, respectively. Thin blue lines are for 100 galaxies and red ones for 200 galaxies. Galaxy gasses range from M to M and are uniformly sampled (lighter colors are for lighter galaxies). The black dashed line corresponds to the pure HBI case. A clear trend for the vorticity to increase with time is seen in this figure, suggesting that modes are present and at least partially trapped, leading to the volumefilling turbulence seen.
As discussed in §2, a growth in vorticity might also lead to growth in the magnetic field; the possibility that magnetic fields could be turbulently amplified in clusters has been repeatedly raised (e.g., Ruzmaikin et al. (1989); Subramanian et al. (2006); Ryu
et al. (2008); Cho et al. (2009)). Whilst a detailed study is beyond the scope of this paper, we check whether these theoretical expectations are satisfied in our simulations. Fig. 9 shows that the magnetic energy density indeed grows in tandem with vorticity, with more vigorous stirring corresponding to greater field amplification. However, the characteristic growth time appears to be somewhat longer. Note that the simulations were initialized with extremely small magnetic fields: the initial plasma beta , where is typically measured in the ICM. These small initial fields were for computational convenience (since the MHD approximation is satisfied with a trivially small magnetic field), and to ensure that magnetic fields never become dynamically important
4.5 Relative contribution to gas heating
In §4.3 we noted a number of interesting features in the temperature profiles of our stable clusters. They remained stable CC clusters, neither becoming isothermal nor developing cooling catastrophes, as clusters stabilized solely by thermal conduction generally do. Furthermore, the central temperature showed timedependent oscillations, sometimes becoming hotter than gas further out. A temperature inversion would not happen if only thermal conduction was at play. This behooves us to take a closer look at what actually stabilizes the customary thermal runaway. We have already discussed the effect that turbulence can have on thermal conduction, by tangling field lines and countering the HBI. However, turbulence itself can be a source of heating, either via viscous dissipation of turbulence, or turbulence diffusion of high entropy gas into low entropy regions (e.g., Dennis & Chandran (2005), and references therein). Let us examine these in turn.
As long as there is sufficient separation of scales that an inertial range can develop (such that the energy per unit mass per unit time is independent of scale), the heating rate from dissipation of turbulent motions is independent of the nature of viscosity. In particular, it is unimportant if our numerical viscosity is different from the actual physical viscosity in the ICM. The heating rate per unit volume due to dissipation of such motions is (Dennis & Chandran, 2005):
(18) 
where is the dominant lengthscale, is the energy density in turbulence, and is the eddyturnover time on the dominant lengthscale. To estimate , we can note that vorticity has units of , and that our scaled vorticity in Fig 9 is . This implies
(19) 
Consistently, note that the vorticity in Fig 9 indeed takes Gyr to rise to its asymptotic value. This implies that the heating time for turbulent dissipation of motions is:
(20) 
where is the thermal energy density, and we have defined the turbulent Mach number (note that our quoted velocities are in 3D). While there are factors of order unity uncertainty, it is clear that the mild subsonic motions we explore are a negligible source of heating via viscous dissipation (and consistent with other estimates; Dennis & Chandran (2005)). This also implies that dynamical friction heating due to galaxy motions (ElZant et al., 2004; Kim et al., 2005; Kim, 2007; Conroy & Ostriker, 2008; Birnboim & Dekel, 2010) is not the source of heating which averts the cooling catastrophe in these simulations.
On the other hand, turbulent heat diffusion is not negligible. One can estimate its contribution from a simple mixing length prescription as in equation 11; this shows that it can be at least comparable to and may exceed the thermal conduction contribution. However, the coefficient of turbulent diffusivity, , is only approximate and subject to order unity corrections. Since we have full knowledge of the density, velocity, temperature and magnetic fields in our simulations, we can attempt to directly compute the heating contributions from thermal conduction and turbulent heat convection. In particular, at a radius we can calculate the inward heat flux due to conduction:
(21)  
where is a unit vector pointing in the direction of the magnetic field and is the SpitzerBraginskii conduction coefficient given by erg scmK, as well as the convective heat flux (Parrish et al., 2008):
(22) 
where is the spatial average of quantity in the shell and is the local deviation of that quantity from its average; generally the second term is dominant. We can then compare these to the total rate of energy loss within radius due to radiative cooling. We can also compute the volumetric heating rate due to these two processes, via , , although these are of course much noisier quantities.
It is useful perhaps to begin by considering a case where the properties of turbulence are well known: the ’strong’ driven turbulence case of RO10, which has volume filling turbulence by construction, and rms velocities of . Conductive (solid line) and convective (dashed line) heating to cooling ratios as a function of time for the ICM within 100 (50) kpc from the cluster center are shown in the upper left (right) panel of Fig 10. It is clear that conduction only contributes of the heat necessary to overcome cooling and convective heat flow is an important part of the energy budget; indeed, in the central regions turbulent advection of heat is the dominant heating process (note that the cluster is not in complete equilibrium, so the sum of the two ratios is not necessarily unity). The convective heat flow shows timedependent fluctuations, as might be expected. In the bottom panels, we show the volumetric heating to cooling ratios as a function of radius, for conduction (left panel), and turbulent convection (right panel). The curves are plotted every 100 Myr; progressively lighter colors denote later times. The divergence of heat fluxes is a much noisier quantity, as reflected in the plots. Nonetheless, it is clear from the plots that convective heating dominates near the center, whilst conductive heating dominates further out. This reminiscent of stable hybrid AGN+conduction heating models (Ruszkowski & Begelman, 2002) where the AGN heats the cluster center and conduction is important further out.
In Fig. 11, we show the same plots, but for the case where turbulence is due to stirring by galaxies. All results presented in this figure are for the case of 200 galaxies, each with M; this is stable against a cooling catastrophe. As before, conduction is only a fraction of the energy budget. However, in this case convective heating shows dramatic oscillations as a function of time; the amplitude of the oscillations near the center is much larger than in the driven turbulence case . The reason for this is that the dominant lengthscales of motion are comparable or larger than the depicted radii, as might be expected if modes are excited (since most of the energy in modes are in the largest lengthscales, comparable to the trapping radius). For instance, from §4.4, a typical lengthscale on which vorticity is excited is . Since fluctuations in the velocity field span larger scales than the ones under interest, our calculation of will show strong time dependence (however, the calculations of the conductive heat flux are of course still valid. Note that has to be unity on average, since the cluster is stabilized against a cooling catastrophe). As noted earlier, Poisson fluctuations in the number of galaxies in the core will also drive timedependent fluctuations in the velocity field. The gas is sloshing in the potential well; we observe this directly too in the simulations, as the gas pressure maximum wanders in time from the center of the potential well. Nevertheless, despite the breakdown of equation 22 in a rigorous sense, it is clear from the amplitude of fluctuations in the bottom panels of Fig 11 that (as in the driven turbulence case) conductive heating increases outwards in radius, while convective heating is more important near the center. In particular, the dominance of convective heating near the center, and its positive and negative fluctuations, allow us to understand the fluctuating temperature profiles seen in Fig 5. Since conduction is only a part of the energy budget, there is no reason for the stabilized temperature profile to approach isothermality. Furthermore, the reason why the central temperature gradient can occasionally become inverted (with the center hotter than its surroundings) is clear: if a highentropy fluid element is compressed at the center, this will result in higher central temperatures. While thermal conduction seeks to make the fluid isothermal (since heat flows down the temperature gradient), turbulent heat diffusion seeks to make the fluid isentropic (since heat flows down the entropy gradient). In this sense, the subsonic turbulence induced by galaxies results in only mild convection, since as seen in Fig 7, the gas remains convectively stable with entropy increasing monotonically outward. Whilst we have not directly calculated the diffusion of metals directly, this also suggests that metal mixing to larger radii will be somewhat enhanced (so that metals will have a broader distribution than the galaxies), but not greatly so. Indeed, a mixinglength theory calculation of metal dispersal via turbulent diffusion by Rebusco et al. (2005), who assumes levels of turbulence very similar to those we have simulated, shows excellent agreement with observations.
Could turbulent heat diffusion alone stabilize a thermal runaway? We tested this hypothesis by running purely hydrodynamic simulations of our most vigorous stirring case (200 galaxies of mass ), without thermal conduction. The results are shown in Fig. 12. The cluster rapidly undergoes a cooling catastrophe (consistent with the results of Kim (2007)), demonstrating that thermal conduction is an essential element in stabilizing the cluster. Note that this model already represents the most extreme choice in parameter space of the amount of stirring possible by galaxies. Also in line with these results, purely hydrodynamic simulations of ’sloshing’ZuHone et al. (2010) show that the cooling catastrophe can be delayed by not disrupted. Besides providing the dominant source of heating in the outer regions, thermal conduction also reduces stabilizing buoyancy forces (as discussed in §2) and thus enables more rapid, efficient mixing and turbulent heat diffusion. Passive scalars such as metals are more efficiently advected in the presence of conduction (Sharma et al., 2009a), and the same is likely true of the advection of entropy. Thus, intriguingly, whilst neither process alone can stabilize the cluster, the interplay between turbulence and conduction does permit stability: turbulence enables conduction, and conduction enables turbulence.
5 Conclusions
Using threedimensional MHD simulations, we have studied the effect of anisotropic thermal conduction and stirring motions due to galaxies orbiting in the cluster potential on the effective cooling rate in cluster cool cores. Such galaxies excite mild subsonic turbulence with . We find that a combination of thermal conduction and turbulent heat transport can stabilize the cluster, for realistic parameter choices consistent with gravitational lensing observations of substructure in clusters. Unlike much previous work, there is no subgrid physics in our simulations: we do not invoke subgrid prescriptions for the topology of the magnetic field (which affects the effective thermal conductivity), the magnitude and volumefilling factor of turbulence, which is calculated directly from the gravity/hydro solver (unlike previous work (Ruszkowski & Oh, 2010; Parrish et al., 2010) in which volumefilling turbulence is inserted by hand), or turbulent heat diffusion (which is directly simulated). We have also simulated a cluster with significantly higher central density than in RO10, and still found it to be thermally stable. Other salient points include:

In order for galaxies to excite volumefilling turbulence, rather than have turbulence confined to galactic wakes, they must excite modes, which requires that , where is the BruntVäisälä frequency appropriate when thermal conduction timescales are rapid (equation 2). On the other hand, overwhelming the stabilizing buoyancy forces to randomize the magnetic field requires that . These two requirements can be simultaneously satisfied since for Kolomogorov turbulence; hence, the lowfrequency, large scale modes can be trapped, while the highfrequency, small scale modes overcome the HBI.

We observed strong growth in vorticity, which is a good tracer of the growth of modes. We also observed turbulent amplification of Bfields in tandem with vorticity.

Thermal conduction provided about of the heating budget, with the rest coming from turbulent heat diffusion. Viscous dissipation of turbulent motions (and hence dynamical friction heating) is negligible. Turbulent heat diffusion tends to be more important in the center of the cluster, while conduction plays a greater role further out. The predominance of turbulent heat diffusion in the center—which is powered by motions on large scales—implies that it exhibits oscillations about the equilibrium temperature profile, and can occasionally exhibit small temperature inversions as high entropy fluid elements are compressed near the center. However, conduction plays a crucial part of the story; our most extreme stirring case still suffered a cooling catastrophe if thermal conduction was omitted. Besides supplying heat further out in the cluster, conduction also reduces stabilizing buoyancy forces and enables more efficient turbulent heat diffusion. It appears that turbulence enables conduction to operate, as well as viceversa. The details of the interplay between turbulence and conduction, as well as the diffusion of metals in our stirring simulations, are interesting topics for future work.
In this paper, we have focussed on a timesteady source of turbulence—stirring by galaxy motions—but we stress that other intermittent sources of turbulence such as mergers or AGN outbursts, can also contribute. Indeed, a sudden rise in heat transport processes such as conduction and turbulent heat diffusion due to an increase in turbulence could effect a CC to NCC transition (Guo & Oh, 2009; Ruszkowski & Oh, 2010; Parrish et al., 2010). Other processes which could reorient field lines in galaxy cluster include rising bubbles, which could amplify and straighten magnetic fields in their wake (Ruszkowski et al., 2007; Guo et al., 2008; Bogdanović et al., 2009). In the future, observations of Faraday rotation by SKA (Bogdanovic et al., 2010) or magnetic draping around galaxies orbiting the cluster center (Pfrommer & Dursi, 2010), could probe the topology of magnetic field lines and test these ideas. Finally, these ideas about the interplay between between the thermal conduction, the HBI and turbulence in the inner regions of the cluster also apply with equal force to the interplay between conduction, the MTI and turbulence in the outer regions of the cluster, which we present elsewhere (Ruszkowski et al 2010).
Acknowledgments
The software used in this work was in part developed by the DOEsupported ASC/Alliance Center for Astrophysical Thermonuclear Flashes at the University of Chicago. MR thanks Jeremy Hallum for his invaluable help with maintaining the computing cluster at the Michigan Academic Computing Center where most of the computations were performed. We are indebted to Dongwook Lee, the author of the MHD module in FLASH for letting us use a proprietary threedimensional set of MHD modules. MR thanks Justin Nieusma for technical assistance with performing the simulations. We thank Ian Parrish, Eliot Quataert, Elena Rasia, Prateek Sharma, MinSu Shin, Maxim Markevitch, John ZuHone, Paul Nulsen, Aneta Siemiginowska, Christine Jones, Larry David, Bill Forman, Renato Dupke, Milos Milosavljević, John ZuHone, Matt Kunz, and Alex Schekochihin for discussions. SPO acknowledges support by NASA grant NNG06GH95G, and NSF grant 0908480. MR and SPO thank Institute of Astronomy, Cambridge, UK and Max Planck Institute for Astrophysics, Garching, Germany for their hospitality.
Footnotes
 Although a WKBJ analysis formally breaks down in this regime, a subsequent numerical study (Lufkin et al., 1995) showed that many of the linear theory results are still valid.
 Note, however, that this analogy is imperfect, since , which leads to a nonlinear coupling in the equations, whereas no such relation exists between and .
 In general, estimates based on rotation measure (RM) lead to stronger magnetic fields, while those based on synchrotron and inverse Compton (IC) analysis give weaker fields. However, RM methods may overestimate fields if singlescale magnetic field correlation length is used (Newman et al., 2002) or when the smallscale fluctuations in density and magnetic field are correlated in a turbulent medium (Beck et al., 2003). Moreover, these estimates depend on whether radio sources used to probe the field strength are embeded in the ICM, with smaller values inferred when background sources rather than embedded ones are used (Carilli & Taylor, 2002).
 From cosmological simulations, Benson (2005) finds that radial and tangential velocities can be correlated (at least at the time of merger), a detail we ignore.
 Although note that all but one of the simulations in RO10 were adiabatic simulations; by contrast, all the simulations presented here simultaneously include radiative cooling.
 Thus, as least in these simulations, the HBI is stabilized by the stirring motions and not by magnetic tension.
References
 Ascasibar Y., Markevitch M., 2006, ApJ, 650, 102
 Balbus S. A., 2000, ApJ, 534, 420
 Balbus S. A., Soker N., 1990, ApJ, 357, 353
 Beck R., Shukurov A., Sokoloff D., Wielebinski R., 2003, A&A, 411, 99
 Benson A. J., 2005, MNRAS, 358, 551
 Binney J., Tabor G., 1995, MNRAS, 276, 663
 Binney J., Tremaine S., 2008, Galactic Dynamics: Second Edition. Princeton University Press
 Birnboim Y., Dekel A., 2010, ArXiv eprints
 Biviano A., Katgert P., 2004, A&A, 424, 779
 Bogdanovic T., Reynolds C., Massey R., 2010, ArXiv eprints
 Bogdanović T., Reynolds C. S., Balbus S. A., Parrish I. J., 2009, ApJ, 704, 211
 Bregman J. N., David L. P., 1988, ApJ, 326, 639
 Carilli C. L., Taylor G. B., 2002, ARA&A, 40, 319
 Cavagnolo K. W., Donahue M., Voit G. M., Sun M., 2009, ApJS, 182, 12
 Cho J., Lazarian A., Honein A., Knaepen B., Kassinos S., Moin P., 2003, ApJ, 589, L77
 Cho J., Vishniac E. T., Beresnyak A., Lazarian A., Ryu D., 2009, ApJ, 693, 1449
 Churazov E., Forman W., Jones C., Sunyaev R., Böhringer H., 2004, MNRAS, 347, 29
 Churazov E., Sunyaev R., Forman W., Böhringer H., 2002, MNRAS, 332, 729
 Conroy C., Ostriker J. P., 2008, ApJ, 681, 151
 Dennis T. J., Chandran B. D. G., 2005, ApJ, 622, 205
 Dolag K., Borgani S., Murante G., Springel V., 2009, MNRAS, 399, 497
 ElZant A. A., Kim W.T., Kamionkowski M., 2004, MNRAS, 354, 169
 Enßlin T. A., Vogt C., 2006, A&A, 453, 447
 Evrard A. E., 1990, ApJ, 363, 349
 Fabian A. C., Sanders J. S., Allen S. W., Crawford C. S., Iwasawa K., Johnstone R. M., Schmidt R. W., Taylor G. B., 2003, MNRAS, 344, L43
 Faltenbacher A., Kravtsov A. V., Nagai D., Gottlöber S., 2005, MNRAS, 358, 139
 Goldshmidt O., Rephaeli Y., 1993, ApJ, 411, 518
 Guo F., Oh S. P., 2008, MNRAS, 384, 251
 Guo F., Oh S. P., 2009, MNRAS, 400, 1992
 Guo F., Oh S. P., Ruszkowski M., 2008, ApJ, 688, 859
 Hernquist L., 1993, ApJS, 86, 389
 Hoeft M., Mücket J. P., Gottlöber S., 2004, ApJ, 602, 162
 Hwang H. S., Lee M. G., 2008, ApJ, 676, 218
 Iapichino L., Niemeyer J. C., 2008, MNRAS, 388, 1089
 Johnstone R. M., Allen S. W., Fabian A. C., Sanders J. S., 2002, MNRAS, 336, 299
 Kazantzidis S., Magorrian J., Moore B., 2004, ApJ, 601, 37
 Kim W., 2007, ApJ, 667, L5
 Kim W., ElZant A. A., Kamionkowski M., 2005, ApJ, 632, 157
 Kim W., Narayan R., 2003a, ApJ, 596, L139
 Kim W.T., Narayan R., 2003b, ApJ, 596, 889
 Kunz M. W., Schekochihin A. A., Cowley S. C., Binney J. J., Sanders J. S., 2010, ArXiv eprints
 Lee D., Deane A. E., Federrath C., 2009, in N. V. Pogorelov, E. Audit, P. Colella, & G. P. Zank ed., Astronomical Society of the Pacific Conference Series Vol. 406 of Astronomical Society of the Pacific Conference Series, A New Multidimensional Unsplit MHD Solver in FLASH3. pp 243–+
 Lin Y.T., Mohr J. J., Stanford S. A., 2004, ApJ, 610, 745
 Lufkin E. A., Balbus S. A., Hawley J. F., 1995, ApJ, 446, 529
 Markevitch M., Gonzalez A. H., David L., Vikhlinin A., Murray S., Forman W., Jones C., Tucker W., 2002, ApJ, 567, L27
 McCarthy I. G., Babul A., Bower R. G., Balogh M. L., 2008, MNRAS, 386, 1309
 McNamara B. R., Nulsen P. E. J., 2007, ARA&A, 45, 117
 Nagai D., Kravtsov A. V., 2005, ApJ, 618, 557
 Nagai D., Kravtsov A. V., Kosowsky A., 2003, ApJ, 587, 524
 Narayan R., Medvedev M. V., 2001, ApJ, 562, L129
 Natarajan P., Kneib J., Smail I., Treu T., Ellis R., Moran S., Limousin M., Czoske O., 2009, ApJ, 693, 970
 Navarro J. F., Frenk C. S., White S. D. M., 1997, ApJ, 490, 493 (NFW)
 Newman W. I., Newman A. L., Rephaeli Y., 2002, ApJ, 575, 755
 Norman M. L., Bryan G. L., 1999, in H.J. Röser & K. Meisenheimer ed., The Radio Galaxy Messier 87 Vol. 530 of Lecture Notes in Physics, Berlin Springer Verlag, Cluster Turbulence. pp 106–+
 Parrish I. J., Quataert E., Sharma P., 2009, ApJ, 703, 96
 Parrish I. J., Quataert E., Sharma P., 2010, ApJ, 712, L194
 Parrish I. J., Stone J. M., 2005, ApJ, 633, 334
 Parrish I. J., Stone J. M., Lemaster N., 2008, ApJ, 688, 905
 Petrosian V., 2001, ApJ, 557, 560
 Pfrommer C., Dursi J. L., 2010, Nature Physics, 6, 520
 Quataert E., 2008, ApJ, 673, 758
 Rebusco P., Churazov E., Böhringer H., Forman W., 2005, MNRAS, 359, 1041
 Ruszkowski M., Begelman M. C., 2002, ApJ, 581, 223
 Ruszkowski M., Brüggen M., Begelman M. C., 2004a, ApJ, 611, 158
 Ruszkowski M., Brüggen M., Begelman M. C., 2004b, ApJ, 615, 675
 Ruszkowski M., Enßlin T. A., Brüggen M., Heinz S., Pfrommer C., 2007, MNRAS, 378, 662
 Ruszkowski M., Oh S. P., 2010, ApJ, 713, 1332
 Ruzmaikin A., Sokolov D., Shukurov A., 1989, MNRAS, 241, 1
 Ryu D., Kang H., Cho J., Das S., 2008, Science, 320, 909
 Sakelliou I., Acreman D. M., Hardcastle M. J., Merrifield M. R., Ponman T. J., Stevens I. R., 2005, MNRAS, 360, 1069
 Sanders J. S., Fabian A. C., Smith R. K., Peterson J. R., 2009, ArXiv eprints, 0911.0763
 Scannapieco E., Brüggen M., 2008, ApJ, 686, 927
 Schekochihin A. A., Cowley S. C., 2007, Turbulence and Magnetic Fields in Astrophysical Plasmas. Springer
 Schuecker P., Finoguenov A., Miniati F., Böhringer H., Briel U. G., 2004, A&A, 426, 387
 Sharma P., Chandran B. D. G., Quataert E., Parrish I. J., 2009a, ApJ, 699, 348
 Sharma P., Chandran B. D. G., Quataert E., Parrish I. J., 2009b, ArXiv eprints
 Sharma P., Colella P., Martin D. F., 2009, ArXiv eprints
 Sharma P., Hammett G. W., 2007, Journal of Computational Physics, 227, 123
 Shin M., Strauss M. A., Oguri M., Inada N., Falco E. E., Broadhurst T., Gunn J. E., 2008, AJ, 136, 44
 Springel V., Di Matteo T., Hernquist L., 2005, MNRAS, 361, 776
 Subramanian K., Shukurov A., Haugen N. E. L., 2006, MNRAS, 366, 1437
 Vazza F., Brunetti G., Gheller C., Brunino R., 2010, New Astronomy, 15, 695
 Vazza F., Brunetti G., Kritsuk A., Wagner R., Gheller C., Norman M., 2009, A&A, 504, 33
 Vogt C., Enßlin T. A., 2005, A&A, 434, 67
 Voigt L. M., Fabian A. C., 2004, MNRAS, 347, 1130
 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
 ZuHone J. A., Markevitch M., Johnson R. E., 2010, ApJ, 717, 908