Reversing cooling flows with AGN jets: shock waves, rarefaction waves, and trailing outflows

TL;DR

This paper uses idealized simulations to explore how AGN jets heat galaxy cluster cores through shock waves and outflows, revealing a cyclical process of cooling and heating that impacts cluster evolution.

Contribution

It provides a detailed physical model of AGN jet interactions with cool-core clusters, highlighting shock, rarefaction, and outflow processes that regulate cooling flows.

Findings

AGN jets drive bow shocks that reverse inflows and heat the core.

Rarefaction waves transport heated energy outward, preventing immediate cooling.

Trailing outflows uplifted by AGN bubbles inhibit gas accumulation, leading to cyclical cooling and heating.

Abstract

The cooling flow problem is one of the central problems in galaxy clusters, and active galactic nucleus (AGN) feedback is considered to play a key role in offsetting cooling. However, how AGN jets heat and suppress cooling flows remains highly debated. Using an idealized simulation of a cool-core cluster, we study the development of central cooling catastrophe and how a subsequent powerful AGN jet event averts cooling flows, with a focus on complex gasdynamical processes involved. We find that the jet drives a bow shock, which reverses cooling inflows and overheats inner cool core regions. The shocked gas moves outward in a rarefaction wave, which rarefies the dense core and adiabatically transports a significant fraction of heated energy to outer regions. As the rarefaction wave propagates away, inflows resume in the cluster core, but a trailing outflow is uplifted by the AGN bubble, preventing gas accumulation and catastrophic cooling in central regions. Inflows and trailing outflows constitute meridional circulations in the cluster core. At later times, trailing outflows fall back to the cluster centre, triggering central cooling catastrophe and potentially a new generation of AGN feedback. We thus envisage a picture of cool cluster cores going through cycles of cooling-induced contraction and AGN-induced expansion. This picture naturally predicts an anti-correlation between the gas fraction (or X-ray luminosity) of cool cores and the central gas entropy, which may be tested by X-ray observations.

Equations

$$\frac{d\rho }{dt}+\rho \nabla \cdot \mathbf{v}=0,$$

$$\rho \frac{d\mathbf{v}}{dt}=-\nabla P-\rho \nabla \Phi ,$$

$$\frac{\partial e}{\partial t}+\nabla \cdot (e\mathbf{v})=-P\nabla \cdot \mathbf{v}-\mathcal{C},$$

$$T=\frac{\mu m_{\mu }P}{k_{\mathrm{B}}\rho }=\frac{\mu P}{{\mu }_{\mathrm{e}}{k}_{\mathrm{B}}{n}_{\mathrm{e}}},$$

$$T(r)=(T_{\mathrm{in}}^5+T_{\mathrm{V}}^5{)}^{1/5},$$

$$T_{\mathrm{V}}(r)=9.68t_{\mathrm{cool}}(r)t_{\mathrm{out}}(r)\text{ keV,}$$

$$t_{\mathrm{cool}}(r)=\frac{0.1+(r/r_1{)}^{1.03}}{1+(r/r_1{)}^{1.03}}\text{ ,}$$

$$t_{\mathrm{out}}(r)=\frac{1}{[1+(r/r_2{)}^a{]}^{b/a}}\text{ ,}$$

$$\Phi ={\Phi }_{\text{DM}}+{\Phi }_{\ast }+{\Phi }_{\text{BH}}\text{.}$$

$${\rho }_{\text{DM}}(r)=\frac{M_0/2\pi }{r(r+r_{\text{s}}{)}^2}\text{,}$$

$${\Phi }_{\text{DM}}=-\frac{2G{M}_0}{r_{\text{s}}}\frac{\text{ln}(1+r/r_{\text{s}})}{r/r_{\text{s}}}\text{,}$$

$${\rho }_{\ast }=\frac{M_{\ast }a}{2\pi r}\frac{1}{(r+a{)}^3}\text{,}$$

$${\Phi }_{\ast }=-\frac{G{M}_{\ast }}{r+a}\text{.}$$

$${\Phi }_{\mathrm{BH}}=-\frac{G{M}_{\mathrm{BH}}}{r-r_{\mathrm{g}}}\text{,}$$

$$P_{\mathrm{jet}}=P_{\mathrm{kin}}+P_{\mathrm{th}},$$

