Cosmic-Ray Feedback Heating of the Intracluster Medium

The gravitational potential and initial conditions for the temperature and density distributions of the gas resemble those adopted by Yang & Reynolds (2016a). In brief, the cluster atmosphere is initially close to hydrostatic equilibrium and its density profile is similar to that corresponding to the Perseus cluster.

We include tangled magnetic fields that are generated using the method similar to that described in Ruszkowski et al. (2007). We assume that in Fourier space the field has the following form,

$$\begin{eqnarray}&&B\propto {k}^{-11/6}\exp \left[-{\left(\displaystyle \frac{k}{{k}_{\mathrm{in}}}\right)}^{4}\right],\end{eqnarray} \tag{ 1 }$$

where k_in = 10²(2π/L), with L = 1 Mpc the size of the computational domain. We perform an inverse Fourier transform to generate real-space magnetic fields and, following Wiener et al. (2013), we rescale the field such that $B\propto {\rho }_{o}^{0.3}$, where ρ_o is the ICM density. This ensures that the magnetic pressure is approximately proportional to the gas pressure. In order to generate a divergence-free field, we Fourier transform the field and perform divergence cleaning as in Ruszkowski et al. (2007). This procedure is repeated until a divergence-free field proportional to ${\rho }_{o}^{0.3}$ is obtained. The final field is normalized such that plasma β ∼ 10². We also impose small isobaric perturbations δρ/ρ on top of the average gas density profiles. Following Gaspari et al. (2012), these fluctuations are approximately characterized by white noise spectrum with the amplitude of 0.1. The resulting ICM gas density distribution is given by $\rho ={\rho }_{o}\max (0.8,1+\delta \rho /\rho )$.

We use adaptive mesh refinement to refine the domain up to the maximum resolution of 1.95 kpc. Refinement is triggered by temperature gradients. We employ diode boundary conditions (the gas is only allowed to flow out of the domain; code variables have vanishing gradient at the boundary) but note that the choice of boundary conditions is not critical as the domain is much larger than the size of the central parts of the cool core.

The black hole feedback model adopted here is based on the "chaotic cold accretion" model (Gaspari et al. 2012, 2013; Li et al. 2015) and closely follows that used by Yang & Reynolds (2016a). In this model, the cooling gas is removed from the hot phase of the ICM when its temperature drops below T = 5 × 10⁵ K. The cold gas is then converted to passive particles that follow the fluid and are allowed to accrete onto the central black hole, triggering feedback. The AGN energy is supplied back to the ICM via bipolar precessing jets.

Compared to the feedback model used by Yang & Reynolds (2016a), the main difference is that here we also include MHD and CR physics and consequently the energy injected by the AGN jets is supplied in kinetic and CR form. We consider jets dominated by the CR component and assume that a fraction of f_cr = 0.8 of the energy of the jet fluid is in the form of CRs. In this model, the accretion rate onto the center is approximated as ${\dot{M}}_{\mathrm{in}}={M}_{\mathrm{cold}}/{t}_{\mathrm{ff}}$, where M_cold is the mass of the gas with T < 5 × 10⁵K and t_ff is the free-fall time at the accretion radius r_accre. The AGN injects mass, momentum, and energy via bipolar jets according to

$$\begin{eqnarray*}{\dot{M}}_{\mathrm{out}} & = & \eta {\dot{M}}_{\mathrm{in}}\\ | \dot{P}| & = & \sqrt{2\eta \epsilon (1-{f}_{\mathrm{cr}})}{\dot{M}}_{\mathrm{in}}\\ {\dot{E}}_{\mathrm{kin}} & = & \epsilon (1-{f}_{\mathrm{cr}}){\dot{M}}_{\mathrm{in}}{c}^{2}\\ {\dot{E}}_{\mathrm{cr}} & = & \epsilon {f}_{\mathrm{cr}}{\dot{M}}_{\mathrm{in}}{c}^{2}.\end{eqnarray*}$$

The adopted model parameters are the jet mass loading factor $\eta =1$, feedback efficiency = 10⁻³, accretion timescale t_ff = 5 Myr, accretion radius r_accre = 5.85 kpc, precession period of the jet t_prec = 10 Myr, and a narrow precession angle of 15°. The feedback energy is injected in a cylinder of 5 kpc in radius and 4 kpc in height.

We solve the MHD equations including CR advection, dynamical coupling between CR and the thermal gas, CR streaming along the magnetic field lines and the associated heating of the gas by the CR, heating of the ICM by Coulomb and hadronic interactions, and radiative cooling:

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

