© R. S. Beckmann et al. 2022

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_⊙ yr⁻¹ (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. 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 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⁴ 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, 2014, 2016; Aleksić et al. 2012; Brunetti et al. 2012, 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 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

(1)

is shorter than the free-fall time

(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_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. 2019, 2022).

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

(3)

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.

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

(4)

In the presence of magnetic fields, this becomes

(5)

where β = 8πP_therm/B² is the plasma beta and P_therm is the thermal pressure. Comparison of the two shows that 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 ℒ_th with respect to temperature T,

   (6)

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

2. the ratio of CR pressure P_CR and thermal pressure P_therm:

   (7)

3. and the ratio

   (8)

   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

(9)

where b is the magnetic field unit vector, and 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.

2.2. 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 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), 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 ρ(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.

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 non-cored 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 small-scale instabilities outside the cluster core.

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

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),

(10)

(11)

(12)

(13)

(14)

(15)

where ρ is the mass density; u is the velocity; e = 0.5ρu² + e_e + e_i + e_CR + B²/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)), F_d = −D_CRb(b.∇)e_CR is the anisotropic diffusion flux (where b is the magnetic field unit vector), F_st = f_b, stu_st(e_CR + P_CR) is the streaming advection term, and −u_st.∇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 u_st = −sign(u_A.∇e_CR)u_A, with 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²⁹ cm² s⁻¹, 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 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 thermo-dynamical equilibrium (LTE) effects into account. The term F_cond = ( − κ_eb(b.∇)T_e) is the anisotropic conductive flux, where

(16)

is the Spitzer conductivity for electrons with n_e = ρ/(μ_em_p) (and n_i = ρ/(μ_im_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

(17)

with the mean free path of electrons

(18)

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 ℓ_Te = 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

(19)

The electron energy equation has an additional term Q_e ↔ i that transfers heat between electrons and ions and is responsible for restoring LTE,

(20)

with the equilibrium timescale

(21)

The radiative loss term for the total energy density ℒ_rad = ℒ_th + ℒ_CR + ℋ_{CR → th} contains the gas radiative loss term ℒ_th (see Sect. 3.5), and the CR radiative loss term due to Coulomb and hadronic collisions (Guo & Oh 2008) ℒ_CR = −7.51 × 10⁻¹⁶(n_e/cm⁻³)(e_CR/erg cm⁻³) erg s⁻¹ reprocessed into the thermal pool at a rate ℋ_{CR → th} = 2.63 × 10⁻¹⁶(n_e/cm⁻³)(e_CR/erg cm⁻³) erg s⁻¹.

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).

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 Δ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⁶ 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⁻⁴ 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 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

(22)

where r_s = r₂₀₀/c is the scale radius, r_core = 13 kpc is the core radius, ρ_s = ρ_cδ₂₀₀ is the density scaling of the profile, with the rescaling factor δ₂₀₀ = 200Ω_mc³/(3f(c)), where f(c) = ln(1 + c)−c/(1 + c), 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¹⁴ M_⊙, the parameters in the NFW profile were as follows: The concentration parameter was chosen to be c₂₀₀ = 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₂₀₀ = 1917 kpc, and a virial velocity of u₂₀₀ = 1344 km s⁻¹. 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⁻⁴.

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¹⁰ 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.

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 rcore assuming an BH with a mass of MBH = 1.6 × 1010 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).

The cluster profile was truncated at the virial radius. Outside of r₂₀₀, the density scales as

(23)

where ρ(r₂₀₀) is the density at r₂₀₀ 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

(24)

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.

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 ∇.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₀ = 20 μG and was then scaled with the initial density profile as B(r) = B₀(ρ(r)/ρ₀)^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⁴ K, with values extended below 10⁴ 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

(25)

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).

with the upper limit and fit taken from Leccardi & Molendi (2008), an the lower limit being the average value of data 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₂₀₀, 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⁻³ and a temperature T < 10⁴ K. This led to a stellar mass resolution of m_*, min = n_H m_pΔx³/X_H = 3.89 × 10⁵ 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 , 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 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.

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_BH = 1.65 × 10¹⁰ M_⊙, following the M_BH − M₅₀₀ 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

(26)

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 ϕ. , where f_grad is a dimensionless pre-factor to scale the magnitude of the acceleration, Δt is the time step of the simulation, and

(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_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,

(28)

where

(29)

is the Bondi-Hoyle-Lyttleton accretion rate, and

(30)

is the Eddington accretion rate. G is the gravitational constant, , , and 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. 2014, 2021, 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.

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

(31)

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⁵ km s⁻¹. 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_CRL_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.

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.

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⁹ to 10¹⁰ M_⊙ at all times. The AGN luminosity varies in the range 10⁴³ − 10⁴⁶ erg s⁻¹ on duty cycles of about 100 Myr (Fig. 6).

Fig. 5. Time evolution of the total gas mass of all gas that has T < 106 K (Mwarm + Mcold, top panel), the warm gas mass (Mwarm, 104 K < T < 106 K, middle panel), and the cold gas mass (Mcold, T < 104 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.

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.

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¹¹ 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⁶² erg for MHDc in comparison to 4.4 × 10⁶¹ 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 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⁶¹ erg for CRc_dsh, and 5.2 × 10⁶⁰ 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¹⁰ 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 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⁴ − 10⁶ K. In comparison to the other simulation, which contains 10⁸ − 10⁹ M_⊙ of gas in this phase, the strong CR simulation CRc_dsh consistently produces more than 10⁹ 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⁴ 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⁻³. When it cools isobarically, gas evolves along a line with slope T ∝ n⁻¹ 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).

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/Vtot, where Vtot 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).

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⁷ K and n < 0.3 cm⁻³ (grey box). As the cooling rate ℒ_th depends on the gas density, this warm and diffuse gas will cool much 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 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⁵ K), r₈₀, 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.

