HOW AGN JETS HEAT THE INTRACLUSTER MEDIUM—INSIGHTS FROM HYDRODYNAMIC SIMULATIONS

H.-Y. Karen Yang; Christopher S. Reynolds

doi:10.3847/0004-637X/829/2/90

1. INTRODUCTION

Galaxy clusters with short central cooling times, or cool-core (CC) clusters, are expected to host massive gas inflows that are absent observationally, i.e., the "cooling-flow problem" (Fabian 1994). The classical cooling-flow model predicts typical mass deposition rates of hundreds to a thousand solar masses per year. However, observationally there is a dearth of gas below ∼1 keV (e.g., Peterson et al. 2003), and the star formation rates within observed CCs are suppressed by roughly an order of magnitude (McNamara & O'Connell 1989; O'Dea et al. 2008; McDonald et al. 2011b; Hoffer et al. 2012). Various heating mechanisms have been proposed in order to circumvent catastrophic cooling (McNamara & Nulsen 2007 for a detailed review). However, processes that do not involve feedback (e.g., Burns et al. 2008), such as thermal conduction (Narayan & Medvedev 2001; Zakamska & Narayan 2003; Voigt & Fabian 2004), often require fine-tuning (Stewart et al. 1984; Bregman & David 1988; Soker 2003) and likely only provide part of the heating, if any.

It is now widely accepted that feedback from active galactic nuclei (AGNs) at cluster centers holds the key to solving the cooling-flow problem. This is strongly motivated by the prevalence of AGN jet-inflated X-ray cavities/bubbles within CC clusters (Best et al. 2007; Mittal et al. 2009; Bîrzan et al. 2012) and the correlation between the cavity power (which measures the mechanical power of AGN jets) and the X-ray luminosity within the cores (Rafferty et al. 2006; Dunn & Fabian 2008). Theoretically, models that incorporate self-regulated energy injections from the central supermassive black holes (SMBHs) successfully prevent run-away cooling (Cattaneo & Teyssier 2007; Sijacki et al. 2007; Dubois et al. 2010; Yang et al. 2012b). More recent simulations (Gaspari et al. 2012; Li & Bryan 2014; Gaspari 2015; Prasad et al. 2015) that are based on SMBH accretion of cold gas, i.e., the "cold-mode" feedback (Pizzolato & Soker 2005), combined with momentum-driven AGN jets, further reproduce the positive temperature gradients (Allen et al. 2001b; Cavagnolo et al. 2008; Hudson et al. 2010), the signatures of multiphase gas of the observed CCs (e.g., Edge 2001; Salomé & Combes 2003; Fabian et al. 2006; McDonald et al. 2011a), and the slope of the soft X-ray spectrum (e.g., Peterson et al. 2003; Sanders et al. 2010). Recent simulations by Yang & Reynolds (2016) that simultaneously include cooling, AGN feedback, and anisotropic conduction, show that even for unimpeded conductivity along magnetic field lines, conductive heating is likely significant only for the most massive clusters, while AGN are still the dominant sources of heating for most clusters. All the above highlight the importance of AGN feedback within CC clusters.

Despite the successes, the details of how the AGN jets heat the intracluster medium (ICM) are still elusive. Many channels of AGN heating are proposed, including cavity heating (Churazov et al. 2001; Brüggen 2003), shock heating (Fabian et al. 2003; Nulsen et al. 2005; Randall et al. 2015), dissipation of sound waves (Ruszkowski et al. 2004; Fabian et al. 2005), mixing by convection or turbulence (David et al. 2001; Kim & Narayan 2003), turbulent dissipation (Dennis & Chandran 2005; Zhuravleva et al. 2014), cosmic rays (Guo & Oh 2008; Pfrommer 2013), and mixing with hot thermal gas within the bubbles (Hillel & Soker 2016). Some of the mechanisms only heat the ICM locally, e.g., cavity heating and bubble mixing should apply in the wakes and surrounding the bubbles, respectively. Moreover, it has been challenging for many processes to produce radial profiles of heating that match those of radiative cooling (Brighenti & Mathews 2002; Fabian et al. 2005; Mathews et al. 2006). Therefore, it is most likely that multiple processes are involved to maintain the intricate balance between cooling and heating in the CCs.

However, two big questions remain: (1) what are the dominating heating mechanisms? and (2) how do AGN jets heat the ICM isotropically, given that the jets are intrinsically directional (Vernaleo & Reynolds 2006)? In this work, we attempt to address these questions by taking a closer look at one of the most successful models to date, namely, three-dimensional (3D) hydrodynamic simulations with cold-mode accretion and momentum-driven feedback. We utilize tracer particles and track specific terms in the hydrodynamic equations as guidance in order to better understand the heating and cooling processes as well as energy transport. We quantify the relative importance of all relevant processes, including radiative cooling, advection, adiabatic compression/expansion, shock heating, and mixing with hot bubble gas. We find that, though the AGN jets primarily heat regions surrounding the bubbles by mixing and weak shocks, the heat is isotropized by the process of "gentle circulation" within the cluster core. Since such simulations have shown broad consistencies with observations, understanding how AGN heating works in the simulations may provide valuable insights into AGN heating in the real universe.

The structure of the paper is as follows. We summarize the essential elements of our simulation and analyses in Section 2. In Section 3.1, we describe the overall cluster evolution and integrated properties. Then we present analyses from the Lagragian and Eulerian perspectives in Sections 3.2 and 3.3, respectively. The relative contributions from relevant processes are quantified in Section 3.4. We discuss the implications of our results on turbulent heating in Section 4.1 and limitations of our simulation in Section 4.2. Finally, we conclude our findings in Section 5.

2. METHODOLOGY

We carry out a 3D hydrodynamic simulation including radiative cooling and self-regulated AGN feedback in an idealized Perseus-like cluster using the adaptive-mesh-refinement (AMR) code FLASH (Fryxell et al. 2000; Dubey et al. 2008). The simulation setup is identical to the AGN-only simulation (Model A) in Yang & Reynolds (2016; except that we neglect magnetic fields), to which we refer the readers for details. In the following, we give a brief summary and emphasize methods that are specific to this study.

The standard hydrodynamic equations are solved:

$\begin{eqnarray}&&\displaystyle \frac{\partial \rho }{\partial t}+{\rm{\nabla }}\cdot (\rho {\boldsymbol{v}})=0,\end{eqnarray} \tag{ 1 }$

$\begin{eqnarray}&&\displaystyle \frac{\partial \rho {\boldsymbol{v}}}{\partial t}+{\rm{\nabla }}\cdot (\rho {\boldsymbol{v}}{\boldsymbol{v}})+{\rm{\nabla }}P=\rho {\boldsymbol{g}},\end{eqnarray} \tag{ 2 }$