$$\begin{eqnarray}&&\displaystyle \frac{\partial \rho {{\boldsymbol{u}}}_{{\rm{g}}}}{\partial t}+{\rm{\nabla }}\cdot \left(\rho {{\boldsymbol{u}}}_{{\rm{g}}}{{\boldsymbol{u}}}_{{\rm{g}}}-\displaystyle \frac{{\boldsymbol{BB}}}{4\pi }\right)+{\boldsymbol{\nabla }}{p}_{\mathrm{tot}}=\rho {\boldsymbol{g}}+{\dot{p}}_{{\rm{j}}},\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\boldsymbol{B}}}{\partial t}-{\boldsymbol{\nabla }}\times ({{\boldsymbol{u}}}_{{\rm{g}}}\times {\boldsymbol{B}})=0,\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial e}{\partial t}+{\boldsymbol{\nabla }}\cdot \left[(e+{p}_{\mathrm{tot}}){{\boldsymbol{u}}}_{{\rm{g}}}-\displaystyle \frac{{\boldsymbol{B}}({\boldsymbol{B}}\cdot {{\boldsymbol{u}}}_{{\rm{g}}})}{4\pi }\right]=\rho {{\boldsymbol{u}}}_{{\rm{g}}}\cdot {\boldsymbol{g}}\\ &&\quad -{\rm{\nabla }}\cdot {{\boldsymbol{F}}}_{{\rm{c}}}-{ \mathcal C }+{{ \mathcal H }}_{{\rm{c}}}+{{ \mathcal H }}_{{\rm{j}}},\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial {e}_{{\rm{c}}}}{\partial t}+{\boldsymbol{\nabla }}\cdot ({e}_{{\rm{c}}}{{\boldsymbol{u}}}_{{\rm{g}}})=-{p}_{{\rm{c}}}{\boldsymbol{\nabla }}\cdot {{\boldsymbol{u}}}_{{\rm{g}}}-{\rm{\nabla }}\cdot {{\boldsymbol{F}}}_{{\rm{c}}}+{{ \mathcal C }}_{{\rm{c}}}+{{ \mathcal H }}_{{\rm{j}}},\end{eqnarray} \tag{ 6 }$$

where ρ is the gas density; ${\boldsymbol{u}}$_g is the gas velocity; ${\boldsymbol{B}}$ is the magnetic field; ${\boldsymbol{g}}$ is the gravitational field; ${\dot{\rho }}_{{\rm{j}}}$ is the rate of injection of thermal gas via jet; ${\dot{p}}_{{\rm{j}}}$ is the rate of momentum injection associated with the AGN; e_c is the specific CR energy density; $e=0.5\rho {{\boldsymbol{u}}}_{{\rm{g}}}^{2}+{e}_{{\rm{g}}}+{e}_{{\rm{c}}}+{B}^{2}/8\pi$ is the total energy density; ${ \mathcal C }$ is the radiative cooling energy loss rate per unit volume; ${\boldsymbol{F}}$_c is the CR flux due to streaming relative to the gas; ${{ \mathcal H }}_{{\rm{c}}}$ is the rate of change of the total specific energy due to streaming instability heating of the gas and Coulomb and hadronic CR losses; ${{ \mathcal C }}_{{\rm{c}}}$ is the CR cooling rate due to the streaming instability, Coulomb, and hadronic CR losses; and ${{ \mathcal H }}_{{\rm{j}}}$ represents heating due to the AGN. The total pressure is ${p}_{\mathrm{tot}}=({\gamma }_{{\rm{g}}}-1){e}_{{\rm{g}}}+({\gamma }_{{\rm{c}}}-1){e}_{{\rm{c}}}+{B}^{2}/8\pi$, where e_g and e_c are the specific thermal energy density of the gas, γ_g = 5/3 is the adiabatic index for ideal gas, and γ_c = 4/3 is the effective adiabatic index of CR fluid.

Radiative cooling is included using the Sutherland & Dopita cooling function (Sutherland & Dopita 1993). In order to speed up the computations, we employ the subcycling method (Anninos et al. 1997; Proga et al. 2003) when the local cooling time becomes shorter than the hydrodynamical time step.

We solve the above equations using the adaptive mesh refinement MHD code FLASH4.2 (Fryxell et al. 2000; Dubey et al. 2008). We employ the directionally unsplit staggered mesh solver (Lee & Deane 2009; Lee 2013). This solver is based on a finite-volume, high-order Godunov scheme and utilizes a constrained transport method to enforce divergence-free magnetic fields. We use the third-order MHD scheme and HLLD Riemann solver.

We include the heating of the ICM by CRs and transport of CRs with respect to the gas. Details of the CR physics module can be found in Yang et al. (2012) and Ruszkowski et al. (2017), where we discuss simulations of the Fermi bubbles and CR-driven galactic winds, respectively. We now summarize key CR physics processes described in that paper and discuss extensions of the CR module specific to the modeling of the ICM presented here.

2.3.1. Streaming of CRs

