SIMULATING THE COOLING FLOW OF COOL-CORE CLUSTERS

Our simulations start with an isolated idealized galaxy cluster roughly in hydrostatic equilibrium. We adopt the Perseus Cluster, a massive relaxed CC cluster that has been well observed, as our model template. Initially, the gas is taken to be spherically symmetric and the azimuthally averaged temperature profile is taken from analytic fits to the observations of the Perseus Cluster (Churazov et al. 2004):

$$\begin{equation} T = 7 \; \frac{1 + (r/71)^3}{2.3 + (r/71)^3} \; \rm {keV}, \end{equation} \tag{ 1 }$$

where r is the radius to the center of the cluster in kpc. This temperature profile is flat for r > 0.5 Mpc, while Perseus, like many clusters, has a declining profile at large radius (e.g., Vikhlinin et al. 2006). Since we are interested in the evolution in the center (r < 0.2 Mpc), we do not expect this difference to be significant. The electron density profile (Mathews et al. 2006) is taken to be

$$\begin{equation} n_e(r) = \frac{0.0192}{1 + \left(\frac{r}{18}\right)^3} + \frac{0.046}{\left[1 + \left(\frac{r}{57}\right)^2\right]^{1.8}} + \frac{0.0048}{\left[1 + \left(\frac{r}{200}\right)^2\right]^{1.1}}\ \rm {cm}^{-3}. \end{equation} \tag{ 2 }$$

The power-law index of the last term is slightly steepened compared to Mathews et al. (2006) so that at small radii (up to a few hundred kpc) the gas density profile does not significantly deviate from that obtained from the X-ray observation, but at large radii, where the fit is poorly constrained, the density drops as r^−2.2, which is more consistent with cosmological simulations and fits better with the NFW dark matter profile (Navarro et al. 1996). We compute the initial pressure from the density and temperature assuming the ideal gas law with γ = 5/3.

The gravitational potential is the combination of the self-gravity of the ICM (which does not dominate but is included for completeness) and three static components: the dark matter halo, the stellar mass of the brightest cluster galaxy (BCG), and the SMBH in its center. The dark matter halo follows an NFW distribution that matches the observed gas density and temperature of Perseus:

$$\begin{equation} \rho ^{{\rm NFW}} (r) = \frac{\rho ^{{\rm NFW}}_0}{\big(\frac{r}{r_s}\big)\big(1 + \frac{r}{r_s}\big)^2} \;\;, \end{equation} \tag{ 3 }$$

where ρ^NFW₀ = 7.5 × 10¹⁴ M_☉ Mpc⁻³, and r_s = 0.494 Mpc is the scale radius. Since the gas density and temperature are only observed up to a few hundred kpc and are uncertain at large radii, the NFW parameters can vary depending on the outer radius we choose when fitting the model. We experimented with slightly different NFW parameters and the results were not affected.

The stellar mass profile of the BCG is

$$\begin{equation} M_*(r) = \frac{r^2}{G} \left[ \left({r^{0.5975} \over 3.206 \times 10^{-7}}\right)^s + \left({r^{1.849} \over 1.861\times 10^{-6}}\right)^s \right]^{-1/s} \end{equation} \tag{ 4 }$$

in cgs units with s = 0.9 and r in kpc (Mathews et al. 2006). The SMBH is treated as a point mass at the very center of the cluster with M_SMBH = 3.4 × 10⁸ M_☉ (Wilman et al. 2005).

The self-gravity of the ICM is computed at each step, but it does not have much impact in our major runs because the total mass of the cool gas in the cluster center is small compared to the BCG and SMBH. The mass of the accumulated gas in the center of the core only becomes significant at late time in our low-resolution simulations where the cooling flow is stronger (see Section 4.5 for a discussion). The gravitational acceleration from different components is plotted in Figure 1. As can be seen from that figure, the radius of influence of the SMBH (where its acceleration dominates) occurs at approximately 50 pc.

Zoom In Zoom Out Reset image size

In order to give the gas some initial angular momentum, we set the cluster initially rotating slowly around the z-axis in our simulation. The angular momentum of a galaxy cluster can be characterized by its spin parameter λ, which was first introduced by Peebles (1971) and later interpreted also as λ = ω/ω_sup (Gottlöber & Yepes 2007), where ω_sup is the angular velocity that would provide rotational support to the system. As a simplification, we use the same rotational velocity v at all radii and ω ∼ v/r_s. Thus, the spin parameter λ can be expressed as λ = v/v_c, where $v_c = \sqrt{{G M_s}/{r_s}}$ with M_s being the total mass inside the scale radius r_s. Typical values of λ for a cluster are around 0.05 (e.g., Bullock et al. 2001). For our model cluster, we take a constant rotational velocity of v = 50 km s⁻¹, which is consistent with cluster simulations.

We also give each cell an initial random velocity that could potentially seed small-scale instabilities. The random velocity obeys a Gaussian distribution with a standard deviation of 200 km s⁻¹ for each of the three components, a typical value for the turbulent motion of gas in galaxy clusters (Ruszkowski & Oh 2010).

To set up the initial cluster configuration, we employ an iterative technique: starting with the root grid, we use the equations above to set the density, temperature, and velocities. We then apply the refinement criteria discussed in the next section to obtain additional levels of refinement, applying the initial conditions to each level and reapplying the initial conditions. The result is a grid hierarchy which is self-consistently refined and contains high-resolution initial data.

