© The Authors 2025

1. Introduction

The cooling flow problem has been one of the major drivers, pushing towards a better understanding of the key physical processes that the plasma hosted in the cores of galaxy clusters is undergoing. In the present work, we review the status of the cooling flow problem in galaxy clusters, the results obtained with the Magneticum simulations, and how the interplay between dynamical state and active galactic nuclei (AGN) feedback shapes the population of cool-core galaxy clusters in the modern universe.

As described by Fabian (2002), the intra-cluster medium gas (ICM) hosted in the central regions of galaxy clusters exhibits short adiabatic cooling times of less than 1 Gyr and mass deposition rates of ∼[100 − 1000] M_⊙/yr, which should theoretically lead to a runaway ‘catastrophic cooling’ situation with significant amounts of cooling gas and star formation. However, spectra from the XMM-Newton Reflection Grating Spectrometer (RGS) detected increasingly less emission at lower temperatures than the cooling flow model would predict, with the coldest gas detected at around ∼6 ⋅ 10⁶ K/0.5 keV, corresponding to Fe_XVII emission (Peterson et al. 2003; Sanders et al. 2010).

Similarly, star formation is detected in bright cluster galaxies (BCGs), but at levels significantly lower than those produced by the mass deposition rates of the cooling flow model. For example, Fraser-McKelvie et al. (2014) report that most (99 ± 0.6%) of BCGs in a large sample of 245 BCGs at low redshift (z < 0.1) have star formation rates (SFR) below ∼10 M_⊙/yr, with a few exceptions such as Cygnus A and Perseus. This result is confirmed by McDonald et al. (2018), based on a sample of 107 galaxy clusters with z < 0.5 (but mostly z < 0.3) with an efficiency of SFR between 1 − 10% of that predicted by the cooling flow mass deposition rates.

This situation, at odds with the theoretical predictions, is known as the ‘cooling flow problem’. To solve these discrepancies, it is widely accepted that a heating source must be present, preventing further cooling of the ICM. However, the potential solutions have to comply with a number of requirements. In the first place, it has to prevent cooling in the whole core region, not only in the most central parts. Also, it has to operate through 2 orders of magnitude in mass, from the smallest galaxy groups with mass 10¹³ M_⊙ to the most massive galaxy clusters with mass 10¹⁵ M_⊙. Finally, it has to be fine-tuned to effectively quench cooling flows without causing obvious overheating (Fabian 2003).

One alternative is stellar feedback. However, simulations have shown that there is an excess of cold gas and star formation in the core regions, even if star formation, metal enrichment, and stellar feedback are also included (Nagai et al. 2007; Borgani & Kravtsov 2011).

Another option is thermal conductivity, as a mechanism to transport heat from outside the core region into it, given the negative temperature gradients seen in the radial direction of the core regions. However, Dolag et al. (2004) showed that even a 1/3 value of the Spitzer conductivity (Spitzer 1962), which would be motivated by randomly orientated magnetic fields, does not significantly alter star formation because most of the cooling and star formation takes place at high redshift, when the temperature of the gas in the ICM is low and therefore thermal conductivity is inefficient. However, it can be an important factor to thermally stabilise the core region and promote mixing once the structure formation process reaches the scales of galaxy clusters and the ICM has temperatures in the order of keVs.

Currently, the most widely accepted mechanism to prevent cooling flows in galaxy clusters is feedback from the central AGN hosted in the brightest cluster galaxy (BCG). This concept works in a straightforward manner at high redshift (z > 1), when most of the AGNs are in quasar mode, with luminosities in the range of 10⁴⁴ − 10⁴⁵ erg/s, comparable to the ICM X-ray luminosity in the cooling region, and jets at the scale of ∼100 kpc. However, in the modern universe, most of the AGNs hosted in BCGs are in radio mode and have low luminosities of < 10⁴³ erg/s, clearly below the ICM X-ray luminosities of the cooling region (Russell et al. 2013; Fujita et al. 2014).

On the other hand, the radio emission of AGNs in the modern universe is usually coincident with cavities whose enthalpy is in the order of ∼10⁵¹ − 10⁶¹ erg, which, when divided by typical times of the system, such as the buoyant rise time, the refill time, or the sound crossing time, yields cavity powers in the range of 10⁴² − 10⁴⁵ erg/s, comparable with the ICM X-ray luminosity in the cooling region (Rafferty et al. 2006). To explain the cavity powers in the absence of radiatively efficient AGNs, Churazov et al. (2001), Churazov et al. (2005) proposed mechanical feedback, which thermalises at almost 100% efficiency, inflating bubbles that then can rise buoyantly. However, the remaining question is how the energy is isotropically transported to the entire core region.

The first implementations of active galactic nucleus (AGN) feedback in cosmological hydrodynamic simulations, by Di Matteo et al. (2005) and Springel et al. (2005), demonstrated its effectiveness in reducing the high redshift star formation and reproducing the M_BH − σ relation, with a radiative efficiency set to the canonical value of ϵ_r = 0.1, corresponding to a standard Shakura-Sunyaev disc (Shakura & Sunyaev 1973), and a relatively low feedback efficiency value of ϵ_f = 0.05. Refinements of the model by Sijacki et al. (2007) and Fabjan et al. (2010) aimed at further suppressing star formation at low redshift by introducing a radio mode with higher feedback efficiency at low accretion rates (Ṁ_BH/Ṁ_Edd < 10⁻²), as suggested by Churazov et al. (2005), to suppress cooling flows in the modern universe, where most AGNs are in radio mode and the accretion rates inferred by observations are low.

However, the simulations of Fabjan et al. (2010), which covered a wide range of masses, with M₂₀₀ in the range [0.52 − 18.51]⋅10¹⁴ M_⊙, showed that the problems at the scale of massive clusters seen in simulations are not solved by increasing the feedback efficiency in radio mode to ϵ_f = 0.2 or even ϵ_f = 0.8. In particular, the stellar mass fractions are still about three times higher than observed. On the other hand, the increased feedback efficiencies introduce differences at the scale of galaxy groups, namely reduced gas fractions and an excess of entropy in the central regions.

The overheating issues seen in simulations at the scale of galaxy groups have also been pointed out by Gaspari et al. (2013), who argued that groups are not scaled-down versions of clusters and require lower feedback efficiencies and gentler feedback to avoid ‘catastrophic heating’. Moreover, recent comparisons between simulations and large samples of galaxy groups from the first eROSITA All-Sky survey (eRASS1) by Bahar et al. (2024) show an excess of entropy in the core regions of simulated low-mass galaxy groups. At this point, it is also important to highlight that, whereas the energetics of AGN feedback in radio mode can halt the development of cooling flows in the modern universe, as inferred from cavities with enthalpy in the order of ∼10⁵¹ − 10⁶¹ ergs, the energy required to transform a cool-core cluster into a hot-core cluster is orders of magnitude higher, ∼[1 − 4]⋅10⁶¹ ergs, far beyond the most energetic AGN bursts observed in the modern universe (McCarthy et al. 2008).

On the other hand, galaxy cluster mergers are the most energetic events since the Big Bang, releasing gravitational binding energies in the order of ∼10⁶⁴ erg (Sarazin 2002). Actually, cool-core clusters are rarely found among dynamically active systems with strong evidence of merge activity. For example, one indicator of merge activity is the projected offset between the BCG and X-ray centroid, which is under 40 kpc for all low entropy systems of the Chandra ACCEPT sample (Hoffer et al. 2012) and below < 0.02 for the strongest cool-cores of the LoCuSS sample (Sanderson et al. 2009).

Moreover, Chen et al. (2007) reported a trend whereby the fraction of cool-core clusters decreases towards the most massive, dynamically young systems. This trend was also reproduced qualitatively by the simulations of Burns et al. (2008) and Planelles & Quilis (2009).

However, Burns et al. (2008) showed that although hot-core clusters are characterised by major merges at the early phases of cluster assembly (up to z < 0.5), which destroy the nascent cores, the situation is different for cool-core clusters, which undergo a smoother assembly process at early epochs and become resilient to late mergers. Also, Poole et al. (2008) studied two-body cluster mergers in simulations and showed that merges between compact cool-cores (CCC) with mass ratios of 1:3 and 1:1 appear disturbed for a period of 4 Gyr, although they recover the CCC state 3 Gyr after the relaxation time. On the other hand, Rasia et al. (2015) show that late mergers can also destroy cool-core systems if AGN feedback is included, whereas the simulations of Burns et al. (2008) and Poole et al. (2008) included only stellar feedback, not AGN feedback.

These results suggest that cluster growth via accretion and merge processes, does have a fundamental role in shaping not only the initial cool-core population but also its evolution. However, preheating via both stellar and AGN feedback is required to soften the nascent cool-cores and make them susceptible to disruption by early and late mergers. The preheating requirement is also based on previous simulations by Motl et al. (2004), which did not include any preheating and produced a cool core in almost all halos.

However, reproducing the fraction of systems hosting a cool core has been quite challenging (Barnes et al. 2018), due to observational biases and the cool-core classification criteria. From the observational side, Andrade-Santos et al. (2017) calculate that cool-core systems are over-represented in observations by a factor of 2.1−2.7 in X-ray selected samples due to the Malmquist bias and obtained a fraction of 44% cool-core systems in the HIFLUGCS X-ray selected sample using the scaled concentration parameter, versus 28% in the early Planck survey (ESZ), with a selection based on the thermal Sunyaev-Zeldovich effect (SZ), which approximates mass selection.

Moreover, the fraction of cool-core systems is highly dependent on the indicator and threshold used, which makes it difficult to compare different studies, especially between observations and simulations. For example, on the side of simulations, Burns et al. (2008) and Planelles & Quilis (2009) obtained 16% of systems hosting a cool-core, defining them as systems where the central temperature drops by at least 20% in comparison with the viral temperature, but Rasia et al. (2015) obtain 38% of cool-core systems using as definition the ratio of pseudo entropy in the 0.05−0.15 R_200c range and a threshold of 0.55. And on the side of observations, Andrade-Santos et al. (2017) obtain 28% of cool-core systems using the scaled concentration ratio in the 0.15−1.0 R_500c range and a threshold of 0.075, but 36% of cool-core systems using the concentration ratio in the 40−400 kpc range and a threshold of 0.4.

Therefore, in this work, we aim at using the same cool-core classification criteria for observations and simulations to make a more direct comparison. Additionally, we resort to a large observational data sample by combining the recent eROSITA field equatorial deep survey data (Bahar et al. 2022; Chiu et al. 2022) with the Planck SZ-selected sample from Lovisari et al. (2020). This can be compared to a statistically relevant cluster sample from the large volume simulation Box2b/hr of the Magneticum Pathfinder simulations (Dolag 2015), allowing to study the mass dependency in cool-core fractions reported by Chen et al. (2007), from the scale of galaxy groups up to the most massive clusters, and to relate this trend to the interplay between dynamical state and AGN feedback. Additionally, we also study which aspects of the AGN feedback model and its numerical implementation relate to the excess of entropy at the scale of galaxy groups reported by Fabjan et al. (2010) and Bahar et al. (2024), as well as the excess of star formation at the scale of galaxy clusters reported by Fabjan et al. (2010).

2. The simulations

In this work, we compare observational data versus the Magneticum Pathfinder simulations (Dolag 2015), which are based on the parallel cosmological Tree Particle-Mesh (PM) Smoothed-particle Hydrodynamics (SPH) code P-Gadget3, an extended version of P-Gadget2 (Springel 2005), with several new improvements as described by Beck et al. (2016): A bias-corrected, sixth order Wendland kernel with 295 neighbours (Dehnen & Aly 2012), an improved locally adaptive treatment for artificial viscosity (Dolag et al. 2005; Cullen & Dehnen 2010), and inclusion of locally adaptive artificial conduction (Wadsley et al. 2008; Price 2008) among others. The cosmological model is based on the flat ΛCDM WMAP7 cosmology (Komatsu et al. 2011), with a dimensionless Hubble constant of h = 0.704, a total matter density of Ω_m = 0.272, a baryonic density of Ω_b = 0.0456 (baryon fraction 16.76%), a spectral index of primordial scalar fluctuations of n_s = 0.963, and an amplitude of the angular power spectrum of σ₈ = 0.809.

The most relevant baryonic processes are implemented in different modules, including isotropic thermal conduction (Dolag et al. 2004) at the level of 1/20 of the Spitzer conductivity (Spitzer 1962), radiative cooling accounting for the local metallicity (Wiersma et al. 2009), and heating from a uniform but time-dependent UV/X-ray component modelling the background radiation from quasars and galaxies (Haardt & Madau 2001). Star formation is modelled by following a hybrid multiphase star formation model that also provides feedback in the form of galactic winds at a velocity of 350 km/s, produced by 10% of the massive stars triggering Type II supernova explosions that release 10⁵¹ erg (Springel & Hernquist 2003). Chemical enrichment is also included as in Tornatore et al. (2004, 2007), assuming the Chabrier (2003) initial mass function and following stellar evolution models (supernova Type Ia, Type II, and stars in AGB phase). Thermal AGN feedback is also incorporated, following an updated implementation of the Di Matteo et al. (2005) and Springel et al. (2005) model, which includes a radio mode with higher feedback efficiency as described by Fabjan et al. (2010), and some new minor changes by Hirschmann et al. (2014), such as not enforcing pinning of black hole particles to the minimum of the potential within the smoothing length and a black hole (BH) seeding criterion based on stellar mass instead of dark matter mass.