Propagation of CRs in the magnetized ICM can be described in the framework of the self-confinement model. In this picture, CRs scatter on waves excited by the streaming instability (Kulsrud & Pearce 1969; Wentzel 1974; Zweibel 2013). In a state of marginally stable anisotropy, the CRs stream at the Alfvén speed down their pressure gradients. However, the waves excited by the streaming instability can be damped by various mechanisms, e.g., by turbulent or Landau damping. When this happens, CRs can stream at speeds exceeding the Alfvén speed. The effective streaming speed increases with the strength of the damping mechanism. The streaming flux is given by ${\boldsymbol{F}}$_cr = (e_cr + p_cr)${\boldsymbol{u}}$_s, where ${{\boldsymbol{u}}}_{{\rm{s}}}=-\mathrm{sgn}(\hat{{\boldsymbol{b}}}\cdot {\rm{\nabla }}{e}_{\mathrm{cr}})f{{\boldsymbol{u}}}_{{\rm{A}}}$ is the streaming velocity, ${\boldsymbol{u}}$_A is the Alfvén velocity, and f is the streaming speed boost factor.

As demonstrated by Wiener et al. (2013), the effective streaming speed in the ICM can significantly exceed the Alfvén speed in the cluster outskirts. For conditions representative of the cluster cool cores, damping mechanisms can lead to moderately super-Alfvénic speeds for the following reasons. Wiener et al. (2013) consider turbulent and nonlinear Landau damping mechanisms. In the turbulent damping case, the effective streaming speed is

$$\begin{eqnarray}&&{u}_{{\rm{s}}}={u}_{{\rm{A}}}\left(1+0.08\displaystyle \frac{{B}_{10\,\mu {\rm{G}}}^{1/2}{n}_{{\rm{i}},-2}^{1/2}}{{L}_{\mathrm{mhd},\ 10}^{1/2}{n}_{{\rm{c}},-9}}{\gamma }_{3}^{n-3.5}{10}^{2(n-4.6)}\right),\end{eqnarray} \tag{ 7 }$$

where ${n}_{{\rm{i}},-2}={n}_{{\rm{i}}}/({10}^{-2}\,{\mathrm{cm}}^{-3})$ is the ion number density, ${n}_{{\rm{c}},-9}={n}_{{\rm{c}}}/({10}^{-9}\,{\mathrm{cm}}^{-3})$ is the CR number density, ${L}_{\mathrm{mhd},10}={L}_{\mathrm{mhd}}/(10\,\mathrm{kpc})$ is the length scale at which turbulence is driven at the Alfvén speed u_A, γ₃ = γ/3 is the average CR Lorentz factor, and n > 4 is the slope of the CR distribution function in momentum (approximately n = 4.6). In the nonlinear Landau damping case, the effective streaming speed is

$$\begin{eqnarray}&&{u}_{{\rm{s}}}={u}_{{\rm{A}}}\left(1+0.03\displaystyle \frac{{n}_{{\rm{i}},-2}^{3/4}{T}_{5\,\mathrm{keV}}^{1/4}{10}^{n-4.6}{\gamma }_{3}^{(n-3)/2}}{{B}_{10\mu {\rm{G}}}{L}_{\mathrm{cr},10}^{1/2}{n}_{{\rm{c}},-9}^{1/2}}\right),\end{eqnarray} \tag{ 8 }$$

where ${L}_{\mathrm{cr},10}={L}_{\mathrm{cr}}/(10\,\mathrm{kpc})$ is the characteristic length scale of the fluctuations in the CR distribution and ${T}_{5\mathrm{keV}}=T/(5\,\mathrm{keV}$) is the ICM temperature. For the conditions representative of cool cores, in both of these cases, CR streaming is not typically super-Alfvénic. However, the damping rate Γ may be further boosted by linear Landau damping, leading to Γ_Landau/Γ_turb ∼ β^1/2, where β is the plasma β ∼ 10² parameter in the ICM (E. Zweibel 2017, in preparation). When this process is included, the second term in Equation (7) needs to be multiplied by β^1/2. For plausible cool core parameters, the CR number density is

$$\begin{eqnarray}&&{n}_{{\rm{c}}}=3\times {10}^{-9}\displaystyle \frac{n-4}{n-3}{q}_{-2}{n}_{{\rm{i}},-2}{T}_{5\mathrm{keV}}{E}_{\min ,\mathrm{GeV}}^{-1},\end{eqnarray} \tag{ 9 }$$

where q is the ratio of CR pressure to the ICM pressure and ${E}_{\min ,\mathrm{GeV}}$ is the low-energy cutoff in CR momentum distribution. Given the uncertainty in β, L_mhd, and n_c, it is plausible that the effective CR streaming speed could be moderately super-Alfvénic, i.e., boosted by a factor of order unity beyond the Alfvén speed. Therefore, in addition to Alfvénic streaming, we also consider super-Alfvénic streaming for f = 1.71, where the streaming speed has been estimated using the fiducial parameters in Equations (7) and (9) assuming Landau damping and β = 100.