We refine a cell whenever any of the following three criteria are met.

• 1.  
  The gas mass criterion—a cell is refined if the gas mass in any cell exceeds 0.2 times the gas mass in one cell of the root grid, or more precisely, when

  $$\begin{equation} m_{{\rm cell}} > 0.2 \;\rho _{\rm mean} \left(\frac{L}{N_{\rm root}}\right)^3 \times 2^{\alpha l}, \end{equation} \tag{ 5 }$$

  where ρ_mean is the mean density of the universe at the initial time, L = 16 Mpc is the comoving box size, l is the level of refinement, and N_root is the number of cells on the root grid in each dimension. We use N_root = 256 for our standard simulation. This defines a critical gas mass that is always refined, which we then modify (with the factor 2^αl) depending on the level of refinement. If α = 0, the refinement is Lagrangian and the limit for m_cell does not change with refinement level l; in this case, a given mass clump will be resolved by the same number of cells as it collapses to smaller sizes. A negative value of α makes the refinement "super-Lagrangian" and the maximum value of m_cell decreases with l. The decrease is faster as α becomes more negative. In our simulation, we find that choosing a more negative α can better resolve the early evolution of the cooling flow (see Section 4.5 for a discussion). We use α = −1.2 for our standard run, which reduces the maximum mass of a cell by a factor of 3.8 × 10⁻⁶ going from l = 0 to l = 15. In the standard run, the level of refinement in the initial conditions is l = 8, which gives a maximum cell mass (i.e., roughly the mass resolution in the initial conditions) of about 8 × 10⁵ M_☉.
• 2.  
  The cooling criterion—we also apply refinement whenever the ratio of cooling time to sound-crossing time over the cell becomes too small (i.e., t_cool/t_cross(Δx) < β). The isobaric cooling time for an optically thin plasma is given by

  $$\begin{equation} t_{\rm cool} = \frac{\frac{5}{2} n k_b T}{n^2 \Lambda (T)}, \end{equation} \tag{ 6 }$$

  where k_b is the Boltzmann's constant, n = ρ/μm_H is the particle number density, and Λ(T) is the cooling function (Schure et al. 2009). The sound-crossing time over a cell is t_cross(Δx) = Δx/c_s, where Δx is the cell size and $c_s=\sqrt{\gamma P/\rho }$ is the sound speed. We use this criterion because when the cooling time is shorter than the sound-crossing time, the gas drops out of pressure equilibrium, and the cooling becomes isochoric. We use the limit β = 6, which is a somewhat arbitrarily chosen value larger than 1 so that we can fully resolve this transition. We have experimented with higher value of β = 12 and β = 20 for this ratio and the results do not change.
• 3.  
  The Jeans length criterion—a cell is refined whenever its size is larger than 1/4 of the Jeans length, to ensure that we properly resolve any gravitational instabilities (Truelove et al. 1997).In our simulations, the baryon mass criterion is important at early stages, well before the cooling catastrophe occurs, while the cooling criterion becomes important later as the gas density grows in the center, and the radiative cooling becomes more significant. The Jeans criterion is never important for our standard runs.

All simulations in this paper include radiative cooling. The radiative cooling function we use is adapted from the SPEX package (Schure et al. 2009), assuming equilibrium cooling. Non-equilibrium effects will play a role in cooler gas, but as this only occurs in the very center of our simulations, we do not expect this to play an important role. We adopt a metallicity of one-half solar for our model cluster (Schmidt et al. 2002). The cooling is truncated at temperatures lower than 10⁴ K because we are not interested in following the evolution of gas below this temperature.

Our standard run, which will be the primary simulation analyzed in this paper, has N_root = 256, with a comoving box size of L = 16 Mpc, and a maximum refinement level of l_max = 15. The minimum physical size of a grid cell in this simulation is given by $\Delta x_{{\rm min}} = L / (N_{\rm root} 2^{l_{\rm max}})$, which is about 2.0 pc.

To test the dependence of the results on the resolution, we also performed simulations with N_root = 64, 128, and 256. We also carry out simulations with the same N_root but different values of α. These runs will be discussed in Section 4.5. See Table 1 for the parameters of all runs.

We also performed several runs with isotropic thermal conduction suppressed from the standard Spitzer value by different factors of f_cond in order to test its influence on the evolution of the cooling catastrophe (see Section 4.3 for details). Two different suppression values were tested (f_cond = 0.1 and f_cond = 0.3). We did not include star formation or any feedback mechanism in our simulations.

In order to facilitate analysis of the simulation, Enzo generates an output when a new refinement level is reached for the first time. This produces outputs throughout the collapsing phase up until the point when the maximum refinement level l_max is achieved. To better resolve the evolution after this, starting from a few Myr before l_max is reached, the simulation writes one output every 1/3 Myr and continues to do so afterward. We run the simulation until the estimated energy from AGN feedback is more than strong enough to offset cooling, which occurs approximately 10–20 Myr after reaching l_max.

For completeness, we note that these simulations were run using comoving coordinates, but over a sufficiently brief span (Δz < 0.2) that expansion did not play a role. We adopted cosmological parameters corresponding to H₀ = 50 km s⁻¹ Mpc⁻¹, Ω_λ = 0, and Ω_m = 1, but note that the cosmological parameters are irrelevant to the physics inside the self-gravitating cluster and thus have no significant influence on its evolution or on the initial condition (they just serve to set the internal code units). The simulation starts at z = 1, which is chosen to give enough time for the cluster to evolve.