As shown in previous studies, the physics model implemented in the Magneticum simulations leads to an overall successful reproduction of the basic galaxy cluster and group properties. Among many other properties, the simulations reproduce the observable X-ray luminosity-relation (Biffi et al. 2013) of clusters, the pressure profile of the ICM (Gupta et al. 2017), and the chemical composition (Dolag et al. 2017; Biffi et al. 2018a) of the ICM. The simulations also reproduce the high concentration observed in fossil groups (Ragagnin et al. 2019) as well as the gas properties within galaxy clusters (Angelinelli et al. 2022) and between galaxy clusters (Biffi et al. 2022). This accordance also extends to the group regime when compared to eROSITA observations (see Bahar et al. 2024; Marini et al. 2024). The Magneticum Pathfinder simulations are offered in a variety of resolutions and volumes, which are publicly available on the Cosmological Web Portal¹ (Ragagnin et al. 2017).

In this work, we use data from Box2b/hr, which is the largest high-resolution (hr) simulation, with a side length of 640 h⁻¹ cMpc. The simulation is resembled by 2 ⋅ 2880³ particles with a mass of 6.9 × 10⁸ h⁻¹ M_⊙ for dark matter particles, 1.4 × 10⁸ h⁻¹ M_⊙ for gas particles, and 3.5 × 10⁷ h⁻¹ M_⊙ for star particles. The gravitational softening length is set to 3.75 kpc/h for dark matter and gas particles and to 2 kpc/h for star particles. We use the output of Box2b/hr at z = 0.25, which corresponds to the average redshift of the observational sample we want to compare with. The large volume of this simulation provides a total of 2027 halos in the galaxy cluster regime with M_500c ≥ 10¹⁴ M_⊙ and 11 656 halos in the galaxy group regime with 0.3 × 10¹⁴ M_⊙ ≤ M_500c < 10¹⁴ M_⊙ for which we can investigate global properties and internal structures.

In Fig. 1 we show X-ray surface brightness (SB) maps of 3 cool-core and 3 hot-core galaxy clusters selected from the 3 most massive bins, to give an impression of the morphologies and internal structures of the simulated galaxy clusters. The cool and hot-core classification criteria is described in Sect. 4.1 and the energy band is XMM-eFEDS as described in Sect. 3. While the hot-core clusters are clearly in an intermediate stage of a merger, the cool-core clusters are showing the effects of the AGN thermal feedback, which is visible as cavity-like features in the SB maps and, although of vastly different origin, seems to effectively resemble some morphological properties of the observed AGN jet-driven cavities in galaxy clusters.

Fig. 1. X-ray SB images of simulated cluster examples selected within the 3 highest mass bins. The upper row panels correspond to cool-core clusters, and the lower row panels to hot-core clusters.

3. Comparing observational data with simulations

Throughout this work, we use detailed cooling functions to obtain luminosity and emission-weighted averages. We generated the cooling functions using PyAtomDB v0.10.10 (Foster & Heuer 2020), which is based on AtomDB v3.0.9 (Foster et al. 2018) and the APEC model (Smith et al. 2001).

This can be obtained by calculating the spectrum, including all radiative processes (line emission, continuum free-free emission from bremsstrahlung, and pseudo-continuum from weak lines) using the PyAtomDB collisional ionisation equilibrium (CIE) module, and then integrating over the desired band to get the total power. This is done for each metallicity and temperature pair in a grid, which can then be used to interpolate the cooling function for each gas particle from the simulation.

Once the cooling function has been interpolated to the temperature and metallicity of each gas particle from the simulation, the total luminosity L_tot of a simulated cluster within the region of interest (e.g. r < R_500c) is computed for all particles within that region through

$$\begin{array}{c}\hfill L_{\mathrm{tot}}={\displaystyle \sum _{\mathrm{particles}}{L}_p=\sum _{\mathrm{particles}}\mathrm{\Lambda }(T_p,Z_p)n_{\mathrm{elec},p}{N}_{\mathrm{ion},p},}\end{array}$$(1)

where Λ(T_p, Z_p) is the interpolated cooling function using the metallicity Z_p and temperature T_p of each gas particle, n_elec, p is the electron number density, and N_ion, p is the total number of ions of the individual gas particles, which can be obtained from the atomic masses and total mass of each metal species assuming full ionization. Then the emission-weighted average temperature T_spec can be calculated via

$$\begin{array}{c}\hfill T_{\mathrm{spec}}=\frac{1}{L_{\mathrm{tot}}}{\displaystyle \sum _{\mathrm{particles}}(T_p·L_p).}\end{array}$$(2)

Notice that although Mazzotta et al. (2004) and Rasia et al. (2004) reported that emission-weighted temperature tends to be biased towards higher temperatures, they also explained that this was due to the simplification of accounting only for bremsstrahlung in the cooling functions used by cosmological simulations at the time, which is the dominant process only above 3 keV. However, our cooling functions do incorporate line emission; therefore, they should not be biased towards higher temperatures.

We consider four different bands for this work; for two of them, we don’t need to consider rest-frame band corrections because they are used to compare the luminosities versus the energy injection from the central AGN in the simulation. These are bolometric in the [0.01−100] keV band and soft in the [0.1−2.4] keV band.

The other two bands are used to compare the simulated data versus observations from clusters with a redshift of z < 0.3; therefore, we had to consider rest-frame band corrections to be consistent with the observations; these are:

• XMM-eFEDS: Using a rest-frame band of [0.55−10.62] keV, corresponding to an average redshift of z = 0.2 from the selected clusters of the Lovisari et al. (2020) and Bahar et al. (2022) samples originally in the observed bands [0.5−7.0] keV and [0.3−10.0] keV respectively.

• ACCEPT: Using a rest-frame band of [0.8−8.0] keV, corresponding to an average redshift of z = 0.14 from the selected clusters of the Cavagnolo et al. (2009) sample originally in the observed band [0.7−7.0] keV.

To obtain relative abundances, we used Asplund et al. (2009) for the XMM-eFEDS cooling function when comparing with data from the Lovisari et al. (2020) and Bahar et al. (2022) samples. On the other hand, we used relative abundances from Anders & Grevesse (1989) when comparing with data from the Chandra ACCEPT sample, as used by Cavagnolo et al. (2009).

As the simulation incorporates star formation physics, we have to carefully consider resolution elements that are describing the multiphase star-forming and interstellar medium (ISM) gas. Therefore, and following Borgani et al. (2004), we explicitly excluded particles with a temperature below 1 ⋅ 10⁵ K as well as gas particles that are star-forming (with a cold fraction above 10%). Additionally, we also excluded particles above 50 keV following Biffi et al. (2012), since these represent the cavities filled with relativistic plasma generated by the AGN feedback model.

We also have to consider a number of limitations, originating in the underlying numerical resolution of the cosmological simulations. One is the gravitational smoothing length, which in the case of Magneticum Box2b/hr is ϵ = 3.75 kpc/h for gas and dark matter particles and ϵ = 2 kpc/h for star particles. Note that this is the Plummer equivalent, and since the simulation uses a spline interpolation instead, the imposed resolution limit is 2.8 times greater (Springel 2005), which results in a gravitational resolution limit of 2.8ϵ = 2.8 ⋅ 3.75 kpc/h = 14.9 kpc.

Additionally, BHs at the centre of galaxies excite feedback to their environment via the implemented sub-grid model. Here, analogous to the gas particles, a sphere containing a fixed number of neighbouring gas particles defines the so-called ‘black hole sphere of influence’ (see discussion in Sect. 7.3). The thermodynamic properties for particles within this scale can be impacted by the numerical implementation in the form of direct thermal feedback from the central BH.

To minimise the impact caused by these numerical limitations, we followed Borgani et al. (2004) and applied a minimum number of 100 particles in every radial bin. Note that this limit is equivalent, on average, to about two times the gravitational resolution limit. Together with the exclusion of particles representing the ISM, this ensures that the results for the ICM profiles are meaningful.

Finally, for all calculations based on the simulated data, we used as the centre the position of the most gravitationally bounded particle, equivalent to the deepest point of the gravitational potential, located at the centre of the BCG. Although observational papers typically use the X-ray peak or the X-ray centroid, we argue that our choice is the most adequate to study cooling flows since colder, denser clumps are expected to precipitate in the deeper, central regions of the gravitational potential (Voit et al. 2015).

4. Cool and hot-core population statistics

4.1. Cool-core indicators

In order to classify cool and hot-core systems, many studies have resorted to properties in the most central regions, such as the central density, central cooling time, central entropy, and central cuspiness. However, an issue with these classification schemes is how well the innermost radius can be resolved in both observations and simulations.

In observations, the minimum possible scales that can be resolved depend on the angular resolution of the telescope and the redshift of the target object. In general, the average value of central properties, such as entropy or cooling time, decreases with the scale of the innermost radius that can be resolved, as reported by Panagoulia et al. (2014), Hogan et al. (2017), and Sanders et al. (2018). Also, consistent comparisons with simulations can get difficult due to the intrinsic limitations of the simulations and their resolution, as described in Sect. 3.

These shortcomings can be prevented by using indicators that cover the whole core region, typically r < 0.15 ⋅ R_500c, as we use throughout this work. For example, Santos et al. (2008) and Maughan et al. (2012) proposed the concentration parameter or core flux ratio, defined as the ratio between the bolometric flux in the core region 0.15 ⋅ R_500c and the total flux within R_500c. However, as Santos et al. (2010) explained, the value of the concentration parameter depends on the redshift due to the K-correction, which can be different for the core emission if it is softer (colder) than the total emission.

Alternatively, Rasia et al. (2015) propose the ratio of pseudo-entropy between the inner region (r < 0.05 ⋅ R₁₈₀) and the outer region (0.05 ⋅ R₁₈₀ < r < 0.15 ⋅ R₁₈₀). However, the main observational papers, including the large samples used for this work, don’t provide entropy measurements directly.

For these reasons, we use a direct comparison between the total temperature of the cluster, including the core region (r < R_500c, i.e. T_X, 500), and the temperature of the cluster, excluding the core region (0.15 ⋅ R_500c < r < R_500c, i.e. T_X, 500cex), and assess the presence of a cool-core cluster by the ratio of these two temperatures (T_ratio, 500), as shown by Equation (3):