CR streaming is incorporated using the method of Sharma et al. (2009). Because the term $-{\rm{\nabla }}\cdot {{\boldsymbol{F}}}_{\mathrm{cr}}$ varies infinitely fast due to the discontinuity in the streaming flux near CR energy local extrema, it leads to a prohibitively small simulation time step. In order to remove the singularity and speed up computations, we regularize the streaming flux by ${{\boldsymbol{F}}}_{{\rm{c}}}\,=-({e}_{{\rm{c}}}+{p}_{{\rm{c}}}){{\boldsymbol{u}}}_{{\rm{A}}}\tanh ({h}_{{\rm{c}}}\hat{{\boldsymbol{b}}}\cdot {\rm{\nabla }}{e}_{{\rm{c}}}/{e}_{{\rm{c}}})$, where h_c is a free (regularization) parameter. In the calculations presented in this paper, we adopt h_c = 100 kpc. This approach for modeling CR transport is appropriate for the cases where wave damping occurs via turbulent, Landau damping or where no damping is included. The super-Alfvénic case (SCHTs) corresponds to Landau damping.

2.3.2. ICM Heating by CRs

As the CRs stream, they also experience an effective drag force. Consequently, CRs lose energy and the gas is heated due to the Alfvén wave heating at the rate of

$$\begin{eqnarray}&&{{ \mathcal H }}_{\mathrm{cr},\mathrm{stream}}=-{{\boldsymbol{u}}}_{{\rm{A}}}\cdot {\rm{\nabla }}{p}_{\mathrm{cr}}.\end{eqnarray} \tag{ 10 }$$

In addition to the heating of the ICM associated with the streaming instability, CRs also heat the gas via Coulomb and hadronic interactions. We approximate the effects of CR cooling due to Coulomb and hadronic losses due to pion production via (Yoast-Hull et al. 2013),

$$\begin{eqnarray}&&{C}_{\mathrm{cr},{\rm{c}}}=-4.93\times {10}^{-19}\displaystyle \frac{n-4}{n-3}\displaystyle \frac{{e}_{{\rm{c}}}\rho }{{E}_{\min }}\displaystyle \frac{\rho }{{\mu }_{{\rm{e}}}{m}_{{\rm{p}}}}\,\mathrm{erg}\,{\mathrm{cm}}^{-3}\,{{\rm{s}}}^{-1},\end{eqnarray} \tag{ 11 }$$

and due to hadronic losses via

$$\begin{eqnarray}&&{C}_{\mathrm{cr},{\rm{h}}}=-8.56\times {10}^{-19}\displaystyle \frac{n-4}{n-3}\displaystyle \frac{{e}_{{\rm{c}}}\rho }{{E}_{\min }}\displaystyle \frac{\rho }{{\mu }_{{\rm{p}}}{m}_{{\rm{p}}}}\,\mathrm{erg}\,{\mathrm{cm}}^{-3}\,{{\rm{s}}}^{-1},\end{eqnarray} \tag{ 12 }$$

where ${E}_{\min }=1\,\mathrm{GeV}$ is the minimum energy of CRs, and μ_e and μ_p are the mean molecular weights per electron and proton, respectively. In the simulations, we assume n = 4.5 and mean proton Lorentz γ = 3. While all of the CR energy loss due to Coulomb collisions is transferred to the gas, only ∼1/6 of the CR energy loss due to pion production is used to heat the gas, and the remainder is removed as gamma-ray emission and neutrinos. Consequently, the rate of change of the total specific energy density of the gas, which includes the thermal and CR specific energy densities, is ${{ \mathcal H }}_{\mathrm{cr}}=(5/6){C}_{\mathrm{cr},{\rm{h}}}/\rho \lt 0$ and the CR specific energy density loss rate is ${{ \mathcal C }}_{{\rm{c}}}=({C}_{\mathrm{cr},{\rm{c}}}+{C}_{\mathrm{cr},{\rm{h}}})/\rho$.

The list of the performed runs is shown in Table 1. Figure 1 presents cross-sections through the cluster center showing the distribution of the specific CR energy density. Clockwise from upper left, these slices correspond to the following cases: (i) hadronic and Coulomb heating but no transport processes (CHT0), (ii) CR streaming and streaming heating (ST1), (iii) CR streaming and heating due to streaming, hadronic, and Coulomb processes (SCHT1), and (iv) same as the last panel but for super-Alfvénic streaming (SCHTs). All snapshots were taken at 3 Gyr. This figure demonstrates that CR transport processes affect the morphology of the radio-emitting plasma and effectively redistribute CRs. The redistribution of CR energy is efficient, especially in the central cooling region of r ∼ 200 kpc, despite the fact that the jet is pointed in an approximately constant direction. As expected, the widening of the CR distribution is most significant when the CR transport is the fastest, i.e., super-Alfvénic. As described in detail below, in the simulations including CR streaming, the ICM generally exhibits larger variations due to more intermittent AGN feedback. This means that the atmosphere can experience both periods of relative calm and more stormy conditions. Recent Perseus data from Hitomi is consistent with relatively low level of turbulence in this cluster (Hitomi Collaboration et al. 2016). It is plausible that the dynamical state of the Perseus cluster currently corresponds to the relatively low-turbulence state captured in Figure 1 in cases including transport processes (see also Li et al. 2016). Alternatively, turbulent motions in the cluster atmosphere could be reduced due to viscosity. We also point out that the iron line shifts corresponding to large gas velocities induced by the AGN at the center of the cool core may be partially diluted by slower moving gas away from the center. This may give an impression of relative calm in the ICM even if fast gas motions are present. This dilution effect has been seen in mock Hitomi simulations that show line shifts consistent with the data (B. Morsony 2017, private communication). We defer to a future publication the study of the iron emission line profiles and observational predictions for the planned Hitomi replacement and the X-ray Surveyor missions.