We show in Figure 2 the evolution of the gas density of our standard simulation with N_root = 256. Throughout the cluster core (r ≲ 30 kpc) the density slowly increases and the temperature decreases over the first few hundred million years. The density in the center (r ≲ 100 pc) grows more compared to the outer part of the core region r ∼ 0.3–10 kpc, where the temperature shows a plateau. After about 300 Myr, which corresponds to the cooling time of the gas in the initial conditions, the evolution rapidly accelerates, and the density increases by several orders of magnitude within a million years. Note that the outputs are not evenly spaced in time (output times are marked and color-coded in the upper-right panel of Figure 2). This marks the onset of the cooling catastrophe and occurs in the central region of the cluster. We refer to the outer boundary of the region where this catastrophe happens as the "transition radius." This is different from the traditional cooling radius which is defined as the radius at which t_cool is equal to the Hubble time (and which occurs near 100 kpc). The huge jump of the central density from ∼10⁻²³ g cm⁻³ to ∼10⁻²¹ g cm⁻³ happens within only a short period of time, about 3 Myr, which also corresponds well with the drastic decrease in the central temperature seen in the lower-left panel of Figure 2. After the cooling catastrophe, the transition radius steadily grows in time, although points exterior to the transition radius evolve only very slowly.

Zoom In Zoom Out Reset image size

The pressure profile (Figure 2, lower-right panel) shows relatively little variation over time, at least in the outer parts. As more gas flows to the cluster center, the pressure slowly builds up, although the slope in the intermediate region (between 100 pc and 10 kpc) is quite shallow. As the cooling catastrophe occurs, the pressure in the very center rises quickly. A small region with a slightly positive pressure gradient forms later right at the inner edge of the transition radius (at r ∼ 50 pc), creating a slight pressure "hole" (or inversion). This pressure hole does act to drive some gas inflow; however, the pressure profile is remarkably constant at the transition radius, while the density and temperature change by orders of magnitude. Both the size and the depth of the pressure hole are sensitive to the resolution of the simulation, as we will discuss in Section 4.5, becoming less prominent as the resolution increases. After the cooling catastrophe, the pressure profile is nearly time-independent.

It is clear from Figure 2 that the cluster can be divided into three regimes in radius. In the outer region (r  >  20 kpc), where the NFW dark matter dominates the gravitational potential, the gas properties stay relatively constant during the simulation. Inside the core, but outside the transition radius (0.1 kpc ≲ r ≲ 20 kpc), where the gravitational potential is dominated by the BCG, the density and pressure increase slowly while the temperature profile is nearly flat. Finally, in the very center of the cluster (r ≲ 0.1 kpc), within the radius of influence of the SMBH where the gravity is dominated by the SMBH (∼50 pc; see Figure 1), the cooling catastrophe first occurs at small radius. Then the transition radius moves slowly outward with time, from ∼10 pc when the collapse first happens to about ∼100 pc, about 10 Myr later.

Figure 3 shows the evolution of the radial inflow velocity of the gas. Initially, it is zero, and in the absence of cooling remains that way except in the outer region (r > 100 kpc), where the cluster is not in perfect hydrostatic equilibrium.¹ As radiative cooling proceeds, a slow but steady inflow develops. This is slow (∼20 km s⁻¹) and nearly constant over a large range in radius, from 20 kpc to a few hundred pc. Note that this is not a steady-state cooling flow because a constant velocity implies an accretion rate that rises with radius—we discuss this in more detail below. As the cooling catastrophe takes hold (recall that the outputs are not evenly spaced in time), the inflow velocity in the central region rises up to a few hundred km s⁻¹ inside 100 pc, and then drops after an accretion disk forms.

Zoom In Zoom Out Reset image size

So far, we have focused on one-dimensional analysis of the cluster evolution. In fact, the evolution is remarkably symmetrical, despite the fact that we seed the initial conditions with substantial small-scale velocity motions which quickly generate density and temperature perturbations. We show projections of the density, (X-ray-weighted) temperature, pressure, and X-ray luminosity of the central 16.6 kpc of the cluster in Figure 4. These are shown 296 Myr after the initial time, shortly (∼16 Myr) after the cooling catastrophe occurred. Clearly, the initial seeds have damped and the profiles are remarkably uniform. A close examination of the temperature image reveals faint spiral structures due to the slight rotation imprinted in the initial conditions. This demonstrates that the flow does not exhibit local thermal instabilities as predicted by linear perturbation theory (Malagoli et al. 1987; Balbus & Soker 1989).

Zoom In Zoom Out Reset image size

The primary instability (or catastrophe) that grows is a global one, as discussed above. This can be seen as the increase in density in the very center of the projected image (and the decrease in the temperature). The pressure map is quite smooth. The X-ray map in Figure 4 shows a small rise in the central emissivity in this region, but it can be appreciated that this would be difficult to observe. We return to the observability of this lower temperature gas in Section 4.2.

We stop the simulation about 20 Myr after the cooling catastrophe first occurs because the estimated energy output from AGN feedback rapidly exceeds the energy loss through radiation. We estimate that the feedback energy from the AGN is

$$\begin{equation} \dot{E}_{\rm AGN}=\epsilon \dot{M}_{400} c^2 , \end{equation} \tag{ 7 }$$