Fig. 8. Time evolution of the radial extent of the warm gas from the cluster centre, measured as the radius r80 that contains 80% of the total gas mass with temperatures T < 2 × 105 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.

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_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).

Fig. 9. Time series for the total warm (104 K < T < 106 K) and cold (T < 104 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.

This is in contrast with the results presented for a galaxy group by Wang et al. (2020), where the case without streaming advection and streaming heating produced a strong run-away 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⁻¹ cm⁻³) to 4.2 Gyr (n_e = 10⁻² cm⁻³) based on the cooling rate from Guo & Oh (2008) in the hot gas, and can potentially be as short as 20 Myr in the warm gas (n_e = 2 cm⁻³). 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.

Fig. 10. Average amount of warm gas Mwarm 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 Mwarm 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.

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

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

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²/D_CR, which leads to L = D_CR/u_A = 974(u_A/km s⁻¹)⁻¹ kpc ≈ 35 kpc for our typical Alfvén velocity of u_A = 30 km s⁻¹. 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.

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

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_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 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.

Fig. 13. Mass-weighted distribution of plasma β and CR pressure ratio η for CRc_dsh_weak (top) and CRc_dsh (bottom) for cold (T < 104 K, blue), warm (104 < T < 106 K, orange), and hot (T > 106 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, PCR > Pmag, 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).

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 drops and η = P_CR/P_therm increases further, leading to the high values of η in cold gas shown in Fig. 13.

Fig. 14. Volume-weighted kernel density estimate for cooling time tcool 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.

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 of the cooling function α_T = ∂log(ℒ_th)/∂log(T), and the fraction . 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.

Fig. 15. Mass-weighted distribution of χ vs. the slope of the cooling function αT for cold (T < 104 K, blue), warm (104 < T < 106 K, orange), and hot (T > 106 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.

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.

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. Volume-weighted kernel density estimated of χ vs. tAGN, 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.

The radial distribution of η and χ (Fig. 18, hot gas only) shows that the cluster core is dominated by CR pressure out to ∼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⁶ K).

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.

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⁻⁴ 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 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 γ-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),

(32)

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 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_π⁰ is the pion mass, and is a shape parameter.

(33)

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 ℬ(a, b) is the beta-function.

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 cool-core clusters with 6 × 10¹⁴ M_⊙ < M₂₀₀ < 10¹⁵ M_⊙ from Ackermann et al. (2014).

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 × 1014 M⊙ < M200 < 1015 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.

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, 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 ≤ 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 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 ≤ 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 Δ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.

Acknowledgments

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

Appendix A: Cooling function

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⁻² < n < 10² cm⁻³. As a consequence, the same is true for its slope α_T. In general, for massive galaxy clusters, α_T > 2 only around temperatures of 10⁴ K.

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

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 sub-grid drag forces due to gas and particles have been switched off to test their impact on the dynamics of the BH.

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 × 10⁶¹ erg s⁻¹ for CRc_dsh_1 to 3.5 × 10⁶¹ erg s⁻¹ 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.

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

Appendix C: Cosmic-ray parameters

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⁸ 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²⁹ cm⁻² s⁻¹, based on recent work by Chan et al. (2019) (see also Hopkins et al. 2021), to constrain 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²⁷ cm⁻² s⁻¹ (CRc_dsh_Dcr1E28) and D_CR = 10²⁸ cm⁻² s⁻¹ (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.

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 DCR nor a higher fb, st change the cooling flow properties of the clusters beyond stochastic variation.

Appendix D: Resolution study

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

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.

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.

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

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/Vtot, where Vtot 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.

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.

All Tables

All Figures

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 non-cored 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 small-scale instabilities outside the cluster core.

In the text

	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 rcore assuming an BH with a mass of MBH = 1.6 × 1010 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).
In the text

	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).
In the text

	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.
In the text

Fig. 5. Time evolution of the total gas mass of all gas that has T < 106 K (Mwarm + Mcold, top panel), the warm gas mass (Mwarm, 104 K < T < 106 K, middle panel), and the cold gas mass (Mcold, T < 104 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.

In the text

	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.
In the text

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/Vtot, where Vtot 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).

In the text

	Fig. 8. Time evolution of the radial extent of the warm gas from the cluster centre, measured as the radius r80 that contains 80% of the total gas mass with temperatures T < 2 × 105 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.
In the text

Fig. 9. Time series for the total warm (104 K < T < 106 K) and cold (T < 104 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.

In the text

Fig. 10. Average amount of warm gas Mwarm 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 Mwarm 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.

In the text

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

In the text

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

In the text

Fig. 13. Mass-weighted distribution of plasma β and CR pressure ratio η for CRc_dsh_weak (top) and CRc_dsh (bottom) for cold (T < 104 K, blue), warm (104 < T < 106 K, orange), and hot (T > 106 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, PCR > Pmag, 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).

In the text

	Fig. 14. Volume-weighted kernel density estimate for cooling time tcool 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.
In the text

Fig. 15. Mass-weighted distribution of χ vs. the slope of the cooling function αT for cold (T < 104 K, blue), warm (104 < T < 106 K, orange), and hot (T > 106 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.

In the text

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.

In the text

	Fig. 17. Volume-weighted kernel density estimated of χ vs. tAGN, 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.
In the text

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.

In the text

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 × 1014 M⊙ < M200 < 1015 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.

In the text

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

	Fig. B.1. Time evolution of the total warm 104 < T < 106 K and cold T < 104 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.
In the text

	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.
In the text

	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 DCR nor a higher fb, st change the cooling flow properties of the clusters beyond stochastic variation.
In the text

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.

In the text

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

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/Vtot, where Vtot 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.

In the text