$$\begin{array}{c}\hfill T_{\mathrm{ratio},500}=\frac{T_{\mathrm{X},500}}{T_{\mathrm{X},500\mathrm{cex}}}\{\begin{array}{c}\ge 1\to \mathrm{hot}-\mathrm{core}\hfill \\ <1\to \mathrm{cool}-\mathrm{core}.\hfill \end{array}\end{array}$$(3)

These two quantities, T_X, 500 and T_X, 500cex, are typically provided by observational papers and can also be calculated from simulation data, which facilitates the comparison. Notice that for observations these quantities are obtained by de-projection; therefore, for the simulation we used a sphere and spherical shell, respectively.

This indicator should be free of resolution issues, since we consider the whole core region by including or excluding it in the temperature measurement. Also, it does not depend on the K-Correction and, therefore, should be constant over redshift. Moreover, since we don’t use absolute values but a ratio, the effect of biases in both the observational measurements and emission-weighted average temperatures obtained from the simulation data should be minimised.

Furthermore, it does not require ad-hoc thresholds depending on the population and has a very direct and simple physical interpretation: If T_ratio, 500 < 1, then the core region is cooler than the average temperature of the cluster, and if T_ratio, 500 ≥ 1, then the core region is hotter or has the same average temperature as the rest of the cluster. For this reason we set the temperature ratio threshold to unity for the main part of this study, although in Appendix A we investigate the effect of changing it.

Finally, notice that although we propose to separate clusters into two broad categories of cool and hot-core clusters using the criteria indicated in Equation (3), this does not imply in itself that we assume the distribution of cool and hot-core clusters to be bimodal. Actually, several observational studies ruled out a bimodal distribution, such as Santos et al. (2010), Rossetti et al. (2017), Yuan & Han (2020), and Ghirardini et al. (2022), as well as studies based on simulations, such as Rasia et al. (2015) and Barnes et al. (2018). The lack of a bimodal distribution of cool and hot-core clusters suggests a rather long timescale in the transition from cool to hot-core systems.

4.2. Cool-core fractions from observational data

The next important aspect is to consider how representative an observational sample of galaxy clusters can be. X-ray selected samples are affected by the Malmquist bias because, at a fixed mass, cool-core clusters are brighter than hot-core clusters. In this context, Andrade-Santos et al. (2017) report that, for the same mass, cool-core clusters are 1.6−1.8 times more luminous than hot-core clusters, which resulted in an over-representation of cool-core clusters by a factor of 2.1−2.7 in their X-ray selected sample.

There are two ways to prevent this problem. On one hand, it is possible to resort to samples selected via the thermal Sunyaev–Zeldovich effect (SZ), which is proportional to the thermal gas pressure integrated along the line of sight and therefore approximates to a mass selection as explained by Lovisari et al. (2020). Still, the SZ signal can be boosted by shock fronts in disturbed systems. However, Planck measures the SZ signal at scales larger than R_500c and therefore should not be significantly affected by small-scale physics such as shocks, as pointed out by Rossetti et al. (2017).

On the other hand, the alternative to SZ-selected samples is to still use X-ray selected samples, but reach enough exposure so that the flux limit can be lowered and all objects above a certain redshift and above a certain mass limit can be detected. For example, the full eROSITA survey is expected to be complete for masses greater than 1 ⋅ 10¹⁴ M_⊙ for redshifts below z < 1 and for masses greater than 0.7 ⋅ 10¹³ M_⊙ for redshifts below z < 0.3 (Comparat et al. 2020).

For these reasons, we use two observational samples to obtain cool and hot-core population statistics. On one hand, we have the sample from Lovisari et al. (2020), consisting of 120 clusters from the early Planck survey (ESZ) observed with XMM-Newton. On the other hand, we can also use the eROSITA field equatorial deep survey (eFEDS), which comprises 542 groups and clusters and has a similar exposure as the full eROSITA All-Sky (Liu et al. 2022). In particular, we used the temperature and luminosity measurements from Bahar et al. (2022), matched to the corresponding weak-lensing masses from Chiu et al. (2022).

Therefore, we combined the samples from Lovisari et al. (2020) and Bahar et al. (2022), both of them providing T_X, tot and T_X, cex measurements. We selected clusters with redshift below z < 0.3, to match the last snapshot of Magneticum Box2b/hr (z = 0.25), and to be below the redshift limit for which eROSITA X-ray selected samples are expected to be complete for masses greater than 0.7 ⋅ 10¹³ M_⊙ (Comparat et al. 2020). The redshift cut produces 115 clusters from the Bahar et al. (2022) sample with an average redshift of z = 0.22 and 93 clusters from the Lovisari et al. (2020) sample with an average redshift of z = 0.17, which combined results in 172 clusters with an average redshift of z = 0.20, quite close to z = 0.25, which is the redshift of the Magnetium Box2b/hr snapshot used for this work.

Then, we divided the clusters in the combined sample with M_500c > 0.7 ⋅ 10¹³ M_⊙ into 5 bins with approximately an equal number of clusters in each bin. Additionally, we produced a bin for the groups from eFEDS with M_500c < 0.7 ⋅ 10¹³ M_⊙. For each bin, we calculated directly the number of cool-core systems using the definition given by Equation (3) and the corresponding fraction of cool-core systems. Additionally, we estimated the cool-core fraction uncertainty via bootstrapping as the standard deviation in the cool-core fraction from 1000 new random samples created from the original cluster sample for each bin, with replacement. The results are shown as black data points in Fig. 2, and the values are also reported in Table 1.

Fig. 2. Observed and simulated cool-core fractions. The gold bars correspond to the simulation, for which the temperature was obtained with emissivity weights in the XMM-eFEDS band. The grey bars for the low-mid mass range correspond to the eROSITA field equatorial deep survey eFEDS data from Bahar et al. (2022) and Chiu et al. (2022). (*Also contains 1 mid-mass cluster from the Lovisari et al. 2020 sample). The black bars for the high mass range correspond to the Planck SZ-selected sample observed with XMM-Newton from Lovisari et al. (2020). (**Also contains 10 high mass clusters from the Bahar et al. 2022 sample). The vertical dashed line corresponds to M500c = 0.7 ⋅ 1013 M⊙, above which the eFEDS survey is expected to be complete for redshifts below z < 0.3 (Comparat et al. 2020).

4.3. Cool-core fractions from Magneticum Box2b

Once we obtained an unbiased observational reference for the cool-core fractions, we wanted to estimate the cool-core fractions from the last snapshot of Magneticum Box2b/hr (z = 0.25) in a similar way so that they are comparable to the observational data. As explained in Sect. 3, we used emission-weighted averages for each cluster to obtain T_X, tot in the r < R_500c region, T_X, cex in the 0.15 ⋅ R_500c < r < R_500c region, and the corresponding T_ratio as described by Equation (3).

Then we binned the clusters, starting with the most massive ones. In the high mass range, there are 100 clusters between M_500c = [10.68 − 3.69]⋅10¹⁴ M_⊙, which we split into 2 bins of 50 clusters each. Below M_500c = 3.69 ⋅ 10¹⁴ M_⊙, there are many more clusters, which allows us to construct bins with width 1 ⋅ 10¹⁴ M_⊙ down to M_500c = 1 ⋅ 10¹⁴ M_⊙ in the middle mass range and with width 0.1 ⋅ 10¹⁴ M_⊙ down to M_500c = 0.3 ⋅ 10¹⁴ M_⊙ in the low mass range.

As done for the observational data, we calculated directly the number and fraction of cool-core systems using the definition given by Equation (3) and estimate the cool-core fraction uncertainty via bootstrapping, as the standard deviation in the cool-core fraction from 1000 new random samples created from the original cluster sample for each bin, with replacement. The results are shown as orange data points in Fig. 2, and the values are reported in Table 2.

4.4. Comparison and interpretation of cool-core fractions

In Fig. 2 we compare the fraction of cool-core systems obtained from the observational data, by combining the samples of Lovisari et al. (2020) and Bahar et al. (2022) for clusters below z < 0.3 with the fraction of cool-core systems obtained from the last snapshot of Magneticum Box2b/hr at z = 0.25. Simulation and observations coincide within the error bars. Both show a characteristic curve for the fraction of cool-core systems, which peaks around a mass scale of M_500c ≈ 10¹⁴ M_⊙ and decreases towards lower mass groups as well as towards the very massive systems in both cases. The position of the peak in the simulation data, which lies in the range M_500c = [1.00 − 1.69]⋅10¹⁴ M_⊙ is quite compatible (within the error bars) with the peak position found in observational data, which is around M_500c = [0.75 − 1.13]⋅10¹⁴ M_⊙.

The transition from cool-core systems to hot-core systems is a very complex interplay between the cooling of the ICM, the additional energy source from AGN, heat transport, and mixing, together with the merger process (see discussions in Rasia et al. 2015; Barnes et al. 2019). Therefore, it is not straightforward to interpret the shape of the obtained cool-core fractions. This is complicated by the large difference in timescales between the different processes, where the current BH accretion rate and the associated energy imprint reflect a timescale much shorter than the timescale of transition, while the overall timescale of a merger most likely reflects a much larger timescale. With this in mind, we propose an interpretation that reflects this shape through the interplay of two distinct factors, as illustrated in Fig. 3.

Fig. 3. Driving factors behind the cool-core fraction characteristic curve. The left and central panels correspond to results from the Magneticum simulation (Box2b), where cool-core clusters are shown in blue and hot-core clusters in orange. Solid lines indicate moving medians, and dashed lines indicate 16% and 84% percentiles (1σ). Left panel: Ratio between the energy injection from the central AGN feedback and the bolometric luminosity in the [0.01−100] keV band for gas particles inside the core region. Central panel: Number of mergers undergone by the central black hole of the brightest cluster galaxy (BCG). Right panel: Sketch to illustrate the concept of how the driving factors combine to produce the characteristic shape of the cool-core fraction curve.

On one hand, the relative impact of the central AGN feedback increases generally when going from more massive, bigger clusters towards less massive, smaller groups. This can be seen in the left panel of Fig. 3, which shows the ratio between the feedback power injected by the central AGN and the bolometric luminosity of the core region computed from the simulation. Therefore, the probability of transitioning from cool-core to hot-core is increased at the scale of smaller galaxy groups due to the extra energy available in the system.

The other factor is the impact of merger activity. From the simulation, we can use the number of mergers that the BH in the central galaxy of the cluster underwent, shown in the central panel of Fig. 3, as a proxy for the number of mergers of the overall cluster. This shows a clear trend increasing towards the large mass end, thereby increasing the probability of transitioning from cool-core to hot-core at the scale of massive galaxy clusters.

Therefore, we propose that the combination of these two factors can produce a characteristic curve, peaking at around the middle mass range (M_500c ≈ 1 ⋅ 10¹⁴ M_⊙), where both factors minimise the probability of a transition from cool-core to hot-core systems, thus increasing the fraction of cool-core systems, as sketched in the right panel of Fig. 3. Although the simulated cool-core fractions shown in Fig. 2 coincide with the observations within the error bars, there are indications of a small, systematic bias of the cool-core fraction in the simulation to be smaller across the entire mass range. For the sample of massive systems from Lovisari et al. (2020), the so-called hydrostatic mass bias could potentially affect the comparison between observational and simulated galaxy clusters. However, accounting for such bias, the observed data points would shift by 10%−20% towards larger masses (Biffi et al. 2016), which would actually increase the difference between simulated and observed cool-core fractions for the massive systems. Moreover, this should not affect the eFEDS data at low mass, since these masses have been obtained via weak gravitational lensing (Chiu et al. 2022).

At the scale of smaller groups, the simulated cool-core fractions decline sharply in comparison with the observations. Possible observational explanations for this discrepancy can be lack of data or undetected hot-core systems in the lower mass group regime (see discussion in Marini et al. 2024). On the other hand, the increased impact of AGN feedback in small mass groups, as seen in the simulation, is consistent with the excess of entropy in the core region (r < 0.15 ⋅ R_500c) reported in the first eROSITA All-Sky survey (Bahar et al. 2024), which accounted for selection effects.

These findings suggest that the details in the numerical implementations of the coupling of the AGN energy to the surrounding ICM, could affect the predicted cool-core fractions from the simulation, especially for small groups, since at this mass range the injected energy can significantly exceed the cooling losses, as can be seen in the left panel of Fig. 3. We look into further insights and possible ways to improve the numerical implementation of AGN feedback in cosmological simulations in Sect. 7.2.

Finally, in Appendix A, we study how the cool-core fractions depend on the choice of the threshold adopted for the temperature ratio, which we set to unity for the main part of the study. Changing the temperature ratio (lowering it) provides a gauge for the mid-to-strong cool-core population, which is more sensitive to the fine-tuning of the simulation and also is less robust in statistical terms due to the reduced population of strong cool cores but still provides a valuable insight to interpret the results overall.

4.5. Dynamical analysis

Having a set of simulated groups and clusters which show good agreement with the observed cool-core fractions, we can now investigate the impact that merger activity has on cool-core and hot-core groups and clusters in the simulation and how that contributes to the observed mass trend. To gauge the impact of recent merge activity, we consider the relation between the kinetic (K_{c, component}), internal (U_{c, component}), and gravitational potential (W_{c, component}) energies, which we computed as a sum over the contributions of the individual particles. Here the sub-index c refers to the core region (r < R_c = 0.15 ⋅ R_500c) and the sub-index component refers to either total, which means all particles (stellar, dark matter, gas, and BH) or only gas particles. We used the potential and internal energy provided directly by the snapshot output for each particle, but for the kinetic energy, we defined the velocities w.r.t. the velocity of the centre of mass for all particles inside R_200c, which is expected to be close to the viral radius of the system.

We define the freedom ratio as the ratio between the sum of kinetic and internal energy divided by the potential energy. We calculated two different freedom ratios, one including all particles F_c, total = (K_c, all + U_c, all)/|W_c, all| and one including only gas particles F_c, gas. Notice that this is not the classical ratio from the Virial theorem, which can only be obtained by either accounting for all particles in the system (not a specific region) or by accounting for boundary conditions around the region of interest (e.g. external potential, surface pressure), which are difficult to estimate (Davis et al. 2011; Klypin et al. 2016). On the other hand, the core freedom ratio can be easily obtained and has a clear interpretation in terms of the average probability of particles in the core region to move to the outer regions of the halo since, for collisionless, non-gas particles, the freedom ratio represents the square of the particle velocity divided by the escape velocity.

Additionally, we computed the kinetic energy fraction, defined as the ratio between the kinetic and the total (kinetic + internal) energy f_{c, kinetic, gas} = K_c, gas/(K_c, gas + U_c, gas). We show the results in Fig. 4. The left panel shows the core freedom ratio for all particles, which allows us to gauge the effect of merge (dynamical) activity independent from the AGN feedback. We see that the core freedom ratio increases towards the most massive clusters, indicating an increasing presence of high-energy particles injected from the merge activity. This confirms that more massive systems are also dynamically young systems (see also discussion in Chen et al. 2007). However, there is no difference between cool-core and hot-core systems regarding the total composition inside the cores of groups and clusters.

Fig. 4. Dynamical analysis. Cool-core clusters are shown in blue, and hot-core clusters in orange. Solid lines indicate moving medians, and dashed lines indicate 16% and 84% percentiles (2σ) from the Magneticum simulation (Box2b/hr). Left panel: Freedom ratio for all particles inside the core region (Fc, total). Central panel: Freedom ratio for gas particles inside the core region (Fc, gas). Right panel: Kinetic energy fraction fc, gas, kinetic for gas particles inside the core region. The black bar indicates the first and only measurement of kinetic energy fraction in galaxy clusters, of 10% measured for the Perseus cluster by Hitomi (Hitomi Collaboration 2018).

This situation changes when only considering the gas particles within the core region shown in the central panel of Fig. 4. While the gas particles show a very similar trend with the mass of the system, there is a clear offset of hot-core clusters having larger total energies than cool-core systems. Interestingly, especially in low-mass systems, the hot-core clusters show a large scatter towards high total energy. This is mainly due to a larger internal energy fraction and again points towards the AGN feedback driving the transition from cool-core to hot-core for lower mass systems.

On the other hand, this offset is of similar strength at all masses, from which we can conclude that the reduced cool-core fraction at higher masses does not directly come from increased merger activity. A clarification of the situation can be seen in the right panel of Fig. 4, where we show the kinetic energy fraction decreases towards the higher mass end, thus indicating that the thermalisation of kinetic energy injected by merge processes is more efficient towards higher mass clusters, thus establishing the first step to increase the core entropy. Then the implementation of thermal conduction in the Magneticum simulations, which has a strong dependency on temperature following the Spitzer (1962) description, provides a more efficient energy transport at larger masses due to the tight connection between cluster mass and temperature as explained in Sect. 4.6. Therefore, the combination of these factors can effectively reduce the cool-core fractions towards massive clusters.

However, the prediction for kinetic energy fractions from the Magneticum simulations is difficult to verify since current X-ray telescopes cannot measure the kinetic energy support (bulk velocity and velocity dispersion), although future missions like Athenea are designed for this (Roncarelli et al. 2018). So far, there are only the Hitomi observations from Perseus (M_500c ∼ 6 ⋅ 10¹⁴ M_⊙), from which (Hitomi Collaboration 2018) derived a ratio between the kinetic and thermal energy in a range of [11−13]%, accounting for geometry corrections to the velocity dispersion as suggested by ZuHone et al. (2012). This is equivalent to a kinetic energy fraction of f_{kinetic, gas} = 0.11 ± 0.08 using our definition of f_{c, kinetic, gas} and as shown by the black bar in the right panel of Fig. 4, which is slightly lower than the predictions of the simulations, although overlapping within the error bars. This is in line with previous findings that simulations over-predict the amount of kinetic energy in the ICM when compared to observations of relaxed galaxy clusters (Eckert et al. 2019). Soon, XRISM will largely increase such measures for massive clusters; however, extending this to the group regime to investigate the mass trend found in the Magneticum simulations will only be possible with possible future instruments like Athena.

4.6. Effect of thermal conductivity

As mentioned in the previous section, the kinetic energy injected by the merge processes has to be first thermalised but also mixed through the core of the ICM. In this sense, the Magneticum simulations provide an implementation of isotropic physical thermal conductivity following the Spitzer (1962) description:

$$\begin{array}{c}\hfill \kappa =4.6·{10}^{13}{\left(\frac{T_{\mathrm{e}}}{{10}^8\mathrm{K}}\right)}^{5/2}\frac{40}{ln\mathrm{\Lambda }}\frac{\mathrm{erg}}{\mathrm{s}\mathrm{cm}\mathrm{K}},\end{array}$$(4)

where T_e is the electron temperature and Λ the Coulomb logarithm set to 37.8, which is an appropriate value for core regions of galaxy clusters (see Arth et al. 2014). Although the effective conductivity coefficient is reduced to 1/20th of the total value obtained from Eq. (4) to account for suppression of the conductivity due to magnetic fields, still the effective value of the conductivity coefficient can be significantly larger for massive clusters due to the strong temperature dependence (T^5/2) and the tight connection between cluster mass and temperature according to scaling relations (T ∼ M^2/3, Böhringer et al. 2012).

In this sense, we show in Fig. 5 the emission-weighted average temperature of the cool-core and hot-core groups and clusters (left panel) and the corresponding value of the effective Spitzer conductivity coefficient for the core temperature, normalised to the value for a system at 1 keV (κ_norm = κ[T_500c]/κ[1 keV], right panel). We see that the median effective value of the conductivity can be almost one order of magnitude higher for the most massive clusters in comparison with smaller mass clusters and groups, thus increasing the mixing of the energy injected by mergers through the core regions and facilitating the conversion of cool-core clusters to hot-core clusters towards the higher mass end.

Fig. 5. Thermal conductivity. Cool-core clusters are shown in blue, and hot-core clusters in orange. Solid lines indicate moving medians, and dashed lines indicate 16% and 84% percentiles (2σ) from the Magneticum simulation (Box2b/hr). Left panel: Temperature of the core region obtained with emissivity weights in the XMM-eFEDS band. Right panel: Effective Spitzer conductivity coefficient for the core temperature, normalised to the value for a system at 1 keV.

5. Radial profiles

Since the simulated cool-core fractions are consistent with the observed ones, we can now investigate the temperature, density, and entropy profiles. To ensure that profiles are resolved properly, we restricted here to clusters with M_500c > 10¹⁴ M_⊙ and divide the 2027 clusters in this sample into 3 broad mass ranges: low: 10¹⁴ M_⊙ < M_500c ≤ 2.69 × 10¹⁴ M_⊙, medium: 2.69 × 10¹⁴ M_⊙ < M_500c ≤ 4.88 × 10¹⁴ M_⊙ and high: 4.88 × 10¹⁴ M_⊙ < M_500c ≤ 9.02 × 10¹⁴ M_⊙, corresponding to the 3 most massive bins of the characteristic cool-core fraction curve presented in Sect. 4.3. For computing the simulated profiles, we used a series of radial annuli in increments of 0.01 ⋅ R_500c, but for the innermost annulus, we imposed a minimum physical radius of 15 kpc above the gravitational resolution limit and a minimum of 100 particles, as described in Sect. 3.

We compare the simulated profiles with median profiles obtained from the Chandra ACCEPT sample (Cavagnolo et al. 2009). For this, we first categorised the Chandra ACCEPT clusters into the previously mentioned mass bins using the corresponding cross-matched masses from the M2C Galaxy Cluster Database. This results in 35, 51, and 50 clusters in the low, median, and high mass bins, respectively, where the observed clusters in these 3 bins have a median redshift of z = 0.08, z = 0.14, and z = 0.20, respectively. The total resulting sample consists of 136 clusters with an average redshift of z = 0.14. To reduce the noise for the constructed median profiles, we drew 100 random instances of each profile, assuming a Gaussian distribution for the errors, and construct a collection of random instances from all clusters included in each mass bin, from which we calculated the overall median and 1-sigma percentiles for each radial bin. In this way, we account for the dispersion in the population but also take into account the measurement errors. Additionally, we also consider a correction of +20% in the masses retrieved from the M2C Galaxy Cluster Database to account for the hydrostatic mass bias (Biffi et al. 2016) and construct two sets of median profiles with the original and corrected observational masses.

The Chandra ACCEPT profiles don’t typically extend up to R_500c for the low-redshift clusters, as this would be outside of the field of view (FoV) of the Chandra Advanced CCD Imaging Spectrometer (ACIS). Therefore, we computed the temperature ratio in two radial ranges corresponding to r < 0.1 ⋅ R_500c and 0.1 ⋅ R_500c < r < 0.2 ⋅ R_500c to classify the clusters from the Chandra ACCEPT sample as cool-core or hot-core.

Provided that we compare the simulated profiles with the Chandra ACCEPT sample, we followed in general the same procedures and definitions for the individual quantities described by Cavagnolo et al. (2009). However, note that our centring procedure always chooses the position of the most gravitationally bounded particle, equivalent to the deepest point of the gravitational potential, whereas the centring procedure used for ACCEPT uses the X-ray peak as the default centre but switches to the X-ray centroid if it is separated by more than 70 kpc from the peak (Cavagnolo et al. 2008).

Since hot-core clusters are typically disturbed, the ACCEPT centring procedure switches to the X-Ray centroid, whereas cool-core clusters are typically relaxed, and it defaults to the X-Ray peak. This procedure results in two very distinct types of profiles for cool-core and hot-core clusters and can be behind the bimodal distribution reported for ACCEPT. On the other hand, our centring procedure produces more similar profiles for cool-core and hot-core clusters but can better track the development of cooling flows, which are expected to precipitate in the deepest regions of the potential.

5.1. Radial temperature profiles

For the temperature profile, we used emission-weighted average temperatures as before, where the emissivity weights are adapted to the energy range of the observations. Following the approach used in the observations, we computed the projected radial temperature profiles out to R_500c in cylindrical shells as done in Cavagnolo et al. (2009), where we used a depth in the z direction corresponding to R_200c.

The comparison of the constructed profiles from the observations with the results from the simulations is shown in Fig. 6. In general, there is a good agreement within the intrinsic error bars between them. Here the simulations reproduce the same characteristic shapes, which differ for cool-core and hot-core systems. Especially the simulated cool-core systems show the same characteristic temperature drop towards the centre to roughly half of the maximum temperature, while the hot-core systems show an almost isothermal core, albeit in the low- and mid-mass range, they still have traces of a peak or turnover similar to the cool-core systems.

Fig. 6. Projected X-ray temperature profiles. Corresponding to the [0.7−7.0] keV band and normalised by the mean X-ray temperature in the range [0.1 − 0.2] R500c. The error bars correspond to the median Chandra ACCEPT sample profiles and ±1σ intervals using the original masses from the M2C Galaxy Cluster Database (black) and a +20% hydrostatic mass bias correction (orange). The blue line and shaded areas correspond to the Magneticum simulation (Box2b) median profiles and ±1σ intervals. Columns are sorted by mass range, as indicated on top of each panel (increasing mass range from left to right). The upper row panels correspond to cool-core clusters, and the lower row panels to hot-core clusters.

The hot-core temperature profiles from the ACCEPT sample show a stronger isothermal core; however, this roots in the previously mentioned centring procedure used for ACCEPT sample that selects the X-ray centroid as the centre if it is separated by more than 70 kpc from the peak, which is typically the case for hot-core clusters. The X-ray centroid is usually located in between the dominant galaxies and corresponds to a rather shallow part of the gravitational potential, resulting in a lack of temperature gradients (isothermal) structure. On the other hand, our simulated clusters are centring at the deepest point of the gravitational potential, which enhances temperature gradients and the accumulation of denser cold clumps. In Appendix B we show how the temperature profiles would look like when following the ACCEPT centring procedure.

Additionally, the hot-core clusters of the low- and mid-mass range appear to be slightly overheated in the innermost regions, which exhibit a ‘hot bubble’ as a remnant of the AGN feedback implementation in the simulations. This hot bubble is preceded by a peak or turnover similar to the cool-core systems and indicates that the AGN feedback energy is not transported effectively away from the immediate surroundings of the AGN where the feedback is initially injected. In the case of Magneticum, this problem is partially addressed by the physical conductivity model; however, it strongly depends on the overall temperature of the system, which decreases towards low-mass systems as explained in Sect. 4.6, and moreover the current model applies a strong suppression down to 1/20 of the Spitzer value assuming turbulent magnetic fields (see Arth et al. 2014) which we plan to review in future works. Also, the implemented AGN feedback could be too strong, especially in low mass systems, as already reported by Vogelsberger et al. (2013), who argued that the residual black hole accretion, when there is no star-forming gas in the vicinity of the black hole, can create artificial hot bubbles around the black hole as an artefact of the imperfections of the sub-grid black hole accretion and feedback model.

Finally, we clearly see that for the cool-core clusters, the median temperature profiles from the ACCEPT sample peak at a larger radius than the simulated clusters. This is due to the significant presence in the ACCEPT sample of clusters with increasing temperatures outwards beyond the core region, such as Abell 1835, Abell 2142, and PKS 0745−191.

5.2. Radial density profiles

While the differences in the cool-core and hot-core temperature profiles are not fully independent from the cool-core criteria used, it is worthwhile to check also the gas density profile and their systematic differences between the two cluster populations. This is especially important, as the density is the fundamental thermodynamic property of the ICM, which is directly linked with the development of cooling flows, star formation, and black hole accretion.

As the density profiles from the Chandra ACCEPT sample are de-projected profiles (Cavagnolo et al. 2009), we used here directly spherical, radial density profiles from the clusters selected in the Magneticum Box2b/hr. In addition, to minimise the impact of sub-structures, which would be typically masked in the observations, we obtained the mean density in each radial bin by summing the mass of all gas particles and dividing by the total volume of the spherical radial bin shell instead of averaging the individual densities and assuming full ionisation of the gas to obtain the electron number density.

The results of this comparison are shown in Fig. 7. Interestingly, the simulated clusters that are classified as cool-core systems show systematically higher central densities compared to the hot-core clusters, well in agreement with the observations.

Fig. 7. Electron number density profiles. The error bars correspond to the median Chandra ACCEPT sample profiles and ±1σ intervals using the original masses from the M2C Galaxy Cluster Database (black) and a +20% hydrostatic mass bias correction (orange). The blue line and shaded areas correspond to the Magneticum simulation (Box2b) median profiles and ±1σ intervals. Columns are sorted by mass range, as indicated on top of each panel (increasing mass range from left to right). The upper row panels correspond to cool-core clusters, and the lower row panels to hot-core clusters.

Generally the simulated profiles of both cool and hot-core clusters show increasingly lower densities, towards the high mass range in comparison with the median profiles from the ACCEPT sample. This discrepancy is reduced when considering the observational masses corrected by the hydrostatic mass bias but not entirely removed. A possible explanation for this behaviour can be associated with excessive star formation towards more massive systems, which has not been effectively quenched by the central AGN and has therefore led to lower gas fractions in the baryonic component. We elaborate more on this scenario in Sect. 6.

5.3. Radial entropy profiles

Finally we constructed gas entropy profiles in the same way as for the Chandra ACCEPT sample, by combining the projected temperature profiles with the spherical density profiles (see Cavagnolo et al. 2009). The results are shown in Fig. 8, where we see that the simulated profiles for the cool-core systems typically decline towards the centre with a central entropy of less than 100 kev/cm², while the hot-core systems typically have a flat entropy profile with a central entropy of larger than 100 kev/cm². This compares very well with the observational trend.

Fig. 8. Electron entropy profiles. The error bars correspond to the median Chandra ACCEPT sample profiles and ±1σ intervals using the original masses from the M2C Galaxy Cluster Database (black) and a +20% hydrostatic mass bias correction (orange). The blue line and shaded areas correspond to the Magneticum simulation (Box2b) median profiles and ±1σ intervals. Columns are sorted by mass range, as indicated on top of each panel (increasing mass range from left to right). The upper row panels correspond to cool-core clusters, and the lower row panels to hot-core clusters.

Driven by the trends in the density profiles, the simulated cool and hot-core clusters show increasingly higher entropy profiles, especially for the high mass bins. As for the densities, assuming a 20% hydrostatic mass bias in the observations brings the simulated and observational profiles closer, although not entirely overlapping. In particular, the simulations consistently produce a flatter profile at the intermediate radii before steepening in the core, diverging from the observed power-law-like shape.

6. Gas and stellar mass fractions

It is well established that the AGN feedback not only suppresses the star formation in very massive systems, it also impacts the gas mass fraction in clusters (see for example Planelles et al. 2013). Here, the overall gas fraction within R_500c of the clusters and groups as predicted by the Magneticum simulations agrees with observationally derived fractions from X-ray observations, both in absolute numbers and also in the trend of having lower gas mass fractions in lower mass systems (Angelinelli et al. 2022, 2023).

In contrast, stellar masses typically are still significantly larger than observationally inferred stellar masses, especially at the high mass end, pointing towards an inefficient suppression of star formation by the AGN feedback implementation in simulations for very massive systems. A similar result was reported by Fabjan et al. (2010), who compared the stellar mass fractions at the scale of R_500c and found them 2−3 times higher than the observational measurements from Gonzalez et al. (2007), which also included the contribution from the intra-cluster light (ICL). To investigate this point in more detail, we evaluate the gas mass fraction, the stellar mass fraction, and the total baryonic mass fraction within R_2500c² to compare to Laganá et al. (2013), who presented the gas and stellar mass fractions at R_2500c for a sample of 27 clusters with average redshift z = 0.22, quite close to z = 0.25 of our sample.

Laganá et al. (2013) obtained the gas mass (M_2500, gas) by modelling the X-ray emission with a modified β-model profile (Maughan et al. 2012), which is then integrated up to R_2500c. The total mass (M_{2500, total}) is calculated assuming hydrostatic equilibrium and isothermality. The stellar mass (M_{2500, stellar}) is calculated by converting apparent magnitudes to absolute magnitudes using K-corrections depending on the morphological type and using two different mass-to-light ratios for early-type and late-type galaxies. The corresponding gas, stellar, and baryon mass fractions (f_2500, gas, f_{2500, stellar}, and f_2500, bar, respectively) are obtained as the ratios between the mass of each component and the total mass, and we calculated the errors by standard error propagation of the 1σ uncertainty in quadrature.

From the simulation, we directly obtained the gas, stellar, and total mass of the components contained in the sphere of R_2500c. The results are shown in Fig. 9, where we see that the simulated gas fractions are systematically lower than the observed gas fractions, as we would have expected from the comparison of the radial density profiles. Interestingly, similar to the findings at R_500c, there is a clear mass trend predicted by the simulations, which, however, is not reflected in the observations. On the other hand, the simulated stellar mass fractions within R_2500c are much higher than the observed one. Being 5 times larger than the observational stellar mass fraction at the very massive end, this strongly deviates from the observations more than the overall stellar mass fraction within R_500c. Surprisingly, the total baryonic mass fraction seems to not be biased in comparison with observations, although the observational results present a much wider scatter than the simulated data. Also worth noting is that the simulation does not show any significant difference between cool-core and hot-core systems.

Fig. 9. Stellar, gas, and baryonic mass fractions. Cool-core clusters are shown in blue, and hot-core clusters are shown in orange. Solid lines indicate moving medians, and dashed lines represent 16% and 84% percentiles (2σ) from the Magneticum simulation (Box2b). Observational data is shown in black, with error bars at the 1σ level, based on XMM-Newton, Chandra, and SDSS from Laganá et al. (2013). Left panel: Gas mass fraction inside R2500c. Central panel: Stellar mass fraction inside R2500c. Right panel: Baryonic mass fraction inside R2500c.

While there are still uncertainties in the observations (see detailed discussion in Laganá et al. 2013), both in the gas mass based on the assumption of hydrostatic equilibrium and isothermality, as well as the stellar mass, where undetected ICL could contribute an additional 10% to 40% to the total stellar mass, this would not solve the observed differences between the simulation and observations. The fact that the total baryonic mass fraction within R_2500c aligns with the observations but there are too many stars formed within the central regions of clusters points towards a deficit in the detailed coupling of the AGN feedback within the simulations rather than a vastly different energy injection by the AGN, which would lift more gas to larger distances. This is also in agreement with the fact that the AGN luminosity function, which reflects the general energy available for the AGN feedback in the Magneticum simulations, agrees well with observations (Hirschmann et al. 2014; Biffi et al. 2018b). Therefore, this points towards the coupling of the feedback from the central AGN within the simulations to the surrounding medium, which is not able to fully suppress star formation at the scale of clusters, in line with the previous findings reported by Fabjan et al. (2010).

7. Implications for the AGN feedback model

The results presented so far are consistent with earlier findings (Fabjan et al. 2010), showing that the AGN feedback model implemented in the simulations is significantly suppressing star formation at the scale of groups and clusters. However, it only partially prevents catastrophic cooling and star formation at the scale of massive galaxy clusters, despite increased AGN feedback efficiencies in radio mode. To investigate the possible origin of this problem, we can compare the effective implementation within the simulations to the observed signatures of AGN feedback in galaxy groups and clusters.

7.1. AGN accretion and energy output

There are currently two models for how cool-core systems are formed. One is described as precipitation, where the deposition of cold gas onto the core is driven by thermal instability (TI, McCourt et al. 2012; Voit et al. 2015). In the other one, the raining of cold gas onto the core is driven by chaotic cold accretion (CCA, Gaspari et al. 2018). The onset of these mechanisms is driven by the ratio of cooling time to free-fall time below a typical threshold of t_cool/t_ff ≲ 10 or the ratio of cooling time over the eddy turnover time of t_cool/t_eddy ≲ 1, respectively. In practice, these two criteria are almost identical, and observations indicate that they correspond to a central entropy threshold of ≈35 keV/cm² across a wide range of redshifts (McDonald et al. 2013). This is in line with the entropy profiles of simulated cool-core clusters being decreasing towards the centre, while the entropy profiles of hot-core clusters show more flattening towards the centre at values above 100 keV/cm². In cool-core clusters, we then expect that a feedback loop is established that prevents low-entropy cooling flows from developing further (Churazov et al. 2005).

As most of the AGN treatment in cosmological simulation, the AGN accretion model used by Magneticum (Hirschmann et al. 2014) is based on the Bondi-Hoyle-Lyttleton model (Hoyle & Lyttleton 1939; Bondi & Hoyle 1944; Bondi 1952), with an implementation that follows Di Matteo et al. (2005), Springel et al. (2005), and Fabjan et al. (2010) and has the form

$$\begin{array}{c}\hfill {\dot{M}}_{\mathrm{BH}}=\alpha \frac{4\pi G^2{M}_{\mathrm{BH}}^2{\rho }_{\mathrm{g}}}{{(c_{\mathrm{s}}^2+v^2)}^{3/2}},\end{array}$$(5)

where G is the gravitational constant, M_BH is the mass of the BH, ρ_g is the gas density, c_s is the speed of sound, v is the black hole velocity with respect to the surrounding gas, and α is a dimensionless parameter (boost factor) used by the simulation, typically set to 100 to account for the unresolved increase in density towards the central regions surrounding the black hole. Since c_s ∼ T^1/2 and the entropy is defined as S ∼ Tρ^−2/3, the Ṁ_B predicted by the Bondi-Hoyle-Lyttleton accretion formula inversely depends on the entropy to the power 3/2, implying that aside from the gas velocity, AGN accretion is larger in the simulation when the entropy of the surrounding medium is low.

How well this implementation describes the general accretion onto the central BH in groups and clusters is shown in the left panel of Fig. 10, where we compare observationally inferred Ṁ_BH as a function of halo mass with the results from the simulations. To do so, we resort to the measurements from Fujita et al. (2014), who obtained the Bondi accretion rate of a sample of BCGs in the near universe (z < 0.35) by modelling the temperature and gravitational contribution from the dark matter halo, galaxy, and central black hole, to then integrate the hydrostatic equation in order to derive the gas density down to the Bondi radius. The obtained accretion of the AGNs hosted in the centres of BCGs agrees very well with the ones predicted by the simulations, indicating that the use of the Bondi-Hoyle-Lyttleton model in the simulations is an appropriate approximation. Also clear to see is that in the simulation the accretion onto the central BHs is systematically larger than in hot-core systems.

Fig. 10. Accretion rates and energy injection of the AGN hosted in the centres on BCGs. Cool-core clusters are shown in blue, and hot-core clusters are shown in orange. Solid lines indicate moving medians, and dashed lines represent 16% and 84% percentiles (1σ) from the Magneticum simulation (Box2b/hr). Observational data is shown with error bars at the 1σ level. Left panel: Comparison of the accretion rates from the AGN hosted in the centres on BCGs, with the estimations from Fujita et al. (2014) and Russell et al. (2013). Right panel: Comparison of black hole energy injection with the observational estimations from Rafferty et al. (2006), Russell et al. (2013), Eckert et al. (2021) and Pasini et al. (2022).

In the simulation the accretion rate of the BH is converted into the energy output (feedback power) $P_{\mathrm{feedback},\mathrm{BH}}={\dot{M}}_{\mathrm{BH}}{c}^2{ϵ}_{\mathrm{r}}{ϵ}_{\mathrm{f}}$, by assuming a radiative efficiency (ϵ_r) which determines what fraction of the accreted mass is converted into energy, and a feedback efficiency (ϵ_f) which determines how much energy released by the black hole accretion is deposited in the surrounding medium. In our case, they are chosen to be 0.2 and 0.15 in the quasar mode, while the latter is 4 times larger in the radio mode (Fabjan et al. 2010). The radio mode is assumed whenever the actual accretion rate is below one hundredth of the Eddington accretion rate. Therefore, in radio mode, the total efficiency reaches 0.12, close to the 10% as inferred in galaxy clusters by Churazov et al. (2005) to compensate for the cooling losses of the ICM.

In the right panel of Fig. 10, we compare the power in observed cavities obtained by Rafferty et al. (2006), Russell et al. (2013), and Eckert et al. (2021) within galaxy clusters with the feedback energy released from the central AGN in the simulation. The observed cavity power (P_cavity) is derived by dividing the cavity enthalpy.

$$\begin{array}{c}\hfill E_{\mathrm{cav}}=\frac{{\gamma }_{\mathrm{c}}}{{\gamma }_{\mathrm{c}}-1}{p}_{\mathrm{s}}{V}_{\mathrm{c}},\end{array}$$(6)

by the buoyancy time scale. Here γ_c is the adiabatic index of the material filling the cavity, V_c is the cavity volume, and p_s is the pressure at the cavity surface. As these cavities are filled with relativistic material, γ_c = 4/3 and E_cav = 4V_cp_s. Additionally, in this work, we apply an x2 factor to account for the shock energy as estimated by Rafferty et al. (2006) and confirmed by the simulations of Guo & Mathews (2010). We also add to the comparison the kinetic luminosity sample inferred by Pasini et al. (2022) based on LOFAR 144 MHz radio power observations of the BCG in the eFEDS sample. We can see that the energy injected in the simulation agrees well with the observations for massive clusters; however, at the scale of galaxy groups, the injected energy is significantly larger than the observed values. The simulation also shows a clear trend that cool-core systems inject a larger amount of energy than hot-core systems.

A possible explanation for the divergence at the scale of galaxy groups, could be that the previously mentioned x2 factor to account for the shock energy is underestimated at the scale of groups. However, the fact that this divergence is also visible for the kinetic luminosity sample from Pasini et al. (2022) indicates that the observed cavity powers cannot be significantly underestimated at the scales of galaxy groups since we would expect at least higher radio emission even if the cavities are not detectable.

In addition, we show in the left panel of Fig. 11 the soft band [0.1−2.4] keV luminosity of the ICM (L_500, soft) from the observational samples of Lovisari et al. (2020) and Bahar et al. (2022)³, as well as the cavity powers previously mentioned. We see that the observational cavity powers fluctuate around the ICM luminosity but are not systematically above. On the other hand, the ICM luminosity of the simulation matches the observations, but the energy input by the AGN feedback of the simulation is systematically above the ICM luminosity, matches the cavity powers only at the scales of very massive clusters, and is much larger at low mass clusters and group scales. This bias is consistent with the bias towards lower cool-core fractions in simulations at the scale of galaxy groups, as shown in Fig. 2. Since in the simulation, as well as in observations, almost all AGNs in the centres of groups and clusters are in radio mode in the modern universe, this indicates that the AGN feedback efficiency in groups must be significantly smaller than in clusters in the radio mode regime.

Fig. 11. Soft luminosity vs. original AGN energy injection and corrected efficiency model. Left panel: Alignment of soft band luminosity with the central AGN feedback. The golden and salmon bars correspond to the ICM L500 luminosity in the [0.1−2.4] keV (soft) band from the Bahar et al. (2022) and Lovisari et al. (2020) samples, respectively, whereas the black, magenta, and brown bars correspond to the cavity powers from Rafferty et al. (2006), Russell et al. (2013), and the kinetic luminosity sample from Pasini et al. (2022), respectively. Observational error bars are at the 1σ level. The solid orange-red and grey lines correspond to the luminosity and original AGN feedback model obtained from the simulation, whereas the blue dashed line corresponds to the ‘corrected’ AGN feedback model obtained from the simulation using Eq. (7), and the blue dotted lines are the corresponding 2σ level intervals. Right panel: Total black hole efficiency in radio mode following the new model described by Eq. (8) based on the original accretion rates from the simulation shown in the left panel of Fig. 10. Cool-core clusters are shown in blue, and hot-core clusters are shown in orange. Solid lines indicate moving medians, and dashed lines represent 16% and 84% percentiles (1σ) from the Magneticum simulation (Box2b/hr). For comparison, the total black hole efficiencies in radio and quasar modes are shown with a dashed and dotted line, respectively.

7.2. A new model for AGN feedback efficiency in radio mode

Given that the accretion rates produced by the simulation seem to be well aligned with the observations as shown in the left panel of Fig. 10, we can check if the total feedback efficiency used in radio mode, especially for groups, can be adjusted. For this, we can interpret the observational relation between Bondi power defined as $P_{\mathrm{Bondi}}=0.1c^2{\dot{M}}_{\mathrm{B}}$ and cavity power (P_cavity) presented by Fujita et al. (2014) as implied total feedback efficiencies in radio mode.

Consistent with the computation of the cavity powers from before and also motivated by Fujita et al. (2016), we introduced two adjustments to the P_cavity − P_Bondi relation originally presented by Fujita et al. (2014): to assume that the jet cavities are filled with relativistic cosmic rays, therefore γ_c = 4/3 and E_cav = 4V_cp_s, and to apply a x2 factor to account for the shock energy as estimated by Rafferty et al. (2006) and Guo & Mathews (2010). We now can express the observed cavity power directly as a function of the inferred Bondi accretion rate Ṁ_B where the final parameters are $P_0=13.{60}_{-35.74}^{+55.56}$ and ${\beta }_P=1.{14}_{-0.22}^{+0.09}$:

$$\begin{array}{c}\hfill P_{\mathrm{cavity}}=\left[P_0{\left(\frac{{\dot{M}}_{\mathrm{B}}}{M_{\odot }/\mathrm{yr}}\right)}^{{\beta }_P}\right]{10}^{44}\mathrm{erg}/\mathrm{s}.\end{array}$$(7)

We further divided both sides by Ṁ_Bc² to obtain the implied total efficiency where ${ϵ}_0=0.{024}_{-0.06}^{+0.10}$ and ${\beta }_P-1=0.{14}_{-0.22}^{+0.09}$, thus indicating a weak dependency of the total feedback efficiency with the accretion rate:

$$\begin{array}{c}\hfill {ϵ}_t=\frac{P_{\mathrm{cavity}}}{{\dot{M}}_{\mathrm{B}}{c}^2}={ϵ}_0{\left(\frac{{\dot{M}}_{\mathrm{B}}}{M_{\odot }/\mathrm{yr}}\right)}^{({\beta }_P-1)}.\end{array}$$(8)

We recomputed the feedback power from the simulation using the above total efficiency depending on the BH accretion rate instead of a fixed value of 0.12 as currently used by the simulation in the radio mode, and we obtain that the injected feedback energy in the simulations matches the observed cavity power over the full range of masses, from massive clusters down to groups, as illustrated by the dashed line in the left panel of Fig. 11. Therefore, the reduced feedback efficiency at lower masses can potentially alleviate the excess of feedback at group scales currently seen in the simulation.

Note that the reduced values of the total feedback efficiency obtained by this scaling are at a similar level as the mechanical efficiencies of the best cold accretion model presented by Gaspari et al. (2013), who found values in the order of 5 ⋅ 10⁻³ − 10⁻² for galaxy clusters. The right panel of Fig. 11 shows a comparison of the total black hole efficiencies in radio mode following the new model described by Eq. (8) based on the original accretion rates from the simulation shown in the left panel of Fig. 10.

7.3. Limitations of AGN feedback models in simulations

In the Magneticum simulations, the AGN feedback is always injected purely as thermal energy within the surrounding gas. Following the original prescription from Springel (2005), this region is defined by requiring a fixed number of resolution elements, similar to the smooth particle hydrodynamics (SPH) treatment. The energy inside this so-called ‘black hole sphere of influence’ (H) is distributed to the neighbouring gas particles by the same kernel weighting scheme as used for the hydrodynamic treatment.

The mixing of this energy received by the particles with the rest of the ICM is mainly promoted by turbulent motions within the cluster core, as well as through thermal conductivity, which is implemented in an isotropic manner at 1/20 of the Spitzer level. However, buoyancy of the heated gas in the centre is largely suppressed by the central, spherical injection of the AGN feedback energy as well as by the smoothed gravitational potential under the gravitational resolution limit of ∼15 kpc.

To verify if the implemented AGN feedback model of the simulation somehow effectively mimics the situation in galaxy clusters and groups, we can compare this ‘black hole sphere of influence’ with the region where the energy of the radio mode feedback by jets is thermalised in observed galaxy clusters and groups. For this, we obtained the cavity maximum radius r_{cavity, max} = R_{cavity, centre} + 0.5(a + b) from the samples of Rafferty et al. (2006) and Eckert et al. (2021), where a is the projected semi-major axis, b is the projected semi-minor axis, and R_{cavity, centre} is the distance to the cavity centre, and we estimated the corresponding error to be 0.5(b − a).

In Fig. 12 we show r_{cavity, max} against the cavity power (left panel) as well as against M_500c (right panel). For comparison, we show the ‘black hole sphere of influence’ of the central BHs from the simulation against their feedback power. While the simulation shows a completely different scaling in the size vs. power relation at the resolution of the Magneticum Box2b/hr simulation, the two distributions have a significant overlap and similar mean values on both axes. Therefore, the current implementation of the AGN feedback model gives a reasonable, effective description for the resolution of the simulation, despite the different scaling. The same is visible in the size vs. mass relation. Although the scatter and trend in the observations differ largely from the model in the simulations, at massive cluster scale the distributions largely overlap, so that the implemented AGN feedback model gives a reasonable, effective description.

Fig. 12. Size of the sphere of influence of the simulated AGNs in the centres of groups and clusters. Cool-core clusters are shown in blue, and hot-core clusters are shown in orange. Observational data from Rafferty et al. (2006) and Eckert et al. (2021) is shown with error bars at the 1σ level. Left panel: Size vs power. Shown is the maximum size or distance of cavities in observed groups and clusters versus cavity power, whereas for the simulation we show the ‘black hole sphere of influence’ versus the AGN energy injection. The dark blue solid line shows the power law fit described by Equation (9), and the dashed lines represent ±1σ confidence intervals. Right panel: Same as the left panel, but plotting versus M500c in the X-axis.

However, the difference increases with increasing feedback power, which is produced by stronger accretion rates in higher gas density environments that typically harbour star formation; therefore, it can be related to the problems in suppressing star formation previously mentioned. Also, towards the group scale, the model in the simulation clearly injects the AGN feedback within a larger region.

Also, the ‘black hole sphere of influence’ scales down with increasing resolution; therefore, we expect this difference to increase as the resolution increases. In addition, alternative AGN feedback implementations used in other cosmological hydrodynamical simulations, even when based on injecting kinetic feedback, do not guarantee that the effective thermalisation of the feedback energy coincides with these observations.

To improve future AGN feedback models, we can quantify the observational relation of cavity reach vs. power. A fit to the observational data from Rafferty et al. (2006) with an orthogonal distance regression (ODR), accounting for the errors in both variables⁴, results in β_r = 0.37 ± 0.05 and r₀ = 4.2 ± 1.3 for a power law relation in the form:

$$\begin{array}{c}\hfill r_{\mathrm{cavity},\max }=\left[r_0{\left(\frac{P_{\mathrm{cavity}}}{{10}^{42}\mathrm{erg}/\mathrm{s}}\right)}^{{\beta }_r}\right]\mathrm{kpc}.\end{array}$$(9)

Using this relation, we can correct the ‘black hole sphere of influence’ used in the simulation to overlap with the observations of cavity reach, as shown by the blue solid and dashed lines in the right panel of Fig. 12, representing the median and ±2σ intervals, respectively. Therefore, this relation can be used to improve future, effective models for AGN feedback in simulations and make them more resolution independent, although the accretion rate in itself still depends on the resolution and requires adjusting the boost factor in Eq. (5) accordingly. Additionally, to fully reproduce the AGN feedback signatures, it is necessary to track the BH spin and inject the feedback along its axis in a bipolar pattern, as shown by the simulations of Sala et al. (2024).

7.4. Interpretation

While the size of the ‘black hole sphere of influence’ is purely given by a technical choice in the implementation and therefore has no physical meaning, the original choice of the feedback efficiency in radio mode was motivated by observations of massive galaxy clusters. Given the tight relation between halo mass and cavity power and the weak dependency between BH accretion rate and total feedback efficiencies presented in Eq. (8), which results in significantly smaller total efficiencies at the scale of groups, it is worthwhile to see if there could be a theoretical argument for such scaling.

The most established theoretical model to produce jets (mechanical feedback) is the Blandford & Znajek (1977) (BZ) model, in which a spinning BH drags the magnetic field lines, producing an electric field by induction that accelerates charged particles. Therefore, in this scenario, the energy extraction efficiency is proportional to the BH spin parameter, which means that slowly spinning BHs should be less efficient at generating jets.

However, according to observational constraints by Reynolds (2021), there is a decreasing tendency in the spin parameter towards higher BH masses typically hosted in the most massive systems. This trend is confirmed by simulations of spin evolution (Sala et al. 2024) and should be expected according to the isotropic principle, since mergers of BHs with spins orientated in random directions should produce higher-mass BHs with an increasingly lower spin parameter. As a consequence, and according to the BZ model, higher-mass black holes should be less efficient at generating jets. This scenario is in tension with observations since stronger jets and cavity powers are actually observed for the most massive systems, as shown by the right panel of Fig. 10.

On the other hand, an alternative route to generate mechanical feedback even for non-spinning black holes can be collisions between particles orbiting the innermost stable orbits (R_ISCO) and free-falling particles. In this sense, theoretical advancements by Frolov (2012) demonstrated that even with a weak uniform magnetic field, the radius of the innermost circular orbit of a non-rotating BH approaches the event horizon thanks to the additional support of a repulsive Lorentz force, which can potentially boost the mass-energy conversion efficiency for the accreted material. The effect is characterised by the ratio b between the Lorentz force and the surface gravity κ at the event horizon (Baker & Frolov 2023) as shown by Equation (10), where q and m are the charge and mass of the particle orbiting R_ISCO, B is the magnetic field at the event horizon, M_BH is the mass of the BH, G is the gravitational constant, and c is the speed of light:

$$\begin{array}{c}\hfill b=\frac{1}{4}\frac{\mathit{qBc}}{\kappa }=\frac{qBG{M}_{\mathrm{BH}}}{m{c}^3}\sim {10}^{15}\left(\frac{q}{e}\right)\left(\frac{m_e}{m}\right)\left(\frac{B}{{10}^4\mathrm{G}}\right)\left(\frac{M_{\mathrm{BH}}}{{10}^9{M}_{\odot }}\right).\end{array}$$(10)

This ratio can indeed become very large even for a weak magnetic field in the order of [1 − 30] G, as measured by the Event Horizon Telescope for M87 (EHT Collaboration 2021). Moreover, Frolov (2012) showed that the energy of a collision between a particle at R_ISCO and a free-falling particle is directly proportional to E_collision ∼ mc²b^1/4.

To quantify the magnetic field, we have to consider that the BHs hosted in the centre of BCGs generally have a low Eddington ratio and luminosity, characterised by the hot advection-dominated accretion flow model (ADAF, Narayan & Yi 1995), which is usually associated with jets, according to observations. In the ADAF model, the magnetic field strength at the Schwarzschild radius is proportional to the accretion rate Ṁ_BH, as shown by Equation (11) (Yuan & Narayan 2014), where α is the viscosity parameter, $\beta =\frac{P_{\mathrm{gas}}}{P_{\mathrm{mag}}}$ the magnetisation parameter, and M_BH the mass of the BH:

$$\begin{array}{c}\hfill B=0.44{(1+\beta )}^{-\frac{1}{2}}{\alpha }^{-\frac{1}{2}}{\left(\frac{M_{\mathrm{BH}}}{{10}^9{M}_{\odot }}\right)}^{-1}{\left(\frac{{\dot{M}}_{\mathrm{BH}}}{M_{\odot }/\mathrm{yr}}\right)}^{\frac{1}{2}}{10}^4\mathrm{G}.\end{array}$$(11)

Then, if we combine Eqs. (10) and (11) and consider the energy extraction efficiency (ϵ) of the Frolov (2012) process as the ratio between the collision energy (E_collision) and the rest mass energy of the free-falling particle (mc²), we obtain that the energy extraction efficiency depends on ${\dot{M}}_{\mathrm{BH}}^{1/8}$ as shown by Equation (12), and it is quite close to the 0.14 exponent obtained from the observational relation shown in Equation (8):

$$\begin{array}{c}\hfill ϵ\sim \frac{E_{\mathrm{collision}}}{m{c}^2}\sim b^{\frac{1}{4}}\sim {\left({\mathit{BM}}_{\mathrm{BH}}\right)}^{\frac{1}{4}}\sim {\dot{M}}_{\mathrm{BH}}^{1/8}.\end{array}$$(12)

This process can provide a channel independent of BH spin to produce jets (mechanical feedback), which can be a key ingredient to explain increasing jet and cavity powers towards the most massive systems and decreasing towards low mass systems to reduce the overall mechanical feedback efficiency at the scale of groups and address the overheating problems described in Sect. 4.4 and also reported by Gaspari et al. (2013) and suggested by the entropy excess at the cores of simulated galaxy groups shown by Bahar et al. (2024).

8. Comparison with the results of other simulations

Comparing the results of this work with other cosmological simulations, we confirm the decreasing trend in cool-core fractions towards higher mass clusters initially reported by Chen et al. (2007) and also confirmed by the simulations of Burns et al. (2008) and Planelles & Quilis (2009). However, Lehle et al. (2024) report that for the TNG-Cluster simulation, the non-cool core fraction increases with mass only when using the central cooling time and entropy as indicators, whereas it decreases when using the central density and cuspiness. Here we argue that, as explained in Sects. 3 and 4.1, the central properties are subject to resolution issues in both observations and simulations, and also, for the case of simulations, are affected by the direct AGN feedback, thus subject to the limitations in the sub-grid model affecting the distribution of the AGN feedback as described in Section 7.3.

Regarding the problems with overheating and excess of entropy in the central regions at the scale of galaxy groups initially reported by Fabjan et al. (2010), and confirmed by Bahar et al. (2024) and this work, Barnes et al. (2017) report a similar problem in the temperature profiles already at the scale of galaxy clusters (median M_500c = 2.1 ⋅ 10¹⁴ M_⊙) for the Cluster-EAGLE simulations (C-Eagle), also confirmed by Altamura et al. (2023) in the entropy profiles at the scales of groups (M_500c = 8.8 ⋅ 10¹² M_⊙) and clusters (M_500c = 2.9 ⋅ 10¹⁴ M_⊙). Similarly, we also see some traces of overheating in our hot-core clusters at the mid-mass range (10¹⁴ M_⊙ < M_500c ≤ 4.88 ⋅ 10¹⁴ M_⊙), although it is possible that this problem is partially mitigated for the Magneticum simulations at the scale of massive galaxy clusters due to the implementation of physical conductivity at 1/20 of the Spitzer value, which nevertheless we plan to review in future works. At the galaxy group regime, our conclusions are aligned with results obtained by the simulations of Gaspari et al. (2013) in that lower mass halos require reduced AGN feedback efficiencies to prevent overheating.

Regarding the distribution of AGN feedback, Le Brun et al. (2014) report that the best feedback model for the OverWhelmingly Large Simulations (OWLS) requires injecting the feedback only in one gas particle per feedback loop and only if the injected energy is above 10⁸ K (8.6 keV) to prevent that the feedback energy gets quickly radiated away if distributed among all neighbouring particles. This contrasts with the approach presented in this work, to distribute the feedback through the thermal channel directly in the region where the signatures of AGN feedback (radio emission, cavities) are observed. Here we point out that a key aspect to prevent catastrophic cooling is to inject the energy into the ICM particles that have not yet been cooled, which can be a more efficient process given the longer cooling times of the ICM. Moreover, although we agree that the AGN feedback in real physical scenarios is generated in the surroundings of the AGN, it is expelled at relativistic speeds and thermalises at a much larger distance corresponding to the observed excavated cavities that overlap with radio emission.

Finally, our conclusion that merge activity does not directly turn cool-cores into hot-cores since it first requires thermalisation, coincides with the conclusions of Poole et al. (2008) based on a suite of idealised mergers with different mass ratios and impact parameters and Hahn et al. (2017) based on the RHAPSODY-G simulations, who report that only direct collisions (low angular momentum mergers) can effectively turn cool-cores into hot-cores. Additionally, we also agree with Rasia et al. (2015) and Hahn et al. (2017) in that thermal conduction is required to effectively redistribute the energy injected by mergers and AGN feedback across the ICM in the core region.

9. Conclusions

We conduct a study to compare the halo mass dependency of cool-core fractions for a large sample of groups and clusters from observational data, combining the samples from Lovisari et al. (2020) and Bahar et al. (2022) at low redshift (z < 0.3) with a large-scale, cosmological hydrodynamical simulation Box2b/hr from the Magneticum set. This covers a mass range of 0.3 × 10¹⁴ M_⊙ < M_500c < 10¹⁵ M_⊙ with a sample size of 201 groups and clusters for the observational dataset and 13 683 for the simulation dataset. For the set of 2027 simulated clusters above a mass of M_500c = 10¹⁴ M_⊙, we compare in detail the radial profiles of density, temperature, and entropy with the results obtained from X-ray observations. We investigate the possible physical mechanism that drives the observed trends of the cool-core fractions in simulations and observations across this mass range. We also compare the assumptions and predictions of the AGN feedback model implemented within the simulations with the observed properties of the central AGN and their inflated radio bubbles in clusters and groups. Our main findings can be summarised as follows:

• In general, the simulation is able to reproduce the overall cool-core fractions at the scale of galaxy clusters (M_500c > 1 ⋅ 10¹⁴ M_⊙) and the mid-to-strong cool-core fractions at the scale of massive clusters (M_500c > 2.7 ⋅ 10¹⁴ M_⊙). When splitting the simulated clusters in cool-core and hot-core systems, the radial profiles of the two classes reproduce the different shapes observed for the two different classes, where hot-core systems show a more cored entropy profile with larger values compared to cool-core systems and hot-core systems show a more isothermal temperature profile in the centre compared to the drop of the temperature of cool-core systems towards the centre. However, some traces of overheating and inefficient energy transport are visible in the temperature profiles of the hot-core clusters at the mid-mass range (10¹⁴ M_⊙ < M_500c ≤ 4.88 ⋅ 10¹⁴ M_⊙), calling for a review of the physical conductivity model currently set at 1/20 of the Spitzer value.

• The cool-core fractions clearly decrease towards high-mass galaxy clusters; this trend is observed in both simulations and observations and confirms the results reported by Chen et al. (2007), Burns et al. (2008), and Planelles & Quilis (2009). However, the relative contribution of AGN feedback to this process also decreases at the high mass end, indicating that an additional factor is required to suppress cool-cores. While within the core, the amount of internal plus kinetic energy compared to the potential energy for the total system, including dark matter, stellar, and gas components, is increasing with mass, there is no visible difference for cool-core and hot-core systems. However, when only considering the gas component, we find a clear separation of cool-core and hot-core systems, with hot-core systems having a larger internal plus kinetic energy compared to cool-core systems as well as higher internal energy fractions, both increasing with mass. This indicates that the thermalisation of the kinetic energy induced by mergers is more efficient for hot-core systems and also increases in the high-mass regime. This factor, combined with the implementation of physical conductivity in the Magneticum simulations, which is also more efficient as the halo mass increases due to its strong temperature dependency, can effectively reduce the cool-core fractions at the scale of massive clusters.

• On the other hand, the cool-core fractions decrease much more sharply in the simulation than in observations for the low-mass galaxy groups. This divergence is equivalent to the excess of entropy in the cores of simulated galaxy groups shown by Bahar et al. (2024). This can be associated with the relatively high impact that the AGN feedback injected in radio mode has in comparison with the luminosity at the scale of galaxy groups. While the simulation produces Bondi accretion rates compatible with the refined estimations by Fujita et al. (2014) for the AGNs hosted in the centres of BCGs, the usage of a relatively large AGN feedback efficiency as found in inferred in massive clusters might be questionable in galaxy groups and might be much lower and in the range as inferred from the cold accretion model presented in Gaspari et al. (2013).

• The AGN feedback efficiencies in the simulations at group scales could be aligned when using the relation between Bondi power and cavity power presented by Fujita et al. (2014). This relation is equivalent to a weak dependency of the total feedback efficiency in radio mode on the accretion rate, with an exponent close to 1/8. Interestingly, the same dependency and exponent are obtained considering an AGN energy injection mechanism based on collisions between free-falling particles and charged particles orbiting the innermost stable orbits in the presence of a weak magnetic field, following the theoretical model of Frolov (2012). This process also provides an AGN energy injection mechanism independent of spin, which decreases towards the supermassive BH hosted in the centres of galaxy clusters, as suggested by the observations of Reynolds (2021) and confirmed by the simulations of Sala et al. (2024).

• Finally, we point out that the current implementation of AGN feedback models in simulations often scales the size of the region where the AGN feedback energy is injected inversely with the local density, which is in contrast to the scaling of the observed size and reach of cavities. This inverted scaling can reduce the effectiveness of the AGN feedback to suppress star formation in high-density environments. Although for the current simulation there is still a significant overlap, which makes such effective models at least work partially, future simulations with increased resolution need to adapt the sub-grid models accordingly. Taking the observational data from Rafferty et al. (2006), we obtain a power law relation of the AGN energy injection and the expected reach that can be used to bring the AGN feedback distribution schema from simulations into agreement with observations.

Acknowledgments

This work was supported by the COMPLEX project from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program grant agreement ERC-2019-AdG 882679. VB acknowledges support by the DFG project nr. 415510302. The calculations for the hydrodynamical simulations were carried out at the Leibniz Supercomputer Center (LRZ) under the project pr83li. We thank the developers of the following software packages that we have used for the data analysis: Python 3.6.15 (Van Rossum & Drake 1995), numpy 1.19.5 (Harris et al. 2020), g3read, g3matcha (Ragagnin et al. 2017), astropy 4.1 (Astropy Collaboration 2013, 2022), scipy 1.5.4 (Virtanen et al. 2020), AtomDB v3.0.9 (Foster et al. 2018), PyAtomDB v0.10.10 (Foster & Heuer 2020), statsmodels 0.12.2 (Seabold & Perktold 2010), asymmetric uncertainty 0.2.1 (Gobat 2022), and matplotlib 3.3.4 (Hunter 2007; Caswell et al. 2023) This research has made use of the M2C Galaxy Cluster Database, constructed as part of the ERC project M2C (The Most Massive Clusters across cosmic time, ERC-Adv grant No. 340519). We thank Adam Foster for his support and insights regarding the APEC model, AtomDB and PyAtomDB; Thomas Pasini for kindly providing the LOFAR kinetic luminosity data; and Lucas Sala for the many useful discussions and feedback regarding accretion flows and AGN feedback.

References

Appendix A: Dependency of cool core fractions on the temperature ratio threshold

In this appendix we show how the cool-core fractions depend on the temperature ratio threshold, set to unity in the main part of the analysis. Although we believe that the choice of T_ratio, 500 < 1, is the most straightforward one to asses the overall cool-core cluster population (including weak, mid and strong cool-core clusters), with a clear physical meaning that the core regions are on average cooler than the cluster average, still it can provide valuable insights to see how the cool-core fractions change depending on this threshold.

To this end we consider two additional thresholds: T_ratio, 500 < 0.98 and T_ratio, 500 < 0.95, which can also be interpreted as the fraction of mid-to-strong and strictly strong cool cores, respectively. Notice that it is not possible to reliably study the population of even stronger cool-cores (e.g. T_ratio, 500 < 0.90) given that we are already separating the data per mass bin, which reduces the amount of objects at the scales of massive clusters for both the simulated and observed samples and limits the possibilities to assess in a statistically robust manner the population of the most rare objects (really strong cool cores at scales of massive clusters).

Fig. A.1 shows the results, where we can see that although the characteristic cool-core fractions curve seen in the simulation data is preserved, the cool-core fractions of mid-to-strong and strictly strong cool cores produced by the simulation are significantly lower than that of the observational data at low and mid mass ranges. On the other hand, the simulation still produces cool core fractions compatible with observations at the high mass range bins M_500c > 2.65 ⋅ 10¹⁴ M_⊙ for all cases (weak, mid, and strong cool-core clusters).

Fig. A.1. Cool core fractions (same as Fig. 2) but using a different temperature ratio threshold. The vertical dashed line corresponds to M500c = 0.7 ⋅ 1013 M⊙, above which the eFEDS survey is expected to be complete for redshifts below z < 0.3 (Comparat et al. 2020). Left panel: Mid-to-strong cool-core fractions determined by the number of clusters for which the temperature ratio between the total temperature, including the core region, and the core-excised temperature is less than 0.98 (Tx500/Tx500, cex < 0.98). Right panel: Strictly strong cool-core fractions determined by the number of clusters for which the temperature ratio between the total temperature, including the core region, and the core-excised temperature is less than 0.95 (Tx500/Tx500, cex < 0.95).

This result extends the conclusions presented in Sec. 4.4, in that the overheating problems caused by the current AGN feedback which suppresses the overall cool-core cluster population at the low mass range, are also visible at the mid mass range when considering the mid and strong cool-core cluster population. Still, both problems should be addressed by the new AGN total feedback efficiency model presented in Sec. 7.2, which produces an effective total AGN efficiency decreasing towards the mid and low mass range.

Appendix B: Temperature profiles following the ACCEPT centering procedure

In this appendix we want to explore how the temperature profiles would look like when following the same centering procedure as the ACCEPT sample described by Cavagnolo et al. (2008). That is, choosing the X-ray peak as the center, unless it is separated by more than 70 kpc from the X-ray centroid, in which case the latter is used as a center.

For this purpose, we used the average rest-frame band XMM-eFEDS described in Sec. 3 to obtain the X-ray peak and centroid based on emissivity weights. Then we followed the ACCEPT centering procedure with the simulated data, using a threshold of $70·\sqrt{3/2}$ kpc, which is the 3-dimensional version of 70 kpc when accounting for projection effects, to switch from X-ray peak to X-ray centroid as a center for the radial profiles.

The results are shown in Fig. B.1, where we see that the hot-core clusters show a stronger isothermal core (with fewer fluctuations) in comparison with the default temperature profiles shown in Fig. 6, which use the most gravitationally bound particle, equivalent to the deepest point of the gravitational potential as the center. The explanation for this behavior is relatively simple:

• The deepest point of the gravitational potential harbors cool, dense gas that falls by precipitation and also hosts the central AGN, which heats the immediate surroundings; hence the temperature profiles swing from cool to hot as we go through the envelope of cool, dense gas into the immediate surroundings of the central AGN.

• On the other hand, the X-ray centroid is distant from the X-ray peak for disturbed clusters (e.g. two dominant galaxies merging), in which case the X-ray centroid is located in between the merging structures, in a rather shallow region of the gravitational potential that does not support strong temperature gradients, harborage of cold dense gas, or stable conditions to host an AGN. All these factors contribute to an isothermal structure at the X-ray centroid of disturbed clusters.

Fig. B.1. Projected X-ray temperature profiles (same as Fig. 6) but using the same centering procedure as the ACCEPT sample described in (Cavagnolo et al. 2008). Columns are sorted by mass range, as indicated on top of each panel (increasing mass range from left to right). The upper row panels correspond to cool-core clusters, and the lower row panels to hot-core clusters.

All Tables

All Figures

	Fig. 1. X-ray SB images of simulated cluster examples selected within the 3 highest mass bins. The upper row panels correspond to cool-core clusters, and the lower row panels to hot-core clusters.
In the text

Fig. 2. Observed and simulated cool-core fractions. The gold bars correspond to the simulation, for which the temperature was obtained with emissivity weights in the XMM-eFEDS band. The grey bars for the low-mid mass range correspond to the eROSITA field equatorial deep survey eFEDS data from Bahar et al. (2022) and Chiu et al. (2022). (*Also contains 1 mid-mass cluster from the Lovisari et al. 2020 sample). The black bars for the high mass range correspond to the Planck SZ-selected sample observed with XMM-Newton from Lovisari et al. (2020). (**Also contains 10 high mass clusters from the Bahar et al. 2022 sample). The vertical dashed line corresponds to M500c = 0.7 ⋅ 1013 M⊙, above which the eFEDS survey is expected to be complete for redshifts below z < 0.3 (Comparat et al. 2020).

In the text

Fig. 3. Driving factors behind the cool-core fraction characteristic curve. The left and central panels correspond to results from the Magneticum simulation (Box2b), where cool-core clusters are shown in blue and hot-core clusters in orange. Solid lines indicate moving medians, and dashed lines indicate 16% and 84% percentiles (1σ). Left panel: Ratio between the energy injection from the central AGN feedback and the bolometric luminosity in the [0.01−100] keV band for gas particles inside the core region. Central panel: Number of mergers undergone by the central black hole of the brightest cluster galaxy (BCG). Right panel: Sketch to illustrate the concept of how the driving factors combine to produce the characteristic shape of the cool-core fraction curve.

In the text

Fig. 4. Dynamical analysis. Cool-core clusters are shown in blue, and hot-core clusters in orange. Solid lines indicate moving medians, and dashed lines indicate 16% and 84% percentiles (2σ) from the Magneticum simulation (Box2b/hr). Left panel: Freedom ratio for all particles inside the core region (Fc, total). Central panel: Freedom ratio for gas particles inside the core region (Fc, gas). Right panel: Kinetic energy fraction fc, gas, kinetic for gas particles inside the core region. The black bar indicates the first and only measurement of kinetic energy fraction in galaxy clusters, of 10% measured for the Perseus cluster by Hitomi (Hitomi Collaboration 2018).

In the text

Fig. 5. Thermal conductivity. Cool-core clusters are shown in blue, and hot-core clusters in orange. Solid lines indicate moving medians, and dashed lines indicate 16% and 84% percentiles (2σ) from the Magneticum simulation (Box2b/hr). Left panel: Temperature of the core region obtained with emissivity weights in the XMM-eFEDS band. Right panel: Effective Spitzer conductivity coefficient for the core temperature, normalised to the value for a system at 1 keV.

In the text

Fig. 6. Projected X-ray temperature profiles. Corresponding to the [0.7−7.0] keV band and normalised by the mean X-ray temperature in the range [0.1 − 0.2] R500c. The error bars correspond to the median Chandra ACCEPT sample profiles and ±1σ intervals using the original masses from the M2C Galaxy Cluster Database (black) and a +20% hydrostatic mass bias correction (orange). The blue line and shaded areas correspond to the Magneticum simulation (Box2b) median profiles and ±1σ intervals. Columns are sorted by mass range, as indicated on top of each panel (increasing mass range from left to right). The upper row panels correspond to cool-core clusters, and the lower row panels to hot-core clusters.

In the text

Fig. 7. Electron number density profiles. The error bars correspond to the median Chandra ACCEPT sample profiles and ±1σ intervals using the original masses from the M2C Galaxy Cluster Database (black) and a +20% hydrostatic mass bias correction (orange). The blue line and shaded areas correspond to the Magneticum simulation (Box2b) median profiles and ±1σ intervals. Columns are sorted by mass range, as indicated on top of each panel (increasing mass range from left to right). The upper row panels correspond to cool-core clusters, and the lower row panels to hot-core clusters.

In the text

Fig. 8. Electron entropy profiles. The error bars correspond to the median Chandra ACCEPT sample profiles and ±1σ intervals using the original masses from the M2C Galaxy Cluster Database (black) and a +20% hydrostatic mass bias correction (orange). The blue line and shaded areas correspond to the Magneticum simulation (Box2b) median profiles and ±1σ intervals. Columns are sorted by mass range, as indicated on top of each panel (increasing mass range from left to right). The upper row panels correspond to cool-core clusters, and the lower row panels to hot-core clusters.

In the text

Fig. 9. Stellar, gas, and baryonic mass fractions. Cool-core clusters are shown in blue, and hot-core clusters are shown in orange. Solid lines indicate moving medians, and dashed lines represent 16% and 84% percentiles (2σ) from the Magneticum simulation (Box2b). Observational data is shown in black, with error bars at the 1σ level, based on XMM-Newton, Chandra, and SDSS from Laganá et al. (2013). Left panel: Gas mass fraction inside R2500c. Central panel: Stellar mass fraction inside R2500c. Right panel: Baryonic mass fraction inside R2500c.

In the text

Fig. 10. Accretion rates and energy injection of the AGN hosted in the centres on BCGs. Cool-core clusters are shown in blue, and hot-core clusters are shown in orange. Solid lines indicate moving medians, and dashed lines represent 16% and 84% percentiles (1σ) from the Magneticum simulation (Box2b/hr). Observational data is shown with error bars at the 1σ level. Left panel: Comparison of the accretion rates from the AGN hosted in the centres on BCGs, with the estimations from Fujita et al. (2014) and Russell et al. (2013). Right panel: Comparison of black hole energy injection with the observational estimations from Rafferty et al. (2006), Russell et al. (2013), Eckert et al. (2021) and Pasini et al. (2022).

In the text

Fig. 11. Soft luminosity vs. original AGN energy injection and corrected efficiency model. Left panel: Alignment of soft band luminosity with the central AGN feedback. The golden and salmon bars correspond to the ICM L500 luminosity in the [0.1−2.4] keV (soft) band from the Bahar et al. (2022) and Lovisari et al. (2020) samples, respectively, whereas the black, magenta, and brown bars correspond to the cavity powers from Rafferty et al. (2006), Russell et al. (2013), and the kinetic luminosity sample from Pasini et al. (2022), respectively. Observational error bars are at the 1σ level. The solid orange-red and grey lines correspond to the luminosity and original AGN feedback model obtained from the simulation, whereas the blue dashed line corresponds to the ‘corrected’ AGN feedback model obtained from the simulation using Eq. (7), and the blue dotted lines are the corresponding 2σ level intervals. Right panel: Total black hole efficiency in radio mode following the new model described by Eq. (8) based on the original accretion rates from the simulation shown in the left panel of Fig. 10. Cool-core clusters are shown in blue, and hot-core clusters are shown in orange. Solid lines indicate moving medians, and dashed lines represent 16% and 84% percentiles (1σ) from the Magneticum simulation (Box2b/hr). For comparison, the total black hole efficiencies in radio and quasar modes are shown with a dashed and dotted line, respectively.

In the text

Fig. 12. Size of the sphere of influence of the simulated AGNs in the centres of groups and clusters. Cool-core clusters are shown in blue, and hot-core clusters are shown in orange. Observational data from Rafferty et al. (2006) and Eckert et al. (2021) is shown with error bars at the 1σ level. Left panel: Size vs power. Shown is the maximum size or distance of cavities in observed groups and clusters versus cavity power, whereas for the simulation we show the ‘black hole sphere of influence’ versus the AGN energy injection. The dark blue solid line shows the power law fit described by Equation (9), and the dashed lines represent ±1σ confidence intervals. Right panel: Same as the left panel, but plotting versus M500c in the X-axis.

In the text

Fig. A.1. Cool core fractions (same as Fig. 2) but using a different temperature ratio threshold. The vertical dashed line corresponds to M500c = 0.7 ⋅ 1013 M⊙, above which the eFEDS survey is expected to be complete for redshifts below z < 0.3 (Comparat et al. 2020). Left panel: Mid-to-strong cool-core fractions determined by the number of clusters for which the temperature ratio between the total temperature, including the core region, and the core-excised temperature is less than 0.98 (Tx500/Tx500, cex < 0.98). Right panel: Strictly strong cool-core fractions determined by the number of clusters for which the temperature ratio between the total temperature, including the core region, and the core-excised temperature is less than 0.95 (Tx500/Tx500, cex < 0.95).

In the text

	Fig. B.1. Projected X-ray temperature profiles (same as Fig. 6) but using the same centering procedure as the ACCEPT sample described in (Cavagnolo et al. 2008). Columns are sorted by mass range, as indicated on top of each panel (increasing mass range from left to right). The upper row panels correspond to cool-core clusters, and the lower row panels to hot-core clusters.
In the text