where $\dot{M}_{400}$ is the mass accretion rate measured at r = 400 pc, which is just outside the region dominated by the potential well of the SMBH, and is an efficiency parameter discussed in more detail below. We stop the simulation shortly after this feedback rate exceeds L_X50, which is the total calculated X-ray luminosity from the cooling gas inside r < 50 kpc.

The parameter is the efficiency relating this mass accretion rate to the feedback energy that heats the ICM, and is really the product of three somewhat poorly constrained efficiencies. The first is the fraction of $\dot{M}_{400}$, which actually accretes onto the black hole. Because we measure this value quite close to the SMBH-dominated region (at 400 pc), we argue that a significant fraction (at least 10%) of this mass will accrete onto the black hole. The second unknown efficiency is that associated with the conversion of accreted mass into energy, which we take to be approximately 10%. Finally, because the accretion rate is likely to be sub-Eddington, this feedback will be the so-called radio mode, and we assume that a large fraction of the feedback energy ends up heating the ICM. For the product of these three efficiencies, we take a conservative value of = 0.01, although we recognize that the final value is uncertain.

As shown in Figure 5, as the gas inflow increases, $\dot{E}_{{\rm AGN}}$ increases and exceeds L_X50 after the cooling catastrophe. At the last step, $\dot{E}_{{\rm AGN}}\sim 2.6 \times 10^{44} \;{\rm erg\ s^{-1}}$ and L_X50 ≈ 1.5 × 10⁴⁴ erg s⁻¹, which has increased only about 50% compared to the initial L_X50. We stop our simulations at that point because AGN feedback, which is not included, would be important and we only want to focus on the onset of the cooling flow. The total X-ray luminosity of the cluster is 6.8 × 10⁴⁴ erg s⁻¹ at the end of our simulation.

Zoom In Zoom Out Reset image size

To better understand the physics behind the cluster evolution, we examine several relevant timescales in Figure 6. The dynamical time of a gaseous system can be defined as

$$\begin{equation} t_{\rm dyn} (r) = \frac{\pi }{2} \frac{r^{3/2}}{\sqrt{G M_{{\rm dyn}}(r)}}, \end{equation} \tag{ 8 }$$

where M_dyn(r) is the enclosed mass of the system.

Zoom In Zoom Out Reset image size

The system undergoes gravitational collapse when t_dyn is shorter than the sound-crossing time

$$\begin{equation} t_{\rm sound} (r) = \frac{r}{c_s}, \end{equation} \tag{ 9 }$$

where r is the radius of the region. Figure 6 shows the cooling, dynamical, and sound-crossing timescales before (left panel) and after (right panel) the cooling catastrophe happens. Initially, the cluster is stable with t_cool  >  t_dyn  >  t_sound. As the gas condenses and cools nearly isobarically in the center, the cooling time t_cool decreases while t_sound slowly increases (both at fixed radius). When t_cool drops below t_dyn, which first occurs around r ≲ 50 pc, the gas undergoes a full-fledged cooling catastrophe. At the same time, t_sound exceeds t_dyn, making the core region gravitationally unstable and causing the gas to collapse to the cluster center under gravity. Meanwhile, the gas can no longer maintain pressure equilibrium since t_cool < t_sound. Also shown in Figure 6 are the thermal conduction and the SNIa heating timescales which will be discussed in Sections 4.3 and 4.4, respectively.

This analysis leaves open two questions. First, why does the gas temperature remain nearly constant in the intermediate region, between the transition radius (at a hundred pc) and the start of the CC (at r ∼ 20 kpc)? Second, what supports the gas inside the transition radius?

The first question arises because of the presence of the large temperature plateau noted earlier, which persists despite the fact that we evolve the system for many cooling times at these radii. Clearly, the cooling must be balanced by some form of heating, and since the gas motions are subsonic, the only real candidate is compression heating. We estimate the compressional heating timescale as

$$\begin{equation} t_{\rm compress} (r) = \frac{1}{(\gamma - 1)\nabla \cdot v}, \end{equation} \tag{ 10 }$$

where γ = 5/3 and $\nabla \cdot v = \frac{1}{r^2} \frac{\partial (r^2 v_r)}{\partial r}$. For simplicity, we only include the radial term, focusing on the role that the cooling flow plays in establishing the temperature plateau (we note that non-radial terms could decrease the heating time in the center, although manual inspection of the velocity indicates that the effect is minor). We plot in Figure 7 the ratio of t_compress over t_cool. At early times, over a large range of radius, t_compress ≈ t_cool. The inflow velocities that drive this heating are slight, a few tens of km s⁻¹ over much of this range, as seen in Figure 3, rising at small radius as the global cooling catastrophe sets in.

Zoom In Zoom Out Reset image size

This brings us to the second question, the evolution of the gas within the transition radius. As the cluster evolves, the cooling rate becomes higher in the center, requiring a larger inflow velocity to provide sufficient compressional heating. When this required velocity exceeds the sound speed, or if the gas becomes rotationally supported, compressional heating can no longer balance cooling and t_compress becomes larger than t_cool inside the transition radius. If the inflow velocity grew to the sound speed, as would occur for a purely radial evolution, the transition radius would be identified as a sonic point, and inside the gas would freely fall toward the black hole. We found this to occur in our lower resolution simulations (discussed in more detail below); however, in our best-resolved, standard simulation, we find that the cold gas inside the transition radius forms a rotationally supported disk. This is shown in Figure 8, which shows an estimate of the rotational velocity of the gas (computed by dividing the magnitude of the total specific angular momentum of the gas in a shell by the radius of that shell), compared to the Keplerian velocity. At late times, inside the transition radius, the gas becomes rotationally supported.

