# Cosmic rays and thermal instability in self-regulating cooling flows of massive galaxy clusters

Ricarda Beckmann, Yohan Dubois, Alisson Pellissier, Valeria Olivares, Fiorella Polles, Oliver Hahn, Pierre Guillard, Matthew Lehnert

► To cite this version:

Ricarda Beckmann, Yohan Dubois, Alisson Pellissier, Valeria Olivares, Fiorella Polles, et al. Cosmic rays and thermal instability in self-regulating cooling flows of massive galaxy clusters. Astronomy and Astrophysics - A&A, 2022, 665, pp.A129. 10.1051/0004-6361/202142527. hal-03649424

HAL Id: hal-03649424

https://hal.science/hal-03649424

Submitted on 11 Apr 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Cosmic rays and thermal instability in self-regulating cooling flows of massive galaxy clusters

Ricarda S. Beckmann<sup>1,2,●</sup>, Yohan Dubois<sup>1,●</sup>, Alisson Pellissier<sup>3</sup>, Valeria Olivares<sup>4,5</sup>, Fiorella L. Polles<sup>4,6</sup>, Oliver Hahn<sup>7,8</sup>, Pierre Guillard<sup>1,9</sup>, and Matthew D. Lehnert<sup>10</sup>

Received 26 October 2021 / Accepted 23 May 2022

## ABSTRACT

One of the key physical processes that helps prevent strong cooling flows in galaxy clusters is the continued energy input from the central active galactic nucleus (AGN) of the cluster. However, it remains unclear how this energy is thermalised so that it can effectively prevent global thermal instability. One possible option is that a fraction of the AGN energy is converted into cosmic rays (CRs), which provide non-thermal pressure support, and can retain energy even as thermal energy is radiated away. By means of magneto-hydrodynamical simulations, we investigate how CR injected by the AGN jet influence cooling flows of a massive galaxy cluster. We conclude that converting a fraction of the AGN luminosity as low as 10% into CR energy prevents cooling flows on timescales of billion years, without significant changes in the structure of the multi-phase intra-cluster medium. CR-dominated jets, by contrast, lead to the formation of an extended, warm central nebula that is supported by CR pressure. We report that the presence of CRs is not able to suppress the onset of thermal instability in massive galaxy clusters, but CR-dominated jets do significantly change the continued evolution of gas as it continues to cool from isobaric to isochoric. The CR redistribution in the cluster is dominated by advection rather than diffusion or streaming, but the heating by CR streaming helps maintain gas in the hot and warm phase. Observationally, self-regulating, CR-dominated jets produce a γ-ray flux in excess of current observational limits, but low CR fractions in the jet are not ruled out.

**Key words.** galaxies: clusters: intracluster medium – methods: numerical – cosmic rays – galaxies: magnetic fields – instabilities – galaxies: jets

## 1. Introduction

From X-ray observations, the hot gas that fills galaxy clusters, the so-called intra-cluster medium (ICM), is cooling rapidly and is expected to produce a cooling flow of about 100–1000 $M_{\odot} \text{yr}^{-1}$ (Fabian 1994). However, observed star formation rates (SFRs) in clusters reach only about 1–10% of the inferred cooling rates (McDonald et al. 2018). Clusters must therefore contain a heating source that prevents over-cooling and slows star formation down while still allowing for the formation of the observed warm $H\alpha$ emission nebulae, as well as cold molecular gas and their extended filamentary morphologies (e.g. Salomé et al. 2011; Olivares et al. 2019, 2022). Research into the contribution of different heating mechanisms, from feedback driven by active galactic nuclei (AGN; e.g. Sijacki & Springel 2006; Dubois et al. 2010; Gaspari et al. 2011; Li et al. 2015; Voit et al. 2017; Beckmann et al. 2019; Ehlert et al. 2022) via

thermal conduction (e.g. Bogdanović et al. 2009; Parrish et al. 2009; Ruszkowski et al. 2011; Kannan et al. 2016) and turbulence (e.g. Kunz et al. 2011; Voit 2018) to cosmic rays (CRs; Sijacki et al. 2008; Ruszkowski et al. 2017; Ehlert et al. 2018, e.g.), has shown that long-term self-regulation of cooling flows is a complex thermodynamic problem that most likely relies on the interplay between different mechanisms that are active at the same time.

Cosmic rays in particular have the potential to significantly change the properties of the multi-phase gas in the cluster center, as they provide non-thermal pressure support and can retain energy even as thermal energy is radiated away. High-energy CRs are accelerated in strong shocks (e.g. Bell 1978) that occur in and around AGN jets (Blandford & Ostriker 1978; Blandford 2000; Matthews et al. 2020).

On large scales, it has been suggested that variations in CR transport physics are responsible for the observed dichotomy

<footnote n="1">¹ Institut d'Astrophysique de Paris, CNRS, Sorbonne Université, UMR7095, 98bis bd Arago, 75014 Paris, France
e-mail: ricarda.beckmann@iap.fr</footnote>

<footnote n="2">² Institute of Astronomy and Kavli Institute for Cosmology, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK</footnote>

<footnote n="3">³ Laboratoire Lagrange, Observatoire de la Côte d'Azur, CNRS, Université Côte d'Azur, UMR7293, Boulevard de l'Observatoire, CS 34229, 06304 Nice, France</footnote>

<footnote n="4">⁴ LERMA, Observatoire de Paris, PSL Research Univ., CNRS, Sorbonne Université, UMR8112, 75014 Paris, France</footnote>

<footnote n="5">⁵ Department of Physics and Astronomy, University of Kentucky, 505 Rose Street, Lexington, KY 40506, USA</footnote>

<footnote n="6">⁶ SOFIA Science Center, USRA, NASA Ames Research Center, M.S. N232-12, Moffett Field, CA 94035, USA</footnote>

<footnote n="7">⁷ Department of Astrophysics, University of Vienna, Türkenschanzstrasse 17, 1180 Vienna, Austria</footnote>

<footnote n="8">⁸ Department of Mathematics, University of Vienna, Oskar-Morgenstern-Platz 1, 1090 Vienna, Austria</footnote>

<footnote n="9">⁹ Institut Universitaire de France, Ministère de l'Éducation Nationale, de l'Enseignement Supérieur et de la Recherche, 1 rue Descartes, 75231 Paris Cedex 05, France</footnote>

<footnote n="10">¹⁰ Centre de Recherche Astrophysique de Lyon, ENS de Lyon, Université Lyon 1, CNRS, UMR5574, 69230 Saint-Genis-Laval, France</footnote>
of galaxy cluster profiles (Guo & Oh 2008; Enßlin et al. 2011), which appear as either cool-core or non-cool-core clusters. CR have the potential to offset cooling flows by providing a source of stable heating (Fujita & Ohira 2011; Jacob & Pfrommer 2017a), but these steady-state solutions require a CR pressure beyond current upper observational limits for clusters with large cooling radii or high SFR (Jacob & Pfrommer 2017b).

On small scales, Ruszkowski et al. (2018) showed that CR could provide the missing heating source that maintains observed Hα filaments in galaxy clusters at their observed temperature of $10^4$ K. Using linear stability analysis, Kempski & Quataert (2020) showed that the presence of CRs modifies local thermal instability criteria, while Butsky & Quinn (2018) studied how the presence of CRs influences the morphology of condensed gas. Huang et al. (2022) showed that the presence of CRs reduces the density contrast of thermally unstable gas.

Early simulations of isolated galaxy clusters by Sijacki et al. (2008) showed that CRs injected by AGN, even in the absence of CR transport and heating mechanisms beyond advection, allow clusters to have cool cores without strong cooling flows merely by providing non-thermal pressure support. More complete recent simulations that include CR transport processes (including diffusion and streaming), CR heating, and a variety of different CR injection mechanisms from the AGN jet (Ruszkowski et al. 2017; Wang et al. 2020; Su et al. 2021) to volume-filling injection (Su et al. 2020), showed that CRs with streaming-instability heating improve cluster self-regulation and allow for the formation of cold gas without a run-away cooling catastrophe. Su et al. (2021) showed that a large amount of CR injection leads to over-quenching, in which cluster cooling is completely suppressed over long periods of time, in contrast to observations. Ji et al. (2020) used cosmological simulations to show that the impact of CRs is not only limited to the cluster centre, but also produces warm volume filling gas at large radii from the centre.

One piece of observational evidence for CR playing a part in galaxy cluster cooling cycles comes from the observations of Faranoff-Riley I (FRI) jet lobes. To be in pressure equilibrium with the surrounding ICM, FRI jet lobes need a dominant non-thermal pressure component (Morganti et al. 1988; Worrall 2009; Croston et al. 2008), with CRs being the most likely candidate (Croston et al. 2018). Another way to observationally search for CRs is via the γ-ray emission produced by the decay of pions. So far, only upper limits have been detected for such pionic γ rays from the ICM of galaxy clusters. In the Coma cluster, combined γ-ray and radio observations place the maximum ratio of CR pressure to thermal pressure, $\eta = P_{CR}/P_{therm} < 0.1$ for both hadronic and reacceleration models (Ackermann et al. 2010, 2014, 2016; Aleksić et al. 2012; Brunetti et al. 2012, 2017). For the Perseus cluster, constraints are even tighter, with $\eta < 0.02$ if the CR density profile is the same as the gas pressure profile. However, if the CR density profile is overall flatter than the gas pressure profile, that is, if CRs preferentially escape from the cluster centre, $\eta < 0.2$ is supported by the data (Ahnen et al. 2016). Finally, CRs can be detected in galaxy clusters via the Sunyaev-Zeldovich effect (Abdulla et al. 2019; Ehlert et al. 2019).

We investigate how CRs injected by the AGN jet in the cluster centre influence the onset of thermal instability in the multi-phase ICM and the cooling flow of a massive galaxy cluster. The paper is structured as follows: A theoretical discussion of the expected impact of CRs on cluster cooling flows is presented in Sect. 2, while a comparative analysis using simulations

(parameters shown in Sect. 3) is presented in Sect. 4. The impact of CRs on the multi-phase gas and local thermal instability is shown in Sect. 6. An observational comparison is presented in Sect. 7, followed by our conclusions (Sect. 8). The robustness of the results with respect to small perturbations is investigated in Appendix B.

## 2. Thermal instability in clusters

### 2.1. Review

To avoid massive cooling flows, clusters need to be in global hydrostatic equilibrium. However, clusters are not described by perfectly uniform density profiles, but instead are turbulent, with local variations in density, temperature, and velocity. This leads to pockets of local thermal instability in which the hot ICM is able to cool and condense.

In the absence of magnetic fields, gas becomes locally unstable when the local cooling time

$$t_{\text{cool}} = \frac{2 n k_{\text{B}} T}{3 n_e n_i \Lambda} \qquad (1)$$

is shorter than the free-fall time

$$t_{\text{ff}} = \sqrt{\frac{2r}{g}}, \qquad (2)$$

where n is the gas number density, $n_e$ and $n_i$ are the electron and ion number density, respectively, $k_B$ is the Boltzmann constant, T is the gas temperature, Λ is the cooling function, g is the local gravitational acceleration, and r is the radius from the cluster centre.

Using linear perturbation theory, McCourt et al. (2012) showed that local thermal instabilities occur when $f_{ff} = t_{\text{cool}}/t_{\text{ff}} \le 1$ if there is a global heating function that ensures global thermal equilibrium. Simulations have generally found that $f_{ff}$ can be as high as $f_{ff} = 10-25$ for spherically averaged profiles (Sharma et al. 2012; McCourt et al. 2012; Beckmann et al. 2019). This is due to the dispersion of entropy perturbations (Voit 2021), which lead to a globally elevated $t_{ff}$ despite the fact that $f_{ff} = 1$ locally (Pal Choudhury et al. 2019; Voit 2021). Observationally, it has been confirmed that dense gas is located at the minima of the radial profile of $f_{ff}$ (Hogan et al. 2017; Pulido et al. 2018; Olivares et al. 2019, 2022).

Another approach to determining local thermal instability is to follow Gaspari et al. (2018) and consider the local eddy turnover time,

$$t_{\text{eddy}} = 2\pi \frac{r^{2/3} L^{1/3}}{\sigma_v}, \qquad (3)$$

in comparison to the cooling time, where $\sigma_v$ is the velocity dispersion at the injection scale L. This criterion is based on the assumption that the kinematics of the warm ionized gas correlates linearly with those of the hot medium, following the development of a multi-phase cascade of turbulent energy. However, there have recently been insights from observations (Li et al. 2020) and simulations (Wang et al. 2021; Mohapatra et al. 2022) that this assumption might not hold in galaxy clusters.

In the presence of non-thermal energies beyond turbulence, such as when magnetic fields or CR are taken into account, the thermal instability changes as the non-thermal energies can provide additional pressure support against collapse. CR energy also is not as susceptible to efficient radiative cooling as thermal energy.
Studying the impact of magnetic fields on thermal instability, Ji et al. (2017) used idealised numerical experiments similar to those used in McCourt et al. (2012) to quantify the impact of magnetic fields on thermal instability. Rather than define a clear threshold for $f_{ff}$, Ji et al. (2017) took a more continuous approach and measured the magnitude of the density fluctuations, $\delta\rho/\rho$ as a function of $f_{ff}$. For the hydro-only case, they reported

$$ \left. \frac{\delta \rho}{\rho} \right|_{\text{hydro}} = 0.1 f_{ff}^{-1}. \qquad (4) $$

In the presence of magnetic fields, this becomes

$$ \left. \frac{\delta \rho}{\rho} \right|_{\text{MHD}} = 3\beta^{-1/2} f_{ff}^{-1}, \qquad (5) $$

where $\beta = 8\pi P_{\text{therm}}/B^2$ is the plasma beta and $P_{\text{therm}}$ is the thermal pressure. Comparison of the two shows that $\left. \frac{\delta\rho}{\rho} \right|_{\text{MHD}} > \left. \frac{\delta\rho}{\rho} \right|_{\text{hydro}}$ for fields as weak as $\beta \le 1000$. Therefore even weak magnetic fields produce stronger density fluctuations in the gas.

Adding CRs is somewhat more complicated, as they bring a range of energy transfer mechanisms that have competing effects on thermal instability. Kempski & Quataert (2020) conducted a linear stability analysis of a plasma in the presence of CRs and concluded that the unstable regime depends on the relative importance of different CR transport processes, such as diffusion and streaming, and on the importance of CR heating. Results depend on three parameters:

1. the slope of the cooling function $\mathcal{L}_{\text{th}}$ with respect to temperature $T$,

$$ \alpha_{\text{T}} = \frac{\partial \log \mathcal{L}_{\text{th}}}{\partial \log T}, \qquad (6) $$

which for the simulations presented here is shown in Appendix A;

2. the ratio of CR pressure $P_{\text{CR}}$ and thermal pressure $P_{\text{therm}}$:

$$ \eta = \frac{P_{\text{CR}}}{P_{\text{therm}}}; \qquad (7) $$

3. and the ratio

