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

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 \gammaray~ flux in excess of current observational limits, but low CR fractions in the jet are not ruled out.


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 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α emission nebulae, as well as cold molecular gas and their extended filamentary morphologies (e.g.Salomé et al. 2011;Olivares et al. 2019Olivares et al. , 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 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, η = P CR /P therm < 0.1 for both hadronic and reacceleration models (Ackermann et al. 2010(Ackermann et al. , 2014(Ackermann et al. , 2016;;Aleksić et al. 2012;Brunetti et al. 2012Brunetti et al. , 2017)).For the Perseus cluster, constraints are even tighter, with η < 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, η < 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 multiphase 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.

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 is shorter than the free-fall time 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 cool /t ff ≤ 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. 2019Olivares et al. , 2022)).
Another approach to determining local thermal instability is to follow Gaspari et al. (2018) and consider the local eddy turnover time, in comparison to the cooling time, where σ 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.
A129, page 2 of 23 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, δρ/ρ as a function of f ff .For the hydro-only case, they reported In the presence of magnetic fields, this becomes where β = 8πP therm /B 2 is the plasma beta and P therm is the thermal pressure.Comparison of the two shows that δρ ρ MHD > δρ ρ hydro for fields as weak as β ≤ 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 L th with respect to temperature T , which for the simulations presented here is shown in Appendix A; 2. the ratio of CR pressure P CR and thermal pressure P therm : 3. and the ratio where D CR is the CR diffusivity, and u A is the Alfvén velocity, χ, 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 χ.When χ > 1, thermal instabilities develop with an amplitude proportional to 2 − α T , that is, thermal instabilities grow strongly when α T < 2. In the opposite regime, where χ ≤ 1, CR diffusion suppresses thermal instability below a maximum length-scale, which can be expressed as an effective field length of the form where b is the magnetic field unit vector, and k is the wave number unit vector.The thermal instability continues to occur on larger scales.λ CRF is at a maximum when η = 1, and it falls off for higher and lower values.

CRF, max [kpc]
Fig. 1.Radial profile of χη (top panel), η crit , the minimum η required for the strong diffusion regime (middle panel), and the maximum λ 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 noncored profile (dashed lines).Grey horizontal lines demarcate the two regimes in χ (top panel) and gas that are CR-pressure dominated vs thermal pressure dominated (middle panel).η crit , the minimum value required for χ < 1 increases with radius.By contrast, CR pressure profiles are expected to decrease with radius.It is therefore unlikely that χ < 1 at large radii, which means that CRs are unable to stabilise smallscale instabilities outside the cluster core.