Zoom In Zoom Out Reset image size

In Figure 9, we show the projected density and density-weighted temperature in the central 330 pc for a slice of gas with a z-thickness of about 16.6 pc. This z-projection clearly shows the disk (x- and y-projections—not shown here—clearly demonstrate that this is a thin disk), which has a radius of about 50 pc. In fact, in this particular run, we find an inner disk and an outer polar ring, which shows some sign of forming denser fragments; typical densities in the disk are of order 10³ cm⁻³, but note that the disk is not well resolved and so the detailed disk structure shown here should be treated with caution.

Zoom In Zoom Out Reset image size

To better understand what determines the structure of the gas and why there are three regimes seen in the gas density, temperature, and pressure profiles in Figure 2, we carry out an approximate analytic analysis assuming hydrostatic equilibrium, which is valid before the cooling catastrophe happens (or at radii larger than the transition radius) when the inflow velocity v_r ≪ c_s.

In hydrostatic equilibrium, the gas pressure P(r) balances the gravitational acceleration g(r):

$$\begin{equation} \frac{{\it dP}(r)}{dr} = - \rho (r) g(r), \end{equation} \tag{ 11 }$$

where P = ρk_bT/μm_H and ρ(r) is the gas density. This by itself is not sufficient to determine the density and temperature profiles uniquely; however, if we also use the fact, as suggested by Figure 7, that compressional heating balances cooling, we can make more progress. We write this as

$$\begin{equation} \rho g v_r \approx n^2 \Lambda (T), \end{equation} \tag{ 12 }$$

where Λ(T) is the cooling rate and v_r is the radial velocity. We have assumed that the work done on a fluid element as it flows inward is balanced entirely by radiative cooling. Simplifying further, if we assume that Λ(T) ≈ Λ₀ is independent of temperature, as is appropriate for gas around a few keV, in the cool region (see Figure 2), and that the inflow velocity is nearly constant (see Figure 3), then we can obtain an expression for the gas density: n = μm_Hgv/Λ₀, and so ρ∝g.

We expect this result to hold primarily in the intermediate region of the flow, from the transition radius at about 100 pc to about 20 kpc. This region is dominated by the potential of the BCG, which is reasonably well fit over this range with a power law g(r) ∼ r^0.75. As predicted by the above argument, this is an excellent approximation to the density profile seen in Figure 2. If we combine this result with the equation of hydrostatic equilibrium, we predict that the temperature profile should be T(r)∝r^0.25, which is also in reasonable agreement with the nearly flat temperature profile we find in Figure 2 over this radial range.

To see how changing the form of the gravitational potential in the cluster core can affect the structure of the cooling gas, we performed a test simulation with only the NFW dark matter but without the BCG and the SMBH and compared it with the standard run with the same number of cells on the root grid (N_root  =  128). The cooling catastrophe starts inside r ∼ 1 kpc, much larger than the standard run. This is because the temperature decreases as T ∼ r toward the center throughout the core under the NFW gravitational potential and no temperature plateau forms.

One problem with the classic cooling flow model is that the expected amount of cool gas (i.e., gas with T ≲ T_vir/3, or T ≲ 2–3 keV for massive clusters) is much higher than that observed in X-ray observations of CC clusters. In our simulations, due to the small size of the region where the cooling catastrophe actually happens, the amount of cool gas is small and the X-ray luminosity does not change much in the core (see Figure 5), which would explain the lack of cool gas observed in the X-ray. This can also be seen from the simulated X-ray luminosity (Figure 4). Note that feedback is still required since this will not last without energy input: the amount of cool gas will grow with time and approaches the classic cooling flow prediction if we let the simulation run longer (see Section 4.7).

Since the observed temperature profile of galaxy clusters is usually obtained by fitting the X-ray spectra assuming a single temperature for the gas, it is hard for us to directly compare our temperature profile to those observed (since along a given line of sight, there exists gas with a range of temperatures). However, Fabian et al. (2006) fitted the X-ray spectra of a long Chandra observation of Perseus with a multi-temperature model, allowing them to compute the amount of mass in each temperature range. We compare the mass distribution from our simulation with their results in Figure 10, as well as with the prediction from the classic cooling flow model. We find that the steep drop-off in mass between 1 and 2 keV in our simulation is in good agreement with the observations and is much lower than the classic cooling flow prediction.

Zoom In Zoom Out Reset image size

Note that the recovery at ∼0.5 keV in the observed mass distribution is due to the observed filaments (this also probably contributes to the 1 keV bin; see Fabian et al. 2006), which may be caused by the interaction between the cooling gas and AGN feedback that is not included in our simulations. The heating from AGN can potentially result in an increased amount of gas at ∼0.5 keV, and we will study the impact of AGN heating in future work. Additionally, there is an excess at around 3 keV in our simulations and a lack of gas at around 5 keV compared to observations. This is likely due to the initial condition we choose in our simulations: we start from a CC configuration that matches Perseus' current configuration. We suspect that if we started from slightly different conditions—say an initial temperature plateau that was marginally higher, the fit could be improved. In any case, we are not trying to fit Perseus in detail, but simply to make the point that the simulation naturally produces a dearth of low-temperature (T ≲ 2 keV) gas.