Zoom In Zoom Out Reset image size

Table 1.  List of Simulations

Run Name	Coulomb/Hadronic Heating	Streaming Heating	Transport Speed
CHT0	yes	no	0
ST1	no	yes	vA
SCHT1	yes	yes	vA
SCHTs	yes	yes	$1.71{v}_{{\rm{A}}}$

Download table as:  ASCIITypeset image

As expected, the dispersal of CRs throughout the core is more pronounced at later times since the onset of feedback and when the speed of CR transport is faster. Interestingly, observations of M87 with LOFAR reveal a sharp radio emission boundary that does not seem to depend sensitively on radio frequency (de Gasperin et al. 2012), i.e., it appears that the boundary corresponds to the physical extent of CRs. At late times, no such boundary is seen in the simulations. However, such boundary in the spatial distribution of CRs could be explained by large-scale sloshing motions that order magnetic fields on large scales and prevent the leakage of CRs to large distances by suppressing cross-field CR transport. The simulations of ZuHone et al. (2013) show that sloshing motions induced by the substructure in the cluster can generate tangential magnetic fields. Such fields could slow down the radial transport of CRs away from the core. Alternatively, weaker or less collimated AGN feedback could prevent the bubbles from overshooting the critical radius at which their internal entropy equals that of the ambient ICM. In such a case, we would expect CRs to exist predominantly within such a critical radius. We defer the exploration of these possibilities to a future publication and point out that there exist counterexamples to the morphological appearance of M87. In Abell 262 (Clarke et al. 2009) and A2597 (Clarke et al. 2005), the radio emission at lower frequencies extends to larger distances from the cluster center.

The pressure support due to CRs is quantified in Figure 2. Pressure support is defined as the ratio of the pressure provided by CRs to the sum of the thermal and CR pressures. In order to exclude CR-filled bubbles that are cooling very inefficiently, this quantity is set to 10⁻² whenever the local cooling time exceeds the Hubble time (this is done only at the post-processing stage, i.e., we do not modify CR pressure in the simulations). All panels show the evolution of the profiles of the pressure support. Dark lines corresponds to 20% of CR contribution to the total pressure (i.e., CR-to-thermal pressure ratio of 25%). In the case excluding CR transport (left panel), CR interaction with the ambient medium is inhibited. Even though hadronic and Coulomb heating processes are included in this case, the CR heating of the ambient ICM is ineffective because CRs do not easily come in contact with the thermal ICM. In order to compensate for this inefficiency in the ICM heating by CRs, more CRs are injected into the ICM. This is a runaway process in which CRs tend to account for a progressively larger fraction of the total pressure support. At the end of the simulation, the CR pressure support in ∼50 kpc builds up to a significant level, creating tension with the observations (Jacob & Pfrommer 2016b).

Zoom In Zoom Out Reset image size

We note in passing that the magnetic fields did not play a significant role in suppressing the mixing of CRs and the ICM for the parameter sets corresponding to the non-transport cases presented here. However, we also performed additional runs (not presented in this paper) for higher CR injection fraction and lower AGN efficiency. We observed that, in the case for higher CR injection fraction (f_cr = 0.9), significant accumulation of CRs in the cool core occurred earlier when magnetic fields were included compared to the case without magnetic fields. This is consistent with the magnetic fields playing a role in suppressing the mixing of CRs with the ICM—on average, the AGN injected more energy to compensate for the inefficiency of the heating caused by the slower mixing. The relative inefficiency of the mixing was likely caused by a slower jet in this case as the kinetic energy of the jet was reduced from 20% to 10% of the injected AGN energy. This is consistent with the findings of Yang & Reynolds (2016a), who considered AGN energy in purely kinetic form and attributed heating to efficient mixing. However, this higher f_cr case resulted in large variations in the profiles of thermodynamical quantities (or close-to-isothermal temperature profiles) in the cases that included transport and, thus, we rejected this set of models. We also performed simulations for f_cr = 0.9 and lower AGN feedback efficiency = 3 × 10⁻⁴. In this case, the accumulation of CRs in the cool core in the non-transport case was the largest of all considered cases. While the profiles of thermodynamical quantities were in general consistent with observations and did not exhibit dramatic variability, all runs for f_cr = 0.9 and = 3 × 10⁻⁴ resulted in excessive mass deposition rates and, thus, we rejected this set of models.