$\begin{eqnarray}&&\displaystyle \frac{\partial e}{\partial t}+{\rm{\nabla }}\cdot [(e+P){\boldsymbol{v}}]=\rho {\boldsymbol{v}}\cdot {\boldsymbol{g}}+{ \mathcal H }-{ \mathcal C },\end{eqnarray} \tag{ 3 }$

where ${\boldsymbol{g}}$ is the gravitational field, ρ, ${\boldsymbol{v}}$ , $e=0.5\rho {v}^{2}+{e}_{{\rm{i}}}$ , and $P=(\gamma -1){e}_{{\rm{i}}}$ are the density, velocity, total energy (including kinetic and internal energy), and thermal pressure of the gas. An equation of state for ideal gas with γ = 5/3 is assumed. ${ \mathcal H }$ and ${ \mathcal C }$ are the heating and cooling terms, respectively. For radiative cooling, ${ \mathcal C }={n}_{{\rm{e}}}^{2}{\rm{\Lambda }}(T)$ , where ${n}_{{\rm{e}}}=\rho /{\mu }_{{\rm{c}}}{m}_{{\rm{p}}}$ is the electron number density and Λ(T) is the cooling function. In order to focus on hydrodynamic effects alone, magnetic fields and other transport effects (e.g., viscosity, conduction) are not included and will be considered in future work.

The simulation domain is 1 Mpc on a side and is adaptively refined on steep temperature gradients up to an AMR refinement level of 7, which corresponds to a resolution element of 1.95 kpc. The diode boundary condition is used, which is similar to the outflow boundary condition but prohibits the inflow of gas into the simulation domain. Radiative cooling is computed using the tabulated table of Sutherland & Dopita (1993) assuming one-third solar metallicity. The gas profiles of the cluster are initialized using empirical fits to the observed Perseus cluster assuming hydrostatic equilibrium within a fixed NFW (Navarro et al. 1996) gravitational potential. The central SMBH is assumed to be fed by cold gas in its vicinity, i.e., the "cold-mode" accretion (e.g., Gaspari et al. 2012, 2015; Li et al. 2015; Voit et al. 2015). The cold gas (T < 5 × 10⁵ K) is dropped out from the hot phase and replaced by passively evolving tracer particles. The assumptions and consequences of this treatment are discussed in detail in Yang & Reynolds (2016). The SMBH accretion rate is estimated by ${\dot{M}}_{{\rm{BH}}}={M}_{{\rm{cold}}}/{t}_{{\rm{ff}}}$ , where M_cold is the total amount of cold gas within an accretion radius of r_accre = 4 kpc and t_ff = 5 Myr is the approximate free-fall time at r_accre for Perseus. The accretion rate is then used to compute the mass, momentum, and energy returned by the bipolar AGN jets along the z-direction following these equations:

$\begin{eqnarray}\dot{M} & = & \eta \,{\dot{M}}_{\mathrm{BH}},\\ | \dot{P}| & = & \sqrt{2\eta \epsilon }{\dot{M}}_{\mathrm{BH}}c\\ \dot{E} & = & \epsilon \,{\dot{M}}_{\mathrm{BH}}{c}^{2},\end{eqnarray} \tag{ 4 }$

where η is the mass loading factor and the injected energy is purely kinetic. We adopt $\eta =1$ and a feedback efficiency of = 0.001 for the current simulation. The feedback is applied to a cylinder with a radius of 2.5 kpc and height of 4 kpc. A small angle (15°) jet precession with a period of 10 Myr is assumed in order to improve coupling between the AGN feedback energy and the ICM. We note that the SMBH accretion and feedback prescriptions adopted here have been extensively investigated (e.g., Gaspari et al. 2012; Yang et al. 2012a; Li et al. 2015; Prasad et al. 2015; Yang & Reynolds 2016). The chosen parameters have successfully reproduced the general properties of the CCs, such as self-regulation, positive temperature gradient, and the signatures of multiphase gas within the cluster cores. Therefore, the simulation is representative and suitable for studying how the heating and cooling balance each other in the self-regulated feedback process.

We utilize passively evolving tracer particles in order to follow the thermodynamic history of the ICM. The tracer particles are initialized at the beginning of the simulations within a sphere of radius r = 100 kpc with particle spacings of 2 kpc. The particles record the hydrodynamic variables of the fluid along their trajectories. The particle data is dumped into output files at a frequency of 0.005 Gyr.

In order to distinguish between AGN jets and the ICM, a tracer fluid is used to identify the AGN jets. For each grid cell, the mass fractions of the jets and the gas are denoted as f_jet and f_gas = 1 − f_jet, respectively. The value of f_jet is initially set to zero everywhere in the simulation domain. When the jets are launched, f_jet = 1 is set within the jet nozzle. With subsequent advection and mixing, regions influenced by the jets could then have values between 0 and 1. It follows that a mass-weighted quantity X is defined as $\int {f}_{{\rm{gas}}}\rho {{Xd}}^{3}x/\int {f}_{{\rm{gas}}}\rho \,{d}^{3}x$ .

One of the candidates for providing heat to the ICM is heating by shocks. As we will see in Section 3.1, the shocks found in our simulation have typical Mach numbers below 2. We hereby briefly summarize how we estimate the ICM heating by weak shocks. We identify the shocks based on negative velocity divergence and jumps in gas pressure

$\begin{eqnarray}&&\displaystyle \frac{{\rm{\Delta }}P}{P}=\displaystyle \frac{{P}_{2}-{P}_{1}}{{P}_{1}}=\displaystyle \frac{2\gamma }{\gamma +1}y,\end{eqnarray} \tag{ 5 }$

where P₁ and P₂ are the preshock and postshock pressure, respectively. The quantity y is related to the Mach number by

$\begin{eqnarray}&&y=\displaystyle \frac{{\rho }_{1}{v}_{{\rm{s}}}^{2}}{\gamma {P}_{1}}-1={{ \mathcal M }}^{2}-1,\end{eqnarray} \tag{ 6 }$

where ρ₁ is the preshock density and v_s is the shock speed. For reference, ${ \mathcal M }=1.1,1.5$ and 2.0 correspond to ${\rm{\Delta }}P/P\sim 0.26,1.56$ and 3.75, respectively. The jump in specific entropy for a weak shock is

$\begin{eqnarray}&&{ds}\simeq \displaystyle \frac{2\gamma {k}_{{\rm{B}}}}{3{(\gamma +1)}^{2}\mu {m}_{{\rm{H}}}}{y}^{3},\end{eqnarray} \tag{ 7 }$