$$ \chi = \frac{D_{\text{CR}}}{t_{\text{cool}} \eta u_{\text{A}}^2}, \qquad (8) $$

where $D_{\text{CR}}$ is the CR diffusivity, and $u_{\text{A}}$ is the Alfvén velocity, $\chi$, which measures the relative strength of the cosmic ray heating rate to the radiative thermal cooling rate, encapsulates the complex interplay between CR diffusion, cooling, and the local energy distribution, which all play a part in determining the local thermal (in)stability of gas in the presence of CRs.

In general, the role of CRs in thermal instability depends on the value of $\chi$. When $\chi > 1$, thermal instabilities develop with an amplitude proportional to $2 - \alpha_{\text{T}}$, that is, thermal instabilities grow strongly when $\alpha_{\text{T}} < 2$. In the opposite regime, where $\chi \le 1$, CR diffusion suppresses thermal instability below a maximum length-scale, which can be expressed as an effective field length of the form

$$ \lambda_{\text{CRF}} = \begin{cases} 2\pi|\mathbf{b} \cdot \tilde{\mathbf{k}}| \sqrt{\eta D_{\text{CR}} t_{\text{cool}}} & \eta \le 1, \\ 2\pi|\mathbf{b} \cdot \tilde{\mathbf{k}}| \sqrt{\frac{D_{\text{CR}} t_{\text{cool}}}{\eta}} & \eta > 1, \end{cases} \qquad (9) $$

where $\mathbf{b}$ is the magnetic field unit vector, and $\tilde{\mathbf{k}}$ is the wave number unit vector. The thermal instability continues to occur on larger scales. $\lambda_{\text{CRF}}$ is at a maximum when $\eta = 1$, and it falls off for higher and lower values.

<figure><img src="image_2.png" /><figcaption>Fig. 1. Radial profile of $\chi\eta$ (top panel), $\eta_{\text{crit}}$, the minimum $\eta$ required for the strong diffusion regime (middle panel), and the maximum $\lambda_{\text{CRF}}$ (bottom panel) for the galaxy clusters of different masses in hydrostatic equilibrium. Different colours stand for different halo masses as indicated in the top inset, either with a cored (solid lines) or with a non-cored profile (dashed lines). Grey horizontal lines demarcate the two regimes in $\chi$ (top panel) and gas that are CR-pressure dominated vs thermal pressure dominated (middle panel). $\eta_{\text{crit}}$, the minimum value required for $\chi < 1$ increases with radius. By contrast, CR pressure profiles are expected to decrease with radius. It is therefore unlikely that $\chi < 1$ at large radii, which means that CRs are unable to stabilise small-scale instabilities outside the cluster core.</figcaption></figure>

## 2.2. Expected behaviour for galaxy clusters

In this section, we predict $\eta_{\text{crit}}$, the value of $\eta$ required to have $\chi = 1$ for galaxy clusters of different masses in order to discuss the potential impact of CRs on thermal instability within different regions of the cluster. In galaxy clusters, temperature, density, and magnetic field strength vary strongly with radius, as does CR pressure. As a consequence, the $\eta$ required to have $\chi = 1$ is expected to be a function of radius. In turn, this means that the impact of CRs on thermal instability will vary throughout the cluster. We would expect CRs to be unable to suppress thermal instability in massive galaxy clusters if $\chi > 1$, as hot gas in galaxy clusters cools predominantly through free-free emission, which has $\alpha_{\text{T}} = 0.5 < 2$. By contrast, if $\chi < 1$ for a significant volume of the cluster, there could be a maximum length scale, the CR Field length $\lambda_{\text{CRF}}$ (see Eq. (9)), below which CRs are able to suppress thermal instability. Instability would be expected to continue on scales larger than $\lambda_{\text{CRF}}$. As both $\chi$ and $\lambda_{\text{CRF}}$ explicitly depend on $\eta$, the impact of CRs on thermal instability in galaxy clusters depends very strongly on the radial distribution of CRs, as expected.

In Fig. 1 we show example profiles of $\chi\eta$ for galaxy clusters of different masses in hydrostatic equilibrium. To test the
impact of different distributions of CRs on the cluster, we computed $\chi\eta$ and then determined $\eta_{crit}$, which will give $\chi < 1$ if $\eta > \eta_{crit}$ (or $\chi > 1$ if $\eta < \eta_{crit}$). Profiles to compute $\chi\eta$ were computed in the same way as the initial conditions for the simulations presented in this paper (see Sect. 3.3). Like for the initial conditions, density profiles are based on a cored Navarro-Frenk-White (NFW) profile (Navarro et al. 1997). The concentration parameter, gas fraction, and black hole mass are scaled with cluster mass according to Mattio et al. (2007), Andreon et al. (2017), and Bogdán et al. (2018), respectively. The central magnetic field was set to 10 µG (e.g. Govoni & Feretti 2004), and was scaled with $\rho(r)^{2/3}$. For each cluster, a cored (solid) and a non-cored (dashed) profile is shown. The solid orange lines in Fig. 1 correspond to the profiles we used as initial conditions.

Figure 1 shows that the profile of $\chi\eta$ varies strongly with galaxy cluster mass, thermal profile (cored or non-cored), and radius. More massive clusters need higher CR pressure to reach $\chi < 1$, which requires $\eta > \eta_{crit}$ (middle panel) as $\chi \propto \frac{1}{\eta}$. The required value of $\eta_{crit}$ increases with radius for all cluster masses.

One source of CR, injection via supernovae and AGN jets, is concentrated in the cluster centre and would cause $\eta$ to decrease with radius. From such injections, it is likely that only the cluster centre, or no part of the cluster, has $\chi < 1$, where CRs can potentially stabilise against thermal instability. Another potential source of CR injection are large-scale structure formation shocks on megaparsec scales (Ryu et al. 2003; Pfrommer et al. 2006; Vazza et al. 2009). However, energy ratios of CR energy and thermal energy for these circum-cluster shocks are low because strong shocks are rare in and around galaxy clusters (Vazza et al. 2009), and the injection scales are very far from the cluster centre. It is therefore unlikely that these CRs make a significant difference to $\eta$ in the cluster centre. Outside regions where $\chi < 1$, even in the presence of CRs, thermally unstable gas will remain, so that even in the presence of CRs as $\alpha_T = 0.5 < 2$ for the free-free emission that is the dominant cooling channel for hot cluster gas.

Inside this potentially CR-stabilised core, CRs can suppress thermal instability on scales below $\lambda_{CRF}$ (see Eq. (9)). The bottom panel of Fig. 1 shows that the maximum length scale below which CRs can suppress thermal instability in a massive cluster is about 10 kpc or more, which occurs if $\eta = 1$. $\lambda_{CRF}$ decreases for higher and lower values of $\eta$, and for dimensions that are not parallel with the magnetic field lines. For a Perseus-like cluster, 10 kpc of $\lambda_{CRF}$ is comparable to the observed length of filaments, but much larger than their 10 pc width (Conselice et al. 2001; Fabian et al. 2016). This lends credibility to the idea that observed filaments might be aligned with magnetic field lines.

As $\eta_{crit} > 1$ for massive galaxy clusters at all radii, there will be gas that is dominated by CR pressure, but in which CRs will be unable to suppress thermal instability, that is, gas for which $1 < \eta < \eta_{crit}$. However, CRs influence not only the local thermal instability, but also the thermodynamic evolution of the gas after the onset of thermal instability. Gas with $\eta > 1$ cools isochorically rather than isobarically. As gas cools, $P_{therm}$ decreases, thus, $\eta = P_{CR}/P_{therm}$ increases even if $P_{CR}$ is constant.

Overall, we conclude that CRs might be able to suppress thermal instability on relevant length scales in the cluster core, but are unlikely to be able to do so outside this central region even for gas dominated by CR pressure. Understanding where and how CRs are able to change the thermal instability criterion in massive galaxy clusters, and how the evolution of this newly condensed gas differs in the presence of CRs from the CR-free case is the aim of this paper.

# 3. Simulation

## 3.1. Magneto-hydrodynamics with cosmic rays

This paper presents a set of hydrodynamical and magneto-hydrodynamical (MHD) simulations of isolated galaxy clusters with and without CRs. The cluster simulations were run with the adaptive mesh refinement code RAMSES (Teyssier 2002) that solves for the MHD (Fromang et al. 2006) with CRs and two ion-electron temperatures (Dubois & Commerçon 2016),

$$ \frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \mathbf{u}) = 0, \qquad (10) $$

$$ \frac{\partial \rho \mathbf{u}}{\partial t} + \nabla \left( \rho \mathbf{u} \mathbf{u} + P_{\text{tot}} + \frac{B \mathbf{B}}{4\pi} \right) = 0, \qquad (11) $$

$$ \frac{\partial e}{\partial t} + \nabla \cdot \left( (e + P_{\text{tot}}) \mathbf{u} - \frac{\mathbf{B}(\mathbf{B} \cdot \mathbf{u})}{4\pi} \right) = \mathcal{L}_{\text{rad}} - \nabla \cdot \mathbf{F}_{\text{st}} - \nabla \cdot \mathbf{F}_{\text{d}} - \nabla \cdot \mathbf{F}_{\text{cond}}, \qquad (12) $$

$$ \frac{\partial \mathbf{B}}{\partial t} - \nabla \times (\mathbf{u} \times \mathbf{B}) = 0, \qquad (13) $$

$$ \frac{\partial e_{\text{CR}}}{\partial t} + \nabla \cdot (e_{\text{CR}} \mathbf{u}) = -P_{\text{CR}} \nabla \cdot \mathbf{u} + \mathcal{L}_{\text{CR}} - \nabla \cdot \mathbf{F}_{\text{st}} - \nabla \cdot \mathbf{F}_{\text{d}} - \mathbf{u}_{\text{st}} \cdot \nabla P_{\text{CR}}, \qquad (14) $$

$$ \frac{\partial e_e}{\partial t} + \nabla \cdot (e_e \mathbf{u}) = -P_e \nabla \cdot \mathbf{u} - \nabla \cdot \mathbf{F}_{\text{cond}} + Q_{e \leftrightarrow i} + \mathcal{L}_{\text{th}}, \qquad (15) $$

where $\rho$ is the mass density; $\mathbf{u}$ is the velocity; $e = 0.5 \rho \mathbf{u}^2 + e_e + e_i + e_{CR} + \mathbf{B}^2/8\pi$ is the total energy density; $e_e, e_i$, and $e_{CR}$ are the electron, ion, and CR energy densities, respectively; $P_{\text{tot}} = P_e + P_i + P_{CR}$ is the total gas pressure; $P_e, P_i$, and $P_{CR}$ are the electron, ion, and CR pressures, respectively; and $\mathbf{B}$ is the magnetic field. Each pressure component is related to its corresponding energy density assuming an adiabatic equation of state with adiabatic index $\gamma_e = \gamma_i = \gamma = 5/3$ and $\gamma_{CR} = 4/3$.

In the CR energy equation (Eq. (14)), $F_d = -D_{CR}b(\mathbf{b} \cdot \nabla)e_{CR}$ is the anisotropic diffusion flux (where $b$ is the magnetic field unit vector), $F_{\text{st}} = f_{b,\text{st}}\mathbf{u}_{\text{st}}(e_{CR} + P_{CR})$ is the streaming advection term, and $-\mathbf{u}_{\text{st}} \cdot \nabla P_{CR}$ is the loss term due to Alfvén wave damping through streaming of the CR energy (directly transferred into the thermal ion pool of energy). The streaming velocity is $\mathbf{u}_{\text{st}} = -\text{sign}(\mathbf{u}_A \cdot \nabla e_{CR})\mathbf{u}_A$, with $\mathbf{u}_A = \mathbf{B}/\sqrt{4\pi\rho}$ the Alfvén velocity, and where $f_{b,\text{st}}$ is a streaming velocity boost term that depends on the strength of the turbulent and non-linear Landau wave damping in the hot ICM (Wiener et al. 2013). This is typically about unity for galaxy clusters (Ruszkowski et al. 2017). For the simulations presented here, $f_{b,\text{st}} = 1$. We set the CR diffusion coefficient to $D_{CR} = 10^{29} \text{ cm}^2 \text{ s}^{-1}$, which is the typical value obtained in simulated galaxies to match the observed $\gamma$-ray flux (Chan et al. 2019; Hopkins et al. 2021), or obtained with detailed models of CR propagation in the Milky Way (e.g. Trotta et al. 2011). To ensure that our results do not sensitively depend on this choice, we test other values of $D_{CR}$ and $f_{b,\text{st}}$ in Appendix C.

The evolution of the CR energy density can be more accurately described by taking the first two moments of the Vlasov equation, which augments the description of the CR evolution with a partial differential equation on the time evolution of the CR flux (as in Jiang & Oh 2018, and in the more complete derivation of Thomas & Pfrommer 2019, who also evolved the energy density of Alfvén waves with their damping processes). In the limit of the steady state for the CR flux (i.e. when the
speed of light tends towards infinity), the two-moment method is equivalent to the standard approach considered here (and commonly employed in the literature; e.g. Ruszkowski et al. 2017; Girichidis et al. 2018; Holguin et al. 2019; Dashyan & Dubois 2020; Buck et al. 2020; Semenov et al. 2021). While these two-moment methods are appealing and show different behaviours in some tailored test problems (see e.g. the CR over-density test case on top of a strong background CR density of Thomas & Pfrommer 2019), they usually rely on a reduced value of the speed of light for the simulations to remain feasible, which needs to remain sufficiently high to guarantee convergence to the correct solution. Nonetheless, Chan et al. (2019) have shown in the context of Milky Way-like simulations that the two methods lead to essentially equivalent results.

Electron energy was treated in a separate equation (Eq. (15)) from the total energy (Eq. (12)) in order to take local thermodynamical equilibrium (LTE) effects into account. The term $F_{\text{cond}} = (-\kappa_e \mathbf{b}(\mathbf{b} \cdot \nabla)T_e)$ is the anisotropic conductive flux, where

$$ \kappa_e = f_{\text{sat}}\kappa_{\text{sp}} = f_{\text{sat}}n_e k_B D_{\text{cond}} \quad (16) $$

is the Spitzer conductivity for electrons with $n_e = \rho/(\mu_e m_p)$ (and $n_i = \rho/(\mu_i m_p)$) the electron (ion) number density, $\mu_e$ and $\mu_i$ the electron and ion mean molecular weights, $m_p$ the proton mass, $k_B$ the Boltzmann constant, and $D_{\text{cond}}$ the thermal diffusivity, which is equal to

$$ D_{\text{cond}} = 1.31 \lambda_e \left( \frac{k_B T_e}{m_e} \right)^{\frac{1}{2}} \\ \approx 8.3 \times 10^{30} \left( \frac{T_e}{10^8 \text{ K}} \right)^{\frac{5}{2}} \left( \frac{n_e}{10^{-2} \text{ cm}^{-3}} \right)^{-1} \text{cm}^2 \text{s}^{-1}, \quad (17) $$

with the mean free path of electrons