$$Z_0=\{\begin{array}{ll}1& \text{ if }r\le 100\text{ kpc,}\\ 0& \text{ if }r>100\text{ kpc.}\end{array}$$

Full text

††pagerange: 000–000††pubyear: 0000

1 Introduction

Presumably located at knots of the cosmic web, galaxy clusters are the largest gravitationally bound systems in the Universe. The baryon component of galaxy clusters is dominated by the hot, diffuse intracluster medium (ICM) at keV temperatures, which emits copiously in X-rays. One of the central problems in galaxy clusters is the so-called cooling flow problem, which states that gas inflows resulting from radiative cooling with typical mass inflow rates of hundreds to a thousand solar masses per year are inconsistent with current multi-wavelength observations (see Peterson & Fabian 2006 for a review). It is now widely accepted that one or more heating sources are present in galaxy clusters, offsetting radiative cooling and suppressing the development of cooling flows.

One of the key physical mechanisms towards solving the cooling flow problem is radio-mode AGN feedback, which operates as radio-emitting jets apparently emanating from central supermassive black holes (SMBHs) penetrate through and heat the ICM (McNamara & Nulsen 2012; Fabian 2012; Heckman & Best 2014; Soker 2016). The jet-ICM interaction inflates kpc-sized active galactic nucleus bubbles (AGN bubbles) of relativistic plasma, which have been frequently observed through radio synchrotron emissions and in X-ray images of galaxy clusters as surface brightness depressions, i.e., “X-ray cavities” (e.g., Boehringer et al. 1993; Fabian et al. 2002; Bîrzan et al. 2004; Croston et al. 2011). Radio-mode AGN feedback also induces shocks (e.g., Randall et al. 2011; Randall et al. 2015) and sound waves (e.g., Fabian et al. 2003; Fabian et al. 2017) in the ICM. Observations suggest that radio-mode AGN feedback is usually energetically sufficient to offset radiative cooling from cool cores of galaxy clusters (Bîrzan et al. 2004; Rafferty et al. 2006).

The key question then is if radio-mode AGN feedback can effectively suppress cooling flows during at least a major fraction of the lifetimes of galaxy clusters, which has been extensively investigated with hydrodynamic simulations (e.g., Ruszkowski & Begelman 2002; Brighenti & Mathews 2006; Vernaleo & Reynolds 2006; Cattaneo & Teyssier 2007; Guo & Oh 2008; Dubois et al. 2011). Recent simulations which implement momentum-driven AGN jets or outflows self-consistently triggered by cold gas accretion have achieved great success in averting catastrophic cooling while simultaneously keeping observed cool-core profiles in a time-averaged sense (e.g., Gaspari et al. 2011; Gaspari et al. 2012; Li & Bryan 2014; Li et al. 2015; Prasad et al. 2015; Yang & Reynolds 2016b). However, to heat the ICM isotropically, the jets in these simulations are often assumed to be precessing with a typical precession angle of $\sim 10^{\circ}$-$25^{\circ}$. It is unclear if this assumption represents reality, although observational signatures of jet reorientation have been found in some clusters (Babul et al., 2013). Furthermore, AGN jets in clusters are presumably powered by SMBH’s hot accretion flows (Yuan & Narayan, 2014), which have never been directly simulated in cluster simulations due to the extremely large dynamical range involved (see relevant discussions in Gaspari & Sa̧dowski 2017). Therefore, there is no guarantee that the sub-grid treatment of SMBH accretion, a key part of the AGN feedback loop, is correct in current simulations.

A related central question is how AGN jets deliver energy to the ICM and avert cooling flows, which has been actively debated in literature. Proposed channels of AGN heating include the pdV work of expanding AGN bubbles (Ruszkowski & Begelman 2002; Guo et al. 2008), shock heating (Brüggen et al. 2007a; Randall et al. 2015; Li et al. 2016), viscous dissipation of sound waves (Ruszkowski et al. 2004; Fabian et al. 2017), mixing of the ICM with jet plasma (Hillel & Soker 2016; Hillel & Soker 2017), turbulence dissipation (Zhuravleva et al. 2014; Zhuravleva et al. 2016), and cosmic ray heating through the dissipation of self-triggered hydromagnetic waves (Guo & Oh 2008; Jacob & Pfrommer 2017; Ruszkowski et al. 2017). The pdV work does not directly increase the ICM’s entropy, but it increases the ICM’s thermal energy, which is subsequently converted to its kinetic and potential energies during the ICM motions and may eventually dissipate into heat and entropy. It is likely that all of these heating channels take place to some extent in jet-mode AGN feedback, but it has been not easy to perform a clean analysis on their relative importance, which may also depend on jet properties (Tang & Churazov, 2017).

In this paper, we use a specific idealized simulation of a typical cool-core cluster to investigate in detail how a powerful AGN jet event triggered by the central cooling catastrophe averts and even reverse cooling flows. Unlike many recent simulations focusing on the jet’s heating effects and the ICM’s long-term balance between heating and cooling, we particularly focus on complex gasdynamical processes during and after the AGN feedback event, including the bow (forward) shock, rarefaction wave, trailing outflow uplifted by the AGN bubble, meridional circulations in the cluster core. We investigate their dynamical and energetic roles in heating and reversing cooling inflows and in averting catastrophic cooling. Such a detailed analysis on a specific AGN jet event would be complementary to other studies that aim to investigate the long-term ability of AGN feedback in solving the cooling flow problem.

The rest of the paper is structured as follows. We describe our model and the numerical setup in Section 2, and present our results in Section 3. In Section 3.1, we study the overall evolution of the ICM’s temperature, density and entropy profiles before and after the AGN jet event. We then analyze the bow shock and rarefaction wave induced by the AGN event in Section 3.2. We study how the AGN event heats the ICM in our simulation in Section 3.3, and the role of the rarefaction wave in energy redistribution in the ICM in Section 3.4. We then study trailing outflows uplifted by the AGN bubble and cool-core circulations along meridional planes in Section 3.5, and describe an emerging picture of cool cluster cores going through cycles of cooling-induced contraction and AGN-induced expansion in Section 3.6. We finally summarize and discuss our results in Section 4.

2 Methodology

2.1 Equations and Numerical Setup

We study the thermal and hydrodynamic evolution of the ICM gas under the influence of gravity, radiative cooling, and AGN feedback. The basic hydrodynamic equations governing the ICM evolution may be written as

[TABLE]

[TABLE]

[TABLE]

where $d/dt\equiv\partial/\partial t+{\bf v}\cdot\nabla$ is the Lagrangian time derivative, $\Phi$ is the gravitational potential, and $\rho$, ${\bf v}$, $e$, $P$ are the density, velocity, energy density, pressure of the ICM gas respectively. $\mathcal{C}$ is the energy loss rate per unit volume due to radiative cooling. The above hydrodynamic equation set is closed by the relation $P=(\gamma-1)e$, where $\gamma=5/3$.

In this paper, we specifically study the evolution of the ICM in a typical cool-core galaxy cluster subject to cooling and an AGN jet episode. The most basic picture here contains a rotational symmetry with respect to the jet axis. We thus solve Equations (1) - (3) in $(R,z)$ cylindrical coordinates using our own two-dimensional Eulerian code similar to ZEUS 2D (Stone & Norman, 1992). The code has been extensively used in many previous studies, e.g., Guo & Mathews (2010b), Guo & Mathews (2010a), Guo & Mathews (2011), Guo & Mathews (2012), Guo et al. (2012), Guo (2015), Guo (2016). Starting from the origin, the computational grid consists of $800$ equally spaced zones in both coordinates out to $200$ kpc plus additional $400$ logarithmically-spaced zones out to $2$ Mpc. The spatial resolution within the inner 200 kpc is 0.25 kpc, while the large outer boundary ($2$ Mpc) minimizes numerical impacts of the outer boundary conditions on the cool cluster core region. We adopt reflective boundary conditions at the inner boundaries. At the outer boundaries, we adopt reflective boundary conditions as well, which prohibit mass exchange between active and ghost zones (Stone & Norman, 1992), facilitating the analysis of the ICM energy evolution in Section 3.3.

The radiative cooling rate in Equation (3) may be written as $\mathcal{C}=n_{\rm i}n_{\rm e}\Lambda(T,Z)$, where $n_{\rm e}$ is the electron number density, $n_{\rm i}$ is the ion number density, and the cooling function $\Lambda(T,Z)$ is adopted from Sutherland & Dopita (1993) with a fixed metallicity $Z=0.4$. The cooling time often becomes very short when the gas cools to low temperatures. In order to accelerate the computations, we enforce a minimum gas temperature floor $T_{\rm floor}=10^{5}$ K, and employ the sub-cycling method for radiative cooling (Anninos et al. 1997; Proga et al. 2003; Ruszkowski et al. 2017) when the local cooling time becomes shorter than the hydrodynamical timestep. Our simulations do not have enough spatial resolutions and physical processes to correctly follow the evolution of cold gas.

In the collisionally-ionized hot ICM, the ion number density $n_{\rm i}$ may be roughly approximated to be $n_{\rm i}=1.1n_{\rm H}$, where $n_{\rm H}$ is the proton number density. Therefore the mean molecular weight per particle is $\mu=0.61$, and that per electron is $\mu_{\rm e}=5\mu/(2+\mu)=1.17$. According to the ideal gas law, the gas temperature is related to the gas pressure and density via:

[TABLE]

where $k_{\rm B}$ is Boltzmann’s constant and $m_{\mu}$ is the atomic mass unit.

2.2 Cluster Setup

For initial conditions, we construct an idealized cool core galaxy cluster in hydrostatic equilibrium. For concreteness, the cluster profiles are initialized according to the archetypal cool core cluster Abell 1795, which has been well observed by both Chandra and XMM-Newton (Tamura et al., 2001; Ettori et al., 2002; Vikhlinin et al., 2006). The initial cluster setup is very similar to that in Guo & Mathews (2010b) and Guo & Mathews (2014), where the readers are referred to for more details.

We first setup the radial temperature profile of the ICM using an analytic fit to the deprojected $3$-dimensional temperature distribution of A1795 derived from Chandra observations:

[TABLE]

where $T_{\rm in}=2.5$ keV is the observed central temperature of A1795, and $T_{\rm V}$ is the best-fit temperature profile to Chandra data of A1795 proposed by Vikhlinin et al. (2006):

[TABLE]

where

[TABLE]

describes that the temperature declines inward in the inner cool-core region, and

[TABLE]

represents the radially declining outer region. Here the parameters are $r_{1}=77$ kpc, $r_{2}=550$ kpc, $a=1.63$ and $b=0.9$ (Vikhlinin et al., 2006). Our initial gas temperature profile $T(r)$ is shown in the bottom panel of Figure 1, indicating clearly that it provides a reasonably good fit to Chandra X-ray data of A1795 covering from the inner few kpc to $\sim 1$ Mpc.

The initial gas density profile is derived from the assumption of hydrostatic equilibrium and a fixed gravitational potential contributed by three components: the dark matter halo ($\Phi_{\text{DM}}$), the central galaxy ($\Phi_{*}$), and the central SMBH ($\Phi_{\text{BH}}$):

[TABLE]

We adopt a Navarro-Frenk-White (NFW) profile (Navarro et al., 1997) for the dark matter density distribution:

[TABLE]

where $r_{\text{s}}$ is the standard scale radius and $M_{0}$ is a characteristic mass. As in Guo & Mathews (2014), we take $M_{0}=4.2\times 10^{14}M_{\odot}$, and $r_{\text{s}}=430$ kpc for A1795, which correspond to a virial mass $M_{\rm vir}=7.06\times 10^{14}M_{\odot}$ and a viral radius $r_{\rm vir}=1.8$ Mpc at the redshift of A1795 $z=0.0632$. Here $r_{\rm vir}$ is defined as the radius within which the mean dark matter density equals to $200\rho_{\rm crit}(z)$, where $\rho_{\rm crit}(z)=3H(z)^{2}/8\pi G$ is the critical density of the universe at the redshift of the system. The corresponding dark matter gravitational potential is:

[TABLE]

where $G$ is the gravitational constant.

We take the stellar mass density of the central galaxy to have a Hernquist profile (Hernquist, 1990):

[TABLE]

where $M_{*}$ is the total stellar mass and $a=R_{\rm e}/1.8153$ is a scale length. Here $R_{\rm e}$ is the radius of the isophote enclosing half the galaxy’s light. For the Brightest Cluster Galaxy (BCG) in A1795, $M_{*}=5.35\times 10^{11}M_{\odot}$ and $R_{\rm e}=40.3$ kpc (Chandran & Rasera 2007; Guo & Mathews 2014). We note that the adopted value of $R_{\rm e}$ is larger than typical BCG sizes, which may underestimate the BCG’s contribution to the gravitational potential:

[TABLE]

The central SMBH’s gravitational potential may be written as

[TABLE]

where $M_{\rm BH}$ is the SMBH’s mass, $r_{\rm g}=2GM_{\rm BH}/c^{2}$ is the Schwarzschild radius, and the $1/(r-r_{\rm g})$ mimics the effects of general relativity (Paczyńsky & Wiita 1980; Quataert & Narayan 2000). For the central SMBH in A1795, we take $M_{\rm BH}=1.66\times 10^{9}M_{\odot}$ (Chandran & Rasera, 2007). While we include the gravitational acceleration contributed by the central SMBH in our simulations, it is negligibly small compared to that contributed by the dark matter halo and the central galaxy in all computational zones.

The initial cluster profiles described here are spherically symmetric. We derive the initial gas density profile from the initial temperature profile and the gravitational potential, assuming hydrostatic equilibrium. We normalize the density profile by the central electron number density $n_{\rm 0}$ (more specifically, the value of $n_{\rm e}$ at $r=9.8$ pc), which is used as a free parameter to fit the observed radial gas density distribution. In our simulations, the value of $n_{0}$ is chosen to be $n_{0}=0.48$ cm*-3*, which leads to an initial density profile that provides a reasonably good fit to observations, as shown in the top panel of Figure 1.

2.3 The Numerical Experiment and Jet Modeling

With the initial setup described in the previous two subsections, we performed a numerical experiment to study the development of cooling flows and more importantly, how a subsequently-triggered powerful AGN outburst transfers energy to the ICM, averting cooling inflows. AGN heating in galaxy clusters is usually considered as a feedback mechanism, which relates the onset of AGN outbursts directly to the ICM cooling and inflow. Such an AGN feedback mechanism is required to stably maintain cool-core cluster profiles, as shown in the global stability analysis of Guo et al. (2008). Recent state-of-the-art hydrodynamical simulations in this field also favor the self-regulated AGN feedback mechanism, where AGN outbursts are triggered by the so-called “cold-mode accretion” of cold gas dropped out from the hot ICM (Gaspari et al. 2011; Gaspari et al. 2012; Li et al. 2015; Yang & Reynolds 2016a; Yang & Reynolds 2016b; Li et al. 2016). These simulations usually start with a spherically-symmetric hydrostatic ICM, and cold gas often first appears right after the central cooling catastrophe happens.

Our numerical experiment is in line with recent progress. We first let the initially hydrostatic cluster evolve under the influence of radiative cooling and gravity, but without AGN outbursts. Once the central cooling catastrophe happens (after around the central gas cooling time; Guo & Mathews 2014), AGN feedback is manually triggered. To be specific, we turn on an AGN jet in our simulated hemisphere of the cluster once the total mass of cold gas with $T<T_{\rm crit}=5\times 10^{5}$ K within the central $1$ kpc becomes larger than $M_{\rm crit}=10^{-3}M_{\odot}$. Once the cooling catastrophe starts, it develops very quickly, and thus the exact values of $T_{\rm crit}$ and $M_{\rm crit}$ do not affect the results in the paper. For simplicity, when it is active, the jet in our simulation is steady and uniform at the jet base, corresponding to a steady accretion flow around the central SMBH. This is different from the “instantaneous” feedback model adopted in recent simulations (e.g., Yang & Reynolds 2016b; Li et al. 2016), where the jet parameters are determined instantaneously by the amount of central cold gas and a fixed “free-fall” timescale. In reality, before being accreted by the central SMBH, cold gas may circularize around it, forming a torus or accretion flow. The jet properties are determined by complex physical processes taking place in accretion flows. Our simple jet modeling method provides a straightforward way to investigate how a powerful AGN jet episode suppresses and even reverses cooling inflows.

In our simulation, the AGN jet is injected into the ICM at the cluster centre along the $z$-axis with an opening angle of [math] degree. Similar to Brüggen et al. (2007b) and Guo & Mathews (2011), we initialize the jet by applying inflow boundary conditions to a cylindrical nozzle placed in ghost zones, which inject the mass-, momentum-, and thermal energy fluxes into active computational zones. The jet’s kinetic energy is automatically injected into the ICM with the injection of mass and momentum. In the initialization nozzle, the jet is uniform and steady, characterized by $4$ parameters: density $\rho_{\rm jet}$, thermal energy density $e_{\rm jet}$, velocity $v_{\rm jet}$, and radius $R_{\rm jet}$. The jet is turned off after an active duration of $t_{\rm jet}$. In our simulation, the values of these parameters are: $\rho_{\rm jet}=1.61\times 10^{-26}$ g cm*-3* (density contrast $\eta=0.1$ with respect to the central ICM density in the innermost active zone at $t=0$), $e_{\rm jet}=1.91\times 10^{-9}$ erg cm*-3* (energy density contrast $\eta_{\rm e}\equiv e_{\rm jet}/e_{\rm amb}=2$ with respect to the central ICM thermal energy density $e_{amb}$ at $t=0$), $v_{\rm jet}=3.0\times 10^{9}$ cm s*-1*, $R_{\rm jet}=1.5$ kpc ($6$ computational grids along the $R$ direction), and $t_{\rm jet}=5$ Myr. These values are specifically chosen to ensure that the AGN event is powerful enough to reheat the cooling flow back to the initial cluster profile, while driving a bow shock with its aspect ratio (major to minor axis ratio) roughly consistent with observations. From $t=250$ Myr to $300$ Myr as seen in the left panel of Figure 3, the shock’s aspect ratio evolves from $3.0$ to $1.6$. For comparison, the shock’s aspect ratios in two observed powerful AGN outbursts Cygnus A (Wilson et al., 2006) and MS 0735.6+7421 (hereafter MS0735; Gitti et al. 2007) are $2.2$ and $1.5$, respectively (considering the projection effect, the real values may be larger; Mathews & Guo 2010). According to the classification of AGN jets in Guo (2015) and Guo (2016), the jet in our simulation is light, and internally supersonic with an internal Mach number $M_{\rm int}\sim 8.26$.

The jet power in the initialization nozzle may be written as

[TABLE]

where $P_{\rm kin}=\rho_{\rm jet}v_{\rm jet}^{3}\pi R_{\rm jet}^{2}/2$ and $P_{\rm th}=e_{\rm jet}v_{\rm jet}\pi R_{\rm jet}^{2}$ are the kinetic and thermal powers of the jet, respectively. The jet in our simulation is dominated by the kinetic power with $P_{\rm kin}=1.46\times 10^{46}$ erg/s and $P_{\rm th}=3.85\times 10^{44}$ erg/s. The total jet power $P_{\rm jet}=1.5\times 10^{46}$ erg/s is close to that estimated in Cygnus A (Wilson et al., 2006) and MS0735 (McNamara et al., 2009), but less than the maximum jet power $\sim 10^{47}$ erg/s reached in some recent simulations (e.g., Gaspari et al. 2011; Li et al. 2015; Yang & Reynolds 2016b). The total jet energy in our simulation is $E_{\rm jet}=P_{\rm jet}t_{\rm jet}\sim 2.4\times 10^{60}$ erg, which is much less than that in Cygnus A and MS0735. The jet episode in our simulation is powerful and brief, while a less powerful jet with a longer duration would result in a more gentle impact on cluster cores. The dominant channel through which the jet delivers energy to the ICM may be different for these two types of jet activities.

3 Results

In this section, we present the results of our main simulation (thereafter denoted as the “AGN feedback” run) described in detail in Section 2.3. For comparison, we also present some results from a control run, denoted as the “pure cooling” run, where cooling flows are developed unimpeded without turning on AGN feedback throughout the whole run.

3.1 Evolution of ICM Profiles

We first investigate the overall cluster evolution by looking at the spherically-averaged profiles. Figure 2 shows the time evolution of mass-weighted spherically-averaged profiles of gas temperature, density, and specific entropy $S\equiv T/n_{\rm e}^{2/3}$ in the “pure cooling” run (left panels) and our main “AGN feedback” run (right panels). The left panels clearly show the development of cooling flows under the influence of gravity and radiative cooling. The hot ICM initially in hydrostatic equilibrium gradually loses pressure support due to cooling, leading to a gradual decrease in both temperature and entropy, and increase in density. At $t\sim 238$ Myr, a central cooling catastrophe suddenly appears, and the gas temperature in central regions drops very quickly down to the temperature floor in our simulation. After that, cold gas drops out and accumulates in central regions quickly. The mass inflow rate across $r=50$ kpc in this “catastrophic cooling” state of A1795 is about $160$ - $200M_{\odot}/$yr, as seen in Figure 9 and relevant discussions in Section 3.6. This picture on the cooling-flow development in galaxy clusters is consistent with recent high-resolution simulations in Li & Bryan (2012) and Guo & Mathews (2014).

The cluster evolution in our main run is shown in the right panels. During the first $238$ Myr, cooling flows develop gradually in the ICM without AGN feedback, which is the same as in the “pure cooling” run. Once the central cooling catastrophe appears at $t\sim 238$ Myr, an AGN jet activity is manually turned on to mimic AGN feedback triggered by the SMBH’s accretion of cold gas resulting from the central cooling catastrophe. As elaborated in Section 2.3, the jet is specifically set up to be powerful enough to avert cooling flows, while its power is still consistent with powerful AGN outbursts observed in some cool core clusters (e.g., Cygnus A and MS0732).

The right panels of Figure 2 clearly show that the jet stops the central cooling catastrophe, and heats up the cluster core. In the inner regions with $r\lesssim 20$ kpc, the gas is over-heated as the temperature and entropy rise above their initial values at $t=0$. The jet event in our simulation is brief with a duration of $t_{\rm jet}=5$ Myr. After the jet event, the cluster re-adjusts its profiles, and cools later on back to roughly its initial distribution at $t\sim 500$ Myr. In the subsections below, we focus on studying in detail the energetical and dynamical processes as the jet averts and even reverses the cooling inflow.

3.2 Shock Waves and Cool Core Expansion

What happens to the cooling inflow as the jet is injected into the ICM? Figure 3 shows the temporal evolution of the inner cluster region in our simulation. The left and right panels show the gas density and velocity distributions, respectively, at $t=237.6$, $250$, $300$, $400$, and $500$ Myr (from top to bottom). The top panels correspond to the ICM state at $t=237.6$ Myr, right before the onset of the central cooling catastrophe. Clearly the ICM gas in the cooling flow is all moving toward the cluster centre, where the gas density is the highest and the cooling time is the shortest. Which naturally explains why the cooling catastrophe first appears in the cluster centre.

The evolution of AGN jets in an ambient medium has been extensively studied with hydrodynamic simulations in the past decades (e.g., Norman et al. 1982; Reynolds et al. 2001; Ruszkowski et al. 2004; Hardcastle & Krause 2013; Guo 2015). As a supersonic jet penetrates through the ambient ICM, it drives a bow shock that encloses the jet and the lobe resulting from its interaction with the ICM. Since the jet moves much faster than the bow shock, it creates a hotspot at its head (working surface) bracketed by the apex of the bow shock and a reverse shock within the jet near its working surface. The jet’s kinetic energy is quickly dissipated into thermal energy through the reverse shock, and the high pressure in the hot spot drives expansion and the bow shock. The outward-propagating bow shock can be clearly seen as a sharp semi-elliptical discontinuity in the second and third rows of Figure 3, corresponding to the times $t=250$, $300$ Myr, respectively. Along the $R$ direction, the shock’s Mach number drops from $\sim 1.26$ at $t=250$ Myr to $1.09$ at $t=300$ Myr. Right behind the shock front, there is a shell of post-shock gas that is moving away from the cluster centre. This is essentially a rarefaction wave following the bow shock. The locations of the bow shock and rarefaction wave are also shown schematically in the middle two panels of Figure 10 (see Sec. 3.6 for more details).

Here we stress important dynamical impacts of the AGN-driven bow shock and the post-shock rarefaction wave on the cool core and the development of cooling flows therein. As the bow shock propagates outward through the ICM, it reverses cooling inflows into outflows, and the shocked gas in the downstream of the shock forms the outward-propagating rarefaction wave, as clearly seen in the second and third rows of Figure 3. In our particular simulation at $t=250$ or $300$ Myr, the typical outflowing velocity in the downstream of the shock is several hundreds km/s, which is much larger than the inflow speed of several tens km/s in the upstream of the shock. AGN-driven bow shocks have been frequently observed in cool-core clusters (Randall et al. 2011; Randall et al. 2015). The post-shock rarefaction wave causes cool core expansion and rarefies the dense cool core region, thus playing an important role in averting the central cooling catastrophe which otherwise develops in a classical cooling flow seen in our “pure cooling” run.

As the shock front propagates far away, the shocked gas in inner regions changes its velocity direction again and moves inward, as seen in the fourth and fifth rows of Figure 3, corresponding to the times $t=400$ and $500$ Myr, respectively. This is consistent with wave nature of the post-shock rarefaction wave. However, a particularly important kinetic feature within the post-shock cool core is the vertical outflow trailing behind the AGN bubble, as clearly seen in the bottom two rows of Figure 3 (also see Fig. 10). These trailing outflows along the jet direction and inflows in other directions form cool-core circulations in meridional planes, which suppresses accumulation of gas in the cluster centre and the subsequent onset of the cooling catastrophe there. We will further investigate trailing outflows and meridional circulations in more detail in the next subsection.

The AGN-driven expansion of cool cluster cores is real in our simulation, which can also be seen in Figure 4. The top panel shows the time evolution of normalized gas masses within three representative spheres of radii $r=50$, $100$, and $200$ kpc. Before turning on the AGN jet, the total gas mass within each spherical region increases gradually with time, as a result of the cooling inflow. Once the jet is turned on at $t\sim 237.6$ Myr, the bow shock and rarefaction wave are produced, and the cluster expands from inner to outer regions. The gas mass enclosed within the inner $50$ kpc region decreases quickly by $\sim 30\%$ within several tens Myrs, while the total gas masses enclosed within $100$ and $200$ kpc regions start to decrease later as it takes time for the bow shock to propagate outward. The peak outflow rate is several thousands $M_{\odot}$/yr, as seen in the bottom panel, which also clearly shows that after the outflow episode, gas flows in again with negative values of $\dot{M}$. However, even at $t=500$ Myr, which is more than $250$ Myr after the onset of AGN feedback, the gas mass enclosed within the inner $50$ kpc is still less than its value at $t\sim 237.6$ Myr (right before the jet event), indicating that the overall cool core expansion induced by the jet event is significant.

3.3 Energy Evolution and Shock Heating

Here we study the energy evolution and heating mechanisms of the ICM during the gasdynamical processes investigated in the previous subsection. We first investigate the global energy evolution of the ICM in Figure 5, which shows the temporal evolution of the kinetic energy $E_{\rm kin}$, the change in the potential energy $\Delta E_{\rm pot}$, and the change in thermal energy $\Delta E_{\rm th}$ of the whole simulated ICM hemisphere in both the “pure cooling” and AGN feedback runs. In the pure cooling run (dotted lines), the decrease of thermal energy is significant during the simulation, which is obviously due to radiative cooling. The cooling-induced inflow causes the cluster contraction, explaining the evident decrease of potential energy, which is gradually converted to thermal energy and kinetic energy. However, the increase of kinetic energy is negligibly small compared to changes in thermal and potential energies during the whole “pure cooling” run.

In the AGN feedback run, an AGN jet is turned on at $t=237.6$ Myr (the left vertical dashed line), and lasts for $5$ Myr, causing sharp increases of the kinetic and thermal energies, as clearly seen in Figure 5. While the initial jet in the initialization nozzle is dominated by kinetic energy ($P_{\rm kin}/P_{\rm jet}\sim 97\%$), the increase of the system’s thermal energy during the jet injection phase is comparable to that in kinetic energy, due to the rapid dissipation of kinetic energy by shocks within the supersonic jet and the bow shock. The kinetic energy of the system starts to decrease when the jet is turned off at $t=242.6$ Myr (the right vertical dashed line), while the increase of the thermal energy by shock dissipation lasts slightly longer until $t\sim 245$ Myr. The potential energy initially decreases with time during the cooling inflow phase, but once the AGN jet is turned on, it becomes relatively flat with time until $t\sim 320$ Myr. The potential energy increases in the post-shock rarefaction wave, which is roughly compensated by its decrease in inflowing outer regions.

To further study jet heating on the cool core region, we show in the top panel of Figure 6 the energy evolution of the ICM gas elements initially located with distances of $r_{0}<100$ kpc to the cluster centre. As described in detail in Section 3.5, we “paint” these gas elements at $t=0$ with the only non-zero concentrations of a tracer fluid: $Z_{0}\equiv\rho^{\prime}/\rho=1$, where $\rho^{\prime}$ is the density of the tracer fluid. As each gas element moves, the value of concentration $Z$ is conserved, and thus we can study the energy evolution of these gas elements in a Lagrangian way. We note, however, that as the jet moves quickly across numerical grids in our simulation, it inevitably mixes with the ambient ICM gas through numerical diffusion (i.e., advection errors), and RT and KH instabilities at the bubble/ICM interface further facilitates the mixing, which results in some grids with $0<Z<1$ (see the top panel of Fig. 8). In the top panel of Figure 6, we show the temporal evolution of the energies associated with the gas elements with both $Z>0.5$ (solid) and $Z>0.1$ (dashed), and their evolutions are very similar (except for the late evolution of the potential energy), suggesting that our method does a reasonably good job in tracing the Lagrangian energy evolution of the ICM.

Now let us look at the Lagrangian evolution of the kinetic energy, change in the potential energy, and change in thermal energy, of the gas elements with $r_{0}<100$ kpc in the top panel of Figure 6. When $t<237.6$ Myr, both the thermal and potential energies decrease with time due to cooling and the ICM contraction. Once the jet is turned on at $t=237.6$ Myr, both the kinetic and thermal energies increase rapidly. The increase phase of the kinetic energy is very brief, and stops at $t\sim 245$ Myr, when shock dissipation surpasses kinetic motions induced by the jet injection. The vertical dotted line corresponds to $t=250$ Myr, when the bow shock starts to propagate across $r=100$ kpc (see the second row of Fig. 3) and the thermal energy starts to decrease due to adiabatic expansion in the rarefaction wave (also see Sec. 3.4). The peak thermal energy after the jet event is much higher than its initial value at $t=0$, indicating that the cool core is over-heated by the jet event.

How does the jet heat the cool core gas in our simulation? Possible heating mechanisms include mixing, pdV work, turbulence, and shock. The jet-ICM mixing region is located within the low-density AGN bubble (cavity), and the ICM heating due to mixing should be appreciably less than the whole thermal energy contained in the AGN bubble (cavity). At $t=250$ Myr (the second row of Fig. 3), the AGN bubble may be approximated by a cylinder with radius $R_{b}\sim 10$ kpc and height $h_{b}\sim 80$ kpc, and the internal thermal energy density $e_{b}$ varies between $2-5\times 10^{-10}$ erg cm*-3*. Taking $e_{b}\sim 4\times 10^{-10}$ erg cm*-3*, the whole thermal energy within the AGN bubble is $E_{b}\sim e_{b}\pi R_{b}^{2}h_{b}\sim 3\times 10^{59}$ erg, which is only about $1/6$ of the sharp increase of the ICM thermal energy $\Delta E_{\rm th}\sim 1.8\times 10^{60}$ erg from $t=237.6$ Myr to $250$ Myr. Similarly, the pdV work done during the inflation of the AGN bubble is $E_{\rm pdV}\sim P_{b}\pi R_{b}^{2}h_{b}\sim(\gamma-1)E_{b}\sim 2\times 10^{59}$ erg, which is only $\sim 10\%$ of $\Delta E_{\rm th}$. Here $P_{b}\sim(\gamma-1)e_{b}$ is the bubble pressure. We further argue that turbulence heating is unlikely to dominate (also see Reynolds et al. 2015), as the velocity fields in Fig. 3 do not show significant levels of turbulence and no explicit viscosity other than artificial viscosity is included to dissipate turbulence in our simulation. Furthermore, the dissipation of the jet’s kinetic energy is very rapid.

We therefore conclude that the ICM heating in our simulation is mainly due to shock heating, which is consistent with the rather rapid dissipation of the jet’s kinetic energy seen in Figure 5 and the top panel of Figure 6. Figure 7 shows the time evolution of the total energy $E_{\rm shock}$ dissipated by shocks within time $t$ in all active zones with a pressure jump: $P_{\rm i}/P_{\rm i+1}>1.6$ (solid), $1.1$ (dotted), $1.01$ (short-dashed), and $1.001$ (long-dashed), where $P_{\rm i}$ is the pressure in each zone and $P_{\rm i+1}$ is the pressure in the adjacent zone along either the positive $z$ or $R$ direction. As in the ZEUS code (Stone & Norman 1992; also see Li et al. 2016), shock dissipation is dealt with a von Neumann-Richtmyer artificial viscosity and $E_{\rm shock}$ can thus be simply evaluated with the dissipation of kinetic energy by the artificial viscosity in the corresponding zones. We note that $E_{\rm shock}$ includes shock dissipation within both the supersonic jet and ICM, although the reverse shock near the jet terminus ($P_{\rm i}/P_{\rm i+1}<1$) has already been excluded by the above pressure jump conditions. For very low values of $P_{\rm i}/P_{\rm i+1}$, $E_{\rm shock}$ may also include dissipation of complex flows and waves (Ryu et al., 2003). Nonetheless, Figure 7 suggests that shock heating may indeed be significant, and happens mostly during and shortly after the jet event. For example, at $t=250$ Myr when the bow shock only reaches $R\sim 30$ kpc along the $R$ direction, shock dissipation with $P_{\rm i}/P_{\rm i+1}>1.6$ has almost finished, and $\sim 87\%$ of energy dissipation by shocks with $P_{\rm i}/P_{\rm i+1}>1.1$ has been completed. As clearly seen in the second and third rows of Figure 3, the most evident shock in the ICM is the bow shock, which presumably dominates over-heating of the cool core. Efficient shock heating during powerful AGN outbursts has also been found recently by Yang & Reynolds (2016a) and Li et al. (2016).

3.4 Energy Transport by the Rarefaction Wave

As studied in the previous subsection (Fig. 5 and the top panel of Fig. 6), shocks driven by the powerful AGN jet dissipate kinetic energy rapidly and over-heat the cool cluster core in our simulation. This is consistent with earlier studies by Fujita & Suzuki (2005) and Mathews et al. (2006), who found that AGN-driven shocks tend to dissipate kinetic energy rapidly in inner cool core regions. Then, how does the jet transfer its energy to the whole cluster core and even outer cluster regions? Bow shocks are followed by rarefaction waves, which cause the cluster expansion. Here we study energy transport by the rarefaction wave in our simulation.

The top panel of Fig. 6 shows that the thermal energy of the gas elements with $r_{0}<100$ kpc decreases significantly by an amount of $2.2\times 10^{60}$ erg from $t=0.25$ Gyr to $0.32$ Gyr, which is clearly due to adiabatic losses by the core expansion during the rarefaction wave phase. Presumably a significant fraction of the lost energy is transported to outer regions. To confirm this, we follow the energy evolution of the gas elements initially located within $100<r_{0}<200$ kpc at $t=0$ in the middle panel of Fig. 6, which clearly shows that the thermal energy of these gas elements increases significantly from $t=0.25$ Gyr (the left vertical dotted line) to $0.31$ Gyr (the right vertical dotted line) by an amount of $1.3\times 10^{60}$ erg. In the rarefaction wave, some thermal energy is also converted to potential energy, while another fraction is lost to radiation.

At $t\sim 0.31$ Gyr, the shock front reaches $r\sim 200$ kpc, and the thermal energy of the gas elements with $100<r_{0}<200$ kpc decreases significantly from $t=0.31$ Gyr to $0.38$ Gyr by an amount of $3.3\times 10^{60}$ erg due to the gas expansion. Most of the lost energy is transported to the gas elements with $r_{0}<100$ kpc (see the corresponding thermal energy increase during the same epoch in the top panel) and $200<r_{0}<300$ kpc (see the thermal energy increase between two dotted lines in the bottom panel of Fig. 6). At $t\sim 0.38$ Gyr, the thermal energy of the cool core gas with $r_{0}<100$ kpc is close to its initial value at $t=0$, indicating that the net energy gain during the first $380$ Myr through shock heating, adiabatic expansion, and the later contraction roughly balances the energy loss through radiation. As also seen in the bottom panel of Fig. 6, the bow shock and rarefaction wave reach $r\sim 300$ kpc at $t\sim 0.38$ Gyr, and the thermal energy of the gas elements with $200<r_{0}<300$ kpc starts to decrease significantly as it is adiabatically transported to the inner ($100<r_{0}<200$ kpc) and outer ($r_{0}>300$ kpc) regions. We note that, to make the middle and bottom panels, we performed two additional simulations that trace the Lagrangian evolution of the gas elements initially located in the regions $100<r_{0}<200$ kpc and $200<r_{0}<300$ kpc, respectively, by assigning the only nonzero values of $Z_{0}=1$ to each corresponding region.

In summary, the rarefaction wave which follows the bow shock transports thermal energy from the over-heated cluster core to outer regions, and thus plays a key role in the energy redistribution in galaxy clusters during AGN feedback events.

3.5 Trailing Outflows and Meridional Circulations

In this subsection, we study in more detail trailing outflows of AGN bubbles and meridional circulations in cool cluster cores found in Section 3.2. The ICM velocity distributions in the bottom two rows of Figure 3 clearly show coherent gas outflows behind the AGN bubble along the $z$ axis at $t=400$ and $500$ Myr, extending from the cluster centre to the bottom of the AGN bubble ($\sim 150$ kpc in our run). In the present paper, we refer these coherent outflows as “intracluster trailing outflows”, as they are apparently trailing flows of AGN bubbles in their wakes. We note that in our simulation the AGN bubble is significantly disrupted by Rayleigh-Taylor (RT) and Kelvin-Helmholtz (KH) instabilities (Reynolds et al., 2005), unlike some old intact X-ray cavities detected in galaxy clusters (e.g., the outer northwestern X-ray cavity in Perseus; see Figure 3 in Fabian et al. 2000). The suppression of these interface instabilities during the development of AGN bubbles may be due to shear viscosity (Kaiser et al. 2005; Reynolds et al. 2005; Guo et al. 2012; Guo 2015) or magnetic tension (Kaiser et al. 2005; Jones & De Young 2005; Ruszkowski et al. 2007), which is not included in our simulation. However, the formation of trailing outflows should be robust, as also seen in simulations of AGN bubbles in viscous media (see Fig. 8 of Guo 2015).

Do coherent outflowing motions trailing behind the AGN bubble represent real outflows? To answer this question, we rerun our simulation by additionally following the evolution of a tracer fluid $\rho^{\prime}$ which follows gas motions, i.e., ${d\rho^{\prime}}/{dt}+\rho^{\prime}\nabla\cdot{\bf v}=0$. Combining with the gas mass conservation equation (1), one can easily derive the governing equation for the concentration of the tracer fluid ($Z\equiv\rho^{\prime}/\rho$): $dZ/dt=0$, i.e., the value of $Z$ for each fluid element does not change with time. In our simulation, we choose as the initial condition at $t=0$ the values of $Z$ to be:

[TABLE]

Thus the concentration $Z$ can be used to trace the spatial evolution of the ICM gas initially located within the central sphere of radius $100$ kpc.

Figure 8 shows the time evolution of the concentration distribution ($Z$) of the tracer fluid. It is clear that the trailing outflow is real: some of the gas originally located within $r\leq 100$ kpc is uplifted to $r\sim 100-200$ kpc along the jet direction by the AGN bubble at $t=400$, $500$, $700$ Myr. The trailing outflow along the $z$ direction and inflows in other directions form gas circulations in the meridional planes, and near the bottom of the AGN bubble, meridional vortexes also appear. At $t=700$ Myr, one can clearly see that the trailing flow is falling back to the cluster centre (see the overplotted velocity vectors), suggesting that the outflow does not leave the cluster eventually, but instead it is confined to the cluster.

Trailing outflows play an important role in suppressing or delaying the central cooling catastrophe. As the rarefaction wave propagates away to large radii, gas inflows resume in the inner core region, as clearly seen in Figures 3 and 8 and discussed in Section 3.2. If there are no trailing outflows, the gas density in the cluster centre would rise rapidly and a central cooling catastrophe would then develop. Trailing outflows uplifted by AGN bubbles transport gas from the inner cluster core to outer regions, suppressing gas accumulation at the cluster centre resulting from inflows along other directions and thus averting the central cooling catastrophe, which is delayed until trailing outflows fall back to the cluster centre at later times (see Section 3.6).

Observational evidence for trailing outflows has been found in some cool core clusters where the spatial distributions of heavy elements are anisotropic and preferentially aligned with the large-scale radio and cavity axes (Kirkpatrick et al. 2011; Kirkpatrick & McNamara 2015), suggesting that hot metal-rich outflows are indeed lifted up by AGN jets. Hydrodynamic simulations following both the jet and metallicity evolution in galaxy clusters would be helpful to further explore this scenario (e.g., Barai et al. 2016, which however does not investigate the angular distribution of metals). Along the axes of trailing outflows, gas density is relatively high, as seen in the bottom two rows of Figure 3, and cold gas may drop out due to local thermal instability, forming cold filaments found in some galaxy clusters (e.g., Mittal et al. 2012; Tremblay et al. 2015).

3.6 Emerging Picture of Cyclic Cool Cores

What is the long-term evolution of cool cluster cores after powerful AGN feedback events? After the onset of the AGN jet event, the cool core in our simulation is quickly over-heated by shock dissipation, and rarefied by the rarefaction wave which transports energy to large radii. As the bow shock weakens and propagates further away, the inner core region starts to contract again and form meridional circulations with trailing outflows uplifted by the AGN bubble, as shown in Figure 3. Trailing outflows avoid rapid accumulation of gas at the cluster centre, and thus delay the onset of central cooling catastrophe. However, as shown by the velocity field in the bottom panel of Figure 8, the trailing outflows could not eventually leave the deep potential well of the cluster, and at $t=700$ Myr, they are clearly falling back to the cluster centre along the $z$ axis. The gas density at the cluster centre would then increase rapidly, and at $t=746$ Myr in our simulation, a cooling catastrophe happens again in the cluster centre. The temporal evolution of gas masses within three representative spheres and the mass flow rates across their surfaces are shown in Figure 9, indicating that the mass inflow rate across $r=50$ kpc in the developed cooling flow in our simulated cluster A1795 is about $160$ - $200M_{\odot}/$yr. The mass inflow rate in the cooling flow phase increases with radius, suggesting that gas densities in the cool core increase continuously, as also seen in the left panel of Figure 2.

Thus, the evolution of cool cluster cores seems to be cyclic, and a typical cycle is summarized schematically in Figure 10. What would be observational evidence for cyclic cool cores going through cooling-induced contraction and AGN-induced expansion. Cool cores are under-heated during the cooling flow phase, but over-heated during the AGN outburst phase. Thus in the former phase, gas densities and the X-ray luminosity of the cool core are relatively higher, while the central gas entropy would be relatively lower. For galaxy clusters with the same virial mass, cyclic cool cores suggested by our simulation would thus lead to an anti-correlation between the gas fraction (or X-ray luminosity) of cool cores and the central gas entropy, which is not naturally expected in steady-state models of cool cluster cores.

4 Summary and Discussion

We investigate the cooling flow problem in a typical cool core cluster – Abell 1795, and in particular, focus on how a powerful AGN jet event triggered by the central cooling catastrophe, heats and reverses the cluster core from the “catastrophic cooling” state back to its initial cool-core state. We use a “pure cooling” control run where no AGN jets are triggered for comparison, and investigate in detail complex gasdynamical processes and energy evolution in an “AGN feedback” run, where a powerful AGN jet event is turned on once the central cooling catastrophe starts. The jet parameters are specifically chosen so that the cluster profiles will reverse back to their initial distributions after the jet event.The supersonic jet is active steadily for a duration of $5$ Myr, and its power $P_{\rm jet}=1.5\times 10^{46}$ erg/s, which is dominated by the kinetic energy, is comparable to that estimated in Cygnus A (Wilson et al., 2006) and MS0735 (McNamara et al., 2009), but less than the maximum jet power $\sim 10^{47}$ erg/s reached in some recent simulations (e.g., Gaspari et al. 2011; Li et al. 2015; Yang & Reynolds 2016b). Our main findings are summarized below:

• •

Without AGN feedback, a cooling inflow gradually develops, and cooling catastrophe first appears at the cluster centre (at $t\sim 237.6$ Myr in our simulation), which is consistent with previous studies (e.g., Li & Bryan 2012 and Guo & Mathews 2014).

• •

The AGN jet triggered after the central cooling catastrophe drives a bow shock, which reverses cooling inflows into outflowing motions. A rarefaction wave following the shock front expands the shocked ICM, rarefying the dense cool core and lowering central cooling rates.

• •

The powerful AGN jet event is able to bring the cluster core from the catastrophic cooling state back to the initial cool-core state. In our simulation, ICM heating is dominated by shock heating, while mixing and turbulent heating are less important. Most shock energy is dissipated rapidly in the inner core region, while the rarefaction wave is efficient in transporting and re-distributing energy to outer regions.

• •

As the AGN bubble rises away from the cluster centre, it uplifts trailing outflows in its wake along the original jet direction. Since inflows resume in the core along other directions as the rarefaction wave propagates away, trailing outflows suppress gas accumulation and cooling rates in the cluster centre by transporting gas from the inner core to outer regions. As a result, meridional circulations also form in the cool core.

• •

Trailing outflows are intracluster (i.e., confined to galaxy clusters), and at later times, a fraction of them fall back to the cluster centre, triggering central cooling catastrophe and potentially a new generation of AGN feedback.

• •

We propose a picture of cyclic cool cores that cycle through cooling-induced contraction and AGN-induced expansion. During the former state, cool cores are under-heated, while during the latter, they are over-heated. This picture thus naturally leads to an anti-correlation between the gas fraction (or X-ray luminosity) of cool cores and the central gas entropy, which may be tested by X-ray observations.

We discuss several limitations of our work here. First, our simulations focus on a specific powerful AGN jet event, while in reality AGN feedback events have a variety of powers and durations. The mechanical variability of AGN outbursts is also quite uncertain. If an AGN jet episode is less powerful but lasts for a longer duration, AGN heating and dynamical feedback would be more gentle. In this case, the dominant heating mechanism and the following energy redistribution may be different than in the present paper, which we leave for future studies.

For simplicity, the parameters of the jet in its initialization nozzle are fixed with time in our simulation when the jet is active. While this assumption potentially corresponds to a steady accretion flow onto the central SMBH after cold gas circularizes around it, it does not reflect a long-term feedback mechanism that automatically links gas cooling, SMBH accretion, and AGN outbursts. Therefore, although our simulation indicates that a powerful AGN jet event is able to reverse a galaxy cluster from the catastrophic cooaling state back to the initial cool-core state, we could not yet answer if cooling flows can be averted by AGN feedback over a long, multi-Gyr timescale. Recent simulations tackling this long-term cooling flow problem often invoke jet precession to avoid the formation of a low-density channel along the jet direction that could rapidly transport jet energy out of cool cores (Gaspari et al. 2012; Li et al. 2015; Yang & Reynolds 2016b), while dynamic ICM weather in realistic clusters may also help distribute jet energy to broader angular regions (Mendygral et al. 2011; Mendygral et al. 2012).

In our simulation, gas motions induced by the AGN jet event are quite regular, and the level of turbulence is not significant, which is consistent with recent studies by Reynolds et al. (2015), Yang & Reynolds (2016b) and Li et al. (2016). However, we note that simulations with grid-based codes inevitably introduce numerical diffusion, which may suppress the development of turbulence. A careful investigation on this effect is required to fully understand the importance of turbulence in AGN feedback.

Finally, our simulation only includes pure hydrodynamics, but neglects some more subtle physics, e.g., viscosity, magnetic fields, cosmic rays, etc. Cosmic rays and magnetic fields are present and may even dominate in AGN jets and bubbles (Guo & Oh 2008; Guo & Mathews 2011; Ruszkowski et al. 2017; Gan et al. 2017). Shear viscosity (Kaiser et al. 2005; Reynolds et al. 2005; Guo et al. 2012; Guo 2015) and magnetic tension (Kaiser et al. 2005; Jones & De Young 2005; Ruszkowski et al. 2007) tend to prevent the disruption of AGN bubbles by suppressing the development of RT and KH instabilities, which may indeed happen in some old intact X-ray cavities detected in galaxy clusters. Nonetheless, we expect that our main results on complex gasdynamical processes induced by AGN jet events, including shocks, rarefaction waves, trailing outflows, and meridional circulations, and their importance in solving the cooling flow problem should be robust.

Acknowledgments

FG thanks James Binney, Fabrizio Brighenti, and Eugene Churazov for helpful discussions. We thank the referee for a useful report that helped to improve the manuscript. This work was supported by Chinese Academy of Sciences through the Hundred Talents Program and the Frontier Sciences Key Project (No. QYZDB-SSW-SYS033), and National Natural Science Foundation of China (Grant No. 11643001 and 11633006). YFY gracefully acknowledges support from National Natural Science Foundation of China (Grant No. U1431228, 11133005, 11233003, and 11421303). Some simulations presented in this work were performed using the high performance computing resources in the Core Facility for Advanced Research Computing at Shanghai Astronomical Observatory.