which then could be used to compute the amount of energy per unit volume gained by the gas, ${dH}=\rho {Tds}$ . The volumetric heating rate by weak shocks is then estimated by ${dH}/{dt}$ , where ${dt}=0.005\,{\rm{Gyr}}$ is the time interval between two subsequent output files. For the evaluation of shock heating (see Section 3.4), we only include shocks with ${\rm{\Delta }}P/P\geqslant 0.2$ because they account for most visible entropy jumps in the simulation (see Figure 8). However, we verified that the amount of shock heating is insensitive to this choice because weak shocks around the threshold do not have a dominant contribution to heating. Note that some previous works apply additional criteria for diagnosing shocks (e.g., Ryu et al. 2003; Brüggen et al. 2005; Bryan et al. 2014) and therefore our estimates of shock heating should be considered as upper limits.

3. RESULTS

3.1. Overall Evolution

We describe the overall evolution and integrated properties of the cluster in this section. Since other aspects of self-regulated AGN feedback in clusters have been extensively discussed in the literature (e.g., Sijacki et al. 2007; Gaspari et al. 2011; Yang et al. 2012b; Li & Bryan 2014; Prasad et al. 2015), we focus on characteristics that are relevant for this study.

Figure 1 shows the evolution of AGN jet power and X-ray luminosity within r = 100 kpc, which is about the cooling radius of the cluster (defined here as the radius at which the cooling time is 3 Gyr). At the beginning of the simulation, the X-ray luminosity increases as the cluster contracts due to energy losses by radiative cooling, forming a cooling flow (Fabian 1994; Li & Bryan 2012). The ICM cools and forms clumps of cold gas due to local thermal instabilities (TI) when t_c/t_ff ≲ 10 (Gaspari et al. 2012; McCourt et al. 2012; Sharma et al. 2012; Li & Bryan 2014; Meece et al. 2015). The central SMBH is fed by the cold gas in its vicinity and then subsequent jet activity is triggered at t ∼ 0.3 Gyr, which is approximately the initial cooling time of the cluster. The kinetic energy of the AGN jets is transformed into heat of the ICM (via processes that will be discussed in detail later), raising t_c/t_ff above the threshold for TI and suppressing cold clump formation until the ICM cools again and fueling is resumed. Thereby, a self-regulated AGN feedback cycle is established, maintaining the cluster in a quasi-equilibrium state after t ∼ 0.7 Gyr.² After t ∼ 0.7 Gyr, the X-ray luminosity stays roughly constant (∼1.3 × 10⁴⁵ erg s⁻¹), and the average jet power is about 3 × 10⁴⁵ erg s⁻¹, implying ∼40% of the kinetic energy is transformed into thermal energy within 100 kpc. The average mass accretion rate for the simulated cluster is ∼15% of the inferred mass deposition rate of the observed Perseus cluster (Allen et al. 2001a), much suppressed compared to that of a pure cooling flow and consistent with the range of mass deposition rates (e.g., Hudson et al. 2010) and star formation efficiencies (e.g., McDonald et al. 2011b; Hoffer et al. 2012) for observed CC clusters.

Figure 2 shows slices of gas density, temperature, projected X-ray emissivity, and the jet mass fraction (f_jet). Snapshots are taken at three epochs, t = 0.935, 2.500, and 2.915, which represent three characteristic phases of AGN activity: single event, quiescent, and multiple episodes. In general, the AGN jets inflate underdense bubbles, producing X-ray cavities similar to those observed (Fabian et al. 2000). The cavities are filled with high fractions of jet materials that are hot (∼few tens of keV) due to the initial thermalization of the jet kinetic energy by shocks. When the bubbles rise and expand due to buoyancy, they mix with the ICM (i.e., f_jet decreases) as they are shredded due to Rayleigh–Taylor, Kelvin–Helmholtz, and Richtmyer-Meshkov instabilities. Weak shocks can also be seen in the temperature and X-ray maps. The Mach numbers of the weak shocks ( $1\lesssim { \mathcal M }\lt 2.1$ ) are consistent with those observed in nearby CC clusters (e.g., Fabian et al. 2006; Blanton et al. 2009; Randall et al. 2015).

Since our jets are directional by construction, there is a clear distinction between the ICM properties within and outside the "jet cones," which are defined as the region ∼30° away from the z-axis (see last panel in Figure 2). The buoyantly rising cavities containing partially mixed jet materials primarily lie within the jet cones; this is also where turbulence is generated by hydrodynamic instabilities during the disruption of bubbles (Yang & Reynolds 2016). Much of the ambient ICM (i.e., outside the jet cones) does not directly contact the jet materials, but is influenced by the weak shocks and sound waves. As we will show later, the primary heating mechanisms are very different between the jet cones and the ambient region. In reality, the bubbles may have a more isotropic distribution like those observed in Perseus (Fabian et al. 2000). We will discuss the implications of our results in this aspect in Section 4.2.

Figure 3 shows the profiles of gas density, temperature, and entropy (top to bottom rows) averaged within the entire cluster (left column), the jet cones (middle column), and the ambient region (right column). The thin curves show mass-weighted profiles at different epochs. The black and red thick lines represent the median of the mass-weighted and emission-weighted profiles, respectively. The profiles integrated over the whole cluster are in agreement with observed CC clusters: the positive gradient of the temperature profiles is preserved, and the entropy profiles follow a baseline power-law with a floor (e.g., Cavagnolo et al. 2009; Lakhchaura et al. 2016, but see also Panagoulia et al. 2014). The profiles (both mass-weighted and emission-weighted) look very similar for the entire cluster and the ambient region, implying that the mass and emission are dominated by the denser ambient ICM. The influence of the AGN is more prominent within the jet cones, as expected. When the AGN is active (e.g., t ∼ 1 and 3 Gyr; see top and bottom rows in Figure 2), the gas density is suppressed, and the temperature and entropy are raised. During the quiescent state (e.g., t ∼ 2.5 Gyr; middle row of Figure 2), the profiles are close to the baseline.

Note that even in the jet cones where most of the heating appears to occur, the emission-weighted profiles do not show significant temperature and entropy inversions (except the very inner ∼10 kpc where the effects of the jets at this radius are likely overestimated due to limited resolution). That is because the gas with very high temperature (see second columns in Figure 2) is mainly contained within the underdense bubbles, and is thus subdominant in terms of mass. The n_e² dependence of the X-ray emissivity further attenuates the contribution from the hot gas in terms of emission. Line-of-sight projections also act to bury the signature of the hot gas. Figure 4 shows a projected, emission-weighted temperature map of the same image in the top, second panel in Figure 2. As can be clearly seen, the projected locations of the bubbles only have slightly higher temperature than their surroundings, instead of a significant trace of heating as one would naively expect. This is in good agreement with observed temperature distributions of well-studied cavities (Fabian et al. 2006; Blanton et al. 2009).

**Figure 4.** Zoomed-in image of projected, emission-weighted temperature map at t = 0.935 Gyr (to be compared with the second panel in the top row in Figure 2).
Download figure:
Standard image High-resolution image