$$ \lambda_e = \frac{3^{3/2}(k_B T_e)^2}{4\pi^{1/2}n_e q_e^4 \ln \Lambda} \\ \approx 2.1 \left(\frac{T_e}{10^8 \text{ K}}\right)^2 \left(\frac{n_e}{10^{-2} \text{ cm}^{-3}}\right)^{-1} \text{kpc}, \quad (18) $$

where $m_e$ is the mass of the electron, $q_e$ the charge of the electron, and $\ln \Lambda \approx 40$ the Coulomb logarithm. This conductive flux saturates when the characteristic scale length of the gradient of temperature $\ell_{T_e} = T_e/\nabla T_e$ is comparable to or shorter than the mean free path of electrons (Cowie & McKee 1977). We followed Sarazin (1986) and introduced an effective conductivity that approximates the solution in the unsaturated and saturated regime by

$$ f_{\text{sat}} = \frac{1}{1 + 4.2\lambda_e/\ell_{T_e}}. \quad (19) $$

The electron energy equation has an additional term $Q_{e \leftrightarrow i}$ that transfers heat between electrons and ions and is responsible for restoring LTE,

$$ Q_{e \leftrightarrow i} = \frac{T_i - T_e}{\tau_{\text{eq,ei}} \gamma - 1}, \quad (20) $$

with the equilibrium timescale

$$ \tau_{\text{eq,ei}} = \frac{3m_e m_p}{8(2\pi)^{1/2}n_i q_e^4 \ln \Lambda} \left(\frac{k_B T_e}{m_e}\right)^{\frac{3}{2}} \\ \approx 20 \left(\frac{T_e}{10^8 \text{ K}}\right)^{\frac{3}{2}} \left(\frac{n_i}{10^{-2} \text{ cm}^{-3}}\right)^{-1} \text{Myr}. \quad (21) $$

The radiative loss term for the total energy density $\mathcal{L}_{\text{rad}} = \mathcal{L}_{\text{th}} + \mathcal{L}_{\text{CR}} + \mathcal{H}_{\text{CR}\to\text{th}}$ contains the gas radiative loss term $\mathcal{L}_{\text{th}}$ (see Sect. 3.5), and the CR radiative loss term due to Coulomb and hadronic collisions (Guo & Oh 2008) $\mathcal{L}_{\text{CR}} = -7.51 \times 10^{-16}(n_e/\text{cm}^{-3})(e_{\text{CR}}/\text{erg cm}^{-3})\text{erg s}^{-1}$ reprocessed into the thermal pool at a rate $\mathcal{H}_{\text{CR}\to\text{th}} = 2.63 \times 10^{-16}(n_e/\text{cm}^{-3})(e_{\text{CR}}/\text{erg cm}^{-3})\text{erg s}^{-1}$.

The induction equation was solved with the constrained transport scheme (Teyssier et al. 2006) that guarantees $\nabla \cdot \mathbf{B} = 0$ at machine precision. The MHD system of equations was first solved by zeroing all source terms on the right-hand side of the equations using the MUSCL-Hancock scheme (Fromang et al. 2006) with linear reconstruction of the conservative variable, a minmod total variation diminishing scheme, and the HLLD approximate Riemann solver (Miyoshi & Kusano 2005). For the hydrodynamics-only simulations we used for comparison, the HLLC Riemann solver (Toro 2009) was used instead. For the anisotropic thermal conduction and CR diffusion, their corresponding fluxes were solved with an implicit method (Dubois & Commerçon 2016), and similarly for the advection part of the streaming, where it was treated just as a diffusion-like term (Dubois et al. 2019). As in Dashyan & Dubois (2020), the diffusion solver included a minmod slope limiter on the transverse component of the face-oriented flux that preserved the positivity of the solution (Sharma & Hammett 2007). The solution of the temperature coupling was obtained with an implicit update (Dubois & Commerçon 2016).

## 3.2. Refinement

The simulations were performed in a box size of 8.7 Mpc with a root grid of 64³ cells (refinement level 6) and a corresponding coarse resolution of 136 kpc, which was adaptively refined by eight more levels up to a maximum resolution of $\Delta x = 531$ pc (level 14). Refinement proceeded according to one of three refinement criteria: (i) Cells were refined when their gas mass exceeded [27 089, 8713, 4098, 1621, 461, 152, 59, 12, 12] × 1.47 × $10^6 M_\odot$ for levels 6 to 14. (ii) All cells within a sphere of radius $4\Delta x$ radius from the black hole (BH) were maximally refined at all times. (iii) To avoid the common problem of mass-based refinement criteria, which derefine and disperse low-density jet bubbles, we ensured that the AGN jet and its cavities remained refined by using a dedicated passive scalar. This scalar was introduced to the simulation by setting $f_{\text{jet}} = 1$ for all cells accelerated at the base of the jet (see Sect. 3.6 for details). It was then advected with the flow and mixed into the background medium as the jet evolved. All cells whose value of the passive scalar injected at the AGN jet base exceeded $f_{\text{jet}} = \rho_{\text{scalar}}/\rho_{\text{gas}} > 10^{-4}$ were refined when the local gradient of this passive scalar exceeded 10%. To avoid filling the entire box with jet scalar over time and to ensure that the highest value of $f_{\text{jet}}$ always marked the jet nozzle and recent jet bubbles, the passive jet scalar exponentially decayed with a decay timescale of 10 Myr. The combination of these refinement criteria ensured that dense collapsing regions and regions that were recently affected by AGN feedback, including hot low-density bubbles that would derefine under a purely Lagrangian refinement scheme, remained well refined throughout the simulation.

## 3.3. Initial conditions of gas and dark matter

The cluster was initialised as a cored NFW profile with a total mass of
<figure><img src="image_3.png" /><figcaption>Fig. 2. Ratio of the cooling time to the free-fall time for the initial conditions of the cluster for various values of the core radius $r_{\text{core}}$ assuming an BH with a mass of $M_{\text{BH}} = 1.6 \times 10^{10} M_{\odot}$. Without the mass of the BH, the initial cluster profile dips lower in the centre for the chosen core radius of 13 kpc (red line).</figcaption></figure>

$$ \rho_{\text{NFW}} = \rho_s \frac{r_s^3}{(r+r_{\text{core}})(r+r_s)^2}, \qquad (22) $$

where $r_s = r_{200}/c$ is the scale radius, $r_{\text{core}} = 13$ kpc is the core radius, $\rho_s = \rho_c \delta_{200}$ is the density scaling of the profile, with the rescaling factor $\delta_{200} = 200 \Omega_m c^3 / (3 f(c))$, where $f(c) = \ln(1+c) - c / (1+c)$, and $\rho_c$ is the critical density of the Universe. The total mass of the halo was split into a fixed dark matter (DM) profile, which has the form $\rho_{\text{DM}} = (1 - f_{\text{gas}}) \rho_{\text{NFW}}$, and a gas profile of the form $\rho_{\text{gas}} = f_{\text{gas}} \rho_{\text{NFW}}$, where $f_{\text{gas}}$ is the initial gas fraction. The DM profile was fixed and did not evolve throughout the simulation, whereas the gas profile did evolve under gravity, cooling, and the influence of the AGN jet.

For the cluster presented here, which had a total mass of $8 \times 10^{14} M_{\odot}$, the parameters in the NFW profile were as follows: The concentration parameter was chosen to be $c_{200} = 4.41$ based on the halo mass according to Eq. (9) of Maccio et al. (2007) for relaxed halos. The combined choice of halo mass and concentration parameter leads to a scale radius $r_s = 434$ kpc, a virial radius of $r_{200} = 1917$ kpc, and a virial velocity of $v_{200} = 1344 \text{ km s}^{-1}$. The gas fraction was chosen according to Eq. (5) of Andreon et al. (2017), based on the cluster mass. For the cluster presented here, it is equal to $f_{\text{gas}} = 0.103$. CR pressure was initialised based on the gas pressure profile using a low and constant value of $\eta = 10^{-4}$.

To ensure both local and global thermal stability in the initial conditions, the core radius was chosen such that the ratio of the cooling time $t_{\text{cool}}$ was greater than the free-fall time $t_{\text{ff}}$ in the cluster centre, which for the current parameters leads to a core radius of $r_{\text{core}} = 134$ kpc. As can be seen by comparing the dashed red line ($r_{\text{core}} = 0$ kpc and $M_{\text{BH}} = 0 M_{\odot}$) to the dotted green lines ($r_{\text{core}} \neq 0$ kpc and $M_{\text{BH}} = 1.6 \times 10^{10} M_{\odot}$) in Fig. 2, the gravitational potential of the BH is crucial in producing an upturn in $t_{\text{ff}}/t_{\text{cool}}$ in the cluster centre. Without it, extremely large core radii would be needed to meet the thermal stability criterion.

The cluster profile was truncated at the virial radius. Outside of $r_{200}$, the density scales as

$$ \rho(r > r_{200}) = \rho(r_{200}) \left(\frac{r_{200}}{r}\right)^5, \qquad (23) $$

where $\rho(r_{200})$ is the density at $r_{200}$ according to Eq. (22).

Avoiding a sharp density cutoff at the outer cluster edge has several advantages. It increases stability in the cluster outskirts, avoiding large-scale bulk flows due to poorly resolved hydrostatic equilibrium in areas of rapid density transition, and it considerably speeds up calculations by reducing conductive heat flows because sharp temperature transitions are smoothed. As this study is focused on the cluster centre, the choice of parameters for the cluster edge has a negligible impact on the results of this study, but it makes an important difference to the overall computational resources required, with parameters here chosen to minimise computational requirements.

Ions and electrons were initial in LTE. The initial ion pressure was set to

$$ P_i = P_e \frac{\mu_i}{1 + \mu_e}, \qquad (24) $$

where the mean molecular weights were $\mu_i = 1.22$ and $\mu_e = 1.14$ for a combined gas mean molecular weight of $\mu_{\text{gas}} = 0.59$. Hydrostatic equilibrium in the initial conditions was calculated so that the total (thermal, CR, and magnetic) pressure gradient balanced the gravitational pull.

### 3.4. Initial conditions of the magnetic field

All magnetised simulations were initialised with a three-dimensional randomly tangled magnetic field on a 512³ grid. In order to guarantee $\nabla \cdot \mathbf{B} = 0$, we set up a random magnetic potential vector $\mathbf{A}$ at each cell corner, and reconstructed the $\mathbf{B}$ field at the centre of cell faces by taking the curl of $\mathbf{A}$. This initial magnetic field was then resampled onto the initial adaptive grid using a procedure that enforced $\nabla \cdot \mathbf{B} = 0$ at machine precision even at refined interfaces. For the box size of 8.7 Mpc used here, this means that the magnetic field was initialised with a characteristic coherence length of 17 kpc. As the simulation progressed, this coherence length evolved in accordance with local gas flows.

The magnetic field was initialised with a central magnetic field strength of $B_0 = 20 \mu\text{G}$ and was then scaled with the initial density profile as $B(r) = B_0(\rho(r)/\rho_0)^{2/3}$. Due to the tangled nature of the initial field, numerical magnetic reconnection took place, which somewhat decreased the average magnetic field strengths throughout the simulation.

### 3.5. Cooling, metallicity, and star formation

Radiative cooling of the gas was calculated according to the tabulated values of Sutherland & Dopita (1993) for temperatures above $T > 10^4$ K, with values extended below $10^4$ K using the fitting functions from Rosen & Bregman (1995). For metal cooling, the local metallicity value of each cell was taken into account, with solar abundance ratios of elements assumed throughout.

Metals were treated as a single species and were advected as a passive scalar. The cluster was initialised with a metallicity profile as shown in Fig. 3, which was computed according to

$$ Z = \min \left( 0.45, \max \left( 0.22, 0.15 \left( \frac{r}{r_{200}} \right)^{-0.28} \right) \right), \qquad (25) $$

with the upper limit and fit taken from Leccardi & Molendi (2008), and the lower limit being the average value of data
<figure><img src="image_4.png" /><figcaption>Fig. 3. Metallicity profile used for the cluster, based on observational data by Leccardi & Molendi (2008) and Urban et al. (2017) in units of solar metallicity according to Anders & Grevesse (1989).</figcaption></figure>

reported in Urban et al. (2017). The metallicity profile we used is therefore flat in the cluster centre, falls off to about 0.3 $r_{200}$, and is flat outside this region. Gas was further metal-enriched throughout the simulation due to supernova explosions.

Stars were formed according to combined density and temperature criteria, with star formation permitted in cells with a hydrogen number density of $n_H > 0.1 \text{ H cm}^{-3}$ and a temperature $T < 10^4 \text{ K}$. This led to a stellar mass resolution of $m_{*,\text{min}} = n_H m_p \Delta x^3 / X_H = 3.89 \times 10^5 \text{ M}_\odot$, where $m_p$ is the proton mass, and $X_H = 0.74$ is the fractional abundance of hydrogen. Star formation produces stellar particles, randomly drawn with a Poisson process (Rasera & Teyssier 2006), following a Schmidt law $\dot{\rho}_* = \epsilon_* \rho / t_{\text{ff}}$, where $\rho$ is the gas density, $t_{\text{ff}}$ is the free-fall time of the gas cell, and $\epsilon_* = 0.1$ is the constant efficiency of star formation.

Stellar feedback was included in the form of type II supernovae only, using the energy-momentum model of Kimm et al. (2015). Each stellar particle releases a total energy of $\epsilon_{*,\text{SN}} = m_* \eta_{\text{SN}} 10^{50} \text{ erg M}_\odot^{-1}$ in one feedback event after 10 Myr, where $\eta_{\text{SN}} = 0.2$ corresponds to the mass fraction of the initial mass function for stars ending their life as type II supernovae, and $m_*$ is the stellar particle mass. These explosions also enriched the gas locally with metals with a constant yield of 0.1.

### 3.6. Active galactic nucleus

To model the self-regulating AGN in the centre of the cluster, an BH was placed at the centre of the cluster as part of the initial conditions. This BH had a mass of $M_{\text{BH}} = 1.65 \times 10^{10} \text{ M}_\odot$, following the $M_{\text{BH}}-M_{500}$ relation in Phipps et al. (2019).

#### 3.6.1. Dynamics

The BH was modelled as a RAMSES sink particle, and it was free to move within the gravitational potential of the halo. To

compensate for unresolved dynamical friction, we employed two related sub-grid algorithms: we calculated the dynamical friction due to gas according to Ostriker (1999), using the high-resolution correction from Beckmann et al. (2018) to prevent the dynamical friction sub-grid algorithms from ejecting the BH when local overdensities are resolved, and the dynamical friction due to the stars (Pfister et al. 2019).

As the initial conditions did not include the stellar component of the central galaxy and the DM halo has a strong core, the gravitational potential in the cluster centre was very shallow at the beginning of the simulation. This caused the BH to wander occasionally from the centre of the galaxy cluster. To diminish this effect, we added a gradient descent correction of the form