Thermal conduction can potentially delay or even suppress cooling instabilities. However, the existence of a magnetic field can decrease the efficiency of thermal conduction below the classical Spitzer value κ_S (Spitzer 1956). This suppression can be parameterized by a factor of f_cond ranging from 0.1 to 1 depending on the orientation and topology of the field (e.g., Zakamska & Narayan 2003). The characteristic timescale associated with thermal conduction can be estimated as

$$\begin{equation} t_{{\rm cond}} (r) = \frac{\gamma }{\gamma -1}\frac{Pr^2}{f_{\rm cond} \kappa _{\rm {S}} T} , \end{equation} \tag{ 13 }$$

where κ_S = 1.84 × 10⁻⁵ T^5/2/ln λ erg s⁻¹ cm⁻¹ K^−7/2 with the usual Coulomb logarithm ln λ ≈ 37.

Thermal conduction is more efficient along the magnetic fields than perpendicular to them. In the center of CC galaxy clusters, where the temperature is increasing with radius, gas with even a sub-dominant magnetic field is susceptible to the heat-flux-driven buoyancy instability, or HBI (Quataert 2008; Parrish & Quataert 2008). This instability, in the saturated phase, orients the magnetic field to be perpendicular to the temperature gradient, and f_cond can be reduced to ≲ 0.1 in the radial direction (Parrish et al. 2009). On the other hand, the turbulent motion of gas can result in a more entangled magnetic field, restoring conduction (Ruszkowski & Oh 2010). For simplicity, we assume isotropic conduction in our simulations, but we perform two additional simulations with the suppression factor f_cond = 0.1 and f_cond = 0.3, which can be thought of as corresponding to the HBI-dominant case and the turbulent case. To save computational time, these simulations are done using 128³ root cells (N_root = 128).

We find that thermal conduction does not significantly affect the overall evolution of the cluster, except that it slightly changes the temperature and density profiles when the cooling catastrophe occurs, and it can somewhat change the timing of the cooling catastrophe itself by a few tens of Myr. This result is consistent with what one would expect from Figure 6: the thermal conduction timescale t_cond is shorter than the cooling timescale at r ≳ 100 kpc, where conduction can stabilize the local cooling instability, but at r ≲ a few tens of kpc where the cooling catastrophe first starts to develop, t_cond is longer than t_cool.

To show the influence of thermal conduction near r ∼ 100 kpc, we plot in Figure 11 the temperature profiles from simulations with different suppression factors. When f_cond is larger, the thermal conduction is more efficient, and thus the gas temperature is higher at r ≲ 100 kpc but slightly lower at r ≳ 100 kpc because thermal energy is transported inward from r ≳ 100 kpc where the temperature peaks. Inside the core region, the temperature evolution looks very similar with different f_cond. Note that the outputs in Figure 11 are shown 2 Myr after the cooling catastrophe occurs in order to make a fair comparison in the central regions. While the f_cond = 0 and f_cond = 0.1 runs collapse at almost exactly the same time (t = 295 Myr), the high conduction run actually collapses slightly earlier (about 20 Myr) because it damps the initial temperature perturbations caused by the relaxation at the beginning of the simulation which, in turn, are due to the slight deviation from hydrostatic equilibrium mentioned earlier.

Zoom In Zoom Out Reset image size

The cooling catastrophe in these simulations takes place very close to the center of the BCG, where heating from SNeIa produced by stars in the BCG may have an effect on the formation of cooling flows (e.g., Domainko et al. 2004) and can be important especially in the central regions (e.g., Conroy & Ostriker 2008). We do not include SN heating in our simulations but here we try to analytically estimate their contribution. SNeII can be ignored in the BCG due to its relatively old stellar population. The SNIa rate is ∼0.1 per 10¹⁰ M_☉ per century in galaxy clusters (Sharon et al. 2007). Assuming each SNIa releases 10⁵¹ erg thermal energy into heating up the ICM, the heating rate then is expressed as M_stellar (in M_☉) × 10³⁸ erg yr⁻¹ M⁻¹_☉ and the SNIa heating timescale can be estimated as

$$\begin{equation} t_{\rm SNIa} = \frac{\rho _{\rm gas}\times \frac{3}{2}\frac{k_b T}{\mu m_H}}{\rho _{\rm stellar} \times 10^{38}\ {\rm erg\ yr}^{-1}\ M_{\odot }^{-1}} , \end{equation} \tag{ 14 }$$

where ρ_gas and ρ_stellar are the local gas and the stellar densities, respectively.

This timescale is plotted in Figure 6, which shows that although SNIa can be an important source of heating at r ≲ 1 kpc and might slightly delay the formation of the cooling flow, it would not prevent the cooling catastrophe because t_cool is always shorter than t_SNIa at all radii.

To test the influence of resolution on our results, we perform simulations with N_root = 64 and 128 (keeping α = −1.2) to compare with our standard N_root = 256 run. Figure 12 shows the density and temperature profiles of the cooling gas at a similar stage of the evolution (∼2 Myr after the cooling catastrophe starts). The overall behavior of the solution is seen to be relatively independent of resolution; however, there is a systematic trend for the temperature plateau in the intermediate region (between 100 pc and 20 kpc) to be even flatter as the resolution increases, while the transition radius shrinks. We do keep the maximum refinement level constant (l_max = 15), so the highest resolution achieved increases by a factor of two in each case; however, even for the lowest resolution run, this corresponds to 7 pc, considerably smaller than the ∼100 pc transition radius. Therefore, we argue that it is not the maximum resolution that is the key, but the resolution achieved in the early stages of the collapse.