Expected behaviour for galaxy clusters
In this section, we predict η crit , the value of η required to have χ = 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 η required to have χ = 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 χ > 1, as hot gas in galaxy clusters cools predominantly through free-free emission, which has α T = 0.5 < 2. By contrast, if χ < 1 for a significant volume of the cluster, there could be a maximum length scale, the CR Field length λ CRF (see Eq. ( 9)), below which CRs are able to suppress thermal instability.Instability would be expected to continue on scales larger than λ CRF .As both χ and λ CRF explicitly depend on η, 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 χη for galaxy clusters of different masses in hydrostatic equilibrium.To test the A129, page 3 of 23 impact of different distributions of CRs on the cluster, we computed χη and then determined η crit , which will give χ < 1 if η > η crit (or χ > 1 if η < η crit ).Profiles to compute χη 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 Maccio et al. (2007), Andreon et al. (2017), andBogdán et al. (2018), respectively.The central magnetic field was set to 10 µG (e.g.Govoni & Feretti 2004), and was scaled with ρ(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 χη varies strongly with galaxy cluster mass, thermal profile (cored or non-cored), and radius.More massive clusters need higher CR pressure to reach χ < 1, which requires η > η crit (middle panel) as χ ∝ 1 η .The required value of η 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 η to decrease with radius.From such injections, it is likely that only the cluster centre, or no part of the cluster, has χ < 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 η in the cluster centre.Outside regions where χ < 1, even in the presence of CRs, thermally unstable gas will remain, so that even in the presence of CRs as α 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 λ 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 η = 1.λ CRF decreases for higher and lower values of η, and for dimensions that are not parallel with the magnetic field lines.For a Perseus-like cluster, 10 kpc of λ 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 η 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 < η < η 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 η > 1 cools isochorically rather than isobarically.As gas cools, P therm decreases, thus, η = 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.

Magneto-hydrodynamics with cosmic rays
This paper presents a set of hydrodynamical and magnetohydrodynamical (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), where ρ is the mass density; u is the velocity; e = 0.5ρu 2 + e e + e i + e CR + B 2 /8π is the total energy density; e e , e i , and e CR are the electron, ion, and CR energy densities, respectively; P 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 B is the magnetic field.Each pressure component is related to its corresponding energy density assuming an adiabatic equation of state with adiabatic index γ e = γ i = γ = 5/3 and γ CR = 4/3.In the CR energy equation (Eq.( 14)), ∇)e CR is the anisotropic diffusion flux (where b is the magnetic field unit vector), F st = f b,st u st (e CR + P CR ) is the streaming advection term, and −u st .∇PCR 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 u st = −sign(u A .∇e CR )u A , with u A = B/ 4πρ the Alfvén velocity, and where f b,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,st = 1.We set the CR diffusion coefficient to D CR = 10 29 cm 2 s −1 , which is the typical value obtained in simulated galaxies to match the observed γ-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,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 A129, page 4 of 23 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 twomoment 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 cond = (−κ e b(b.∇)T e ) is the anisotropic conductive flux, where is the Spitzer conductivity for electrons with n e = ρ/(µ e m p ) (and n i = ρ/(µ i m p )) the electron (ion) number density, µ e and µ i the electron and ion mean molecular weights, m p the proton mass, k B the Boltzmann constant, and D cond the thermal diffusivity, which is equal to with the mean free path of electrons where m e is the mass of the electron, q e the charge of the electron, and ln Λ 40 the Coulomb logarithm.This conductive flux saturates when the characteristic scale length of the gradient of temperature T e = T e /∇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 The electron energy equation has an additional term Q e↔i that transfers heat between electrons and ions and is responsible for restoring LTE, with the equilibrium timescale The radiative loss term for the total energy density L rad = L th + L CR + H CR→th contains the gas radiative loss term L th (see Sect. 3.5), and the CR radiative loss term due to Coulomb and hadronic collisions (Guo & Oh 2008) L CR = −7.51× 10 −16 (n e /cm −3 )(e CR /erg cm −3 ) erg s −1 reprocessed into the thermal pool at a rate H CR→th = 2.63 × 10 −16 (n e /cm −3 )(e CR /erg cm −3 ) erg s −1 .
The induction equation was solved with the constrained transport scheme (Teyssier et al. 2006) that guarantees ∇.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).

Refinement
The simulations were performed in a box size of 8.7 Mpc with a root grid of 64 3 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 ∆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 for levels 6 to 14. (ii) All cells within a sphere of radius 4∆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 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 jet = ρ scalar /ρ 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 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 lowdensity bubbles that would derefine under a purely Lagrangian refinement scheme, remained well refined throughout the simulation.