We have demonstrated that kinetic AGN feedback indeed could self-regulate and yield gross ICM properties that are consistent with observations. The quasi-equlibrium state of the cluster over gigayear timescales implies a rough balance between radiative cooling and AGN heating within the cluster core as a whole. In the following sections, we will look into the details of how this balance is established.

3.2. The Lagrangian View

In this section, we show results from analyzing the thermodynamic properties of the ICM tracer particles. In essence, these tracers evolve following the hydrodynamic equations in the Lagrangian form. Particularly, the evolution of entropy of the tracers directly probes the balance between heating and cooling, i.e.,

$\begin{eqnarray}&&\rho T\displaystyle \frac{{ds}}{{dt}}={ \mathcal H }(l)-{ \mathcal C }(l),\end{eqnarray} \tag{ 8 }$

where s is the entropy, and ${ \mathcal H }(l)$ and ${ \mathcal C }(l)$ are the heating and cooling terms along trajectories of the tracer particles. The entropy is important because it tells us whether and where new sources of energy are generated (as opposed to being transported, see the next section). We will show in the following that there is net heating within the jet cones, while cooling dominates over heating for the ambient region. Also, the primary heating mechanisms within the jet cones are shocks as well as mixing between the ICM and the jet materials.

Figure 5 shows the evolution of the projected locations of the ICM tracers. The colors represent the jet mass fraction, f_jet, at the locations of the particles. At t = 0, the particles are uniformly distributed within 100 kpc and have f_jet = 0. Before the AGN is triggered at t ∼ 0.3 Gyr, the cluster contracts due to radiative cooling and thus the sphere containing the particles shrinks radially. Later on, the AGN jets break the symmetry. Particles tracing the ambient region (outside the jet cones, or having low f_jet) keep contracting, similar to a cooling flow. However, note that the initial cooling time is 3 Gyr at r = 100 kpc. Therefore, the fact that the sphere bounding the ambient-gas particle tracers has a radius of r ∼ 40 kpc at t = 3 Gyr instead of sinking to the cluster center implies that the gas inflow is suppressed compared to a pure cooling flow. As we will show below (and also in Section 3.4), though cooling is dominant within the ambient region, the cooling rate is partially compensated by heating from weak shocks.

On the other hand, tracer particles that start from the jet cones or from smaller radii in the ambient region (those have enough time to reach the center by 3 Gyr), are uplifted to larger distances by the rising bubbles. As can be seen from Figure 5, the particle distributions near the jet cones become broader over time. This is because of the combination of a few effects. First, particles that were lifted can slip toward the side of the bubbles and flow inward radially to fill the vacuum behind the cavities. As they flow in, they can be pushed sideways by elliptical shocks generated by subsequent AGN activity. Turbulent diffusion also helps to spread the particle distribution. Note that some of these particles eventually leave the jet cones and become parts of the ambient medium. If the simulation were run longer, they would cool to the center and be lifted by the bubbles again. In other words, there is a gentle circulation within the whole cluster. In Section 3.4, we will show that this circulation plays a crucial role for the transportation of heat.

Next, we study in more detail the thermodynamic histories of the ICM tracers. Figure 6 shows the evolution of all relevant hydrodynamic variables along the trajectories of eight particles that have initial positions within the region [−2 kpc ≤ x ≤ 2 kpc, 8 kpc ≤ y ≤ 12 kpc, 8 kpc ≤ z ≤ 12 kpc]. Hereafter, we denote them as TR10 since these particles are initially located ∼10 kpc from the center right above the equatorial plane. Consistent with the picture illustrated by Figure 5, for the first t ∼ 0.45 Gyr before they mix with the jet materials, TR10 cool and flow in because cooling dominates over heating within the ambient region (i.e., the entropy $K\equiv T/{n}_{{\rm{e}}}^{2/3}$ decreases). A few weak shocks (i.e., jumps in $\rho ,T,K,{e}_{{\rm{int}}},{e}_{{\rm{kin}}}$ , and ΔP/P) heat the gas and raise the entropy, which slows down but does not completely stop the cooling flow. After TR10 start to enter the jet cones (f_jet ≳ 0.01), the evolution becomes more complicated because several processes operate simultaneously. First, the entropy of the traced gas increases whenever f_jet increases (as can be seen clearly from the light blue curve), providing evidence for heating from mixing between the ICM and the jet materials. Second, weak shocks generated by repetitive AGN outbursts also add entropy to the gas. Though each shock only increases the entropy slightly, the shocks occur very frequently (see also Figure 8). Lastly, as the gas is lifted away from the center by buoyantly rising bubbles, its density and internal energies decrease due to adiabatic expansion. All the above effects combined result in the evolution of the temperature as seen in Figure 6. Note that as long as the particles stay within the jet cones (all but the light blue curve during t ∼ 1–2 Gyr), their entropy increases, meaning that on average heating dominates over cooling within the jet cones.

**Figure 6.** Evolution of all relevant variables of eight ICM tracers that are initially located ∼10 kpc from the cluster center right above the equatorial plane (TR10; see the text for a definition). From top to bottom, the panels show the particle coordinates (the r and θ components), the jet mass fraction (f_jet), gas density (ρ), temperature (T), entropy ( $K\equiv T/{n}_{{\rm{e}}}^{2/3}$ ), internal energy (e_int), kinetic energy (e_kin), and the pressure contrast (ΔP/P) as a diagnostic for shocks (Equation (5)).
Download figure:
Standard image High-resolution image

Figure 7 shows the same plot for tracers TR30 that are initially located within the region [−2 kpc ≤ x ≤ 2 kpc, 28 kpc ≤ y ≤ 32 kpc, 8 kpc ≤ z ≤ 12 kpc]. The general trends are similar to TR10, except that they start to mix with the jet fluid at a later time, t ≳ 0.5 Gyr. Again, we find that cooling dominates over heating during the initial period when TR30 are in the ambient region. After they enter the jet cones where heating prevails over cooling, their entropy increases with time almost monotonically. The heating within the jet cones again are done by mixing and shocks. We see significant entropy jumps associated with jumps in f_jet (e.g., red curve at t ≳ 0.9 Gyr, brown and navy lines at t ∼ 1.15 Gyr). The shocks are weaker (ΔP/P < 1.5 or ${ \mathcal M }\lesssim 1.5$ ), but also contribute to heating (e.g., red curve at t ∼ 0.7–0.8 Gyr, where f_jet stays roughly constant but there are small jumps in entropy coincident with shocks; see also Figure 8).

**Figure 8.** Evolution of one of the TR10 tracers (light blue curve in Figure 6). The panels from top to bottom show jumps in f_jet (which traces mixing between the ICM and the jet materials), ΔP/P (which diagnoses shocks) and jumps in entropy, respectively. Cyan and yellow vertical lines mark the time when there is more significant mixing and shock heating, respectively. Both mixing and shocks are important sources of ICM heating within the jet cones.
Download figure:
Standard image High-resolution image