Zoom In Zoom Out Reset image size

Since the point where the cooling catastrophe happens (the transition radius) is sensitive to the temperature and density profiles at earlier stages, we want to resolve well the early evolution of the cooling gas in the center of the cluster. Changing α in the baryon mass refinement criterion (see Section 2.2) can affect the resolution, especially in the early stage inside the cluster core, and has a significant impact on the final results. The two bottom panels of Figure 12 show a comparison of runs with the same number of root cells N_root = 128, but different values for α. Again, the results are shown at a point roughly 2 Myr after the cooling catastrophe starts. At fixed radius, the gas density is lower and the temperature is higher with a more negative α, and the transition radius is smaller. Noticeably, with α = −0.1, the initial cooling region has a size of ∼1 kpc, more than an order of magnitude larger than that in the run with α = −1.2. This is despite the fact that each calculation has the same maximum resolution.

We also note that there is a degeneracy between α and N_root: a higher N_root has an effect similar to a more negative α. This is because they both result in better mass resolution (see Equation (5)) in the center of the cluster as the cooling catastrophe develops. If we only increase the maximum refinement level l_max without changing N_root or α, the result does not change despite a smaller cell size at the highest refinement level. Thus, we argue that it is crucial to have sufficiently high resolution at the early stage of the gas evolution.

Although we have not been able to achieve complete convergence in these calculations, it is clear that higher resolution tends to produce flatter temperature profiles and smaller transition radii.

Another difference between simulations with different resolution is that in low-resolution runs (e.g., with N_root = 64 and α = −0.2), the pressure drops dramatically inside r ≲ 1 kpc when the central gas density and temperature show a sudden change, forming a pressure hole which then grows deeper and larger with time. The pressure hole is deeper when the resolution is lower. The gas inflow velocity is larger than that in the high-resolution runs and becomes supersonic at r ∼ 1 kpc where the pressure gradient is the steepest, forming a sonic point and leading to a cooling catastrophe. Cold gas does not become rotationally supported as in the high-resolution runs but fragments inside the pressure hole via local cooling instability.

Finally, we examine the impact of different methods for solving the hydrodynamics equations. As noted earlier, we use the Zeus method for our simulations because of its fast performance and robust treatment of cold regions. To test the accuracy of its results, we also carried out test runs using the PPM, which is third-order accurate with its high-order spatial interpolation. We used N_root = 64, α = −1.2 and found results which were in good agreement with the Zeus solver.

In this section, we briefly examine how our results change as we vary details in the initial conditions. We start our simulations assuming that the gas is in hydrostatic equilibrium. However, the NFW dark matter potential we use is fitted to the observed gas properties in the inner few hundred kpc region, excluding the central ∼10 kpc where the BCG starts to dominate the potential, and thus the resulting NFW parameters can vary slightly depending on the boundary radius we choose when fitting the parameters. Therefore, the gas is not in perfect equilibrium when the simulation starts, especially in the outskirts and inside ≲ 10 kpc. Examination of the early evolution of the cluster shows that the temperature and density of the gas in the central kpc increases by about 30% and 50%, respectively, after about 20 Myr (∼t_cross at a few kpc). They reach these values quickly, and then show little evolution until cooling becomes important on a few hundred Myr timescale. In order to check if the initial transients are important for the cooling, we also carried out a simulation where we let the system settle down for ∼600 Myr (∼t_dyn at r ∼ 100 kpc) without turning on cooling, so that the gas in the region of interest is in hydrostatic equilibrium when cooling starts. The later evolution of the cluster is the same as the run that starts with the cooling on. This indicates that the transients are not playing a role in cooling; however, it does mean that the true "initial state" of the core region (i.e., the density and temperature profiles after the transients die out) is slightly different from our given initial profiles. This is inevitable in the sense that the initial (observed) profiles are not in hydrostatic equilibrium. To test the general robustness of our results, we have also experimented with slightly different sets of NFW parameters for the dark matter and find that the gas properties at the outskirts are slightly different with different NFW parameters, but the cooling flow evolution in the core region is not affected.

In many galaxy clusters, there is an offset between the X-ray emission center and the BCG, although the offset tends to be smaller in CC clusters (e.g., Sanderson et al. 2009a). To see if the offset would significantly change the results, we have performed one simulation with an initial offset of 20 kpc between the center of the gas and the gravitational potential. We find that the cluster gas settles down and re-centers on the BCG before the cooling catastrophe happens. This is consistent with the fact that t_cool > t_dyn initially, and the cluster relaxes before cooling starts, and therefore the results do not differ from the simulations without the offset.

To test the effect of other changes in our initial conditions, we performed one simulation without the initial random velocities. We find that changing the initial random velocity does not have a significant impact on the evolution of the CC because the initial random velocity is damped before the cooling catastrophe happens.

Since velocity perturbations do not directly perturb gas entropy, to further confirm that small-scale perturbations do not grow outside the transition radius in our simulation, we performed a run with initial density perturbations instead of velocity perturbations. To do this, we multiplied the density in each cell in the initial conditions by a Gaussian factor with a mean of unity and a standard deviation of 10%. Again, we do not see the growth of any local instabilities. This is in agreement with Joung et al. (2011), who found that perturbations in a hydrostatic atmosphere did not cool unless the perturbation was sufficiently nonlinear that the cooling time in the perturbation dropped below the time for the clump to accelerate to the local sound speed (roughly the dynamical time).