$$ \mathbf{x}_{n+1} = \mathbf{x}_n - \tilde{\gamma}_d \mathbf{f}(\mathbf{x}_n) \quad (26) $$

to the BH trajectory according to Pellissier (in prep.). Here, $\mathbf{x}_n$ is the position of the BH at step $n$, and $\mathbf{f}(\mathbf{x}_n) = \nabla \phi$ is the gradient on the gravitational potential $\phi$. $\tilde{\gamma}_d = f_{\text{grad}} \sqrt{\tilde{\gamma}_d \Delta t}$, where $f_{\text{grad}}$ is a dimensionless pre-factor to scale the magnitude of the acceleration, $\Delta t$ is the time step of the simulation, and

$$ \gamma_d = \frac{|(\mathbf{x}_{n+1} - \mathbf{x}_n)^{\text{T}}[\mathbf{f}(\mathbf{x}_{n+1}) - \mathbf{f}(\mathbf{x}_n)]|}{||\mathbf{f}(\mathbf{x}_{n+1}) - \mathbf{f}(\mathbf{x}_n)||^2} \quad (27) $$

following Barzilai & Borwein (1988). This adds a small force along the steepest local gradient descent, helping the BH to remain attached to the centre of the cluster. The contribution of the gradient descent acceleration can be controlled with the free parameter $f_{\text{grad}}$.

#### 3.6.2. Accretion and spin

The BH accretes according to the Bondi-Hoyle-Lyttleton (BHL) accretion rate, capped at a maximum of 1% the Eddington accretion rate,

$$ \dot{M}_{\text{BH}} = \min(\dot{M}_{\text{BHL}}, 0.01 \dot{M}_{\text{Edd}}), \quad (28) $$

where

$$ \dot{M}_{\text{BHL}} = \frac{4\pi(GM_{\text{BH}})^2\bar{\rho}}{(\bar{c}_s^2 + \bar{v}_{\text{BH,rel}}^2)^{3/2}} \quad (29) $$

is the Bondi-Hoyle-Lyttleton accretion rate, and

$$ \dot{M}_{\text{Edd}} = \frac{4\pi G M_{\text{BH}} m_p}{\epsilon_{\text{mad}} c \sigma_{\text{T}}} \quad (30) $$

is the Eddington accretion rate. $G$ is the gravitational constant, $\bar{\rho}$, $\bar{c}_s$, and $\bar{v}_{\text{BH,rel}}$ are the weighted local average density, sound speed, and relative velocity between BH and gas, respectively (see Dubois et al. 2012, for details). $m_p$ is the proton mass, $c$ is the speed of light, $\sigma_{\text{T}}$ is the Thompson cross section, and $\epsilon_{\text{MAD}}(a_{\text{BH}})$ is the spin-dependent feedback efficiency from Dubois et al. (2021) based on MAD simulations by McKinney et al. (2012).

All quantities were measured within a radius $4\Delta x$ of the BH location. Throughout the evolution of the BH, its spin vector $a_{\text{BH}}$ was followed self-consistently through accretion (see Dubois et al. 2014, 2021, for technical details). The BH was initialised as non-spinning, and was spun up or down according to the angular momentum of the accreted gas throughout the simulation. The spin magnitude and orientation were updated assuming a magnetically arrested disc (MAD; McKinney et al. 2012) at all times.
<table><caption>**Table 1.** Simulation properties for simulations with an AGN.</caption><thead><tr><th>Simulation name</th><th>Initial B field</th><th>Jet energy fraction</th><th>f<sub>grad</sub></th><th>Thermal conduction</th><th>CR diffusion</th><th>CR streaming advection</th><th>CR streaming heating</th></tr></thead><tbody><tr><td>HYDRO</td><td>No</td><td>100% kinetic</td><td>0.2</td><td>No</td><td>No</td><td>No</td><td>No</td></tr><tr><td>MHDc</td><td>Yes</td><td>100% kinetic</td><td>0.2</td><td>Yes</td><td>No</td><td>No</td><td>No</td></tr><tr><td>CRc_dsh</td><td>Yes</td><td>10% kinetic, 90% CR</td><td>0.2</td><td>Yes</td><td>Yes</td><td>Yes</td><td>Yes</td></tr><tr><td>CRc_dsh_weakcr</td><td>Yes</td><td>90% kinetic, 10% CR</td><td>0.2</td><td>Yes</td><td>Yes</td><td>Yes</td><td>Yes</td></tr><tr><td>CRc</td><td>Yes</td><td>10% kinetic, 90% CR</td><td>0.2</td><td>Yes</td><td>No</td><td>No</td><td>No</td></tr><tr><td>CRc_d</td><td>Yes</td><td>10% kinetic, 90% CR</td><td>0.2</td><td>Yes</td><td>Yes</td><td>No</td><td>No</td></tr><tr><td>CRc_s</td><td>Yes</td><td>10% kinetic, 90% CR</td><td>0.2</td><td>Yes</td><td>No</td><td>Yes</td><td>No</td></tr><tr><td>CRc_h</td><td>Yes</td><td>10% kinetic, 90% CR</td><td>0.2</td><td>Yes</td><td>No</td><td>No</td><td>Yes</td></tr><tr><td>CRc_ds</td><td>Yes</td><td>10% kinetic, 90% CR</td><td>0.2</td><td>Yes</td><td>Yes</td><td>Yes</td><td>No</td></tr><tr><td>CRc_dh</td><td>Yes</td><td>10% kinetic, 90% CR</td><td>0.2</td><td>Yes</td><td>Yes</td><td>No</td><td>Yes</td></tr><tr><td>CRc_sh</td><td>Yes</td><td>10% kinetic, 90% CR</td><td>0.2</td><td>Yes</td><td>No</td><td>Yes</td><td>Yes</td></tr><tr><td>CRc_dsh_f0.5</td><td>Yes</td><td>10% kinetic, 90% CR</td><td>0.5</td><td>Yes</td><td>Yes</td><td>Yes</td><td>Yes</td></tr><tr><td>CRc_dsh_f1</td><td>Yes</td><td>10% kinetic, 90% CR</td><td>1</td><td>Yes</td><td>Yes</td><td>Yes</td><td>Yes</td></tr><tr><td>CRc_dsh_f1_2*</td><td>Yes</td><td>10% kinetic, 90% CR</td><td>1</td><td>Yes</td><td>Yes</td><td>Yes</td><td>Yes</td></tr><tr><td>CRc_dsh_1</td><td>Yes</td><td>10% kinetic, 90% CR</td><td>0.2</td><td>Yes</td><td>Yes</td><td>Yes</td><td>Yes</td></tr><tr><td>CRc_dsh_2</td><td>Yes</td><td>10% kinetic, 90% CR</td><td>0.2</td><td>Yes</td><td>Yes</td><td>Yes</td><td>Yes</td></tr></tbody></table>

Notes. All simulations have a total cluster mass of $8 \times 10^{14} M_{\odot}$ and an BH with a mass of originally $1.6 \times 10^{10} M_{\odot}$. All MHD simulations treat electrons and ions separately, while the non-MHD simulation (HYDRO) has only one temperature (and no CR). $f_{grad}$ stands for the strength of the gradient descent for the BH (see Sect. 3.6.1). The asterisk in $CRc\_dsh\_f1\_2*$ indicates that the sub-grid drag force from gas and stars is turned off.

### 3.6.3. Feedback

We modelled AGN feedback with jets following the method presented in Dubois et al. (2012). At each time step of the simulation, the AGN has a luminosity of

$$L_{\text{AGN}} = \epsilon_{\text{MAD}}(a_{\text{BH}}) \dot{M}_{\text{BHL}} c^2. \qquad (31)$$

This luminosity was used to accelerate gas within a jet cylinder of width of $2\Delta x_{min}$ and a height of $4\Delta x_{min}$, whose major axis was (anti)aligned with the BH spin vector. Jets were injected with a velocity of $10^5$ km s$^{-1}$. As this spin vector naturally evolved throughout the simulation under the influence of the precipitation of cold gas onto the cluster centre and its chaotic accretion onto the BH (Gaspari et al. 2013), there was no need to add an ad hoc precession to produce several generations of inflated AGN bubbles (Beckmann et al. 2019). At an average density ratio of $\rho_{jet}/\rho_{ICM} = 0.03$, our simulation modelled light jets.

The energy deposited in the jet was split into two channels: kinetic energy $L_{kin} = (1 - f_{CR})L_{AGN}$, and CR energy $L_{CR} = f_{CR}L_{AGN}$ according to a fixed CR energy fraction $f_{CR}$.

### 3.7. Simulation parameters

We present a suite of simulations to explore the impact of individual physical mechanism such as magnetic fields and CRs and their associated choice of parameters in a comparative manner. A summary of the parameters we used for all simulations can be found in Table 1.

## 4. Impact of magnetic fields and cosmic rays on the intracluster medium

Adding magnetic fields and CRs has a significant and persistent impact on the amount and distribution of warm and cold gas in the cluster centre. This is shown in Fig. 4.

In the absence of magnetic fields and CRs, (run HYDRO), a dense, centralised clump forms early on, with little dense

gas found beyond 5 kpc from the cluster centre at any point during the simulation. The timeseries in Fig. 5 shows that the AGN is able to regulate cooling in the cluster for the 3 Gyr of evolution probed here, with the central gas feature containing $10^9$ to $10^{10} M_{\odot}$ at all times. The AGN luminosity varies in the range $10^{43}-10^{46} \text{ erg s}^{-1}$ on duty cycles of about 100 Myr (Fig. 6).

With the addition of magnetic fields and anisotropic thermal conduction, MHDc, the central dense core remains, but locally thermally unstable regions extend farther from the cluster centre, as shown by the filamentary dense gas structure that surrounds the central core in Fig. 4. This increased cooling flow leads to large cold gas reservoirs that can exceed $10^{11} M_{\odot}$ (see Fig. 5), although the AGN injects more than an order of magnitude more cumulative energy throughout the 3 Gyr of evolution than in the HYDRO case (see Fig. 6), at $6.1 \times 10^{62} \text{ erg}$ for MHDc in comparison to $4.4 \times 10^{61} \text{ erg}$ for HYDRO. The largest cold gas reservoir, which builds up from 2.8 Gyr, continues to grow by the end of the simulation, which is a clear sign that the AGN has become overwhelmed and fails to self-regulate gas cooling in the cluster. As was shown in previous simulations, this run-away cooling is rarely brought back under control by the AGN (Li et al. 2015). A comparison of the MHDc and HYDRO shows the effect discussed in Ji et al. (2017) that magnetic fields enhance thermal instability in clusters (see Sect. 2 for a discussion).

The evolution of MHDc also shows that despite the presence of a powerful AGN jet that continuously stirs the cluster centre, anisotropic thermal conduction is unable to prevent cooling flows in massive galaxy clusters. This is due to a strong heat-buoyancy instability (Quataert 2008; Parrish et al. 2009) that preferentially aligns magnetic field lines in a tangential configuration and shuts of the radial heat flows required to offset cooling in the cluster centre. The impact of thermal conduction on cluster cooling flows is discussed in detail in the companion paper (Beckmann et al. 2022).

When CRs are added to an AGN jet, the choice has to be made how much energy is placed in CRs rather than kinetic energy. We tested two values of $f_{CR}$: a fiducial strong CR
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Projected gas density showing the multi-phase structure of the gas in the cluster centre for different physics (rows) and at three different times 1, 2, and 3 Gyr from left to right. The location of the BH is marked in black. The morphologies and the extent of the condensed gas depend strongly on the presence of magnetic fields and CRs.</figcaption></figure>

simulation with $f_{CR} = 0.9$, called CRc_dsh, and a complimentary weak CR simulation with $f_{CR} = 0.1$, called CRc_dsh_weak.

By injecting the majority of the jet energy as CRs, the AGN needs one-third less energy to regulate the cluster cooling flow over 3 Gyr (see Fig. 6). As a result, after also taking the CR injections from SNe into account, which contribute 8.4% and 1.7% of the overall CR energy for CRc_dsh_weak and CRc_dsh, respectively, the cluster under the influence of a CR-dominated jet receives a total CR energy that is roughly five times that of

CRc_dsh_weak over 3 Gyr (at $2.6 \times 10^{61}$ erg for CRc_dsh, and $5.2 \times 10^{60}$ erg for CRc_dsh_weak).

Simulation CRc_dsh_weak in Fig. 4 shows that a small number of CRs does not significantly change the morphology of the dense gas. However, with CRs, the total cold gas mass remains below $10^{10} M_{\odot}$ (Fig. 5) and the AGN continues to self-regulate the cluster for the whole 3 Gyr. While it is still possible that the AGN becomes overwhelmed at a later stage, this is unlikely as the very consistent amount of cold gas in the simulation over
<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Time evolution of the total gas mass of all gas that has $T < 10^6$ K ($M_{\text{warm}} + M_{\text{cold}}$, top panel), the warm gas mass ($M_{\text{warm}}$, $10^4$ K < $T < 10^6$ K, middle panel), and the cold gas mass ($M_{\text{cold}}$, $T < 10^4$ K, bottom panel) in the galaxy cluster under the influence of different physical mechanisms as indicated in the top inset. With weak CRs, a cooling catastrophe is avoided, and with strong CRs, the bulk of the condensed gas is maintained in the warm phase for long periods of time.</figcaption></figure>

more than 2 Gyr and the low variability in the AGN luminosity over the same time span are strong signs for a stable configuration. Overall, $\text{CRc\_dsh\_weak}$ closely resembles HYDRO in its evolution. The total AGN energy injected in both simulations differs by only 5%. This suggests that for a given cluster, a given amount of energy is required for self-regulation, which is robust to small changes in physics (this is discussed further in Appendix B).

The images in Fig. 4 show that the morphology of the condensed gas is much more volume filling and diffuse in the presence of a large CR fraction. The bottom panel of Fig. 5 shows that the bulk of this extended nebula is in the warm phase, that is, it has a temperature in the range $10^4-10^6$ K. In comparison to the other simulation, which contains $10^8-10^9 M_{\odot}$ of gas in this phase, the strong CR simulation $\text{CRc\_dsh}$ consistently produces more than $10^9 M_{\odot}$ of warm gas (middle panel), which is compensated for by a decrease in the cold gas mass (bottom panel). The total gas mass, however, is very similar for $\text{CRc\_dsh}$ and $\text{CRc\_dsh\_weak}$ on average, showing only a variation of 30%

<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Time evolution of the AGN luminosity (top panel) and of the cumulative energy injected by the AGN (bottom panel) for the different simulations with colours as indicated in the top inset. Even a small fraction of CRs in the jet allows the AGN to maintain self-regulation of the cluster over long periods of time.</figcaption></figure>

between the two simulations. Therefore, a comparable amount of gas condenses from the hot phase in both simulations, but a stronger fraction of CRs keeps this gas warm for long periods of time, rather than allowing it to cool to $10^4$ K or lower.