Note that the kinetic energy throughout the cluster evolution remains only at the percent level compared to the internal energy (see Figures 6 and 7). This is consistent with a recent work by Reynolds et al. (2015) who have found that driving of turbulence by AGN in the ambient gas is inefficient. Even if we assume all the kinetic energy goes into turbulent heating (in reality the kinetic energy is also contributed by bulk motions), it is not enough to provide substantial heat to the cluster gas. Therefore, turbulent heating is not the primary source of ICM heating in our simulation. We will discuss this result in more detail in Section 4.1.

We now have a closer look at heating from mixing and shocks within the jet cones, as some of the entropy jumps may not be apparent in Figures 6 and 7. In order to see clear associations among entropy jumps, mixing and shocks, we plotted in Figure 8 the difference in f_jet, ΔP/P and the difference in entropy for one of TR10, namely, the light blue curve in Figure 6. The difference is done between two consecutive output files (i.e., every Δt = 0.005 Gyr). Jumps in f_jet greater than 0.005 are marked in cyan; shocks with ΔP/P > 0.2 or ${ \mathcal M }\gtrsim 1.08$ are painted yellow (these thresholds, though arbitrary, are chosen since they appear to cover most of the visible jumps). This figure shows clearly that both mixing and shocks individually provide entropy to the gas (i.e., entropy jumps that are marked exclusively yellow or cyan). Nevertheless, large entropy jumps preferentially occur when both processes work together (i.e., entropy jumps marked with a blended green color).

Note also that shock heating occurs much more frequently than mixing. Therefore, although heating by each shock is less efficient than each mixing event, as also found by Hillel & Soker (2016), it has been unclear whether the cumulative amount of heating from shocks over the course of the cluster evolution is significant or not. Our simulation is suitable for answering this question because we follow the cluster evolution over billions of years. Furthermore, the strengths and frequencies of shocks are determined by the self-regulated AGN feedback. In Section 3.4, we will present a more quantitative analysis on the relative amounts of heating from different processes.

3.3. The Eulerian View

In this section, we investigate the sources of heating from the Eulerian point of view, as if there were static ICM tracer particles. We will first show that the thermodynamic quantities of the gas at any fixed point in space are on average nearly constant. Together, with results derived from the previous section, we argue that advection of heat as well as adiabatic heating/cooling are important within the ambient region/jet cones.

To study the gas within the ambient region, we plot the relevant hydrodynamic variables at fixed points along the $+y$ -axis in Figure 9. In general, all quantities are remarkably steady, with only small fluctuations caused by weak shocks and sound waves. The only exception is the static tracer at r = 10 kpc (the navy curve). Since it is very close to the cluster center where cooling time is short and also where jets are launched, greater variations are expected. Nevertheless, aside from the large fluctuations caused by mixing (e.g., t ∼ 0.55, 1.90, 2.60 Gyr), the evolution is nearly constant after t ∼ 0.7 Gyr (see also the right column of Figure 3).

**Figure 9.** Evolution of hydrodynamic variables at fixed points in space along the $+y$ -axis. Meanings of variables are the same as in Figure 6.
Download figure:
Standard image High-resolution image

The quasi-steady nature is also found for gas within the jet cones. Figure 10 shows the hydrodynamic variables at fixed points along a sightline 20° away from the $+z$ -axis. Since these locations undergo more events of mixing and stronger shocks, they have greater fluctuations than the ambient gas (see Figure 3). However, except for the excursions caused by mixing and shocks, the gas density, temperature, entropy, and especially internal energy remain essentially unchanged over gigayear timescales.

**Figure 10.** Same as Figure 9, but at fixed points along a line of sight 20° away from the $+z$ -axis. Meanings of variables are the same as in Figure 6.
Download figure:
Standard image High-resolution image

How do we reconcile the nearly constant temperature or internal energy with the fact that there is net heating/cooling within the jet cones/ambient region? We can understand this by recalling the internal energy equation:

$\begin{eqnarray}&&\displaystyle \frac{\partial {e}_{{\rm{i}}}}{\partial t}+{\rm{\nabla }}\cdot [({e}_{{\rm{i}}}+P){\boldsymbol{v}}]-{\boldsymbol{v}}\cdot {\rm{\nabla }}P={ \mathcal H }({\boldsymbol{x}})-{ \mathcal C }({\boldsymbol{x}}),\end{eqnarray} \tag{ 9 }$

where ${ \mathcal H }({\boldsymbol{x}})$ and ${ \mathcal C }({\boldsymbol{x}})$ are functions of space. The ${\rm{\nabla }}\cdot ({e}_{{\rm{i}}}{\boldsymbol{v}})$ term represents transport of the internal energy by the flow via advection or convection. The part of the equation, ${\rm{\nabla }}\cdot (P{\boldsymbol{v}})-{\boldsymbol{v}}\cdot {\rm{\nabla }}P=P({\rm{\nabla }}\cdot {\boldsymbol{v}})$ is the change of internal energy due to adiabatic compression and expansion. From the previous section, we have learned that the right-hand side of the equation is positive/negative within the jet cones/ambient region. In order to make the internal energy roughly constant over time, the net heating/cooling must be balanced by the energy transport and adiabatic energy losses/gains. In the next section, we quantify each term and compare their relative importance.

3.4. Relative Contribution from Each Process

In the previous sections, we identified the processes that are relevant for maintaining the quasi-equilibrium state of the cluster, including radiative cooling, bubble mixing, weak shocks, advective and convective transport, and adiabatic compression and expansion. In the following, we quantify all of these processes and compare their relative importance by computing the volumetric heating and cooling rates.

In order to construct radial profiles of the heating and cooling rates, we divide the domain within 100 kpc into 20 sectors. First, we distinguish the jet cones and the ambient region, and then we further divide them into 10 radial bins, i.e., the ith bin contains region with (i − 1) × 10 kpc ≤ r ≤ i × 10 kpc, where $i\in [1,10]$ . For each process, we compute the heating or cooling "luminosities" (in units of erg s⁻¹) in these sectors. The luminosity of a process is the volume integral of the corresponding term in Equation (9). More specifically, the luminosities due to cooling, advection/convection, adiabatic compression/expansion, and weak shocks are, respectively,