We also performed a simulation without initial rotation. The gas in the very center in this run still eventually becomes rotationally supported because the random initial velocities eventually are amplified due to the conservation of angular momentum and a gas disk forms. In fact, even in our standard run, a smaller disk along the x-axis forms inside the major disk along z-axis, which can be seen in Figure 9. The size of the disk in the run without initial rotation is smaller at early times, and when the cooling catastrophe first occurs in that run, the gas in the very center has yet to become rotationally supported; however, the required inflow velocity to balance cooling has already exceeded the sound speed, and so the flow passes through a sonic point (in the run with initial rotation, rotational support occurs before a sonic point develops).

As noted earlier, the gravitational potential does play an important role in the cooling catastrophe, and runs without a BCG produce significantly different results (see Section 3.2 for more details).

Finally, we carry out one simulation with a non-CC configuration, where we use the same NFW dark matter profile for Perseus but set the initial temperature to be isothermal and compute the initial gas density assuming hydrostatic equilibrium. The initial t_cool in the center is about 2 Gyr. Our simulation shows that after about 2 Gyr, the cooling starts to run away and a cooling flow develops in a way quite similar to what we see in the simulations with initial CC configurations. The temperature plateau is slightly different, which we think has to do with the difference in the initial gas to dark matter ratio. We will return to this point in a later paper.

In this section, we try to place these results in context, first making a link to steady-state cooling flow solutions, and then comparing to other (primarily simulation) work which looked specifically at only the developing cooling flow, and did not include feedback.

The classic cooling flow model (Fabian 1994) predicts a "cooling flow" of 100 s M_☉ yr⁻¹ for rich clusters, assuming that in a steady state, without other heating sources, the gas flows inward at a constant rate to replace the central gas that has cooled down and formed stars. The mass dropout occurs over the central cooling flow region, and was originally assumed to cool and condense out in small clumps via a local cooling instability. This picture has been known to be in disagreement with observations, which usually indicate a star formation rate at least an order of magnitude lower than predicted by the steady-state cooling flow model (Tamura et al. 2001; O'Dea et al. 2008; Rafferty et al. 2008). Our simulations show that a steady state is not reached before the AGN feedback is potentially strong enough to balance cooling, and therefore we argue that this solution is not relevant. However, it is interesting to see if we recover the steady-state result if we run the simulation long enough. In Figure 13, we plot the gas inflow for a simulation with N_root = 64 that runs much longer than our major runs. We find that after less than a hundred Myr, without any heating mechanism, the system approaches a steady state with a roughly constant mass flow of $\dot{M} \sim 300$ M_☉ yr⁻¹ in the cluster core (r < 100 kpc), consistent with the classic cooling flow prediction.

Zoom In Zoom Out Reset image size

Most previous simulation work on CC clusters has focused on the heating process, especially AGN feedback, but they do usually include a pure cooling flow simulation in which heating is turned off (e.g., Croton et al. 2006). Our results are consistent with these results inasmuch as there is overlap. For example, Brighenti & Mathews (2006) examined two-dimensional models and also found that a cooling-only model results in a relatively flat temperature profile (falling only by a factor of 2–3 over a range of 100 in radius). However, these simulations did not have the resolution (∼1 kpc) to follow the cooling catastrophe in detail, and instead used a parameterized mass dropout term in the mass-conservation equation. Similar results were found for a one-dimensional cooling flow in Mathews & Brighenti (2003).

Our results are also consistent with previous theoretical work; for example, Bertschinger (1989) presented one-dimensional steady accretion models normalized to self-similar cooling waves that demonstrated slowly declining—or even increasing—temperature profiles. They even found a sonic radius (similar in radius to our transition radius) that occurred very close to the SMBH.

Simulations that do not include a specific dropout term and did not have high resolution resulted in a much larger cooling region (a larger transition radius in our terminology) and a larger cooling flow immediately after the cooling catastrophe than found in our simulations.

Kaiser & Binney (2003) provide a semianalytic model of the time evolution of the cooling flow assuming hydrostatic equilibrium at every step of the evolution. Their density and temperature evolution agrees well with what we find for r > 1 kpc, when the inflow is subsonic and the gas is roughly in hydrostatic equilibrium. However, a semianalytic model cannot describe what happens inside r ≲ 100 pc where cooling runs away.

Guo et al. (2008) analyze the local and global instabilities in CC clusters and show that AGN feedback can suppress global radial instabilities. An interesting set of recent simulations (McCourt et al. 2011; Sharma et al. 2011) has addressed the role of local cooling instabilities in simulations of a CC, suppressing the cooling catastrophe by instituting a global heating term. They showed that local cooling instabilities can form as long as the timescale for a cooling instability is less than 10 times the dynamical time. This result does not conflict with our finding in this paper that local cooling instabilities do not occur before the global instability (or catastrophe) because of the global heating term instituted in those papers.

Finally, we note that we do not include feedback and so do not compare to simulations that try to explicitly model AGN energy injection (e.g., Omma et al. 2004; Brüggen & Scannapieco 2009; Dubois et al. 2010; Fabjan et al. 2010), although we will examine this in future work. Also, we have focused on the evolution of CC clusters rather than the formation of CC clusters, and so it remains possible that isotropic thermal conduction can play a role in the earlier evolution of non-CC clusters (Voit 2011).