The impact of CR appears clearly in the phase plots shown in Fig. 7. In the HYDRO case, gas cools in a narrow isobaric fashion, except for a brief isochoric phase around $0.5 \text{ cm}^{-3}$. When it cools isobarically, gas evolves along a line with slope $T \propto n^{-1}$ in the plots shown here (solid line), while isochoric cooling causes it to evolve along a vertical line (dashed line). With the addition of magnetic fields, MHDc, the phases in which gas can be found do not significantly change in comparison to the non-magnetised case. With the addition of some CR, $\text{CRc\_dsh\_weak}$, the hot gas distribution narrows again, but the cold gas phase is broader. In general, there is no significant difference between properties stacked across different points in time (blue) versus individual snapshots (reds).

Only when a much larger fraction of the AGN luminosity is injected into CRs ($f_{\text{CR}} = 0.9$) does the distribution of the gas phases change significantly. As discussed in Kempski & Quataert (2020), with the addition of a sufficient number of CRs, the cooling in the cluster becomes predominantly isochoric, with little isobaric evolution. This leads to the broad distribution of gas shown in $\text{CRc\_dsh}$. This isochoric cooling in the presence of CRs was also reported by Butsky et al. (2020) using simulations of stratified cooling boxes. This work also showed that the presence of CR turns cooling isochoric, but efficient CR transport again flattens the slope of $T$ versus $n$.