$\begin{eqnarray}{L}_{{\rm{c}}} & = & -\,{\displaystyle \int }_{{V}_{i}}{n}_{{\rm{e}}}^{2}{\rm{\Lambda }}(T){dV},\\ {L}_{{\rm{[adv,}}{\rm{conv]}}} & = & -\,{\displaystyle \int }_{{V}_{i}}{\rm{\nabla }}\cdot {{\boldsymbol{Q}}}_{{\rm{[adv,}}{\rm{conv]}}}{dV}=-\,{\oint }_{{S}_{i}}{Q}_{{\rm{[adv,conv]}}}{dS},\\ {L}_{{\rm{ad}}} & = & -\,{\oint }_{{S}_{i}}P{\boldsymbol{v}}\cdot d{\boldsymbol{S}}+{\displaystyle \int }_{{V}_{i}}{\boldsymbol{v}}\cdot {\rm{\nabla }}{PdV},\\ {L}_{{\rm{sh}}} & = & {\displaystyle \int }_{{V}_{i}}\displaystyle \frac{{dH}}{{dt}}{dV},\end{eqnarray} \tag{ 10 }$

where V_i and S_i are the volume and surface of each sector and ${dH}/{dt}$ is the heating rate by weak shocks (see end of Section 2). The advective and convective fluxes are defined as

$\begin{eqnarray}{Q}_{{\rm{adv}}} & = & \displaystyle \frac{1}{\gamma -1\,}{k}_{{\rm{B}}}(\langle n\rangle \langle T\rangle \langle {v}_{\perp }\rangle +\langle T\rangle \langle \delta n\delta {v}_{\perp }\rangle ),\\ {Q}_{{\rm{conv}}} & = & \displaystyle \frac{1}{\gamma -1\,}{k}_{{\rm{B}}}(\langle n\rangle \langle \delta {v}_{\perp }\delta T\rangle +\langle {v}_{\perp }\rangle \langle \delta n\delta T\rangle \\ & & +\langle \delta n\delta T\delta {v}_{\perp }\rangle ),\end{eqnarray} \tag{ 11 }$

where ${v}_{\perp }$ is the velocity component perpendicular the surface of integration, angle brackets denote averages over the surfaces, and δ refers to the local deviation from the mean value of a quantity, e.g., $\delta n\equiv n-\langle n\rangle$ . The volumetric heating or cooling rates, or "emissivities" (in units of erg s⁻¹ cm⁻³), are then calculated by dividing the luminosities by the volume of each sector, V_i. The luminosity and emissivity profiles for both the jet cones and the ambient region are obtained in this fashion for every simulation output. Note that all the relevant processes could be directly estimated in the way described above, except for heating from mixing. Since the process of mixing is a complex fluid dynamic problem and no simple models exist for our current application, we will infer the amount of mixing heating from the fact that there is rough balance between heating and cooling within the jet cones.

Figure 11 shows the time-averaged profiles of heating and cooling emissivities for the jet cones (left column) and the ambient region (right column). The profiles are averaged between t = 1 and 3 Gyr in order to study the heating and cooling balance when the cluster is in the quasi-equilibrium state. Within the jet cones, we find that because the gas is relatively underdense, radiative cooling is not the dominant cooling process. Instead, cooling due to transport and adiabatic processes is more significant. By breaking down the contributions further (see the bottom left panel), we find that the energy losses are primarily due to adiabatic expansion, while energy transport by advection and convection is subdominant. In terms of heating, the amount of heat provided by repetitive weak shocks is enough to offset the radiative losses, but not the transport plus adiabatic losses (except for larger radii, i.e., r ≳ 80 kpc). In order to maintain the ICM in quasi-equilibrium (see Figure 10), heating from bubble mixing is needed to balance the substantial energy losses from transport and adiabatic processes. We derived such an inferred profile for bubble mixing and plotted it in the upper left panel. As can be seen from the figure, bubble mixing appears to be a more important heating source than weak shocks within the jet cones, especially within the inner 80 kpc. This is consistent with the expectation that mixing is more efficient than shocks (see Section 3.2). However, here we quantitatively show that the cumulative amount of heat provided by mixing is also greater than weak shocks despite of their high frequency. Note that the emissivities shown here does not exclude the bubbles. However, we repeated the analyses by computing all quantities in Equations (10) and (11) weighted by f_gas and the conclusions are unchanged. This implies that adiabatic expansion is associated with not only the bubbles, but also the uplifted ICM as they move outward. This is consistent with the density evolution of the ICM tracer particles as seen in Figures 6 and 7.

The emissivity components for the ambient region are plotted in the right columns of Figure 11. As have been alluded to in the previous sections, there is net cooling within the ambient region, i.e., the only available heating source is weak shocks, and it is insufficient to overcome radiative cooling. What maintains the ambient gas at a given radius to be quasi-steady (see Figure 9) is energy from transport and adiabatic processes. As shown in the bottom right panel, adiabatic compression and advection are both non-negligible (the former has slightly more contribution), while convection plays a lesser role. By estimating the amount of advective energy transfer through the surfaces of each sector, we find that the energy primarily flows from the interface between the jet cones and the ambient region, as opposed to radial inflow from larger radii. In other words, the ambient region, though having net cooling, is pumped by flows of energy from already-heated gas. The gas further gains energy by adiabatic compression as it flows inward radially as a part of gentle circulation. Note that for the first two radial bins in the upper right panel of Figure 11, there appears to be some net energy losses. This loss of internal energy could be compensated by mixing with the non-negligible amount of jet materials within the inner 20 kpc in the ambient region (see the right columns of Figure 2).

In terms of the overall contributions by different processes within the whole cluster core (i.e., the averaged heating and cooling luminosities integrated within 100 kpc), we find that within the jet cones, L_c:L_mix:(L_adv + L_conv + L_ad):L_sh =(−1.52 × 10⁴⁴):9.14 × 10⁴⁴:(−1.37 × 10⁴⁵):4.08 × 10⁴⁴ ≃ (−0.11):0.67:(−1.0):0.30. That is, cooling is dominated by adiabatic expansion instead of radiative cooling; heating from bubble mixing contributes about twice the amount of heating from shocks within the cooling radius. For the ambient region, L_c:(L_adv + L_conv + L_ad):L_sh = (−1.20 × 10⁴⁵):9.95 × 10⁴⁴:2.05 × 10⁴⁴ ≃ (−1.0):0.83:0.17. In other words, shock heating slows down radiative cooling by ∼17%, while the rest of the gas internal energy comes from advection and adiabatic compression. The overall energetics suggest that in our simulation, most injected energy from the AGN is stored within the bubbles, whereas the energy associated with shocks contain a minor portion. This is consistent with the energy partitions inferred by recent analysis of the Perseus cluster (Zhuravleva et al. 2016).

4. DISCUSSION

4.1. Turbulent Heating

In Section 3.2, we showed that the kinetic energy associated with the tracer particles is always only at the percent level compared to the internal energy (Figures 6 and 7). The kinetic energy contains contributions not only from turbulence but also from shocks and waves. To better understand the energetics, we decompose the velocity fields into the compressible (which traces shocks and waves) and incompressible (which measures turbulence and g-modes) components (Reynolds et al. 2015; Yang & Reynolds 2016). The total kinetic energies within 100 kpc for regions containing the bubbles (f_jet ≥ 0.01) and excluding the bubbles (f_jet < 0.01) are plotted in Figure 12. For comparison, the injected AGN energy multiplied by 0.01 is overplotted.

