CAUSE AND EFFECT OF FEEDBACK: MULTIPHASE GAS IN CLUSTER CORES HEATED BY AGN JETS

Unopposed radiative cooling of the intracluster medium (ICM) would lead to the cooling catastrophe—an inflow of gas toward the center, due to the loss of pressure, at the rate of several hundred M_☉ yr⁻¹ (see Fabian 1994 for a review). The main consequences would be strongly peaked surface brightness maps, central temperatures falling below T_vir/3, and a huge accumulation of cold gas in the cluster core, well over 10¹² M_☉. The fact that none of these cooling flow features are observed (e.g., Peterson et al. 2001, 2003; Tamura et al. 2003; Peterson & Fabian 2006) points out that either the cooling proceeds non-radiatively or it is balanced by some kind of heating.

The most popular and promising heating mechanism is the active galactic nucleus (AGN) feedback. Radio and X-ray data provide overwhelming evidence for the interaction of jets/outflows from the supermassive black holes (SMBHs) with the surrounding ICM (Böhringer et al. 1993; Blanton et al. 2001; Finoguenov & Jones 2001; Jones et al. 2002; Rafferty et al. 2006; McNamara & Nulsen 2007; Finoguenov et al. 2008; Rafferty et al. 2008). The power supplied by the AGN is in the range of 10⁴³–10⁴⁶ erg s⁻¹, adequate to globally offset radiative cooling losses. Further support for the outflow model comes from blueshifted absorption lines in the core of several galaxies: in the optical band (Nesvadba et al. 2008, 2011), UV and X-ray bands (Tombesi et al. 2010; Crenshaw et al. 2003 for a review), and at 21 cm (Morganti et al. 2005, 2007). These lines originate at kiloparsec distances and correspond to velocities of up to several 10⁴ km s⁻¹ and mass outflow rates around a M_☉ yr⁻¹ or more. The geometry is still unclear, but radio images are suggestive of conical outflows with narrow opening angles and possibly consistent with precession.

While the global energy balance may be achieved with the help of the AGN feedback and/or other mechanisms (Ruszkowski & Begelman 2002; Brighenti & Mathews 2003; Guo et al. 2008; Conroy & Ostriker 2008; Ruszkowski & Oh 2011), the issues regarding the coupling of the mechanical AGN energy to the thermal ICM energy, along with the spatial distribution of heating, are not settled. The coupling may occur, for example, via turbulent mixing of AGN bubbles with the ICM (Scannapieco & Brüggen 2008), viscous dissipation of weak shocks (Ruszkowski et al. 2004a, 2004b), and cosmic ray heating (Guo & Oh 2008; Ruszkowski et al. 2008; Sharma et al. 2009). Other mechanisms may also help heat the ICM and redistribute the energy throughout the cool core: dynamical friction (El-Zant et al. 2004), thermal conduction (Narayan & Medvedev 2001; Bogdanović et al. 2009; Ruszkowski & Oh 2010, 2011; Parrish et al. 2010), or sloshing motions (Morsony et al. 2010; ZuHone et al. 2010).