One consequence of this evolution is that with strong CR feedback, gas can be found in a warm and diffuse phase that is absent from any of the other simulation, that is, where $T < 10^7$ K and $n < 0.3 \text{ cm}^{-3}$ (grey box). As the cooling rate $L_{\text{th}}$ depends on the gas density, this warm and diffuse gas will cool much
<figure><img src="image_8.png" /><figcaption>**Fig. 7.** 2D phase plots showing the volume-weighted distribution of gas number density *n* vs. temperature *T* within the central 60 kpc of the cluster for a stacked sample measured every 0.25 Myr between *t* = 0 Gyr and *t* = 3 Gyr. Solid (dashed) grey lines show example lines along which gas evolves when it cools in an isobaric (isochoric) fashion. Histograms along the *x*-axis (y-axis) show 1D histograms for the number density (temperature), weighted by *V*/*V*<sub>tot</sub>, where *V*<sub>tot</sub> is the total volume probed here. Solid lines show the distribution at individual points in time, and the dashed line shows the distribution for the stacked sample. With a CR-dominated jet, gas cooling in the cluster centre changes from isobaric (solid lines) to isochoric (dashed lines).</figcaption></figure>

more slowly than gas at the same temperature but higher density. This leads to the buildup of warm gas in the cluster centre as discussed previously, and, hence, prevents the runaway cooling of the hot phase at *T* ≈ 10⁷–10⁸ K into the cold star-forming phase *T* ≪ 10⁴ K.

We caution that the isochoric cooling discussed here might be partially or entirely due to the limited resolution of our simulations. As shown in Fielding et al. (2020), low resolution produces artificial pressure dips during rapid cooling, which vanish as the structure of cooling gas becomes better resolved. It is possible that the transition to isochoric cooling reported here and in Butsky et al. (2020) is influenced by this lack of resolution, and that cooling would remain isobaric if individual cloudlets with radii of <1 pc (McCourt et al. 2018), and their mixing layers, were fully resolved. This will remain technically impossible in cluster-scale simulations for the foreseeable future. We test

the impact of resolution on cooling flows and phase transitions in our simulations in Appendix D and find little change within the accessible parameter space.

We conclude that if even a fraction as small as 10% of the total jet energy is converted into CRs, the impact of magnetic fields on thermal instability is offset and AGN self-regulation of the cluster becomes more robust. Converting a larger fraction of AGN jet energy into CRs means that less overall jet energy is required to self-regulate the cluster cooling flow, but it also significantly changes the bulk cooling properties of the gas to isochoric. This leads to a broad distribution of warm diffuse gas in the cluster centre.

The images in Fig. 4 show that not only the morphology of condensed gas, but also its spatial distribution changes with increasing the CR fraction. As shown in Fig. 8, the radius that contains 80% of the warm and cold gas ($T < 2 \times 10^5$ K), $r_{80}$,
<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Time evolution of the radial extent of the warm gas from the cluster centre, measured as the radius $r_{80}$ that contains 80% of the total gas mass with temperatures $T < 2 \times 10^5$ K for the different simulations with different colours, as indicated in the top inset. With a CR-dominated jet, the central nebula is significantly more extended.</figcaption></figure>

varies with time but is broadly similar for simulations HYDRO, MHDc, and CRc_dsh_weak. However, for CRc_dsh, the bulk of the gas is found significantly farther from the cluster centre, out to radii as large as 20 kpc or more.

**5. Cosmic-ray transport mechanisms**

In this section, we investigate the role played by the three terms on the right-hand side of Eq. (14), that is, CR diffusion (“d” for the letter standing in the simulation name after “CRc_”), CR redistribution through streaming advection (“s”), and transfer of energy from CR to thermal energy through streaming heating (“h”). For this investigation, we compared simulations with only a single term (e.g. CRc_d with the diffusion term only) to simulations with two terms (e.g. CRc_dh for diffusion and streaming heating) to the full-physics simulation CRc_dsh and a simulation with CR, but without CR transport and related processes, CRc. We note that these reduced simulations have a limited physical validity as all transport mechanisms are physically tightly linked. They were performed here to gain insight into which term in Eq. (14) dominates the impact of CR on the cluster cooling flows. All simulations included for CRs advection by the gas, CR pressure work, and CR radiative losses.

Figure 9 shows that different combinations of transport terms lead to variations in the evolution of the total cold and warm gas mass over time, but the chaotic nature of the simulations makes it difficult to isolate clear trends in comparison to the range of evolution covered by the full-physics simulations CRc_dsh (filled grey area, based on all simulations from Appendix B). The simulation with CRs but without any transport mechanisms (CRc) does not produce a strong cooling flow. Like for all simulations with only partial CR transport, $M_{\text{warm}}$ is enhanced during the early evolution of the cluster, that is, while $t < 1.5$ Gyr. However, the late-time evolution of CRc is consistent with the range of behaviour exhibited by CRc_dsh for both warm and cold gas, without a sign of the run-away build-up of cold gas mass $M_{\text{cold}}$ exhibited for example by MHDc (see Sect. 5).

This is in contrast with the results presented for a galaxy group by Wang et al. (2020), where the case without stream-

<figure><img src="image_10.png" /><figcaption>**Fig. 9.** Time series for the total warm ($10^4$ K $< T < 10^6$ K) and cold ($T < 10^4$ K) gas mass for different CR transport mechanisms in comparison to the range covered by all full-physics simulations based on CRc_dsh presented in Appendix B (grey area). The dotted grey line denotes the mean of the grey area at each point in time. Overall trends are difficult to differentiate when individual simulations are compared due to the chaotic nature of cooling flows.</figcaption></figure>

ing advection and streaming heating produced a strong runaway cooling flow that was not regulated by the AGN, which is only balanced when CR streaming was included. Being a galaxy group, rather than the massive galaxy cluster studied here, the hot gas in Wang et al. (2020) is cooler, which results in shorter cooling times, and the AGN jet is about an order of magnitude less powerful. We note that in the presence of CR streaming heating, the redistribution of CRs within the cluster due to CR streaming advection becomes inefficient, as can be seen by comparing CRc_h and CRc_sh in Fig. 9, which only diverge after 1.5 Gyr of evolution. Therefore, the impact on cooling flows in Wang et al. (2020) is entirely due to streaming heating and not to a redistribution of CR by streaming advection. This suggests that CRs might play a very different role in galaxy groups and clusters of different masses.

Despite the lack of cooling flows in CRc, CR transport mechanisms do influence the long-term evolution of the cluster, as can be seen in the binned time series in Fig. 10. When the average warm gas mass of all simulations that include the (streaming) heating term (top panel, blue line) are compared to all simulations that do not (top panel, orange line), it becomes clear that streaming heating dominates the late-time evolution of the warm gas. Beyond $t > 1.5$ Gyr, simulations with heating have the same amount of warm gas on average as those with all transport terms (grey), whereas simulation without streaming heating produce significantly more warm gas. Therefore, CR streaming heating is able to keep some of the gas hot that would otherwise cool into the warm phase. The CR radiative cooling times due to Coulomb and hadronic collisions ranges from 84 Myr ($n_e = 5 \times 10^{-1}$ cm$^{-3}$) to 4.2 Gyr ($n_e = 10^{-2}$ cm$^{-3}$) based on the cooling rate from Guo & Oh (2008) in the hot gas, and can potentially be as short
<figure><img src="image_11.png" /><figcaption>**Fig. 10.** Average amount of warm gas $M_{\text{warm}}$ over time for all simulations with (blue) in comparison to all simulations without (orange) a given CR transport mechanism. Average masses for full-physics simulations are shown in grey. Shaded areas denote the minimum to maximum $M_{\text{warm}}$ at a given point in time for all simulations in a given bin. Only CR heating by the streaming instability has a significant impact on the evolution of the warm gas in the cluster, and only after $t > 1.5$ Gyr.</figcaption></figure>

as 20 Myr in the warm gas ($n_e = 2 \text{ cm}^{-3}$). For this reason, CR pressure is sustained over long periods of time, so that streaming heating provides a steady long-term source of energy that is able to maintain gas in the warm and hot phase for long periods of time. One notable feature of Fig. 10 is that the average $M_{\text{warm}}$ for full-physics simulations (grey line) continues to increase even after 3 Gyr of evolution. This suggests that the cluster has not yet reached a steady-state configuration. This raises an interesting question as to the eventual fate of this condensed gas. Baring a major reheating event, most of this gas will eventually have to cool to the cold phase, potentially producing a large amount of cold gas mass at later times. However, longer simulations runs by Ruszkowski et al. (2017) showed a continued gradual evolution of clusters under the influence of CRs on timescales of 6 Gyr or more, making a sudden cooling flow in CRc_dsh unlikely.

Streaming advection itself has a weaker influence on cluster evolution. Simulations with the streaming advection term have half the amount of $M_{\text{warm}}$ on average than those without, but still twice that of CRc_dsh. Results that CR streaming and streaming

heating play an important role in shaping long-term galaxy cluster cooling flows agree with those reported by Ruszkowski et al. (2017) and Wang et al. (2020).

The last panel of Fig. 10 shows that there is no significant difference in cooling flows with and without CR diffusion. The lack of importance of CR diffusion seems at odds with the conclusions from Sect. 6, which concluded that most of the gas is in a diffusion-dominated phase (i.e. $\chi > 1$). The important difference are the length-scales considered in the two scenarios. Section 6 is concerned with small-scale diffusion and its influence on locally thermally unstable regions of gas. By contrast, changing the large-scale cooling flow properties of the whole cluster would require diffusion on very large scales to redistribute CR in the cluster centre and flatten CR density profiles. Figure 11 shows that jets with an age of 200 Myr extend to around 50–70 kpc. In comparison, the diffusion length over the same period of time is $L_D = \sqrt{D_{\text{CR}}t} = 8.14$ kpc, which is significantly smaller than the length of the jet. As a result, CR remain confined to a region in and around the jet cone (see Fig. 4), and the overall impact of diffusion on the long-term evolution of the cluster is greatly diminished.

The relative impact of the different CR transport terms is confirmed in Fig. 12, which shows the distribution of gas densities (top panel) and temperatures (bottom panel) in the cluster centre. None of the transport terms has a significant impact on the distribution of gas densities, which is very similar for all stacked samples. However, with CR transport (CRc_dsh), the distribution of temperatures is more peaked in the hot phase. The only simulation that can reproduce this trend is CRc_h, that is, this peak in the temperature distribution is produced by streaming heating. A typical length scale on which heating deposits its energy can be found by comparing the streaming heating timescale $t_{\text{st}} = (\gamma_{\text{CR}} - 1)L/u_{\text{A}}$ with the diffusion timescale $t_D = L^2/D_{\text{CR}}$, which leads to $L = D_{\text{CR}}/u_{\text{A}} = 974(u_{\text{A}}/\text{km s}^{-1})^{-1}$ kpc $\approx 35$ kpc for our typical Alfvén velocity of $u_{\text{A}} = 30 \text{ km s}^{-1}$. This value is approximately comparable with the extent of diffuse warm structures seen in Fig. 4. Increasing $D_{\text{CR}}$ would increase both the efficiency of CR diffusion and the scale over which energy is deposited via streaming heating.

Overall, we conclude that in our simulations, CR transport processes are not required to suppress cooling flows. Even simulations without these processes show the broad distribution of density and temperature and long-term self-regulation of cooling flows associated with a CR-dominated jet in our simulations. However, CR streaming heating, and to some extent, the redistribution of CR within the cluster centre through streaming advection, do influence the late-time evolution of the cluster cooling flow by helping to reduce both the cold and warm gas mass on timescales of 1.5 Gyr and longer.

## 6. Thermal instability and cosmic rays

In this section we study the relative importance of thermal, magnetic, and CR energy in the gas, and assess when and why criteria for thermal (in)stability are met in galaxy clusters for different fractions of CR energy injected in the jet, $f_{\text{CR}}$. This section applies the thermal stability analysis of a thermal-CR composite fluid by Kempski & Quataert (2020), which has been conducted only in the linear regime. Recent work by Butsky et al. (2020) on the late non-linear evolution of thermal instability in the presence of CRs broadly confirmed results from Kempski & Quataert (2020), but the full applicability of the analysis by Kempski & Quataert (2020) in the saturated
<figure><img src="image_12.png" /><figcaption>**Fig. 11.** Temperature-weighted projections of $\eta = P_{cr}/P_{therm}$ for CRc\_dsh\_weak (**left panel**) and CRc\_dsh (**right panel**). Contours show the outline of remnants created by AGN jet ejections that occurred 200 Myr (white), 100 Myr (light grey), and 50 Myr (dark grey) before the snapshot. The black cross shows the current position of the AGN. The CR-dominated jet leads to an isotropic distribution of CRs in the cluster centre through the weak kinetic energy of the jet and the dense central nebula (**right**). For smaller fractions of CRs, the CRs remain confined to the jet cone (**left**).</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Fig. 12.** Volume-weighted probability distribution $\Phi(V/V_{\text{tot}})$ of the number density and temperature within the central 60 kpc of the cluster. Each line shows a stacked sample, measured at 0.25 Myr intervals throughout the simulation. While no CR transport mechanism influences the distribution of densities in the cluster centre, CR heating is required to reproduce the temperature distribution of the full-physics simulation.</figcaption></figure>

non-linear regime found in the ICM in galaxy clusters is yet to be determined.

Figure 13 shows that the main factor determining the relative importance of thermal pressure, CR pressure, and magnetic pressure is the underlying thermal phase of the gas: In general, thermal pressure dominates both magnetic pressure ($\beta > 1$) and CR pressure ($\eta < 1$) in the hot phase, but in the cold and warm phase,

CR pressure dominates thermal pressure ($\eta > 1$), and magnetic pressure contributes significantly or dominates (i.e. $\beta < 100$) in any fraction of CRs injected in the jet. As expected, a stronger injection of CR leads to a higher average CR fraction for cold and warm gas, and the wider spread in density and gas phases for the strong CR case also leads to a broader distribution of $\eta$ and $\beta$. However, when it is condensed out of the hot phase, all dense gas found here cools isochorically (i.e. has $\eta > 1$), even for CRc\_dsh\_weak.

Figure 11 shows that the spatial distribution of CRs in the hot gas is very different for the two cases. Because only 10% of the energy is injected into CRs, the kinetic energy of the jet allows it to remain collimated and extend to much larger radii. As a result, CRs remain confined within a broad cylinder aligned with the jet axis, and $\eta$ is very low in regions far from the jet.

By contrast, when 90% of the jet energy is injected as CR energy, in CRc\_dsh, the volume-filling fraction of CRs is much higher. This is partially due to the lower kinetic energy of the jet. When 90% of $L_{\text{AGN}}$ is injected into CR energy, only 10% are injected as kinetic energy, and CRc\_dsh generally already has a lower total $L_{\text{AGN}}$ than CRc\_dsh\_weak (see Fig. 6). With such high $f_{cr}$, jets carry much less momentum, which limits their ability to propagate to large distances from the cluster centre. Their propagation is further limited by the dense extended gas structure that fills the cluster centre out to 25 kpc or more (see Fig. 4). As a result, gas recently affected by AGN jets (white contours in Fig. 11) mixes more efficiently into the ICM, which leads to a more isotropic distribution of CRs in the cluster centre.

One important impact of stronger CR injection is that a hot phase dominated by CR pressure develops, which is entirely absent in CRc\_dsh\_weak (see Fig. 13), and which will cool isochorically even before transitioning to the warm phase. The distribution of cooling time $t_{\text{cool}}$ vs. $\eta$ for the hot phase alone (Fig. 14) shows that this CR-dominated hot gas has a wide range of short cooling times, but that the gas that is expected to cool and condense next (i.e. the gas with the shortest cooling time) also has the highest concentration of CRs. Therefore, newly cooled gas is well-saturated with CRs, which might explain the heating of filamentary nebulae in cluster centers (Ruszkowski et al. 2018). As gas cools, the thermal pressure
<figure><img src="image_14.png" /></figure>

**Fig. 13.** Mass-weighted distribution of plasma $\beta$ and CR pressure ratio $\eta$ for CRc_dsh_weak (top) and CRc_dsh (bottom) for cold ($T < 10^4$ K, blue), warm ($10^4 < T < 10^6$ K, orange), and hot ($T > 10^6$ K, red) gas, respectively. The small panels in the top row show only a single temperature bin, while all three temperature bins are superimposed in the large bottom panel. Above the dotted line, $P_{CR} > P_{mag}$, while the reverse is true below. The data shown are for gas within the central 100 kpc of the cluster at time $t = 2$ Gyr. All cold and warm gas is CR-pressure dominated, while all (most) of the hot gas is dominated by thermal pressure in CRc_dsh_weak (CRc_dsh).

drops and $\eta = P_{CR}/P_{therm}$ increases further, leading to the high values of $\eta$ in cold gas shown in Fig. 13.

As discussed in Sect. 2, the presence of magnetic fields enhances thermal instability, but the presence of CRs can stabilise small-scale thermal instability. Whether magnetised CR-saturated gas is thermally stable depends on the slope of the cooling function with respect to the temperature $T$, the power

<figure><img src="image_15.png" /></figure>

**Fig. 14.** Volume-weighted kernel density estimate for cooling time $t_{cool}$ vs. $\eta$ for hot gas in CRc_dsh_weak and CRc_dsh at time $t = 2$ Gyr. Contours enclose 100%, 99%, 90%, and 50% of the total volume. Gas with the shortest cooling time has the highest ratio of CR to thermal pressure, and is therefore well saturated with CRs.

of the cooling function $\alpha_T = \partial \log(L_{\text{th}})/\partial \log(T)$, and the fraction $\chi = D_{\text{CR}}/(t_{\text{cool}}\eta u_A^2)$. CRs can suppress the onset of thermal instability if $\chi < 1$, whereas thermally unstable gas remains so if $\chi > 1$ and $\alpha_T < 2$.

As predicted, the bulk of the hot gas in both CRc_dsh and CRc_dsh_weak can be found in the phase where CR make no difference to thermal instability ($\chi > 1, \alpha < 2$, lower right quadrant, Fig. 15). Only when the gas has already cooled to the warm phase, that is, when it has condensed out of the hot phase, does it reache $\alpha_T > 2$, where CR pressure can suppress instabilities and stabilise the gas against further collapse. CRc_dsh_weak has only a small amount of gas in the regime where $\chi < 1$, but in CRc_dsh, a fraction of gas across all three phases (2% of hot gas, 4% of warm gas, and 38% of cold gas, by mass) are found in this regime.

This leaves the question where in the cluster this potentially stabilised gas is located, and whether stabilisation against thermal instability occurs on a sufficiently large scale to influence the cluster cooling flow. The projections in Fig. 16 show that the lowest $\chi$ gas occurs within the AGN jet bubbles, in particular at their upper edges, where bubbles are rising buoyantly. This link between low $\chi$ and jet bubbles is shown more quantitatively in Fig. 17, where the lowest $\chi$ values are found for the gas that has been most recently directly affected by the AGN jet (i.e. has low $t_{\text{AGN}}$). Due to their long cooling times, driven by high temperatures and low densities, AGN jets and bubbles are locally thermally stable even in the absence of CRs (Gaspari et al. 2011; Li et al. 2015; Bourne et al. 2019; Beckmann et al. 2019). Therefore, we conclude that the location of (thermal) instability within massive galaxy clusters is not meaningfully altered in the presence of CRs. Nonetheless, the thermodynamic evolution of gas after it has cooled out of the hot phase changes significantly with the number of CRs present in the gas.

The radial distribution of $\eta$ and $\chi$ (Fig. 18, hot gas only) shows that the cluster core is dominated by CR pressure out to
<figure><img src="image_16.png" /><figcaption>Fig. 15. Mass-weighted distribution of $\chi$ vs. the slope of the cooling function $\alpha_T$ for cold ($T < 10^4$ K, blue), warm ($10^4 < T < 10^6$ K, orange), and hot ($T > 10^6$ K, red) gas. The small panels in the top row show only a single temperature bin, while all three temperature bins are superimposed in the large bottom panel. The top plot shows values for CRc_dsh_weak, and the bottom panel shows them for CRc_dsh. Data are measured within the central 100 kpc of the cluster at time $t = 2$ Gyr. As $\alpha_T < 2$ for all hot gas, CRs are unable to prevent the onset of thermal instability in galaxy clusters.</figcaption></figure>

~10 kpc for CRc_dsh_weak and out to ~30 kpc for CRc_dsh. To quantify the variation in the two quantities depending on the location of the jet axis, we compared the average radial distribution for the whole sphere (solid line) to two cylinders with a diameter of 15 kpc: “on-jet” is aligned with the z-axis of the simulation,

which roughly aligns with the jet-axis of the simulation at this point in time (see the black contours in Fig. 16 for both simulations), while “off-jet” is aligned with the y-axis of the simulation, and therefore oriented perpendicular to the jet axis. All measurements are volume weighted and were only performed for hot gas ($T > 10^6$ K).

Figure 18 shows that for both CRc_dsh_weak and CRc_dsh the radial profile converges to a CR-dominated core ($\eta > 1$) and gas-pressure-dominated cluster outskirts ($\eta < 1$) after 1 Gyr of evolution. The core is larger in CRc_dsh, where it extends to 40 kpc, in comparison to CRc_dsh_weak, where it extends only to 15 kpc. In both simulations, $\eta = 10^{-4}$ at large radii due to the initial conditions of the simulation. The profiles in $\chi$ build up over time; the lowest values are found in the range $5 < r < 20$ kpc in CRc_dsh_weak and $5 < r < 60$ kpc in CRc_dsh. In neither cluster is there a region where $\chi < 1$ on average. The fact that the radial profile of $\eta$ converges after 2–3 Gyr agrees well with Ruszkowski et al. (2017).

When we compare the radially averaged distribution (solid line) to the profiles along the jet (dashed) to perpendicular to the jet (dotted), the distribution for CRc_dsh is almost isotropic. The distributions of $\eta$ and $\chi$ on-jet and off-jet are very similar to the radial average case at all three times we tested here. This is not the case for the jet-dominated CRc_dsh_weak case, where for $10$ kpc $< r < 150$ kpc, the off-jet region has a significantly lower $\eta$, while $\chi$ of the jet is significantly higher than average.

Overall, we conclude that the fraction of AGN luminosity injected into CRs has a significant and lasting impact on the properties of the multi-phase ICM. We report a long-term self-regulation of cooling flows, and convergence in cluster properties over timescales of billions of years for a jet-dominated configuration with only a small fraction of CR ($f_{CR} = 0.1$) and for a CR-dominated feedback mode ($f_{CR} = 0.9$). The former closely resembles the MHD-only case in morphology and distribution of dense gas, but the small contribution of CR pressure to total gas pressure prevents over-cooling and allows the cluster to remain in self-regulation for long timescales. Strong CR feedback leads to a strongly CR-dominated cluster core and to an isotropic volume-filling distribution of warm gas that slows cluster cooling down over long timescales. Using the linear stability analysis of Kempski & Quataert (2020), we find no evidence for regions within the cluster in which CRs are able to prevent existing local thermal instability.

## 7. Comparison to $\gamma$-ray observations

While it has been conclusively shown in this work that CRs can influence the long-term evolution of galaxy clusters, the question remains whether this evolution is supported by observational evidence. One way to observationally constrain the CR content of galaxy clusters is via the $\gamma$-rays emitted by hadronic processes.

To compare our cluster simulations to observations, we computed the expected hadronic $\gamma$-ray photon flux $F_\gamma = \int q_\gamma(r, E_\gamma) dE_\gamma dV$ from our simulations using the source function $q_\gamma$ from Pfrommer & Enßlin (2004),

$$q_{\gamma}(\mathbf{r}, E_{\gamma}) \approx \sigma_{\text{pp}} c n_{\text{N}}(\mathbf{r}) \tilde{n}_{\text{CRp}}(\mathbf{r}) \times \frac{2^{4-\alpha_{\gamma}}}{3 \alpha_{\gamma}} \left(\frac{m_{\pi^0} c^2}{\text{GeV}}\right)^{-\alpha_{\gamma}} \left[ \left( \frac{2 E_{\gamma}}{m_{\pi^0} c^2} \right)^{\delta_{\gamma}} + \left( \frac{2 E_{\gamma}}{m_{\pi^0} c^2} \right)^{-\delta_{\gamma}} \right]^{-\alpha_{\gamma}/\delta_{\gamma}}, \quad (32)$$

where $\sigma_{\text{pp}} = 32(0.96 + e^{4.4-2.4\alpha_{\gamma}})$ mbarn is the effective cross section, c is the speed of light, $n_N = n_e(r)/(1 - 0.5X_{\text{He}})$ is the target particle number density, $n_e(r)$ is the electron number density
<figure><img src="image_17.png" /><figcaption>Fig. 16. Projections of $\chi$ weighted by the jet scalar (see Sect. 3.6) for $CR_c_{dsh\_weak}$ (left panel) and $CR_c_{dsh}$ (right panel). Contours show the outline of remnants created by AGN jet ejections that occurred 300 Myr (white), 200 Myr (light grey), 100 Myr (dark grey), and 50 Myr (black) before the snapshot. The black cross shows the current position of the AGN. The only diffusion-dominated gas (where $\chi < 1$) is found within buoyantly rising AGN bubbles.</figcaption></figure>

taken from our simulation, $X_{He} = 0.24$ is the fraction of helium, $a_{\gamma} = 4(a_{p} - 0.5)/3$ is the $\gamma$-ray spectral index, $m_{\pi^{0}}$ is the pion mass, and $\delta_{\gamma} = 0.14 a_{\gamma}^{-1.6} + 0.44$ is a shape parameter.

$$ \tilde{n}_{CRp}(\mathbf{r}) = \varepsilon_{CR}(\mathbf{r}) \frac{2(a_p - 1)}{m_p c^2} \left(\frac{m_p c^2}{\text{GeV}}\right)^{\alpha_p - 1} \mathcal{B}^{-1}\left(\frac{a_p - 2}{2}, \frac{3 - a_p}{2}\right) \quad (33) $$

is an effective CR number density, where $a_p$ is the slope of the CR spectrum, $\varepsilon_{CR}(\mathbf{r})$ is the kinetic CR energy density, taken from our simulation, $m_p$ is the proton mass, and $\mathcal{B}(a, b)$ is the beta-function.

To compare to the observed $\gamma$-ray flux limits from Ackermann et al. (2014), we integrated over the energy range $E_{\gamma} = [0.5, 200]$ GeV. To calculate the flux, we assumed that our cluster is at the distance of the Perseus cluster, 73.8 Mpc. For the model cited here, $F_{\gamma}$ is at its maximum for $a_p = 2.2$, but $F_{\gamma} \to 0$ as $a_p \to 2$ and $a_p \to 3$. To establish upper limits, we adopted $a_p = 2.2$ as our fiducial value, but also explored the impact of varying $a_p$.

The resulting total $\gamma$-ray energy flux is shown in Fig. 19. We also show for comparison point-mass $\gamma$-ray fluxes for cool-core clusters with $6 \times 10^{14} M_{\odot} < M_{200} < 10^{15} M_{\odot}$ from Ackermann et al. (2014).

Figure 19 shows that $CR_c_{dsh}$ produces a significantly higher $F_{\gamma}$ than $CR_c_{dsh\_weak}$. Emission from $CR_c_{dsh\_weak}$ is beginning to saturate from 2.5 Gyr onwards, whereas emission from $CR_c_{dsh}$ continues to rise. When it is compared with the observed upper limits for $\gamma$-ray fluxes from Ackermann et al. (2014), $CR_c_{dsh}$ quickly exceeds observational constraints by up to an order of magnitude. Emission from $CR_c_{dsh\_weak}$, by comparison, builds up much more slowly and converges at a value that agrees well with the observational upper limits.

We therefore conclude that for a self-regulated AGN, CR-dominated jets with high $f_{CR}$ are not supported observationally for massive galaxy clusters. A small contribution of CR to the AGN luminosity of $f_{CR} \le 0.1$, however, produces $\gamma$-ray emission that is compatible with current observational limits. This is true for a range of different slopes of the CR spectrum

<figure><img src="image_18.png" /><figcaption>Fig. 17. Volume-weighted kernel density estimated of $\chi$ vs. $t_{AGN}$, the time since gas has been directly affected by AGN feedback for all hot gas within a sphere of radius 100 kpc, centred on the cluster centre. Contours enclose 100%, 99%, 90%, and 50% of the volume in the cluster centre. Gas most recently affected by the AGN jet has the lowest $\chi$ values.</figcaption></figure>

$2.03 < a_p < 2.73$ (shaded region, covering the range of $a_p$ such that $F_{\gamma} > 0.25 \times F_{\gamma,\max}$). For $CR_c_{dsh}$ to agree with the observed limits from Ackermann et al. (2014), the CR spectrum would have to be either very shallow or very steep indeed.

In addition to heating and pressure support, CRs impact the gas chemistry, in particular in the cold molecular gas, where secondary electrons produced by CR dissociate CO molecules. Although we did not probe the molecular phase in these
<figure><img src="image_19.png" /><figcaption>**Fig. 18.** Time evolution of $\eta$ (top) and $\chi$ (bottom) profiles for CRc_dsh_weak (left) and CRc_dsh (right). Solid profiles are a volume-weighted average over the whole sphere. At $t = 1.5$ and 3 Gyr, we also show the profile measured along the jet (dashed) and along a cylinder perpendicular to the jet (dotted) for comparison. For comparison, the $t = 3$ Gyr value of the radial averaged profile for the other simulation is shown in red in each plot. In the CR-dominated and kinetic-dominated jets, the cluster core is dominated by CR pressure, but the extent of the CR-dominated core is much larger for a CR-dominated jet.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>**Fig. 19.** Time series of the total $\gamma$-ray flux within 100 kpc from the cluster centre, $F_\gamma$, over time. Point-mass $\gamma$-ray fluxes for all cool-core clusters with $6 \times 10^{14} M_\odot < M_{200} < 10^{15} M_\odot$ from Ackermann et al. (2014) are shown as vertical lines for comparison. A self-regulating CR-dominated jet produces a $\gamma$-ray emission that exceeds observational limits, while small fractions of CR in the jet are compatible observational constraints.</figcaption></figure>

simulations, we anticipate that high $f_{CR}$ values would help maintain a high C+/CO abundance (Mashian et al. 2013). This could

explain the very strong C+ cooling lines observed in brightest cluster galaxies (e.g. Edge et al. 2010; Mittal et al. 2012; Polles 2020).

## 8. Discussion and conclusions

We studied how injection of a fraction of the luminosity of the central AGN of a galaxy cluster as CRs influences the cooling flow of a massive galaxy cluster. We focused on how changing the fraction of AGN jet energy that is injected into CR energy rather than in kinetic energy, $f_{CR}$, changes the evolution of local thermal instability on timescales of billions of years, and what impact the resulting CR pressure and heating have on the multi-phase structure of the ICM. Our results are listed below.

1. Generally, injecting a fraction of the AGN luminosity as CR energy at the jet base aids the AGN in regulating cooling flows on timescales of billions of years. Even values as low as $f_{CR} \le 0.1$, that is, a case in which only 10% of the AGN luminosity is injected as CRs, prevent thermal instability and strong cooling flows.

2. For a CR-dominated jet (here $f_{CR} = 0.9$, i.e. 90% of $L_{AGN}$ is injected as CR energy), the structure of multi-phase gas changes significantly. Rather than cooling very quickly from the hot, diffuse phase to a cold, dense phase, gas spends long periods of time in a warm, diffuse phase and remains at lower densities even at cold temperatures. This is due to the additional pressure support from CRs and to the energy transfer from CRs to thermal energy via streaming instability heating.

3. Within massive galaxy clusters, CRs are distributed mainly via advection and not via diffusion or streaming. For low $f_{CR}$,
CRs are concentrated in and around the AGN jet. For high $f_{CR}$, the low kinetic energy in combination with the build-up of condensed gas in the cluster centre means that jets struggle to propagate to large distances and deposit their energy and CRs almost isotropically within the central 70 kpc of the cluster.

4. With any $f_{CR}$ tested here, hot gas in the cluster remains dominated by thermal pressure over CR pressure and magnetic pressure. Conversely, cold and warm gas is dominated by CR pressure for any $f_{CR}$.

5. Cosmic rays are unable to suppress local thermal instability in the hot gas because the conditions for them to do so, according to the linear perturbation analysis derived in Kempski & Quataert (2020), are only met within the buoyantly rising outer edges of the AGN bubbles, which are thermally stable already due to their long cooling times.

6. Self-regulating CR-dominated jets (i.e. those that have high $f_{CR}$) lead to γ-ray emission beyond current observational upper limits. Values as low as $f_{CR} = 0.1$, however, are still within observational upper limits, and therefore represent the more likely scenario for massive galaxy clusters.

A low fraction of CRs in the jet also agrees well with predictions from particle-in-cell simulations (e.g. Crumley et al. 2019), which show that only about 10% of the kinetic energy can be converted into CRs in strong shocks. Higher CR fractions can only be explained if at injection scales, $f_{CR} < 0.1$, but that by the resolution scale of the simulation (here $\Delta x \sim 500$ pc) sufficient kinetic energy has been transferred to gas internal energy and radiated away such that the remaining energy balance in the jet has shifted significantly. For strongly cooling gas, this scenario is certainly possible, but it does imply that the jet was significantly more powerful at injection. For example, for a jet to change from 20% CRs to 90%, it must have been 4.5 times as powerful at injection and must have lost 77% of its original energy by 500 pc.

Overall, we conclude that if up to 10% of the AGN luminosity is injected as CR energy, cooling flows self-regulate over timescales of billions of years due to the additional pressure support and to the additional heating of the gas via CR streaming for gas that condenses out of the hot phase. In this scenario, multi-phase gas properties look globally similar to the extensively studied non-magnetised CR-free case, without evidence for extensive diffuse, warm gas seen for gas with higher CR pressures.

Acknowledgements. R.B. and Y.D. designed the project. R.B., Y.D. and A.P. contributed to code development. R.B. designed, executed and processed the suite of simulations. R.B. and Y.D. analysed and interpreted results, and wrote the paper. V.O., F.P., M.L., O.H. and P.G. provided observational context and discussion of results. We would like to thank the anonymous referee for their comments. This work was supported by the ANR grant LYRICS (ANR-16-CE31-0011) and was granted access to the HPC resources of CINES under the allocations A0080406955 and A0100406955 made by GENCI. R.B. would like to thank Newnham College, University of Cambridge, for support. A.P. acknowledges funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement No. 679145, project 'COSMO-SIMS'). This work has made use of the Infinity Cluster hosted by Institut d'Astrophysique de Paris. We thank Stéphane Rouberol for smoothly running this cluster for us. Visualisations in this paper were produced using the yt project (Turk et al. 2011).

## References

Abdulla, Z., Carlstrom, J. E., Mantz, A. B., et al. 2019, ApJ, 871, 195

Ackermann, M., Ajello, M., Allafort, A., et al. 2010, ApJ, 717, L71

Ackermann, M., Ajello, M., Albert, A., et al. 2014, ApJ, 787, 18

Ackermann, M., Ajello, M., Albert, A., et al. 2016, ApJ, 819, 149

Ahnen, M. L., Ansoldi, S., Antonelli, L. A., et al. 2016, A&A, 589, A33

Aleksić, J., Alvarez, E. A., Antonelli, L. A., et al. 2012, A&A, 541, A99

Anders, E., & Grevesse, N. 1989, Geochim. Cosmochim. Acta, 53, 197

Andreon, S., Wang, J., Trinchieri, G., Moretti, A., & Serra, A. L. 2017, A&A, 606, A24

Barzilai, J., & Borwein, J. M. 1988, IMA J. Numer. Anal., 8, 141

Beckmann, R. S., Slyz, A., & Devriendt, J. 2018, MNRAS, 478, 995

Beckmann, R. S., Dubois, Y., Guillard, P., et al. 2019, A&A, 631, A60

Beckmann, R. S., Dubois, Y., Pellissier, A., Polles, F. L., & Olivares, V. 2022, A&A, in press, https://doi.org/10.1051/0004-6361/202243873

Bell, A. R. 1978, MNRAS, 182, 147

Blandford, R. D. 2000, Phys. Scr., T85, 191

Blandford, R. D., & Ostriker, J. P. 1978, ApJ, 221, L29

Bogdán, Á., Lovisari, L., Volonteri, M., & Dubois, Y. 2018, ApJ, 852, 131

Bogdanović, T., Reynolds, C. S., Balbus, S. A., & Parrish, I. J. 2009, ApJ, 704, 211

Bourne, M. A., Sijacki, D., & Puchwein, E. 2019, MNRAS, 490, 343

Brunetti, G., Blasi, P., Reimer, O., et al. 2012, MNRAS, 426, 956

Brunetti, G., Zimmer, S., & Zandanel, F. 2017, MNRAS, 472, 1506

Buck, T., Pfrommer, C., Pakmor, R., Grand, R. J. J., & Springel, V. 2020, MNRAS, 497, 1712

Butsky, I. S., & Quinn, T. R. 2018, ApJ, 868, 108

Butsky, I. S., Fielding, D. B., Hayward, C. C., et al. 2020, ApJ, 903, 77

Chan, T. K., Kereš, D., Hopkins, P. F., et al. 2019, MNRAS, 488, 3716

Conselice, C. J., Gallagher, J. S., III, & Wyse, R. F. G. 2001, AJ, 122, 2281

Cowie, L. L., & McKee, C. F. 1977, ApJ, 211, 135

Croston, J. H., Hardcastle, M. J., Birkinshaw, M., Worrall, D. M., & Laing, R. A. 2008, MNRAS, 386, 1709

Croston, J. H., Ineson, J., & Hardcastle, M. J. 2018, MNRAS, 476, 1614

Crumley, P., Caprioli, D., Markoff, S., & Spitkovsky, A. 2019, MNRAS, 485, 5105

Dashyan, G., & Dubois, Y. 2020, A&A, 638, A123

Dubois, Y., & Commerçon, B. 2016, A&A, 585, A138

Dubois, Y., Devriendt, J., Slyz, A., & Teyssier, R. 2010, MNRAS, 409, 985

Dubois, Y., Devriendt, J., Slyz, A., & Teyssier, R. 2012, MNRAS, 420, 2662

Dubois, Y., Volonteri, M., & Silk, J. 2014, MNRAS, 440, 1590

Dubois, Y., Commerçon, B., Marcowith, A., & Brahimi, L. 2019, A&A, 631, A121

Dubois, Y., Beckmann, R., Bournaud, F., et al. 2021, A&A, 651, A109

Edge, A. C., Oonk, J. B. R., Mittal, R., et al. 2010, A&A, 518, L46

Ehlert, K., Weinberger, R., Pfrommer, C., Pakmor, R., & Springel, V. 2018, MNRAS, 481, 2878

Ehlert, K., Pfrommer, C., Weinberger, R., Pakmor, R., & Springel, V. 2019, ApJ, 872, L8

Ehlert, K., Weinberger, R., Pfrommer, C., Pakmor, R., & Springel, V. 2022, MNRAS, submitted [arXiv:2204.01765]

Enßlin, T., Pfrommer, C., Miniati, F., & Subramanian, K. 2011, A&A, 527, A99

Fabian, A. C. 1994, ARA&A, 32, 277

Fabian, A. C., Walker, S. A., Russell, H. R., et al. 2016, MNRAS, 461, 922

Fielding, D. B., Ostriker, E. C., Bryan, G. L., & Jermyn, A. S. 2020, ApJ, 894, L24

Fromang, S., Hennebelle, P., & Teyssier, R. 2006, A&A, 457, 371

Fujita, Y., & Ohira, Y. 2011, ApJ, 738, 182

Gaspari, M., Ruszkowski, M., & Sharma, P. 2011, ApJ, 746, 94

Gaspari, M., Ruszkowski, M., & Oh, S. P. 2013, MNRAS, 432, 3401

Gaspari, M., McDonald, M., Hamer, S. L., et al. 2018, ApJ, 854, 167

Girichidis, P., Naab, T., Hanasz, M., & Walch, S. 2018, MNRAS, 479, 3042

Govoni, F., & Feretti, L. 2004, Int. J. Mod. Phys. D, 13, 1549

Guo, F., & Oh, S. P. 2008, MNRAS, 384, 251

Hogan, M. T., McNamara, B. R., Pulido, F. A., et al. 2017, ApJ, 851, 66

Holguin, F., Ruszkowski, M., Lazarian, A., Farber, R., & Yang, H.-Y. K. 2019, MNRAS, 490, 1271

Hopkins, P. F., Squire, J., Chan, T. K., et al. 2021, MNRAS, 501, 4184

Huang, X., Jiang, Y.-F., & Davis, S. W. 2022, ApJ, 931, 140

Jacob, S., & Pfrommer, C. 2017a, MNRAS, 467, 1449

Jacob, S., & Pfrommer, C. 2017b, MNRAS, 467, 1478

Ji, S., Oh, S. P., & McCourt, M. 2017, MNRAS, 476, 852

Ji, S., Chan, T. K., Hummels, C. B., et al. 2020, MNRAS, 496, 4221

Jiang, Y.-F., & Oh, S. P. 2018, ApJ, 854, 5

Kannan, R., Vogelsberger, M., Pfrommer, C., et al. 2016, ApJ, 837, L18

Kempski, P., & Quataert, E. 2020, MNRAS, 493, 1801

Kimm, T., Cen, R., Devriendt, J., Dubois, Y., & Slyz, A. 2015, MNRAS, 451, 2900

Kunz, M. W., Schekochihin, A. A., Cowley, S. C., Binney, J. J., & Sanders, J. S. 2011, MNRAS, 410, 2446

Leccardi, A., & Molendi, S. 2008, A&A, 487, 461

Li, Y., Bryan, G. L., Ruszkowski, M., et al. 2015, ApJ, 811, 73

Li, Y., Gendron-Marsolais, M.-L., Zhuravleva, I., et al. 2020, ApJ, 889, L1
Maccio, A. V., Dutton, A. A., Van Den Bosch, F. C., et al. 2007, MNRAS, 378, 55

Mashian, N., Sternberg, A., & Loeb, A. 2013, MNRAS, 435, 2407

Matthews, J. H., Bell, A. R., & Blundell, K. M. 2020, New Astron. Rev., 89, 101543

McCourt, M., Sharma, P., Quataert, E., & Parrish, I. J. 2012, MNRAS, 419, 3319

McCourt, M., Oh, S. P., O'Leary, R., & Madigan, A.-M. 2018, MNRAS, 473, 5407

McDonald, M., Gaspari, M., McNamara, B. R., & Tremblay, G. R. 2018, ApJ, 858, 45

McKinney, J. C., Tchekhovskoy, A., & Blandford, R. D. 2012, MNRAS, 423, 3083

Mittal, R., Oonk, J. B. R., Ferland, G. J., et al. 2012, MNRAS, 426, 2957

Miyoshi, T., & Kusano, K. 2005, J. Comput. Phys., 208, 315

Mohapatra, R., Jetti, M., Sharma, P., & Federrath, C. 2022, MNRAS, 510, 2327

Morganti, R., Fanti, R., Gioia, I. M., et al. 1988, A&A, 189, 11

Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 490, 493

Olivares, V., Salome, P., Combes, F., et al. 2019, A&A, 631, A22

Olivares, V., Salome, P., Hamer, S. L., et al. 2022, A&A, in press, https://doi.org/10.1051/0004-6361/202142475

Ostriker, E. C. 1999, ApJ, 513, 252

Pal Choudhury, P., Sharma, P., & Quataert, E. 2019, MNRAS, 488, 3195

Parrish, I. J., Quataert, E., & Sharma, P. 2009, ApJ, 703, 96

Pfister, H., Volonteri, M., Dubois, Y., Dotti, M., & Colpi, M. 2019, MNRAS, 486, 101

Pfrommer, C., & Enßlin, T. A. 2004, A&A, 413, 17

Pfrommer, C., Springel, V., Ensslin, T. A., & Jubelgas, M. 2006, MNRAS, 367, 113

Phipps, F., Bogdán, Á., Lovisari, L., et al. 2019, ApJ, 875, 141

Polles, F. L. 2020, Proc. Int. Astron. Union, 15, 185

Pulido, F. A., McNamara, B. R., Edge, A. C., et al. 2018, ApJ, 853, 177

Quataert, E. 2008, ApJ, 673, 758

Rasera, Y., & Teyssier, R. 2006, A&A, 445, 1

Rosen, A., & Bregman, J. N. 1995, ApJ, 440, 634

Ruszkowski, M., Lee, D., Brüggen, M., Parrish, I., & Oh, S. P. 2011, ApJ, 740, 81

Ruszkowski, M., Yang, H.-Y. K., & Reynolds, C. S. 2017, ApJ, 844, 13

Ruszkowski, M., Yang, H.-Y. Y. K., & Reynolds, C. S. 2018, ApJ, 858, 64

Ryu, D., Kang, H., Hallman, E., & Jones, T. W. 2003, ApJ, 593, 599

Salomé, P., Combes, F., Revaz, Y., et al. 2011, A&A, 531, A85

Sarazin, C. L. 1986, Rev. Mod. Phys., 58, 1

Semenov, V. A., Kravtsov, A. V., & Caprioli, D. 2021, ApJ, 910, 126

Sharma, P., & Hammett, G. W. 2007, J. Comput. Phys., 227, 123

Sharma, P., McCourt, M., Quataert, E., & Parrish, I. J. 2012, MNRAS, 420, 3174

Sijacki, D., & Springel, V. 2006, MNRAS, 366, 397

Sijacki, D., Pfrommer, C., Springel, V., & Enlin, T. A. 2008, MNRAS, 387, 1403

Snodin, A. P., Shukurov, A., Sarson, G. R., Bushby, P. J., & Rodrigues, L. F. S. 2016, MNRAS, 457, 3975

Su, K.-Y., Hopkins, P. F., Hayward, C. C., et al. 2020, MNRAS, 491, 1190

Su, K.-Y., Hopkins, P. F., Bryan, G. L., et al. 2021, MNRAS, 507, 175

Sutherland, R. S., & Dopita, M. A. 1993, ApJS, 88, 253

Teyssier, R. 2002, A&A, 385, 337

Teyssier, R., Fromang, S., & Dormy, E. 2006, J. Comput. Phys., 218, 44

Thomas, T., & Pfrommer, C. 2019, MNRAS, 485, 2977

Toro, E. F. 2009, Riemann Solvers and Numerical Methods for Fluid Dynamics (Berlin, Heidelberg: Springer)

Trotta, R., Jóhannesson, G., Moskalenko, I. V., et al. 2011, ApJ, 729, 106

Turk, M. J., Smith, B. D., Oishi, J. S., et al. 2011, ApJS, 192, 9

Urban, O., Werner, N., Allen, S. W., Simionescu, A., & Mantz, A. 2017, MNRAS, 470, 4583

Vazza, F., Brunetti, G., & Gheller, C. 2009, MNRAS, 395, 1333

Voit, G. M. 2018, ApJ, 868, 102

Voit, G. M. 2021, ApJ, 908, L16

Voit, G. M., Meece, G., Li, Y., et al. 2017, ApJ, 845, 80

Wang, C., Ruszkowski, M., & Yang, H.-Y. K. 2020, MNRAS, 493, 4065

Wang, C., Ruszkowski, M., Pfrommer, C., Oh, S. P., & Yang, H.-Y. K. Y. K. 2021, MNRAS, 504, 898

Wiener, J., Oh, S. P., & Guo, F. 2013, MNRAS, 434, 2209

Worrall, D. M. 2009, A&ARv, 17, 1
# Appendix A: Cooling function

<figure><img src="image_21.png" /><figcaption>**Fig. A.1.** Cooling function $\mathcal{L}_{\text{cool}}$ (top panel) and $\alpha_T$ (bottom panel) used in the simulations presented here, as a function of temperature *T* and metallicity *Z*.</figcaption></figure>

Fig. A.1 shows that the cooling function used in the simulations presented here, based on Sutherland & Dopita (1993) and Rosen & Bregman (1995), is a function of temperature and metallicity, but is entirely independent of number density for $10^{-2} < n < 10^2 \text{ cm}^{-3}$. As a consequence, the same is true for its slope $\alpha_T$. In general, for massive galaxy clusters, $\alpha_T > 2$ only around temperatures of $10^4 \text{ K}$.

<figure><img src="image_22.png" /><figcaption>**Fig. B.1.** Time evolution of the total warm $10^4 < T < 10^6 \text{ K}$ and cold $T < 10^4 \text{ K}$ gas mass for variations of the simulations CRc_dsh. The time evolution of the cold and warm gas mass is robust to perturbations in AGN dynamics.</figcaption></figure>

# Appendix B: Stochastic nature and robustness of results

Cluster cooling cycles are inherently chaotic. The energy input via the AGN jet is driven not by the overall amount of cold gas, but by the consistent supply of cold gas very close to the BH. The zone of influence of the BH is small (at the resolution limit for our simulations). The total injected energy therefore not only depends on global cooling, but also, among others, on the distribution of cold gas within the cluster centre and on the position of the BH. In this section, we test how robust the trends observed in the previous sections are to chaotic perturbations of the system over long periods of time.

To do this, we conducted five variations of simulation CRc_dsh, the details of which can be found in Table 1:

- CRc_dsh_1 and CRc_dsh_2 have identical parameter choices to CRc_dsh, but initial conditions were set up with a different random seed for the initial tangled magnetic field and velocity perturbations.

- CRc_dsh_f0.5 and CRc_dsh_f1 have identical initial conditions to CRc_dsh, but have a higher pre-factor for the gradient descent acceleration. In these simulations, the BH is anchored more securely to the cluster centre.

- CRc_dsh_f1_1 is identical to CRc_dsh_f1 except that subgrid drag forces due to gas and particles have been switched off to test their impact on the dynamics of the BH.

<figure><img src="image_23.png" /><figcaption>**Fig. B.2.** Time series of cumulative AGN energy injected throughout the simulation and distance of the AGN from the cluster centre for variations of the simulations CRc_dsh. The overall energy input required for cluster self-regulation is very similar even for different AGN dynamics.</figcaption></figure>

Fig. B.1 confirms that the global evolution of gas phases is similar for all six simulations. While the turbulent variations in the warm (top panel) and cold (bottom panel) gas mass is readily apparent when the simulations are compared, no overall trends or strong divergence of results emerge over the 3 Gyr of evolution probed here. This conclusion is confirmed by the AGN time series in Fig. B.2. The AGN luminosity, picking up the variations in mass and distribution of the condensed gas, varies strongly but remains within the same overall luminosity range for all six simulations. The overall AGN energy injected into all six
simulations after 3 Gyr of simulation is even more robust. It ranges from $2.2 \times 10^{61}$ erg s$^{-1}$ for $CRc\_dsh\_1$ to $3.5 \times 10^{61}$ erg s$^{-1}$ for $CRc\_dsh\_f1$ in Fig. B.2. We therefore conclude that overall, the conclusions drawn in this paper are robust to reasonable perturbations of the system, despite its overall chaotic nature.

## Appendix C: Cosmic-ray parameters

<figure><img src="image_24.png" /><figcaption>**Fig. C.1.** Time evolution of the cumulative energy injected by the AGN (top) and the warm (middle) and cold (bottom) gas mass for the three variations of $CRc\_dsh$ listed in Table C.1. Grey shaded areas show the statistical variation of $CRc\_dsh$ from Appendix B. Neither a lower $D_{CR}$ nor a higher $f_{b,st}$ change the cooling flow properties of the clusters beyond stochastic variation.</figcaption></figure>

To ensure that our results do not depend sensitively on the parameters chosen as part of the CR model, in this section we conduct a series of variations of our fiducial simulation, $CRc\_dsh$, with different values for $D_{CR}$ and $f_{b,st}$ (see Sect. 3.1), to determine how sensitively our results depend on these two parameters. The CR diffusion coefficient $D_{CR}$ is poorly constrained for astrophysical plasma because it sensitively depends on the small-scale morphology of the magnetic field on length scales of the order of the Larmor radius. This is an estimated $10^8$ orders of magnitude below the resolution of our simulation (Snodin et al. 2016). We therefore relied on a simplified description of diffusion, whose efficiency is captured in the magnitude of $D_{CR}$. For our fiducial simulation, we selected $D_{CR} = 10^{29}$ cm$^{-2}$ s$^{-1}$, based on recent work by Chan et al. (2019) (see also Hopkins et al. 2021), to con-

<table><caption>**Table C.1.** Variation of $CRc\_dsh$ to explore the parameter choices of the CR model. All parameters not listed in this table are identical to $CRc\_dsh$. The simulation with the asterisk is included as a reminder.</caption><thead><tr><td>Simulation name</td><td>log($D_{CR}$) [cm<sup>-2</sup> s<sup>-1</sup>]</td><td>$f_{b,st}$</td></tr></thead><tbody><tr><td>$CRc\_dsh^*$</td><td>29</td><td>1</td></tr><tr><td>$CRc\_dsh\_Dcr28$</td><td>28</td><td>1</td></tr><tr><td>$CRc\_dsh\_Dcr27$</td><td>27</td><td>1</td></tr><tr><td>$CRc\_dsh\_fst1.7$</td><td>29</td><td>1.7</td></tr></tbody></table>

strain $D_{CR}$ through a comparison of simulations to observations. A second free parameter is $f_{b,st}$, which depends on the efficiency of damping processes that are excited by the streaming instability. Landau damping is thought to be the dominant mechanism for typical ICM properties, which allows for (moderate) super-Alfvénic streaming speeds such that $f_{b,st} = 1.7$ (Ruszkowski et al. 2017). To test the robustness of our results with respect to these two parameters, we examined three variations of our fiducial simulation: lower diffusion coefficients of $D_{CR} = 10^{27}$ cm$^{-2}$ s$^{-1}$ ($CRc\_dsh\_Dcr1E28$) and $D_{CR} = 10^{28}$ cm$^{-2}$ s$^{-1}$ ($CRc\_dsh\_Dcr1E27$), and an increased streaming boost factor of $f_{b,st} = 1.7$ ($CRc\_dsh\_fst1.7$). The parameter choices are summarised in Tab. C.1.

Fig. C.1 shows that lowering $D_{CR}$ by up to two orders of magnitude or increasing $f_{b,st}$ to $f_{b,st} = 1.7$ has little effect on the cooling flow in the cluster. The total energies injected by the AGN in each case are very similar, and the time series of $M_{warm}$ and $M_{cold}$ remain within the range for the stochastic variation of $CRc\_dsh$ shown in Appendix B.

## Appendix D: Resolution study

<figure><img src="image_25.png" /><figcaption>**Fig. D.1.** Time evolution of the cumulative energy injected by the AGN (top) and the warm and cold (solid bottom) and warm gas only (dotted bottom) gas mass over time. Grey shaded areas show the statistical variation of $CRc\_dsh$ from Appendix B. The case with the lowest resolution, $CRc\_dsh\_1kpc$, requires significantly more energy to balance cooling flows at both early and late times, and produces predominantly warm gas. All three runs with higher resolution show good convergence in the cumulative AGN energy and total cold and warm gas. A higher resolution leads to more cold gas vs warm gas early on.</figcaption></figure>
<figure><img src="image_26.png" /><figcaption>**Fig. D.2.** Projected gas density showing how the multi-phase structure of the gas changes with resolution.</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>**Fig. D.3.** 2D phase plots showing the volume-weighted distribution of gas number density *n* vs temperature *T* within the central 60 kpc of the cluster for a stacked sample measured every 0.25 Myr between *t* = 0 Gyr and *t* = 3 Gyr. Solid (dashed) grey lines show example lines along which gas evolves when it cools in an isobaric (isochoric) fashion. Histograms along the x-axis (y-axis) show 1D histograms for the number density (temperature), weighted by *V*/*V*<sub>tot</sub>, where *V*<sub>tot</sub> is the total volume probed here. Solid lines show the distribution at individual points in time, and the dashed line shows the distribution for the stacked sample. With increasing resolution, the distribution becomes broader, but the generally isochoric nature of cooling remains.</figcaption></figure>