**Figure 12.** Evolution of the kinetic energies of the compressible (dotted) and incompressible (solid) velocity fields, which are further separated by regions containing the bubbles (f_jet ≥ 0.01; blue) and excluding the bubbles (f_jet < 0.01; red). Dashed and dashed–dotted lines show 1% of the injected energy from the AGN and the total gas thermal energy within 100 kpc, respectively.
Download figure:
Standard image High-resolution image

We find that the turbulent energy contained within the bubbles is the largest among all components, ∼2%–3% of the gas thermal energy. However, all kinetic energies are only about 1% of the gas thermal energy or the injected feedback energy of the AGN, suggesting that the driving of turbulence by AGN feedback is inefficient. This result is broadly consistent with the recent simulations with more idealized AGN energy injections by Reynolds et al. (2015) and Hillel & Soker (2016). As discussed in Yang & Reynolds (2016), the inefficiency of AGN-driven turbulence is due to the facts that the typical frequency of AGN outbursts is higher than that required to excite g-waves and that the g-modes are likely trapped within a small volume within the CC.

Given the small amount of energy available for turbulent cascade, we conclude that turbulent heating is not the primary source of heating in our simulated cluster. This is in tension with recent analyses of the observational data for the Perseus cluster (Zhuravleva et al. 2014, 2016), which suggests that turbulent energy is ∼10% of the gas thermal energy and the estimated timescale for turbulence dissipation is comparable to what is required to compensate radiative cooling losses. More efforts, both theoretically and observationally, are needed to resolve this issue. It remains a possibility that predicting the level of AGN-driven turbulence requires models beyond pure hydrodynamics (Reynolds et al. 2015). Observationally, a glimpse at the turbulent dynamics of the ICM will be given by the ASTRO-H/Hitomi observation of the core of the Perseus cluster obtained prior to the loss of that mission in March 2016 (Takahashi et al. 2016). A more complete observational characterization of ICM turbulence must wait for a future X-ray micro-calorimeter mission. If future works prove both statements true (i.e., observed turbulence within cluster cores is energetically significant but AGN driving is inefficient), then it would imply that turbulence within the cores originates primarily from other mechanisms such as sloshing (e.g., ZuHone et al. 2013), galaxy motions (e.g., Ruszkowski & Oh 2011), and cascades from large-scale turbulence generated by mergers (e.g., Vazza et al. 2013). Note that these are not "feedback" mechanisms as they are associated with the dynamical state of the whole cluster, and therefore it may still be challenging to solve the cooling-flow problem via turbulent heating. Alternatively, it might be that AGN feedback is still the solution to cooling catastrophes (e.g., by mixing and shocks), but the dissipation of turbulence is somehow significantly suppressed due to poorly understood plasma physics on microscopic scales. For example, the fact that dissipation of magnetohydrodynamic turbulence is highly localized both spatially and temporally requires heat to diffuse faster than being radiated away (Zhdankin et al. 2015). In the regime where the ICM is weakly collisional, the Braginskii viscosity only dissipates certain components of the velocity gradients (Schekochihin & Cowley 2007), and thus the available energy for turbulence dissipation depends on the partitions between different (compressible versus incompressible) modes (Miniati 2015). To address these questions is beyond the scope of this paper and will be part of our future work.

4.2. Limitations of Current Simulation

In this work, we intentionally focused on the purely hydrodynamic effects of AGN heating, while magnetic fields and plasma transport processes of the ICM are omitted. The details of mixing between the bubbles and the ambient ICM may be affected by the magnetic fields (e.g., Robinson et al. 2004; Ruszkowski et al. 2008) and viscosity (Reynolds et al. 2005). Also, some of the potential heating mechanisms are not explored in the current study. For example, while sound waves in our simulation simply propagate out of the core, their energy can be dissipated within a conducting and/or viscous medium (e.g., Brüggen et al. 2005; Fabian et al. 2005). Cavity heating (e.g., Churazov et al. 2001) is also not included in our analysis, since the gravitational energy gained by the uplifted gas is transformed into kinetic energy when it refills the wakes of the bubbles. In principle, this kinetic energy could be viscously dissipated into heat, but we do not explicitly include viscosity in the current simulations. These additional mechanisms will change the relative partitions among different heating mechanisms. We will incorporate these processes in our future simulations in order to quantify their importance. Note, however, that both cavity heating and sound-wave heating come from the kinetic energy of the gas (within the wakes and sound waves, respectively), which appears to be subdominant in terms of the total feedback energy (see Figure 12). Therefore, the current study suggests that these two mechanisms may not be energetically dominant.

The current simulations have neglected externally driven turbulence in order to probe perturbations solely due to the AGN. In reality, turbulence generated by other sources (e.g., Ruszkowski & Oh 2011; Vazza et al. 2013; ZuHone et al. 2013) could co-exist with the AGN. The pre-existing turbulent motions are expected to displace the buoyantly rising bubbles (e.g., Morsony et al. 2010) so that the nearly isotropic distributions of the observed bubbles can be better reproduced (Fabian et al. 2006). In this case, the heating due to mixing between the bubbles and the ICM would not be confined within the jet cones, but would be distributed more isotropically. Therefore, instead of large-scale circulation consisting of a meridian outflow and an equatorial inflow (Soker et al. 2015; Hillel & Soker 2016), the ICM within the CCs may resemble a pot of boiling water, with outwardly rising gas parcels directly heated by the AGN and under-heated gas inter-penetrating between bubbles and sinking toward the center (Lim et al. 2008). In any case, our current simulation has shown the success of AGN feedback even when the volume of heating is limited. Including external turbulence should only aid the AGN by the additional transport of heat.

Our simulation, as well as recent simulations of AGN feedback that aim to reproduce the properties of the CCs and TI (e.g., Gaspari et al. 2011; Li & Bryan 2014; Prasad et al. 2015), have been focusing on momentum-driven feedback. Since the kinetic energy of the jets is quickly thermalized by shocks, one implicit assumption of these simulations is that the bubbles are filled with ultra-hot (but non-relativistic) gas. Therefore, heating the ICM by mixing is relatively easy since the bubbles contain mainly thermal energy. In reality, the composition of the radio bubbles is still uncertain. Future high-resolution observations of the Sunyaev–Zel'dovich effect of cavities could potentially constrain the bubble composition (Pfrommer et al. 2005). If relativistic particles dominate the bubble pressure, i.e., if the bubbles are inflated by cosmic-ray dominated jets (Guo & Mathews 2011), heating by mixing would likely be less efficient (Mathews & Brighenti 2008). The details of this scenario, e.g., the mechanisms of heating, the interaction between cosmic-ray filled bubbles and the surrounding ICM, and how these processes depend on plasma microphysics, require further investigation. We will address some of these issues in our forthcoming papers.