The remaining three panels illustrate that the role of transport processes is essential for removing this tension with observations. The second panel shows that including CR streaming and associated streaming heating dramatically reduces CR contribution to the pressure support. This reduction in CR pressure occurs because CRs can now come into contact with the thermal ICM and heat it, thus reducing the CR energy density and associated CR pressure. Similarly, as shown in the third panel, CR pressures are further reduced when, in addition to the processes included in the second panel, we also include CR hadronic and Coulomb losses. These two processes further drain the energy from CRs and heat the thermal gas. Finally, the last panel demonstrates the consequences of including faster (super-Alfvénic) streaming. Compared to the Alfvénic streaming case, in this case CRs are redistributed more efficiently from the very center of the cool core to larger distances still within the cool core. This redistribution of CRs also reduces the CR pressure gradient and, thus, it decreases the amount of CR heating associated with the streaming instability. These two factors lead to the reduction of the CR-to-thermal pressure in the very center of the cool core and a slight increase of this ratio at larger distances within the cool core. Note that this boost in the CR streaming speed only affects the rate of CR transport rather than the Alfvén wave heating. In all cases but the one shown in the leftmost panel, the CR pressure support is very small. We note that the buoyant stability condition in cluster cores (in the presence of anisotropic thermal conduction and anisotropic CR transport) is nk_b∇ T + ∇P_cr > 0, where n is the gas number density and k_b is the Boltzmann constant (Chandran & Rasera 2007). Since P_cr is generally very small in most of the cool core volume in our simulations, convection should not play a major role in redistributing CRs (with the exception of buoyancy associated with young unmixed AGN bubbles and jets.

While predicting detailed observational gamma-ray and radio signatures based on these simulations is beyond the scope of this paper, we point out that the typical levels of CR pressure support that we find in simulations including CR transport are generally consistent with the data. Based on one-dimensional models that include heating by thermal conduction and CRs, Jacob & Pfrommer (2016b) argue that in those cool core clusters that do not host radio mini halos, AGN activity and CR heating are the strongest, and that CRs can provide an adequate level of heating without violating observational radio and gamma-ray constraints. They further argue that primary and secondary CR electron radio emission associated with the AGN outbursts could be difficult to detect due to the small physical extent of the radio emission in this case and the large flux dynamical range of the AGN jet and the halo. This picture is likely to be consistent with the elevated CR pressure support during AGN outbursts that is seen in Figure 2 (e.g., near ∼3 Gyr in the third panel). In Jacob & Pfrommer (2016a), typical values of CR-to-thermal pressure are on the order of 0.1 and vary substantially from object to object and thus presumably depend on the cluster dynamical state. Interestingly, Pfrommer (2013) shows that in the Virgo cluster, in the absence of thermal conduction, adequate CR heating rate can be supplied when CR-to-thermal fraction is around 0.3 while not violating observational data. When streaming is included, the levels of CR pressure support that we observe in our simulations during outbursts are comparable to those suggested by Jacob & Pfrommer (2016a) and could presumably be reduced further if we included thermal conduction. In the case of cool cores that are associated with radio mini halos, Jacob & Pfrommer (2016b) predict that the amount of CR pressure support needed to stably heat the cool core exceeds observational limits and suggest that such objects are expected to be dominated by radiative cooling. This situation could correspond to the periods in between the outbursts seen in Figure 2. Thus, the general properties of our simulations, and in particular the presence of the feedback loop and two classes of cool cores, are broadly consistent with the picture based on the above one-dimensional models. We also note that the simulations that do not include CR transport processes (left panel in Figure 2) do not show intermittent AGN activity and would therefore not be able to account for the transitions between cool cores with and without radio mini halos.

Giacintucci et al. (2017) report that most cool core clusters more massive than ${M}_{500,\mathrm{th}}=6\times {10}^{14}$ M_⊙ host radio mini halos. This appears to be consistent with the suggestions in Jacob & Pfrommer (2016b). Using an approximate conversion between M₂₀₀ and M₅₀₀ (White 2001), it is clear that most of the clusters with radio mini halo detections presented in Jacob & Pfrommer (2016b) also tend to be those with M₅₀₀ above this threshold. Interestingly, Jacob & Pfrommer (2016b) also suggest that it is the clusters without observable radio halos that are CR heated. Most of the clusters in the sample considered by Jacob & Pfrommer (2016b) with ${M}_{500}\lt {M}_{500,\mathrm{th}}$ are not predicted to have radio mini halos and would thus be consistent with being heated by CRs. Finally, we note that here we focus on general trends and defer to a future publication the study of the parameter space of the models (e.g., transport speed, CR injection fraction, feedback efficiency) that meets observational constraints in detail for a range of cool core cluster masses.

In all three cases that include transport processes (panels 2–4 in Figure 2), there is a significant variation in the CR pressure support over time. This behavior of the atmosphere is reflected in Figure 3, which shows the AGN jet power as a function of time. In all four cases but the one shown in the first panel, the black hole feedback is highly variable. Note that despite the large variability, the AGN jet never completely switches off. In all cases, the average jet power corresponds to mass accretion rate of approximately 60 M_⊙ yr⁻¹. This is about 20% of the nominal unimpeded mass deposition rate in the Perseus cluster and is thus in good agreement with observations.

Zoom In Zoom Out Reset image size

The evolution of the X-ray luminosity within the central 100 kpc is shown in Figure 4. The green line corresponds to bolometric brehmsstrahlung luminosity and the black line to the X-ray emission integrated in the 0.5–10 keV range. As the X-ray emission is dominated by the densest central region of the cool core, an increase in the X-ray luminosity implies a larger accretion of gas onto the central supermassive black hole. This boost in the accretion rate consequently implies stronger AGN feedback, and this is why peaks in the X-ray luminosity closely correlate with the times when the jet power increases (see Figure 3). This cyclic behavior of the X-ray luminosity is evident in the cases including CR streaming.

Zoom In Zoom Out Reset image size

The evolution of the profiles of the ratio of heating to radiative cooling is shown in Figure 5. As in the case of the profiles of the CR pressure support shown in Figure 2, in order to exclude regions that are cooling very inefficiently, the heating-to-cooling ratio is set to 10⁻² whenever the local cooling time exceeds the Hubble time (this is done only at the post-processing stage, i.e., we do not modify the heating in the simulations). From left to right, the top row corresponds to the heating due to streaming in the case with (i) streaming heating (ST1), (ii) streaming heating and hadronic and Coulomb heating (SCHT1), and (iii) super-Alfvénic streaming heating and hadronic and Coulomb heating (SCHTs). The bottom row shows the ratio of combined Coulomb and hadronic heating to radiative cooling. Shown from left to right in the bottom row are the following cases: (i) Coulomb and hadronic heating without CR streaming transport (CHT0), (ii) Coulomb and hadronic heating with CR streaming transport (SCHT1), and (iii) same as (ii) but for super-Alfvénic CR streaming transport (SCHTs).

Zoom In Zoom Out Reset image size

Let us begin discussing Figure 5 by focusing on the bottom-left panel. This panel shows the ratio of the combined heating due to Coulomb and hadronic interactions to radiative cooling without including CR transport effects. This panel mirrors what is shown in the left panel in Figure 2 in the sense that the regions characterized by high heating-to-cooling ratios coincide with those corresponding to high CR pressure support. In this case, CR heating is very centrally concentrated as the CRs cannot easily escape to larger distances. Since the coupling of CRs to the gas at larger distances inside the cool core is very weak, gas accretion from those more distant regions onto the center is relatively unopposed, and the jet is constantly turned on injecting CRs. Consequently, CR pressure, and the associated heating, continues to build up in the center for a longer time than in the cases involving CR transport (3 Gyr rather than 2 Gyr; see, e.g., left and third panels in Figure 3, respectively). This process eventually increases the central cooling time, slows down accretion in the very central parts of the cool core, and delays the onset of the next peak in X-ray luminosity (see Figure 4). We therefore suggest that the average variability timescale of the AGN in this case is more influenced by the cooling rate of the gas farther away from the core center, where the gas cooling time is on average longer. Thus, the AGN variability is significantly weaker in this case compared to the cases with streaming included.

The non-streaming case is deceptively similar to the cases considered by Yang & Reynolds (2016a), who simulated AGN feedback using hydrodynamical simulations. The main differences between the non-streaming case presented here and their simulations is that (i) in their model, AGN jets inflate bubbles dominated by thermal energy whereas in our case the injection is dominated by CRs, and (ii) we include magnetic fields. As discussed above, for the feedback efficiency considered here, AGN jets are faster than for lower efficiencies, and the magnetic fields are not very efficient in suppressing mixing. While most of the AGN energy is in CR form, some of the energy in kinetic form quickly thermalizes, and consequently in the non-transport case the ICM is heated by a combination of hadronic and Coulomb processes and the mixing of the thermal gas. On the other hand, in the Yang & Reynolds (2016a) simulations, the heating of the ambient ICM proceeds mostly via the mixing of the thermal AGN jet fluid with the ICM.

We point out that the increase of the ICM entropy in cool cores may be dominated by CR heating rather than by, for example, turbulent dissipation. After the ICM has come into contact with CRs and experienced localized heating, it can expand locally. Such generated gas motions could eventually decay via turbulent dissipation. However, the primary heating mechanism in this case would be CR heating rather than "secondary" turbulent dissipation. However, we also note that the framework we are using does not allow for the dissipation of sound waves by conductive and viscous processes. Although these processes are likely to play an important role, too (see, e.g., Ruszkowski et al. 2004a, 2004b; Fabian et al. 2017), including these processes is beyond the scope of this paper.

Next we discuss simulations that include the effect of CR transport. In terms of CR physics, these models are somewhat similar to those of Sijacki et al. (2008), who include CRs and CR transport via diffusion rather than streaming and ignore magnetic fields. In their case, cooling catastrophe is prevented most likely as a result of more efficient mixing and diffusion of the AGN fluid containing CRs and the subsequent interactions of CRs with the ambient ICM via processes other than streaming heating. Typical patterns in the evolution of the heating-to-cooling ratios shown in Figure 5 are dramatically different when CR streaming is included, i.e., in all panels except for the bottom-left panel. It is evident that including streaming increases temporal variability in the CR heating profiles. This variable behavior also mirrors what is seen in Figure 2, which shows the evolution of the CR pressure support. Importantly, as shown in the top row, each significant AGN outburst results in CR streaming heating rates being comparable to radiative cooling throughout the cool core. When hadronic and Coulomb processes are included, streaming heating rates in the Alfvénic streaming case (upper middle panel) are comparable to those in the super-Alfvénic streaming case (upper right panel). Recall that in the latter case, CR pressure support is larger at larger radii but the CR pressure gradients tend to be weaker due to faster CR dispersal. These two competing factors result in streaming heating rates in this case being comparable to the Alfvénic streaming case.

Compared to the simulation without CR streaming, where most of the heating occurs close to the center, relatively more CR are transported to larger distances. Consequently, in the very central parts of the cool core where the cooling times are shortest, the gas can begin to accrete onto the black hole on a shorter timescale and trigger a shorter AGN outburst. This is again evident in Figure 3 where the jet switches off sooner in the cases with CR transport included. The energy injected by this outburst is again more spatially uniformly distributed when CR transport is included and, thus, subsequent outbursts also occur on shorter timescales governed by the conditions operating closer to the cool core center. In other words, more efficient use of CR energy has a clear impact on the dynamics of the atmosphere, intermittency of the AGN and, as mentioned above, should lead to more diffuse radio emission. We suggest that the reason for relatively steady jet power in the simulations of Yang & Reynolds (2016a) is that their jets were heavier and faster as they did not include CRs, which led to generally more efficient mixing throughout the cool core.

We can also compare the contributions of CR streaming heating and the combined Coulomb and hadronic losses to the total heating budget by comparing the top and bottom panels in the middle and right columns. The top panels show the contribution from the CR streaming case while the bottom ones show the contribution due to the sum of Coulomb and hadronic heating. Interestingly, it is the CR streaming heating that dominates in all cases.

In Figure 6, we show profiles of emission-weighted temperature, emission-weighted density, and emission-weighted entropy (from top to bottom, respectively; ordering of columns is the same as in Figure 1; weighting is computed using the X-ray band extending from 0.5 to 10 keV). Color-coded lines correspond to different times. These profiles are broadly consistent with observations. Most importantly, the cluster cool core is preserved despite AGN activity, i.e., the temperature and density are generally increasing and decreasing functions of radius, respectively. Consequently, the entropy profiles also demonstrate that the AGN feedback is gentle enough to preserve the positive entropy gradient in agreement with observations. Only at the very early times is there departure from this trend (see blue profiles). However, this is not unexpected as the simulations, by construction, start from a state that is out of thermodynamical equilibrium. This means that, at least initially, we do expect larger temperature variations compared to what one could have predicted starting from hydrostatic and thermal equilibrium in the initial state.

Zoom In Zoom Out Reset image size

In the case without streaming, profile variability is the weakest, but as mentioned above, this case leads to systematic and excessive build-up of CR in the cool core. The cases corresponding to Alfvénic streaming (ST1, SCHT1) exhibit only a mild level of variability. The level of variability is somewhat stronger in the super-Alfvénic streaming case (SCHTs). This is consistent with what is seen in the right panels in Figures 3 and 4 that show somewhat more variable jet power and larger variability in the X-ray luminosity of the core (see Section 3.4 for the physical interpretation of this result). As a result of the stronger feedback, temperature profiles in this case can at times approach isothermality and gas densities can be lower.

As mentioned above, we also performed runs for lower feedback efficiency and higher CR injection fraction (not shown here). In those runs, there were significant qualitative difference between the evolution of the temperature profiles in the non-streaming case and all other cases. In the non-streaming case, the temperature systematically decreased over time due to the development of global thermal instability caused by inefficient mixing of CRs with the thermal ICM and thus inefficient heating of the bulk of the ICM. Very low gas temperatures would lead to significant line emission and star formation in excess of what is observed in cool cores.