Initial conditions of gas and dark matter
The cluster was initialised as a cored NFW profile with a total mass of A129, page 5 of 23 where r s = r 200 /c is the scale radius, r core = 13 kpc is the core radius, ρ s = ρ c δ 200 is the density scaling of the profile, with the rescaling factor , and ρ 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 ρ DM = (1 − f gas )ρ NFW , and a gas profile of the form ρ gas = f gas ρ NFW , where f 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 × 10 14 M , 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 u 200 = 1344 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 gas = 0.103.CR pressure was initialised based on the gas pressure profile using a low and constant value of η = 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 cool was greater than the free-fall time t ff in the cluster centre, which for the current parameters leads to a core radius of r core = 134 kpc.As can be seen by comparing the dashed red line (r core = 0 kpc and M BH = 0 M ) to the dotted green lines (r core 0 kpc and M BH = 1.6 × 10 10 M ) in Fig. 2, the gravitational potential of the BH is crucial in producing an upturn in t ff /t 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 where ρ(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 where the mean molecular weights were µ i = 1.22 and µ e = 1.14 for a combined gas mean molecular weight of µ 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.

Initial conditions of the magnetic field
All magnetised simulations were initialised with a threedimensional randomly tangled magnetic field on a 512 3 grid.In order to guarantee ∇.B = 0, we set up a random magnetic potential vector A at each cell corner, and reconstructed the B field at the centre of cell faces by taking the curl of A. This initial magnetic field was then resampled onto the initial adaptive grid using a procedure that enforced ∇.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 µG and was then scaled with the initial density profile as B(r) = B 0 (ρ(r)/ρ 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.

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 with the upper limit and fit taken from Leccardi & Molendi (2008), an the lower limit being the average value of data Fit to Leccardi08: Z= 0(r/r 200 ) 0.28 Z , And89 Leccardi08, Z max = 0.45 Z , And89 Urban17, <Z> = 0.22 Z , And89 Profile used for simulation Leccardi2008 Urban2017 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).
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 H cm −3 and a temperature T < 10 4 K.This led to a stellar mass resolution of m * ,min = n H m p ∆x 3 /X H = 3.89 × 10 5 M , 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 ρ * = * ρ/t ff , where ρ is the gas density, t ff is the free-fall time of the gas cell, and * = 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 e * ,SN = m * η SN 10 50 erg M −1 in one feedback event after 10 Myr, where η 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.

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 BH = 1.65 × 10 10 M , following the M BH −M 500 relation in Phipps et al. (2019).

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 highresolution 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 to the BH trajectory according to Pellissier (in prep.).Here, x n is the position of the BH at step n, and f (x n ) = ∇φ is the gradient on the gravitational potential φ. γd = f grad √ γ d ∆t, where f grad is a dimensionless pre-factor to scale the magnitude of the acceleration, ∆t is the time step of the simulation, and 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 grad .

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, ṀBH = min( ṀBHL , 0.01 ṀEdd ), where is the Bondi-Hoyle-Lyttleton accretion rate, and is the Eddington accretion rate.G is the gravitational constant, ρ, cs , and vBH,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, σ T is the Thompson cross section, and MAD (a 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∆x of the BH location.Throughout the evolution of the BH, its spin vector a BH was followed self-consistently through accretion (see Dubois et al. 2014Dubois et al. , 2021, for technical details), for technical details).The BH was initialised as non-spinning, and wa 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.
A129, page 7 of 23 Notes.All simulations have a total cluster mass of 8 × 10 14 M and an BH with a mass of originally 1.6 × 10 10 M .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.

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 This luminosity was used to accelerate gas within a jet cylinder of width of 2∆x min and a height of 4∆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 ρ jet /ρ 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 .

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.

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 at all times.The AGN luminosity varies in the range 10 43 −10 46 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 (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×10 62 erg for MHDc in comparison to 4.4 × 10 61 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 heatbuoyancy 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 A129, page 8 of 23 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 × 10 61 erg for CRc_dsh, and 5.2 × 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 (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 A129, page 9 of 23 Fig. 5. Time evolution of the total gas mass of all gas that has T < 10 6 K (M warm + M cold , top panel), the warm gas mass (M warm , 10 4 K < T < 10 6 K, middle panel), and the cold gas mass (M 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.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, 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 of gas in this phase, the strong CR simulation CRc_dsh consistently produces more than 10 9 M 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 CRc_dsh and CRc_dsh_weak on average, showing only a variation of 30% 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 cm −3 .When it cools isobarically, gas evolves along a line with slope T ∝ 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, 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 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 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 cm −3 (grey box).As the cooling rate L th depends on the gas density, this warm and diffuse gas will cool much A129, page 10 of 23 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 tot , where V tot 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).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 7 −10 8 K into the cold star-forming phase T 10 4 K.We caution that the isochoric cooling discussed here 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 × 10 5 K), r 80 , A129, page 11 of 23 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 × 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.
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.

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 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 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- 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×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 A129, page 12 of 23 Fig. 10.Average amount of warm gas M 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 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.
as 20 Myr in the warm gas (n e = 2 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 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 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.χ > 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 = √ D 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 st = (γ CR − 1)L/u A with the diffusion timescale t D = L 2 /D CR , which leads to L = D CR /u A = 974(u A /km s −1 ) −1 kpc ≈ 35 kpc for our typical Alfvén velocity of u A = 30 km s −1 .This value is approximately comparable with the extent of diffuse warm structures seen in Fig. 4. Increasing D 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.

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 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 A129, page 13 of 23  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 (β > 1) and CR pressure (η < 1) in the hot phase, but in the cold and warm phase, CR pressure dominates thermal pressure (η > 1), and magnetic pressure contributes significantly or dominates (i.e.β < 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 η and β.However, when it is condensed out of the hot phase, all dense gas found here cools isochorically (i.e. has η > 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 η 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 AGN is injected into CR energy, only 10% are injected as kinetic energy, and CRc_dsh generally already has a lower total L 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 cool vs. η 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 A129, page 14 of 23 Fig. 13.Mass-weighted distribution of plasma β and CR pressure ratio η 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 CRpressure dominated, while all (most) of the hot gas is dominated by thermal pressure in CRc_dsh_weak (CRc_dsh).
drops and η = P CR /P therm increases further, leading to the high values of η 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 CRsaturated gas is thermally stable depends on the slope of the cooling function with respect to the temperature T , the power A ). CRs can suppress the onset of thermal instability if χ < 1, whereas thermally unstable gas remains so if χ > 1 and α 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 (χ > 1, α < 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 α 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 χ < 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 χ gas occurs within the AGN jet bubbles, in particular at their upper edges, where bubbles are rising buoyantly.This link between low χ and jet bubbles is shown more quantitatively in Fig. 17, where the lowest χ values are found for the gas that has been most recently directly affected by the AGN jet (i.e. has low t 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 η and χ (Fig. 18, hot gas only) shows that the cluster core is dominated by CR pressure out to A129, page 15 of 23 ∼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 sim-ulation, 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 (η > 1) and gas-pressure-dominated cluster outskirts (η < 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, η = 10 −4 at large radii due to the initial conditions of the simulation.The profiles in χ 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 χ < 1 on average.The fact that the radial profile of η 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 η and χ 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 η, while χ 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 selfregulation 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.

Comparison to γ-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 γ-rays emitted by hadronic processes.
To compare our cluster simulations to observations, we computed the expected hadronic γ-ray photon flux F γ = q γ (r, E γ ) dE γ dV from our simulations using the source function q γ from Pfrommer & Enßlin (2004), where σ pp = 32(0.96+ e 4.4−2.4αγ ) mbarn is the effective cross section, c is the speed of light, n N = n e (r)/(1 − 0.5X He ) is the target particle number density, n e (r) is the electron number density A129, page 16 of 23 taken from our simulation, X He = 0.24 is the fraction of helium, a γ = 4(α p − 0.5)/3 is the γ-ray spectral index, m π 0 is the pion mass, and δ γ = 0.14 α −1.6 γ + 0.44 is a shape parameter.
is an effective CR number density, where α p is the slope of the CR spectrum, ε CR (r) is the kinetic CR energy density, taken from our simulation, m p is the proton mass, and B(a, b) is the betafunction.
To compare to the observed γ-ray flux limits from Ackermann et al. (2014), we integrated over the energy range E γ = [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 γ is at its maximum for α p = 2.2, but F γ → 0 as α p → 2 and α p → 3. To establish upper limits, we adopted α p = 2.2 as our fiducial value, but also explored the impact of varying α p .
The resulting total γ-ray energy flux is shown in Fig. 19.We also show for comparison point-mass γ-ray fluxes for coolcore clusters with 6 × 10 14 M < M 200 < 10 15 M from Ackermann et al. (2014).
Figure 19 shows that CRc_dsh produces a significantly higher F γ than CRc_dsh_weak.Emission from CRc_dsh_weak is beginning to saturate from 2.5 Gyr onwards, whereas emission from CRc_dsh continues to rise.When it is compared with the observed upper limits for γ-ray fluxes from Ackermann et al. (2014), CRc_dsh quickly exceeds observational constraints by up to an order of magnitude.Emission from CRc_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, CRdominated 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 ≤ 0.1, however, produces γ-ray emission that is compatible with current observational limits.This is true for a range of different slopes of the CR spectrum 2.03 < α p < 2.73 (shaded region, covering the range of α p such that F γ > 0.25×F γ,max ).For CRc_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 A129, page 17 of 23 Fig. 18.Time evolution of η (top) and χ (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.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).

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 multiphase 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 ≤ 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 , A129, page 18 of 23 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 ∆x ∼ 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.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 cm −3 .As a consequence, the same is true for its slope α T .In general, for massive galaxy clusters, α T > 2 only around temperatures of 10 4 K.

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.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), 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 signficiant 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 find a stable configuration, with M warm + M cold continuing to increase after 3 Gyr.
By comparison, the two simulations with a higher resolution, CRc_dsh_265pc and CRc_dsh_132pc, show a good convergence.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 -density space, as is shown in Fig. D.3.However, the general evolution 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 cooling flows in our galaxy clusters are well resolved.The internal structure of the multi-phase gas continues to evolve, however.

Fig. 2 .
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 core assuming an BH with a mass of M BH = 1.6 × 10 10 M .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).

Fig. 4 .
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.

Fig. 6 .
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.

Fig. 9 .
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.

Fig. 11 .
Fig. 11.Temperature-weighted projections of η = 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).

Fig. 12 .
Fig. 12. Volume-weighted probability distribution Φ(V/V 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.

Fig. 14 .
Fig.14.Volume-weighted kernel density estimate for cooling time t cool vs. η 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.

Fig. 15 .
Fig.15.Mass-weighted distribution of χ vs. the slope of the cooling function α 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 α T < 2 for all hot gas, CRs are unable to prevent the onset of thermal instability in galaxy clusters.

Fig. 16 .
Fig. 16.Projections of χ weighted by the jet scalar (see Sect. 3.6) for CRc_dsh_weak (left panel) and CRc_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 χ < 1) is found within buoyantly rising AGN bubbles.

Fig. 17 .
Fig.17.Volume-weighted kernel density estimated of χ 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 χ values.

Fig. 19 .
Fig. 19.Time series of the total γ-ray flux within 100 kpc from the cluster centre, F γ , over time.Point-mass γ-ray fluxes for all cool-core clusters with 6 × 10 14 M < M 200 < 10 15 M from Ackermann et al. (2014) are shown as vertical lines for comparison.A self-regulating CR-dominated jet produces a γ-ray emission that exceeds observational limits, while small fractions of CR in the jet are compatible observational constraints.

Fig
Fig. A.1.Cooling function L cool (top panel) and α T (bottom panel) used in the simulations presented here, as a function of temperature T and metallicity Z.

Fig
Fig. A.1 shows that the cooling function used in the simulations presented here, based onSutherland & Dopita (1993) andRosen & Bregman (1995), is a function of temperature and metallicity, but is entirely independent of number density for 10 −2 < n < 10 2 cm −3 .As a consequence, the same is true for its slope α T .In general, for massive galaxy clusters, α T > 2 only around temperatures of 10 4 K.

Fig. B. 1 .
Fig. B.1.Time evolution of the total warm 10 4 < T < 10 6 K and cold T < 10 4 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.

Fig. B. 2 .
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.

Fig
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 A129, page 21 of 23

Fig. D. 2 .
Fig. D.2.Projected gas density showing how the multi-phase structure of the gas changes with resolution.

Fig
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 tot , where V tot 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.

Table 1 .
Simulation properties for simulations with an AGN.