For most microscopic heating mechanisms the ICM is expected to be locally thermally unstable. In other words, even if the ICM maintains global thermal balance without a cooling flow, some of the gas can experience runaway cooling. Indeed, the observations by Edge (2001) and Salomé et al. (2003) revealed that cooling flows are not totally stopped: the mass of cold gas, mainly dominated by molecular H₂ associated with CO emission, reaches 10⁹–3 × 10¹¹ M_☉ in the majority of cool-core clusters. Cooling may also be inferred from star formation (e.g., O'Dea et al. 2008, typically 1–10 M_☉ yr⁻¹). More recently, a wide survey conducted by McDonald et al. (2010, 2011a, 2011b) utilized Hα emission as a robust tracer of the cold phase. In more than 60% of clusters and groups, extended and/or nuclear cold gas was detected. The morphology of this phase is usually clumpy or filamentary and elongated toward the center (e.g., in the Perseus cluster; Salomé et al. 2006; Fabian et al. 2008).

Cavagnolo et al. (2008), Rafferty et al. (2008), Voit et al. (2008), and McDonald et al. (2011a) have shown the existence of a tight correlation between the amount of cold gas and the ICM entropy or the cooling time. Only clusters with central entropies below 30 keV cm² show evidence for the extended cold gas. This criterion can be expressed in terms of the ratio of the cooling time (t_cool⁷) and the free-fall time (t_ff).

Analytical estimates and Cartesian simulations by McCourt et al. (2011; hereafter M11) pointed out that a spatially extended and inhomogeneous cold phase can form whenever t_cool/t_ff ≲ 1. Previous studies (e.g., Balbus & Soker 1989) did not predict the presence of multiphase gas because they neglected a spatially distributed heating, which provides global thermal stability. Figure 12 in M11, based on observational data, suggests that the critical threshold is actually higher, ≈10. Sharma et al. (2011; hereafter S11) extended the work of M11, considering spherical feedback models. They found that, due to the geometrical compression, the criterion is in fact less stringent: instability occurs when t_cool/t_ff ≲ 10.

The results presented in M11 and S11 are very promising, but the distributed heating models adopted there are idealized. In these models the heating is either added using a prescribed spherically symmetric functional form, or it is constructed to be equal to the shell-averaged radiative cooling loss at every radius. Given the simplicity of these prescriptions, it is necessary to consider more realistic heating models to reliably study the formation of the multiphase medium, along with reproducing several important features, such as shocks, cavities, and turbulence. We perform high-resolution three-dimensional simulations which incorporate ICM heating via collimated AGN jets. This study is also a natural extension of our previous efforts (Gaspari et al. 2009, 2011a, 2011b), which included AGN jet feedback but mainly focused on the dynamics of the hot ICM phase, suppressing the formation of the multiphase gas via cold phase dropout.

Our more realistic AGN feedback simulations confirm that extended multiphase ICM forms only when the instantaneous ratio of the thermal instability and free-fall timescales falls below a critical threshold of ≈10, irrespective of the initial condition. The cluster core attains rough thermal equilibrium after a few cooling times. In addition to occasional extended cold filaments, centrally concentrated cold gas exists at most times. Thus, AGN feedback primarily occurs via the cold phase. Our simulations also show that sometimes the cold gas can be carried out by AGN jets. The evidence for such cold filaments dredged out from the central galaxy is seen in the Perseus cluster (Fabian et al. 2003).

The paper is organized as follows. In Section 2 we present the numerical setup, emphasizing our implementations of cooling and jet feedback. In Section 3 we present the results. Section 4 involves the discussion of the detailed thermal state of the core and Section 5 summarizes the main conclusions.

We use the adaptive mesh refinement code FLASH 4.0 (Fryxell et al. 2000), modified to simulate mechanical AGN feedback in the cooling ICM. The numerical finite-volume method employed is the third-order accurate split piecewise parabolic method, suited to solve the Euler equations including the required heating and cooling terms:

$$\begin{equation} \frac{\partial \rho }{\partial t} + {\bf \nabla }\cdot \left(\rho {\bf v}\right) = \pm {\bf \nabla }\cdot {(\dot{m}_{\rm jet}\,{\bf n}}_{\rm z}) -\dot{\rho }_{\rm bh}, \end{equation} \tag{ 1 }$$

$$\begin{equation} \frac{\partial \rho {\bf v}}{\partial t} + {\bf \nabla }\cdot \left(\rho {\bf v} \otimes {\bf v}\right) + {\bf \nabla }{P} = \rho {\bf g}\, \pm {\bf \nabla }\cdot ({\dot{m}_{\rm jet} \bf {n}_{\rm z}\otimes {\bf v}_{\rm jet}})-\dot{p}_{\rm bh}, \end{equation} \tag{ 2 }$$

$$\begin{eqnarray} \frac{\partial \rho \varepsilon }{\partial t} + {\bf \nabla }\cdot \left[\left(\rho \varepsilon + P\right) {\bf v}\right] &=& \rho {\bf v}\cdot {\bf g} -\: n_{\rm {e}} n_{\rm {i}} \Lambda (T,Z)-\dot{E}_{\rm bh}\nonumber\\ &&\pm {\bf \nabla }\cdot {\left(\frac{1}{2}\dot{m}_{\rm jet} v_{\rm jet}^2 \;{\bf n}_{\rm z}\right)}, \end{eqnarray} \tag{ 3 }$$

$$\begin{equation} P = \left(\gamma -1\right)\rho \left(\varepsilon -\frac{1}{2}v^2\right), \end{equation} \tag{ 4 }$$

where ρ is the gas density, $\bf {v}$ is the velocity, ε is the specific total energy (internal plus kinetic), P is the pressure, and γ = 5/3 is the adiabatic index. The temperature is computed from P and ρ using Equation (4), with an atomic weight of μ ≃ 0.62, appropriate for a totally ionized plasma with 25% He in mass. The externally imposed gravitational acceleration g accounts for the dark matter NFW halo (Navarro et al. 1996) characterized by a virial mass and radius of 10¹⁵ M_☉ and 2.6 Mpc, respectively (concentration parameter ≃ 6.6). Our implementations of cooling, jet flux and black hole sink terms in Equations (1)–(3) deserve careful discussion and are presented in Sections 2.1 and 2.2.

We decided to neglect magnetic fields and anisotropic thermal conduction. M11 showed that anisotropic conduction changes the morphology of the cold gas (more filamentary), but does not affect the criterion for the generation of the multiphase gas due to thermal instabilities. We defer the study of these effects to future work.

2.1. Radiative Cooling

The ICM emits radiation mainly in the X-ray band due to bremsstrahlung (T ≳ 10⁷ K) and due to line emission at lower temperatures (10⁴–10⁷ K). Assuming collisional ionization equilibrium, the cooling source term in Equation (3) can be modeled with a cooling function Λ(T, Z): we use an analytical fit to the Sutherland & Dopita (1993) tabulated values for a fully ionized plasma with metallicity Z = 0.3 Z_☉ (e.g., Tamura et al. 2004), setting the temperature floor at 10⁴ K. The total emissivity of the gas is also proportional to the electron and ion number densities: $\mathcal {L}=n_{\rm {e}} n_{\rm {i}} \Lambda (T,Z)$. For simplicity, we assume that molecular weights are constant, with μ_eμ_i = 1.464.

An important difference between the current implementation of the cooling term and that adopted in our previous models (Gaspari et al. 2011a, 2011b; hereafter G11a and G11b; Brighenti & Mathews 2002) is the absence of the dropout term, which removes the cold gas from the domain. Here we are interested in the formation and growth of thermal instabilities, producing cold clumps and filaments in the ICM.

The presence of the cooling term can significantly slow down the computation. Explicit schemes, such as Runge–Kutta, are usually simple and fast (per time step), but impose severe restrictions on the time step Δt, making it much smaller than the regular hydrodynamical time step. Implicit schemes allow us to use larger time steps, but these methods are less precise. In a splitting framework, i.e., first solving the hydrodynamic partial differential equation and then the cooling source term ordinary differential equation (ODE, assuming fixed density), no implicit or explicit solver is actually needed for the ODE step.

The new temperature T^{n + 1}_j of a cell j at time n + 1 can be computed from the old temperature Tⁿ_j by exactly evaluating the following integral:

$$\begin{equation} \int ^{T_j^{n+1}}_{T_j^n}{\frac{dT^{\prime }}{\Lambda (T^{\prime })}}=-\frac{(\gamma -1)\mu }{k_{\rm B}\mu _{\rm e}\mu _{\rm i}m_{p}}\rho _j\Delta t\:. \end{equation} \tag{ 5 }$$

A practical way to compute the integral is to split the integration interval into two parts: one from Tⁿ_j to T_ref and the other one from T_ref to T^{n + 1}_j, using an arbitrary reference value (see Townsend 2009). This leads to

$$\begin{equation} T^{n+1}_j=Y^{-1}\left[Y\big(T^n_j\big)+\frac{T^n_j}{T_{\rm ref}}\frac{\Lambda (T_{\rm ref})}{\Lambda \big(T^n_j\big)}\frac{\Delta t}{t_{\rm cool}} \right], \end{equation} \tag{ 6 }$$

where the cooling time and Y are defined as

$$\begin{equation} t_{\rm cool}\equiv \frac{e}{n_{e}n_{i}\Lambda \big(T_{j}^{n}\big)}=\frac{k_{\rm B} \mu _{\rm e} \mu _{\rm i}m_{\rm p}}{(\gamma -1)\mu }\frac{T^n_j}{\rho _j\Lambda \big(T_j^n\big)}, \end{equation} \tag{ 7 }$$

$$\begin{equation} Y(T)\equiv \frac{\Lambda (T_{\rm ref})}{T_{\rm ref}}\int ^{T_{\rm ref}}_T\frac{dT^{\prime }}{\Lambda (T^{\prime })}, \end{equation} \tag{ 8 }$$

where e ≡ P/(γ − 1) is the internal energy density. The temperature T^{n + 1}_j can be exactly computed if the reciprocal of the cooling function is analytically integrable. Thus, we decompose the cooling function into a large number of piecewise power laws in logarithmic temperature intervals

$$\begin{equation} \Lambda (T)=\Lambda _k\left(\frac{T}{T_k}\right)^{\alpha _k}\quad \quad \quad \quad \quad T_k\le T\le T_{k+1} \end{equation} \tag{ 9 }$$

and precompute at initialization the constants of integration for each interval. The integral in Equation (8) can be computed analytically using this power-law decomposition of the cooling function. The only effort in retrieving T^{n + 1}_j comes from finding the correct interval in which the argument of Y⁻¹ resides (Equation (6)).

Aside from being exact (in the ODE step), this method is faster than either explicit or implicit methods. The error from the piecewise power-law approximation is negligible, because the number of the temperature intervals can be made very large (matching the tabulated Sutherland & Dopita values), without large memory requirements.

In principle, our solver does not require a time step limiter for stability. However, the coupling between the hydrodynamics and source term integration is still approximate. Thus, we preferred to impose an upper limit on the cooling time step (∼t_cool). We note that this did not dramatically increase the computational requirements. In order to further improve the coupling, we have also implemented a second-order accurate temporal evolution in the splitting scheme, i.e., we alternate the order of the integration of the hydrodynamical (H) and source terms (S) as follows: H-S-S-H.

We have compared explicit calculations, setting a strong limiter (≲ 0.1 t_cool), with runs based on the "exact" solver. Tests were carried out in idealized conditions (e.g., an isolated cooling cylinder or sphere) and simulating a short pure cooling flow evolution. Overall, we obtained equivalent results, e.g., the cold mass and time of collapse.

2.2. AGN Heating

The heating module consists of two key elements: computation of mass accretion onto the central SMBH and the energy injection via collimated bipolar outflows. The quantity of accreted material is estimated as the inflow through a "spherical"⁸ surface very near the center, r ∼ 500 pc. Each cubical cell j containing a piece A_j of the spherical surface contributes $\dot{M}_{{\rm acc}_j}=\rho _j v_{r_j} A_j$ to the accretion rate only if the radial velocity $v_{r_j}$ is negative (inflow). The accreted gas mass, $\dot{M}_{{\rm acc}_j}\Delta t$, is then removed from the same zone⁹ (along with the related momentum and total energy; "bh" sink terms in Equations (1)–(3)).

After computing the total mass accretion rate, we calculate the kinetic power of the bipolar outflows from

$$\begin{equation} P_{\rm jet} = \epsilon \,\dot{M}_{\rm acc}\, c^2, \end{equation} \tag{ 10 }$$

where is the (mechanical) efficiency of the feedback, which is a free parameter of our model. We note that the actual accreted gas should be ideally calculated a few gravitational radii from the black hole. This would require an extreme dynamical range, which is currently not possible due to hardware limitations. However, we can simply absorb any uncertainty in the actual accretion rate into our efficiency parameter . If only a fraction f of the gas entering our sink region falls into the black hole, the real efficiency would need to be rescaled to /f.

The AGN mass, momentum, and energy are injected into the volume through the internal boundaries in the middle of the computational domain. The total cross-sectional area of the nozzle is A_noz = (2Δx)², where Δx is the highest resolution in the domain. Both jets originate at the same interface along the positive and negative z-direction (n_z is the z-vector). The jet flux terms appear on the right-hand side of Equations (1)–(3). This method of including AGN feedback appears to be slightly better than injecting quantities at zone centers, which may lead to artificial evacuation of the regions perpendicular to the jets.

The integrated mass outflow rate, $\dot{M}_{\rm jet}\equiv \dot{m}_{\rm jet} A_{\rm noz}$, is based on the conservation of momentum (M_jetv_jet), leading to

$$\begin{equation} \dot{M}_{\rm jet} = 2\,P_{\rm jet}/v^2_{\rm jet}, \end{equation} \tag{ 11 }$$

where v_jet is the subrelativistic outflow velocity fixed at 5 × 10⁴ km s⁻¹. The velocity of AGN outflows is still poorly constrained and observations suggest values from several 10³ to several 10⁴ km s⁻¹ (e.g., Crenshaw et al. 2003; Tombesi et al. 2010). Our choice of the outflow velocity gives reasonable values for the mass outflow rate of ∼M_☉ yr⁻¹. Note that more massive (and slower) outflows tend to carve a more visible channel in the ICM, while lighter (and faster) jets release most of their energy in the central region, producing strong shocks and big bubbles similar to X-ray data (see G11a,b).

We also implemented wobbling jets, with variable inclination (θ) and azimuthal angle (ϕ). In these computations the location of the nozzle remains fixed in space, but all three velocity components (e.g., v_x = v_jet sinθ cosϕ) vary with time. The injection is still "cylindrical," in the sense that the velocity vectors at the nozzle remain parallel. We consider simple models where jet orientation varies randomly within a specified solid angle. These models are further discussed in Section 3.2.3.

2.3. Initial Conditions

We model our initial conditions on A1795, which is a typical massive relaxed cool-core cluster characterized by the virial mass of ∼10¹⁵ M_☉. In order to compute the initial hydrostatic gas density profile, we take the observed temperature fit (Figure 1, second column, black dashed line) and use the analytical formula for an NFW dark matter gravitational acceleration. The normalization of the density profile is chosen such that the gas fraction is ≈0.15 at the virial radius. The obtained initial conditions correspond to the initial minimum ratio of t_cool/t_ff ∼ 7, where the free-fall time is given by

$$\begin{equation} t_{\rm ff} = \left(\frac{2r}{g(r)}\right)^{1/2}. \end{equation} \tag{ 12 }$$

Zoom In Zoom Out Reset image size

Our computational domain is 1.27 Mpc on a side and we use outflow boundary conditions. We adopt static mesh refinement, imposing the highest zoom at the center of the box/system.¹⁰ The linear size of each (cubical) shell doubles among adjacent levels. The highest refinement level is 10 and every refinement shell has a width of two FLASH blocks, i.e., 2 × 8 zones. Such choice of the refinement pattern does not affect our conclusions regarding cold clumps and other structures. Finally, the finest cell of the box has a size of Δx ∼ 300 pc.

As noted by S11, symmetric initial conditions do not produce thermal instabilities. Therefore, we perturb the initial density and temperature profiles to initiate the formation of cold filaments and blobs. Depending on the model parameters, such features can develop even in the early stages of evolution, before the AGN has become very powerful. We also note that the presence of the instabilities early on facilitates dispersal of the energy injected by the collimated AGN jets.

The fluctuations leading to the thermal instability can be produced by the turbulence associated with cosmological mergers. It is beyond the scope of this work to study such effects. It is thus sufficient to adopt a simple Gaussian random field with a white noise power spectrum of gas density. Using dispersion relation and direct simulation, M11 and S11 demonstrated that thermal instability does not depend on the range of the driving modes in Fourier space. Moreover, in the late stage of the evolution, the density fluctuations are dominated by the AGN jets.

We generate initial conditions using IDL. First, we initialize a complex stochastic field (1024³) in the k-space, with the amplitude of each mode given by

$$\begin{equation} W\propto k^{-1}, \end{equation} \tag{ 13 }$$

where k = (k²_x + k²_y + k²_z)^1/2. The above prescription ensures that the "3D" power spectrum, which is ∝W² k², remains constant, i.e., is described by a white noise spectrum. Second, we combine the above amplitude with a Gaussian random field (e.g., Ruszkowski et al. 2007):

$$\begin{equation} [{\rm Re}(W({\bf k})),\; {\rm Im}(W({\bf k}))] = [G(u_1)W,\;G(u_2)W], \end{equation} \tag{ 14 }$$

where G is a function of a uniform random deviate u₁ or u₂ that returns Gaussian distributed values. This method guarantees that the phases are random. We also impose cutoffs at k_max = 1024/4 and k_min = 3 (in units of 2π/L). The high-k cutoff is introduced in order to avoid spurious numerical effects. Third, we compute an inverse Fourier transform of W(k) to obtain W(x, y, z). We normalize it by the mean of its absolute value 〈|W(x, y, z)|〉. Finally, we superpose the normalized fluctuations onto the initial density distribution ρ_unp:

$$\begin{equation} \rho (x,y,z) = \rho _{\rm unp}(x,y,z)\left(1+\xi \frac{W(x,y,z)}{\langle |W(x,y,z)|\rangle }\right), \end{equation} \tag{ 15 }$$

where ξ = 0.3 is the amplitude of the perturbations. The new temperature distribution is computed from the new density assuming isobaric fluctuations.

We present simulation results for different AGN efficiencies , methods of energy injection (steady and wobbling jets), and different initial values of the minimum t_cool/t_ff in the ICM (TI-ratio; minimum value over all radii). Our previous results (G11a; S11) demonstrated that a suitable choice of the mechanical AGN efficiency may be in the range 5 × 10⁻³–10⁻². This is supported by observational constraints (e.g., Merloni & Heinz 2008). Common observed values for the minimum initial TI-ratio are in the range from 5 to 25 (M11). Motivated by these results, we perform computations with initial TI-ratio 7 or 21, and = 6 × 10⁻³ or 10⁻² (see Table 1).¹¹

Table 1. Model Parameters

Model	Efficiency	Min(tcool/tff)	Notes
Pure CF	...	7, 21	No AGN heating
r7-6em3	6 × 10−3	7	Along ±nz
r7-6em3w	6 × 10−3	7	Random wobbling
r7-6em3c	6 × 10−3	7	Convergence (2Δx)
r7-1em2	10−2	7	Along ±nz
r21-6em3	6 × 10−3	21	Along ±nz
r21-1em2	10−2	21	Along ±nz

Download table as:  ASCIITypeset image

3.1. Pure Cooling Flow

We first analyze a reference model, i.e., a simulation without AGN heating. In this simulation we keep the central gas sink term but neglect feedback. This case leads to the classic "cooling flow" catastrophe. The plasma in the central region of the cluster (∼100–200 kpc) loses energy due to radiative cooling. Consequently, the gas pressure drops, leading to a subsonic spherical inflow. The increase in the central gas density results in higher cooling rates, which, in turn, accelerates the mass accretion rate, leading to runaway cooling and accretion. The results from a cooling flow run for the initial TI-ratio of 7 are presented in Figures 1 and 2.

Zoom In Zoom Out Reset image size

In the first column of Figure 1 we show the total mass, the average mass rate, and the accretion rate as a function of time. The quantities in the first two panels correspond to the gas within the central 20 kpc. Different gas phases are coded according to "cold" (T < 5 × 10⁵ K, blue), "warm" (5 × 10⁵ ⩽ T < 10⁷ K, yellow), and "hot" (T ⩾ 10⁷ K, red). The choice is not arbitrary because the limits of these temperature ranges correspond approximately to different critical slopes in the cooling function Λ(T). Top panel reveals that the total mass of cold gas exceeds 10¹² M_☉ after a few Gyr. Middle panel shows that after 1 Gyr the cold phase continues to accumulate in the core (20 kpc) at a steady average rate of ∼500 M_☉ yr⁻¹ (peaks can reach ∼800 M_☉ yr⁻¹), a blatant discrepancy with observations (Peterson et al. 2001, 2003; Peterson & Fabian 2006). The cold gas never decreases with time, usually overwhelming the hot and warm phases by one and three orders of magnitudes, respectively. The accretion rate in the sink region is closely linked to the cooling rate. However, not all of the cooling gas gets accreted because the inflowing cold gas forms a rotationally supported structure.

The second column displays the radial profiles of the emission-weighted temperature, electron number density, and entropy (for the X-ray-emitting gas with T > 0.3 keV). Curves are plotted every 500 Myr, and the color changes gradually from dark blue to light violet; black dashed lines correspond to the initial conditions. The pure cooling flow manifests itself as a sudden decrease in temperature, lesser than T_vir/3. The projected emission-weighted (T_ew) temperatures falls in particular below the observational data of A1795 (circles, Tamura et al. 2001; triangles, Ettori et al. 2002). The electron number density n_e, linked to the surface brightness, increases rapidly within 70 kpc from the center, with the simulated profiles lying above the observational data points. The entropy (k_bT/n^2/3_e) near 1 kpc is ∼1 keV cm², without showing any sign of the typical observed floor (Cavagnolo et al. 2008).

In Figure 2, we present the maps of X-ray surface brightness, the midplane cut of electron number density (n_e), and the gas phase diagrams with temperature versus n_e (values are color coded by gas mass in each bin). Each map has logarithmic scale and contains the data only for the central core (20 kpc); the snapshot is taken at 1.51 Gyr. The density maps reveal significant cold gas accumulation in the core, accompanied by a strong increase in the central surface brightness. Since in the "warm" temperature range the gas has a short cooling time, the ICM tends to accumulate in the "cold" and "hot" phases. This effect is clearly visible in the phase diagram.

In the cooling flow run distinct cold gas clouds and filaments can be seen predominantly in the early stages of the evolution. At late times (e.g., 1.51 Gyr) the collapse proceeds in a monolithic way, and the atmosphere is globally unstable. Note that the t_cool/t_ff criterion does not apply to a cooling flow, because the theoretical analysis assumes thermal balance (see M11). In order to generate multiphase gas in a cooling flow, a substantial amplitude of perturbations is required, e.g., ≳ 0.1 (see also Pizzolato & Soker 2005; S11). Both initial TI-ratio of 7 and 21 (not shown) produce similar results. The latter case has longer saturation time and slightly smaller mass in the cold phase (the initial total mass in the box is also 8% smaller).

All the features of the pure CF model, discussed above, are in contrast with observations of galaxy clusters (Peterson et al. 2001, 2003; Ettori et al. 2002; Tamura et al. 2003; Peterson & Fabian 2006). This simulation merely serves as a reference case for the runs presented in the subsequent sections. A successful model will instead be characterized by a quasi-equilibrium between heating and cooling, with low cooling rates (≲ 10% of the pure cooing flow), profiles consistent with observations (e.g., emission-weighted radial profiles with reasonable gradients and characteristic length scales). At the same time, the model needs to account for the spatially extended (up to ∼15–20 kpc) distribution of cold blobs and filaments while not violating the constraints on the total mass in these features. Below we discuss models that meet these requirements.

3.2. AGN Feedback: t_cool/t_ff = 7

Before including AGN jet feedback, we verified that using the idealized prescription of S11, i.e., local heating rate given by the cooling rate averaged in spherical shells, produces thermal instabilities and spatially extended multiphase gas. The results were consistent with S11. However, unless cooling and heating are turned off in the central zones (as in M11), this simple prescription led to the massive accumulation of cold gas in the central few kiloparsecs, especially when the gas dropout or outflow through a central boundary was neglected. These findings provide additional motivation to include a realistic feedback mechanism. Below we discuss the results of our simulations that include self-regulated feedback due to collimated AGN jets.

3.2.1.  = 6 × 10⁻³

Mechanical feedback due to AGN jets dramatically changes the entire evolution of the galaxy cluster, compared to the pure cooling case (Section 3.1). In Figures 3 and 4, we present diagnostic plots and maps for the runs that include AGN jets with mechanical efficiency = 6 × 10⁻³.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The first column in Figure 3 shows the same quantities as those shown in the left panels in Figure 1. The gas mass in the cold phase M_c(t), shown in the top panel, reaches ∼10¹¹ M_☉ at the end of the evolution. This is consistent with molecular cold gas observations (Edge 2001; Salomé et al. 2003) and is reduced by over a factor of 10 compared to the pure cooling flow value. The middle panel shows that the average cooling rate $\langle \dot{M}_{\rm c}\rangle$ is also significantly quenched: after an initial transition phase, in which cooling is slightly elevated, the system starts to self-regulate through more powerful outbursts, and eventually reaches a steady rate of ∼20 M_☉ yr⁻¹. This value is about 4% of the value seen in the pure CF simulation, entirely consistent with observations (e.g., Peterson & Fabian 2006). We note that the instantaneous cooling rate (not shown) is quite variable, changing rapidly from positive to negative values due to the action of AGN heating.

As the cold gas mass M_c in the central 20 kpc exceeds the mass in the hot (three times) and warm phases (10³ times), the accretion rate is determined by the cold gas (last panel of first column); most of the time the black hole accretes cold gas at a rate between 1–5 M_☉ yr⁻¹. Cold accretion rate (blue) exceeds hot accretion rate (red) by at least two orders of magnitude. The latter is commonly associated with a Bondi-like mode; since the hot accretion rate is so weak, feedback simulations based on Bondi prescription usually require an artificial boost factor or high mechanical efficiencies ≳ 0.1, even if the Bondi radius is resolved (e.g., G11a; Di Matteo et al. 2008). However, we warn against interpreting our cold gas accretion rates literally, as the actual accretion rate onto the black hole may be lower (due to non-negligible angular momentum of the infalling gas and the possibility of significant black hole outflows, e.g., Blandford & Begelman 1999, or magnetically- and line-driven accretion disk winds).

The mass accretion rate determines the AGN power and mass outflow rate (see Equations (10) and (11)). The properties of one of the jets are shown in the second column in Figure 4. The cumulative energy injection as a function of time is presented in the top panel. At 5 Gyr the total injected energy is ∼10⁶² erg. Even if the entire mass of the gas that falls into the sink region is added to the SMBH, its total mass will still be within reasonable limits (less than several 10⁹ M_☉). Typical values of the jet power (second panel) and mass outflow rate are P_jet ∼ 3– 4 × 10⁴⁴ erg s⁻¹ and $\dot{M}_{\rm jet}\sim 0.5\hbox{--}1$ M_☉ yr⁻¹. The cumulative distribution function (CDF¹²; last panel) of the jet power suggests that the feedback is characterized by a "duty cycle." For example, assuming that the AGN is considered "active" only when the mechanical power exceeds 10⁴⁴ erg s⁻¹, the duty cycle would be about 70% (i.e., AGN active 70% of the time). This quantity may be poorly defined because of its sensitivity to the chosen threshold. In particular, deeper observations may reveal more signs of AGN activity, leading to higher duty cycles. Current observational estimates (e.g., Dunn et al. 2006) place AGN duty cycle at a level ≳ 70%.

The radial profiles (Figure 3, third column) of the emission-weighted temperature, electron number density, and entropy reveal the presence of the feedback heating in two ways. First, the density inside 70 kpc is similar to the initial (observed) conditions, meaning that the inflow due to the cooling is substantially quenched. Second, the temperature and entropy profiles do not fall below T_vir/3 and 10 keV cm², respectively. The profiles maintain a positive gradient, a sign that the cool core has not been destroyed by the outbursts. Note that the real comparison with X-ray observations must be based on emission-weighted temperatures: the pure mass-weighted T profile would appear much more disturbed with higher central values. The electron number density profile (n_e), also linked to the X-ray emissivity, appears consistent with observational data, with a flatter gradient compared to the pure CF run.

In Figure 3 (fourth column), we show the radial profile of t_cool/t_ff and its minimum as a function of time. Although we discuss the relationship between t_cool/t_ff and the presence of the multiphase gas in more detail in Section 4, here we briefly point out that whenever this timescale ratio falls below ≈10, thermal instabilities grow and cold gas clumps condense out of the hot phase. Initially, the cold phase has a filamentary structure; then it turns into small clumps which fall toward the center in about a free-fall time.

The SB_X map and the two-dimensional snapshot of n_e at 1.5 Gyr, shown in Figure 4 (first and second panels), correspond to a moment when the TI-ratio falls below 10. The jet outburst creates two cavities along the z-direction, which on larger scales become larger and move buoyantly outward (not shown in the density cuts). At the same time the perturbations due to the jets seem to lead to the formation of dense cold blobs up to ∼10 kpc away from the center. The presence of multiphase gas is also apparent in the n_e − T phase diagram (Figure 4, third panel), which displays two distinct diagonal islands occupied by the gas in the cluster core, with small variations away from pressure equilibrium.

Because the AGN heating is anisotropic close to the base of the jet and the cooling gas possesses a small amount of angular momentum, the cold gas tends to accrete and condense more along the equatorial region. In a realistic situation this gas cannot be entirely removed in a brief time (as in S11); consequently, a cold rotationally supported torus is formed, which provides most of the fuel for the AGN. The fuel supply from the torus, or the torus itself, can be temporarily disturbed by AGN outbursts, shocks, and the ICM turbulence. The central accumulation of cold gas could be further reduced by including star formation.

Overall, the heating due to realistic AGN jets, with mechanical efficiency ≈6 × 10⁻³, appears to reproduce the main properties of the observed clusters: quenching the cooling flows while preserving the cool cores for several Gyr. At the same time, the model predicts the formation of the multiphase medium in the form of cold blobs and filaments while satisfying the observational constraints on the amount of cold gas (see McDonald et al. 2010, 2011a, 2011b).

3.2.2.  = 10⁻²

We now discuss a run for higher mechanical efficiency ( = 10⁻²) compared to the one considered in the previous section ( = 6 × 10⁻³). We note that these values likely bracket the range of acceptable efficiency levels, as pointed out by previous computations (G11a; S11) and observations (Merloni & Heinz 2008). We use the same set of diagnostic plots as before and we only discuss the differences between these models. The results for = 10⁻² are shown in Figures 5 and 6, which are direct counterparts to Figures 3 and 4, respectively.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The increase in the mechanical efficiency of the jets leads to stronger outflows characterized by a slightly more peaked jet power distribution centered around ∼6 × 10⁴⁴ erg s⁻¹. Because in the adopted feedback scheme $\dot{M}_{\rm jet}$ is proportional to the jet power, the outflows also reach slightly higher mass-loading rates (≳ 0.5 M_☉ yr⁻¹). This means that some jet events carry larger amount of linear momentum and, thus, have more "piercing power" with which they plow through the ICM. Consequently, the jets release more of their energy at relatively larger distances from the cluster center, allowing a bit more cold gas accumulation along the equatorial region. The projected surface brightness map clearly shows the formation of a torus-like structure that was less pronounced in the lower efficiency run.

The initial increase in the average cooling (and accretion) rate is strongly quenched at around 1 Gyr since the beginning of the simulation, which is earlier than in the previous run. After this short phase, the system gradually settles to a quasi-equilibrium state, with cooling rates around 40–50 M_☉ yr⁻¹, twice the level seen in the lower efficiency case. Paradoxically, more powerful jets are less efficient in stopping the cooling flow, due to tunnel carving. The total injected energy is comparable in both runs (∼10⁶² erg), a sign that the self-regulation process is overall very similar, although with a slightly shorter active phase and a lower frequency of the outbursts (see jet power CDF). The total cold gas mass (2 × 10¹¹ M_☉) is slightly larger than in the preceding case, but still only 10% of the pure cooling flow value, which is consistent with the observational constraints.

During the entire evolution of this simulation, the outflows tend to produce bigger bubbles or a deeper channel (see also the n_e cut). The effect of more powerful jets is also to increase the scatter in the phase diagram: the gas in hot phase occupies a wider range of values in the n_e − T plane about the (imaginary) diagonal line representing average pressure.

The minimum TI-ratio falls below 10 only early in the evolution. Extended cold gas, in the form of blobs and filaments, is mainly formed at that epoch (see the n_e cut). At later times, more powerful events raise t_cool/t_ff well over 10, and the ICM becomes hotter. Even if the system settles later around a minimum TI-ratio of 13, extended cold filaments may be occasionally seen at moderate radii due to dredge-up of the cold torus by outflows. The uplift of the cold gas and the expansion of bubbles may further increase density, possibly triggering new thermal instabilities. This picture is also consistent with the observational analysis presented in M11 (see their Figure 12, right panel), suggesting that value 10 marks an approximate transition to the regime where the multiphase gas may appear.

3.2.3.  = 6 × 10⁻³, Wobbling

To test the robustness of our results we included jet wobbling in the run with mechanical efficiency = 6 × 10⁻³. Observational data are consistent with the possibility that the jets may precess or their orientation may change. Since the data are very limited, we keep the model very simple and change jet orientation whenever the accreted mass reaches a threshold of 10⁷ M_☉. For typical mass accretion rate of the order of ∼1 M_☉ yr⁻¹, this corresponds to a reorientation timescale of ∼10 Myr (e.g., Dunn et al. 2006).

In the first simulation with jet wobbling, we assumed that the probability of jet pointing in a given direction is given by a random uniform distribution over 4π steradians (the bipolar jets are always collinear). Initially, the jet orientation changes at a slower rate, ∼100 Myr, because the accretion rate is relatively low. Thus, in this early stage, the evolution is very similar to that for the non-wobbling case with efficiency of = 6 × 10⁻³. However, the energy is better spatially distributed, and thus the cooling rate is more quenched, reaching values of ≲ 10 M_☉ yr⁻¹. As the accretion rate accelerates, and the jet orientation starts to change on timescales of 10 Myr or less, the deposition of energy becomes nearly isotropic. A more uniform distribution of energy makes it more difficult to feed the black hole via the cold torus, which is now strongly disturbed. Rather than being used to drill a channel through the ICM, the AGN jet is no longer able to efficiently deliver the energy to the ICM beyond r > 15 kpc. This enables the gas to cool rapidly at these radii, which leads to the cooling catastrophe in about 1 Gyr. Larger efficiencies could possibly stop the cooling, but the nearly isotropic centrally concentrated energy distribution would destroy the cool core.

In order to remedy this situation, we considered models where the jet orientation is confined to two conical regions of half-opening angles of 75°, 45°, and 25°. As before, the probability distribution of the jet orientations is randomly uniform within these regions. In the first two cases, the final result is similar to the case considered above, but the cooling flow is quenched for an additional time of ≲ 1 Gyr.

For a half-opening angle of 25°, the jet energy distribution becomes narrower and more anisotropic. Apart from larger bubbles and a more frequent jet channel fragmentation, this computation resembles the previous best model r7-6em3: the cooling rate is steady for several Gyr (∼10–30 M_☉ yr⁻¹), with total cold gas mass around 10¹¹M_☉. Cold blobs are produced mainly in the first 1.5 Gyr when TI-ratio oscillates around 10. At later times, the system settles to a hotter state characterized by the TI-ratio of ∼15–16. Jet events are also more variable than in the r7-6em3 model, changing rapidly from 10⁴³ to 10⁴⁵ erg s⁻¹ because more inclined outbursts are able to partially dismantle the central torus, temporarily inhibiting accretion.

In summary, randomly wobbling jets with small and moderate opening angles produce similar results to unidirectional narrow outflows, especially for stopping the cooling flow catastrophe and generating extended multiphase gas. On the contrary, almost spherically symmetric injections prevent proper fueling of the black hole and, in the long term, the feedback cannot quench the cooling flow. Therefore, we warn against implementing subgrid models based on isotropic injection, as is customarily done in smoothed particle hydrodynamic simulations (e.g., Power et al. 2011). The other relevant feature is that the central accumulated cold torus can be frequently destroyed, depending on the wobbling parameters.

3.3. AGN Feedback: t_cool/t_ff = 21

3.3.1.  = 6 × 10⁻³

In order to study the effect of feedback in slower cooling systems, we generated new initial hydrostatic conditions as follows. We decreased the central density preserving the initial profile slope at r ≳ 80 kpc; simultaneously, we increased the central temperature to ensure that the ICM remains in hydrostatic equilibrium. The new initial conditions correspond to the minimum TI-ratio of ∼21. The new initial profiles are plotted in Figure 7, black dashed lines (third column). The entropy is elevated to 30 keV cm² at 10 kpc—still within observational bounds (Cavagnolo et al. 2009)—while the total gas mass is decreased only by 8%. The results from the runs corresponding to these new initial conditions and AGN mechanical efficiency of = 6 × 10⁻³ are shown in Figures 7 and 8. These figures are analogous to those previously discussed.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Overall, the behavior of the system is quite similar to what we observed before. The model is successful in predicting a globally stable ICM core, low cooling rates (∼50 M_☉ yr⁻¹, 10% of pure CF), preservation of the cool core (e.g., entropy, SB_X), typical features of AGN jets (bubbles, shocks, turbulence), and the amount of multiphase medium consistent with observations. The main result is that the local thermal instability is triggered despite the fact that the initial minimum TI-ratio exceeds ∼10. It is only when the instantaneous value of this ratio falls below the critical threshold of ∼10 that the multiphase medium suddenly appears in the simulation (e.g., after 1 Gyr). As before, this leads to increased feedback that brings the system back to the high TI-ratio regime, which in this case is well above ∼20 and higher than in the previously considered systems.

The typical jet power (few 10⁴⁴ erg s⁻¹) and total injected energy (4 × 10⁶¹ erg; single jet) are lower compared to r7-6em3, mainly due to the shallower initial density profile and thus reduced accretion rates. However, the AGN outflows have a relatively stronger impact on the lighter ICM, secularly increasing the TI-ratio back above 20. As a consequence, spatially concentrated cooling, in the form of a rotating torus, is present during the entire evolution. In Section 4 we make an in-depth comparison between the models with different initial t_cool/t_ff ratio.

3.3.2.  = 10⁻²

The total cold gas mass in the model discussed above reaches 2 × 10¹¹ M_☉, which is higher than r7-6em3. Although this is observationally acceptable, some systems show smaller amounts of cold gas mass. Therefore, we experimented with increasing the mechanical efficiency in order to decrease the amount of cold gas (figures not shown). However, as argued in Section 3.2.2, higher efficiencies lead to stronger and more mass-loaded jets that deposit their mechanical energy farther away from the center. If in the initial outbursts the system appears to be more heated due to shocks, the situation reverses later, when the jet carves a deeper channel in the ICM. Thus, more accretion occurs along the equator and a strong cooling flow develops. As seen before, this feature could be alleviated by a randomly wobbling outflow, with a less anisotropic energy deposition.

Overall, the simulations with higher minimum initial TI-ratio 21 can again produce extended multiphase gas whenever the instantaneous t_cool/t_ff falls ≲ 10. In fact, in the early phase of the evolution, when cooling slightly prevails over heating, the TI-ratio slowly decreases, crossing the threshold in about a Gyr. Later, stronger outbursts set the system back to a hotter state and higher TI-ratio, until the cycle restarts again (it may require several Gyr). The main trend for the entire 5 Gyr is therefore concentrated cooling occurring in the core.

Contrary to the more massive models (TI-ratio of 7), the right efficiencies in order to quench the cooling flow and preserve the cool core should not be more than 6 × 10⁻³. A less massive system is in fact more susceptible to mechanical outbursts and requires a more delicate feedback, analagous to the environment of galaxy groups or ellipticals, as has been shown in previous works (G11b, S11).

The fundamental result of our computed models, consistent with observations, is that extended cold filaments and blobs condense out of the hot phase whenever t_cool/t_ff ≲ 10 or, equivalently, when the central entropy becomes ≲ 30 keV cm². The cold clumps then fall toward the center in about a free-fall time, forming a central cold rotating torus and in part boosting the black hole accretion rate. The formation of cold filaments is a result of interplay between thermal instability and gravity. Previous idealized analyses by M11 and S11 enforced strict thermal equilibrium and found that multiphase gas is generated according to the same threshold.¹³

A realistic feedback displays a variegated self-regulated evolution, producing complex observed phenomena, which are not captured by simple idealized prescriptions. In the present work, we tested bipolar AGN jets, tied to the mass accretion rate at very small radii.¹⁴ The asymmetrical input of energy produces strong elliptical shocks and regularly inflates hot bubbles, which rise outward. At the same time the inflow is not entirely stopped: along the equatorial direction the gas, mainly cold, continues to fuel the black hole. Turbulence, entrainment, and mixing dissipate the directional kinetic energy of the jets, further inhibiting the cooling flow near the center. When the jet ram pressure is substantial, the cold gas is also dredged up, as observed in a few cases such as Perseus and Hydra (Hatch et al. 2006; Gitti et al. 2011). The cycle repeats when cooling starts to dominate again and the cold gas condensation ignites more powerful outbursts, setting the system in a slightly hotter state.

Overall, thermal equilibrium is only roughly maintained, on large spatial and temporal scales. For example, the two maps in Figure 9, for best model r7-6em3, display the net thermal energy increase (red) or decrement (blue) per volume and time. The internal energy equation is given by

$$\begin{equation} \frac{\partial e}{\partial t} = -{\bf \nabla }\cdot \left(e\, {\bf v}\right) -P({\bf \nabla }\cdot {\bf v}) -\: n_{\rm {e}} n_{\rm {i}} \Lambda (T), \end{equation} \tag{ 16 }$$

where e is the internal energy density (erg cm⁻³). The maps indicate that "heating" (∂e/∂t > 0) and "cooling" (∂e/∂t < 0) alternate on scales of a few kiloparsecs. Averages of ∂e/∂t in radial shells of width ≳ 4–5 kpc suggest that a rough thermal equilibrium is maintained, although cooling appears slightly dominant in the central 10 kpc, while the jets tend to deposit their energy at larger radii. Equation (16) does not, however, fully capture the irreversible heating associated with turbulent mixing, shocks, and wave damping due to numerical effects. A better direct indicator of the balance between cooling ($\mathcal {L}$) and heating ($\mathcal {H}$) is the variation of the entropy S (∝ln K, where K ≡ k_bT/n^2/3_e), which is insensitive to adiabatic compressions or expansions:

$$\begin{equation} \rho T \frac{dS}{dt} = \mathcal {H} - \mathcal {L}, \end{equation} \tag{ 17 }$$

where d/dt is the Lagrangian derivative.

Zoom In Zoom Out Reset image size

The right panels in Figure 9 show the time variations in radial shells of the "astrophysical" entropy K for the whole ICM, with respect to the initial value K₀. As expected the entropy increases on a secular timescale due the continuous action of the AGN feedback. On smaller temporal scales, e.g., a few typical cooling times (100–200 Myr), the amplitude of the fluctuations in the six thick spherical shells (centered at r = 22, 18, 14, 10, 6, and 2 kpc; shown from the top to bottom rows, respectively) can reach two or three times the mean value, especially in the nuclear region. Although the fluctuations in K can be substantial, its mean saturates at a stable level. This suggests that the cluster core does not dramatically deviate from the state of (quasi) equilibrium between cooling and heating, at least on scales of tens of kiloparsecs and hundreds of Myr. This was also suggested by the low cooling/accretion rates compared to a pure cooling flow. The assumption of strict global thermal equilibrium made in the work of M11 and S11 is therefore very idealized. It is thus particularly remarkable that our criterion for the formation of the multiphase gas is in good agreement with their result.

Unlike in S11, the central kiloparsec region is not excised, allowing the cold gas (either nuclear cooling or fallen blobs) to form an angular momentum-supported torus as observed in the core of M87 (e.g., Ford et al. 1994). Such a disk can be partially dismantled when the jet direction is perturbed by turbulence or by random jet precession (Section 3.2.3). Nevertheless, the central dense disk is a common feature among elliptical galaxies (e.g., Macchetto et al. 1997; Mathews & Brighenti 2003; Verdoes et al. 2006). Centrally concentrated atomic gas is also seen in the survey of McDonald et al. (2010, 2011a). The central disk is certainly a nursery for new stars. Since the star formation efficiency is quite uncertain (a few percent, up to 50%, with an average value of 15%, according to McDonald et al. 2011b), we decided to not model star formation via a sink term. Thus, the central accumulation of cold gas could be further reduced.

The presence of multiphase gas is clearly seen in the n_e − T joint distributions. The cold dense phase (10⁴–10⁵ K), in the upper left corner, is present in large quantities and detached from the hot rarefied gas (>few 10⁷ K), with a lack of gas at intermediate temperatures. A recent Hα survey (McDonald et al. 2010, 2011a) has shown that extended cold filaments are found in ≳ 35% of all clusters. It is striking that some of their maps (see their Figure 4) resemble our computed morphologies, with extended blobs and filaments up to ∼20 kpc: e.g., A0496, Hydra A, A1644, Sersic 159-03, and A1795 (the cluster on which our model is based). Another case with substantial extended Hα emission is the Perseus cluster (Salomé et al. 2006; Fabian et al. 2008). Moreover, the total cold gas mass we obtain, ∼10¹¹ M_☉, is consistent with the constraints of Edge (2001), whose observations indicate clusters harboring up to 3 × 10¹¹M_☉ of molecular gas.

The white noise perturbations in the initial conditions are wiped out over a few cooling times. In fact, in the pure cooling flow computation, multiphase gas is only possible in such a restricted interval, while the secular evolution is dominated by monolithic cooling. On the contrary, AGN feedback drives and sustains the fluctuations in the long-term evolution, although the original seeds initially help to distribute the jet energy better due to the interaction of jet with the inhomogeneous ICM. The thermal state of the cluster core depends mainly on the feedback efficiency: with smaller the frequency of cold filaments condensing out of the hot phase is larger (Figures 3, 5, 7). Nevertheless, higher initial t_cool/t_ff (and lower density) generates extended multiphase gas according to the TI criterion, but thermal equilibrium is reached at later epochs. It is interesting to note that AGN feedback is more effective in lighter systems, therefore requiring lower efficiencies (≲ 6 × 10⁻³), as has been shown for galaxy groups (G11b; S11).

To better quantify the presence of extended multiphase gas, Figure 10 shows a scatter plot of the average distance of cold gas (R_cold; excluding r < 3.5 kpc) as a function of t_cool/t_ff. The points are color coded according to time: from blue (early) to red (late). For the run with the initial t_cool/t_ff = 7 (upper panel), the cold gas is much more spatially extended, especially in the first 2 Gyr. The run with less dense initial conditions (TI-ratio 21) shows a lower frequency of cold blobs at large distances because t_cool/t_ff ≳ 10 most of the time (compare Figures 3 and 7). In both cases, the points clustered around ∼4 kpc correspond to the rotating cold torus. Since this plot shows the current average distance rather than the radius where the cold gas originates, extended filaments appear even when TI-ratio is elevated to ∼12–15 due to AGN heating. This tendency is further enhanced if cold blobs are dredged up by powerful outflowing jets. Moreover, as the TI-ratio is averaged in spherical shells, a few cold zones may contribute to R_cold, without changing the average t_cool/t_ff much. That is, some extended Hα emission is consistent with t_cool/t_ff slightly greater than 10.

Zoom In Zoom Out Reset image size

Figure 11 displays R_cold gas as a function of time. With denser initial conditions (r7-6em3), t_cool/t_ff is less than 10 for several short intervals within 2 Gyr and we expect extended cold gas.¹⁵ Indeed, the average distance of the cold gas is significantly larger than that of the cold torus (≈4 kpc) during this period. Similarly, for the lower initial density run we see extended cold filaments at 1–1.5 Gyr, again when t_cool/t_ff ≲ 10; later times are, however, dominated by the nuclear cold gas. We conclude that the instantaneous value of t_cool/t_ff is a robust indicator of extended cold gas formation, while nuclear cold gas can be present irrespective of the TI-ratio (e.g., bottom panel of Figure 10).

Zoom In Zoom Out Reset image size

It is likely that other physical processes play an important role in shaping the morphology of condensing multiphase gas. Sharma et al. (2010) and M11 showed that anisotropic thermal conduction leads to cold gas stretched along the local magnetic field direction. In order to be able to fully understand the exact topology of cold filaments, this physics should be included in the simulations (with an adequately higher resolution). However, the generation of cold filaments, whenever t_cool/t_ff ≲ 10, is not likely to be affected by the presence of anisotropic thermal conduction.

The other physics that is missing in our simulations involves cosmological accretion and mergers, which drive additional turbulence, inducing more instabilities. Like AGN feedback, major mergers are associated with entropy generation and strong heating of cluster cores, raising the cooling time over the cluster age. Thus, mergers may be responsible for extreme non-cool-core clusters. However, they are also rare with timescales much longer than the cluster core cooling times (e.g., Fakhouri et al. 2010). Therefore, AGN feedback, rather than mergers, is expected to be responsible for balancing cooling in most systems.

Our theoretical models are able to solve the cooling flow problem for several Gyr. The key factor resides in adopting realistic mechanical AGN jets. The bulk of AGN heating is in fact a combination of shock dissipation and gradual thermalization of the mechanical energy (see also Brüggen et al. 2007). Random motions, either due to AGN, multiphase gas, and cosmological turbulence strongly help the deposition of the feedback energy in the cool core, in contrast to idealized symmetric atmospheres (e.g., Vernaleo & Reynolds 2006).

An important issue is whether our results depend sensitively on resolution. The morphology of the cold phase depends on the resolution. We note that finite resolution may prevent the disruption and mixing of the small infalling clouds with the ambient ICM. Thus, limited resolution could in fact mimic some of the effects of the physics that we do not include (e.g., magnetic draping). We performed a test (r7-6em3c) with two-times-lower resolution, but with the same resolution around the jet injection region¹⁶ and the same accretion/depletion radius. We obtained very similar results. In particular, a steady global thermal equilibrium is reached in a few Gyr, with only a slightly larger final cold mass (2 × 10¹¹ M_☉). The evolution of the TI-ratio and condensed cold blobs appears also similar to the higher resolution model (r7-6em3), although the increased numerical dissipation tends to slow down the growth of thermal instabilities at larger radii. We note also that our simulations have a resolution ∼10 times higher compared to the work done in G11a. Although in G11a we adopted a cold mass dropout and a different jet injection method, AGN jets could also stop the cooling flow, in line with observations. Long-term thermalization seems therefore not dramatically affected by resolution or the exact heating prescription. The fine details of the feedback can, however, change with higher resolution, e.g., jets tend to become more frequent (cf. G11a).

Our main conclusions can be summarized as follows.

• 1.  
  Mechanical AGN feedback self-consistently linked to black hole accretion (with efficiency of 5 × 10⁻³–10⁻²) can solve the cooling flow problem, quenching the cooling rate down to ∼10 M_☉ yr⁻¹ (1%–10% of the pure cooling flow value), preserving the cool-core structure in agreement with observations (Peterson et al. 2001, 2003; Ettori et al. 2002; Tamura et al. 2003; Peterson & Fabian 2006; McNamara & Nulsen 2007). Typical outflow powers are 10⁴⁴–10⁴⁵ erg s⁻¹, with velocities ∼10⁴ km s⁻¹ and mass-loading rates ∼M_☉ yr⁻¹ (see Crenshaw et al. 2003; Morganti et al. 2005, 2007; Nesvadba et al. 2008, 2011; Tombesi et al. 2010). Best models are associated with narrow jets with fixed inclinations or fast random wobbling with a half-opening angle of ≲ 25°.
• 2.  
  Multiphase gas, with cold (T ∼ 10⁴ K) and hot (T > 10⁷ K) components, is a direct consequence and source of feedback. The cold gas comes in two phases: extended phase (up to 10–20 kpc of distance form the center) in the form of blobs and filaments generated by thermal instabilities (and occasionally due to dredge-up by jet ram pressure) and a more spatially concentrated one in the form of central rotating torus due to nuclear cooling or accumulated fallen blobs. The Hα survey of McDonald et al. (2010, 2011a) shows both kinds of Hα emitting gas. Nuclear cold reservoirs of gas especially in the form of rotating disks are also a common feature among ellipticals (Ford et al. 1994; Macchetto et al. 1997; Mathews & Brighenti 2003; Verdoes et al. 2006; a star formation sink term and jet wobbling can reduce it considerably though). Total cold gas masses are ∼10¹¹ M_☉ after 5 Gyr, within the observational constraints (≲ 3 × 10¹¹ M_☉; Edge 2001; Salomé et al. 2003).
• 3.  
  The generation of extended cold gas is dictated by an important threshold: whenever the instantaneous ratio of the cooling time to the free-fall time is ≲ 10 thermal instabilities lead to cold gas condensing out of the hot ICM at large radii. On the contrary, when t_cool/t_ff > 10, slow runaway cooling occurs only within the central few kiloparsecs. This threshold seems to be surprisingly robust, as predicted by simple analytical arguments (M11) and simulations with idealized feedback (S11). The criterion is also supported by recent observations (Figure 12 in M11; Cavagnolo et al. 2008; Rafferty et al. 2008; Voit et al. 2008).
• 4.  
  Multiphase ICM can exist in cool-core clusters irrespective of their initial minimum t_cool/t_ff. With a realistic feedback the cluster core always finds a way to self-regulate toward a state of global quasi-thermal equilibrium on large spatial and temporal scales. The system forgets about the initial state in a few cooling times and the eventual evolution depends primarily on the feedback efficiency. Lower results in frequent episodes of extended cold gas. Intense cooling episodes are followed by increased jet activity and the core returns to a higher entropy level. This phase is more prolonged and intense in less massive systems (e.g., initial t_cool/t_ff = 21). The reduced gas inflow leads to the decreasing jet power, allowing a slow spatially concentrated cooling to start again (t_cool/t_ff > 10).

Overall, the main thermal properties of cluster cores are determined by global quasi-thermal equilibrium, with an amazingly detailed complexity of cold filaments/torus, jets, bubbles, and shocks. Additional processes, such as thermal conduction, magnetic fields, and cosmological accretion are required to make more precise comparisons with observations. Nevertheless, we showed that multiphase gas due to thermal instabilities is an essential element of the feedback cycle: an end product of cooling and fuel for the central black hole.

The software used in this work was in part developed by the DOE NNSA-ASC OASCR Flash Center at the University of Chicago. We acknowledge the NASA awards SMD-10-1609, SMD-11-2209 (Pleiades), the CINECA awards HP10BPTM62, HP10BOB5U6 (SP6), and the University of Michigan for the availability of high-performance computing resources and support. M.G. thanks F. Brighenti for useful discussions and P. Temi, as the principal support scientist at the NASA/Ames base. M.R. acknowledges the NSF grant 1008454. Support for P.S. was partially provided by NASA through Chandra Postdoctoral Fellowship grant PF8-90054 awarded by the Chandra X-Ray Center, which is operated by the Smithsonian Astrophysical Observatory for NASA under contract NAS8-03060. We are also thankful for the hospitality of the Kavli Institute for Theoretical Physics (KITP) at UC Santa Barbara, where this work was initiated.