5. CONCLUSIONS

The key to solving the cooling-flow problem of galaxy clusters is most likely linked to AGN feedback. However, how exactly the feedback energy is transformed into heat, in particular, isotropically, remains unclear. In this work, we use a representative hydrodynamic simulation of a CC cluster with self-regulated, momentum-driven AGN feedback to gain insights into the relative contributions from various physical mechanisms, including radiative cooling, shock heating, bubble mixing, fluid motions (advection, convection, and turbulence), and adiabatic processes.

The picture of AGN feedback drawn from our simulation is summarized as follows. The AGN jets inflate bubbles and produce weak shocks, which increase the entropy of the ICM by shocks and mixing with the hot bubble materials. The heat is primarily dumped within the jet cones (≲30° away from the z-axis); however, the Eulerian internal energy of the gas is not dramatically increased because the gas adiabatically expands as it is uplifted by the buoyantly rising bubbles. A gentle circulation over multiple billions of years is established by the episodic AGN activity, where there is net gas outflow within the jet cones, and some of the gas can refill the ambient region and eventually move inward toward the cluster center. Within the ambient region, heating from weak shocks slows down the cooling losses but not completely, and hence forms a "reduced" cooling flow. However, the thermodynamic properties of the ambient ICM stay nearly constant over gigayear timescales because the ambient region is replenished by previously heated gas from the jet cones, and also because the gas is adiabatically compressed as it moves inward.

In terms of X-ray emission, the dominating component is from the reduced cooling flow, whereas the heated gas mixed with the bubbles is relatively underdense and less luminous. This allows the kinetic AGN feedback model to explain the positive temperature gradients of observed CCs without producing signatures of over-heated gas in the vicinity of bubbles. This aspect of our model is qualitatively similar to the one-dimensional circulation model proposed by Mathews et al. (2004), though the mechanisms of heating were not specified in their work. We also confirmed that AGN heating occurs inside-out (Figure 11), which is an essential ingredient to explain the slope of the soft X-ray spectrum in the cold accretion plus momentum-driven feedback scenario (Gaspari 2015).

Among the sources of heat that we explored, we find that both bubble mixing and shock heating are important, especially within the jet cones. Moreover, the amount of heat provided by mixing is greater than by shocks. There has been a long debate in the community about whether bubble mixing or shock heating is the dominant heat source. On one hand, observations that measure the shock strengths and frequencies suggest enough heating from weak shocks to counteract radiative cooling (e.g., Randall et al. 2015). On the other hand, simulations using tracer fluids found more efficient heating by bubble mixing than shocks (e.g., Hillel & Soker 2016). Interestingly, we find that the two seemingly contradictory statements are both true. Indeed, the cumulative amount of heating from repetitive shocks turns out to be comparable to the radiative losses. However, greater energy losses are associated with adiabatic expansion of the gas, which is invisible from the observed X-ray images. These adiabatic cooling losses have to be compensated by mixing heating.

For the simulated cluster, the kinetic energy of the incompressible velocity field (which mainly traces turbulence) is only at the percent level compared to the injected feedback energy (Figure 12), suggesting that turbulent heating is not the dominant source of heating in our simulation. The implications of this result are discussed in Section 4.1.

One of the biggest puzzles of AGN feedback in CCs is, given that AGN jets are directional, how to heat the whole cluster core isotropically. Specifically, previous works have been looking for sources of heating that can balance radiative cooling at nearly every location within the cluster core. Such a requirement poses great challenges to all the proposed heating mechanisms since none of them can easily provide heating with radial and azimuthal profiles similar to that of radiative cooling. This difficulty can be overcome by taking into account energy transport and adiabatic processes, which are natural consequences of buoyantly rising bubbles. The results of our simulation suggest that the AGN jets do not need to deposit energy isotropically. Instead, they provide excess entropy to parts of the core (i.e., the jet cones in our setup) by bubble mixing and weak shocks. The excess entropy is spent both on heating and adiabatic expansion of the gas. Through the process of gentle circulation, the heated gas brings energy into the under-heated parts of the core. This advection, together with adiabatic compression, provide energy that balances the net entropy losses for the under-heated gas. In other words, the fluid dynamics self-adjusts such that it compensates and transports the heat provided by the AGN, and thus no fine-tuning of the heating profiles is needed.

Our model suggests that self-regulated, momentum-driven AGN feedback, averaged over time, does not provide heat that balances radiative cooling at every location within the cluster core, i.e., the ICM is not in a global thermal equilibrium. Instead, there can be a significant variation in ${ \mathcal H }-{ \mathcal C }$ (of magnitude comparable to ${ \mathcal C }$ ) within the core. This result has important implications for studies of TI within cluster cores, which often assume global thermal equilibrium with small local perturbations (McCourt et al. 2012; Sharma et al. 2012; Meece et al. 2015). In particular, the properties of TI are expected to be different for a background medium that is in global thermal equilibrium versus a (reduced) cooling flow (Sharma et al. 2010). Therefore, the former assumption may not be adequate, but the dynamic nature of the ICM must be taken into account in order to properly model TI and multiphase gas within the CCs (e.g., Gaspari et al. 2012; Li & Bryan 2014).

The current simulation has several limitations, including the assumptions of pure hydrodynamics, the absence of external turbulence, and thermal-pressure dominated bubbles (see Section 4.2 for detailed discussions). Future simulations including processes of increased levels of complexity, such as magnetic fields, anisotropic thermal conduction, Braginskii viscosity, external turbulence, and cosmic rays, will further improve our understanding of the intricate balance between heating and cooling within cluster cores.

The authors thank Yuan Li, Noam Soker, Massimo Gaspari, and an anonymous referee for helpful discussion and comments that helped to improve the manuscript. H.Y.K.Y. acknowledges support by NASA through Einstein Postdoctoral Fellowship grant number PF4-150129 awarded by the Chandra X-ray Center, which is operated by the Smithsonian Astrophysical Observatory for NASA under contract NAS8-03060. C.S.R. is grateful for financial support from the Simons Foundation (through a Simons Fellowship in Theoretical Physics) and the National Science Foundation (under grant AST1333514). The simulations presented in this paper were performed on the Deepthought2 cluster, maintained and supported by the Division of Information Technology at the University of Maryland College Park. FLASH was developed largely by the DOE-supported ASC/Alliances Center for Astrophysical Thermonuclear Flashes at the University of Chicago.

HOW AGN JETS HEAT THE INTRACLUSTER MEDIUM—INSIGHTS FROM HYDRODYNAMIC SIMULATIONS

Article metrics

Permissions

Author e-mails

Author affiliations

Author notes

ORCID iDs

Dates

ABSTRACT

1. INTRODUCTION

2. METHODOLOGY