To test whether our results are strongly dependent on resolution, we present three variations of our fiducial simulation CRc_dsh with different maximum resolutions Δx. As a reminder, CRc_dsh has Δx = 531 pc (refinement level 14). In this section, we also present CR_dsh_1kpc, with Δx = 1.06 kpc (level 13), CR_dsh_265pc (Δx = 265 pc, level 15), and CRC_dsh_132pc (Δx = 132 pc, level 16).

Fig. D.1 shows that the results converge reasonably well for a
simulation of this complexity, but only just. CRc_dsh_1kpc, the
run with a lower resolution than the fiducial simulation, shows
that a significant increase in cumulative AGN energy is needed to
balance the cooling flow at early (t < 1 Gyr) and late (t > 1 Gyr)
times. This simulation produces very little cold gas overall (see
the difference between the dashed and dotted line in the bottom
panel of Fig. D.1), and while the onset of warm gas production
is delayed in comparison to CRc_dsh, it does not seem to find a
stable configuration, with M_warm + M_cold continuing to increase
after 3 Gyr.

By comparison, the two simulations with a higher resolu-
tion, CRc_dsh_265pc and CRc_dsh_132pc, show a good conver-
gence. The cumulative AGN energy required for all three is very
similar, as is M_warm + M_cold within the bounds of the expected
stochastic variation of the cluster. The only notable trend is a
reduction in M_warm and a corresponding increase in M_cold with
increasing resolution early on.

This difference in the distribution of gas is also visible in
the density projection in Fig. D.2. A lower resolution leads to a
more diffuse central nebula that is less extended at early times.
Increasing the resolution allows for more clumping of the gas,
which in turn means a broader distribution in temperature – den-
sity space, as is shown in Fig. D.3. However, the general evo-
lution of cooling, that is, isobaric at high and low temperatures,
but a broad isochoric distribution in between, is captured for the
whole range of resolutions we tested here.

We therefore conclude that for Δx < 500 pc, large-scale cool-
ing flows in our galaxy clusters are well resolved. The internal
structure of the multi-phase gas continues to evolve, however.