The FLAMINGO Project: Galaxy clusters in comparison to X-ray observations

Galaxy clusters are important probes for both cosmology and galaxy formation physics. We test the cosmological, hydrodynamical FLAMINGO simulations by comparing to observations of the gaseous properties of clusters measured from X-ray observations. FLAMINGO contains unprecedented numbers of massive galaxy groups ($>10^6$) and clusters ($>10^5$) and includes variations in both cosmology and galaxy formation physics. We predict the evolution of cluster scaling relations as well as radial profiles of the temperature, density, pressure, entropy, and metallicity for different masses and redshifts. We show that the differences between volume-, and X-ray-weighting of particles in the simulations, and between cool-core non cool-core samples, are similar in size as the differences between simulations for which the stellar and AGN feedback has been calibrated to produce significantly different gas fractions. Compared to thermally-driven AGN feedback, kinetic jet feedback calibrated to produce the same gas fraction at $R_{\rm 500c}$ yields a hotter core with higher entropies and lower densities, which translates into a smaller fraction of cool-core clusters. Stronger feedback, calibrated to produce lower gas fractions and hence lower gas densities, results in higher temperatures, entropies, and metallicities, but lower pressures. The scaling relations and thermodynamic profiles show almost no evolution with respect to self-similar expectations, except for the metallicity decreasing with redshift. We find that the temperature, density, pressure, and entropy profiles of clusters in the fiducial FLAMINGO simulation are in excellent agreement with observations, while the metallicities in the core are too high.


INTRODUCTION
The largest collapsed structures in our Universe, galaxy clusters, are an excellent probe of some of the most violent events we can observe, such as mergers with other clusters and feedback from accretion onto the most massive black holes.Furthermore, their abundance and clustering are important probes for cosmology.Because of their large mass, they are, unfortunately, quite rare, and statistical samples are much smaller than for lower-mass objects like galaxies.
★ E-mail: braspenning@strw.leidenuniv.nl The number of such galaxies, called satellites, is a good measure of the mass of the cluster (e.g.Pereira et al. 2018).
Because clusters contain large amounts of hot gas, X-ray observations provide an alternative way to identify clusters.X-ray observations probe the intracluster medium (ICM), unlike longer wavelengths (such as the optical) which only identify galaxies.By imaging the ICM, X-ray observations can be used to construct surface brightness profiles (e.g.Neumann & Arnaud 1999), and the (deprojected) temperature, metallicity, and density profiles can be derived by fitting the spectrum (e.g.Vikhlinin 2006;Sun et al. 2009).These can then be combined to measure entropy and pressure profiles (e.g.Ponman et al. 2003;Arnaud et al. 2010).
The density and temperature profiles can be used to derive the hydrostatic mass of clusters through the equation of hydrostatic equilibrium (e.g.Ettori et al. 2013).However, such masses tend to be biased, as clusters are never truly in hydrostatic equilibrium, and there may be significant non-thermal pressure.This ratio between the measured hydrostatic mass and the true mass, is called the hydro-static bias.Comparing predicted hydrostatic masses with true masses in simulations leads to a hydrostatic bias of ≈ 0.8 − 0.9 (e.g.Le Brun et al. 2014;Biffi et al. 2016;Barnes et al. 2021;Gianfagna et al. 2021;Jennings & Davé 2023).Alternatively, observational comparisons between X-ray inferred masses and masses obtained from weak lensing can be made (e.g.Mahdavi et al. 2013;von der Linden et al. 2014;Henson et al. 2017), giving a slightly larger hydrostatic bias (≈ 0.73) (e.g.Muñoz-Echeverría et al. 2023;Hoekstra et al. 2015).
Even though clusters deviate from perfect hydrostatic equilibrium, they are still the largest collections of particles in the Universe showing behaviour which can, at least in part, be described by such simple physics.As a result, it is of particular interest to compare scaling relations of different physical quantities associated with clusters.Common scaling relations in the literature are the mass-X-ray-luminosity relation (e.g.Lovisari et al. 2015Lovisari et al. , 2020;;Gaspari et al. 2019), masstemperature relation (e.g.Pearson et al. 2017), and the temperature-X-ray-luminosity relation (e.g.Pratt et al. 2009;Migkas et al. 2020).
At the lower mass end, lower temperatures lead to more metal emission lines, boosting the integrated X-ray luminosity over the thermal bremsstrahlung expectation (e.g.Lovisari et al. 2021).However, due to non-gravitational processes such as AGN feedback pushing out gas from the centers of groups, the scaling relations do not show any such boost towards the low mass and low temperature end.The gas fraction within  2500c 1 is observed to be strongly mass dependent, whereas it is almost mass independent in the range  2500c −  500c , suggesting that gas is pushed out of the centers of groups, but does not escape the system (Sun et al. 2009).The mass dependence of the mass-X-ray-luminosity relation also depends on the X-ray band within which the flux is integrated, as some lines may or may not be captured.In particular, it has been shown that the 0.5 − 2.0 keV band is fairly insensitive to the metallicity, but the wider 0.1 − 2.4 keV band, as well as the bolometric luminosity, are sensitive to the metallicity and the assumed element abundances for lower temperatures (< 1 keV) (e.g.Lovisari et al. 2021).
Whereas scaling relations characterize clusters with single quantities such as mass, luminosity, or temperature, radial profiles give insight into the physical processes at play at different distances from the cluster centre.They are of particular interest when comparing simulations with observations, as agreement with observations would be a strong indication that the physical processes implemented in those simulations yield, on the scales represented by galaxy clusters, conditions that are realistic.Direct comparisons are, however, difficult.For example, observed density, temperature, and metallicity profiles are derived from fitting models to the X-ray spectrum of clusters within radial bins, whereas in simulations mass-or volumeweighted quantities are often used (e.g.Lehle et al. 2023;Li et al. 2023;Nelson et al. 2023;Pakmor et al. 2023;Towler et al. 2023), but there are already long standing efforts to create mock observations and measure cluster properties directly from those (e.g.Nagai et al. 2007;Rasia et al. 2012;McCarthy et al. 2017;Robson & Davé 2023).Because the X-ray luminosity, and hence the contribution to the total observed spectrum of a region, scales with the square of the gas density, denser regions contribute more to the quantities inferred from observations than they would in mass-or volume-weighted simulations.For the inferred temperature there is an empirical formulation, the spectroscopic-like temperature, to mimic the temperature derived 1  Δc is the radius within which the mean internal density is Δ times the critical density of the universe.The mass contained within  Δc is denoted  Δc .
Observational X-ray selection tends to be biased towards high surface brightness clusters, especially at higher redshifts (e.g.Eckert et al. 2011).High surface brightness clusters are typically coredominated, and classified as cool-core clusters since their temperature profile shows a central drop.The bias decreases but is still present for Sunyaev-Zel'dovich (SZ) selected samples (e.g.Lin et al. 2015;Rossetti et al. 2017).Furthermore, observational analyses have to deproject the X-ray maps to obtain 3D profiles, which adds significant uncertainty as it relies on the assumed sphericity of the cluster and limited thermal structure (e.g.Bartalucci et al. 2023).In particular, this approach tends to overestimate central temperatures (Lakhchaura et al. 2016).
Due to their rarity, even in simulations the number of galaxy clusters tends to be low, as enormous volumes are required to obtain a large population.Here we analyze the FLAMINGO cosmological hydrodynamical simulations (Schaye et al. 2023;Kugel et al. 2023), which include the largest full physics simulation run to  = 0 to date.FLAMINGO includes different numerical resolutions and volumes up to 2.8 Gpc on a side.The subgrid feedback in the fiducial model has been calibrated to reproduce the observed low-redshift galaxy stellar mass function and cluster gas fractions, while model variations produce different mass functions and/or gas fractions, use jet-like kinetic AGN feedback instead of thermally-driven AGN feedback, or assume different cosmologies.The unprecedented combination of volume and resolution offers a very large number of resolved clusters to compare with observational samples, e.g. more than two million haloes with mass  500c > 10 13 M ⊙ at  = 0. Feedback and cosmology variations allow the study of the relative importance of different effects, and elucidate the physics driving the evolution and observational appearance of these rare objects.
In this paper, we analyse the scaling relations and thermodynamic profiles of clusters in the FLAMINGO simulations and compare them with observations.All runs have the same resolution, with the exception of convergence tests in the Appendix, for model comparisons we use the (1 Gpc) 3 volumes, all other results are based off the (2.8 Gpc) 3 volume.We will show how the results from the simulations depend on sample composition in mass, redshift, and cool-core fraction, as well as algorithmic choices, and how these choices affect the match with observational results.First, we introduce the FLAMINGO simulations, halo selection, and X-ray calculation in Section 2. We then focus on the effect of model variations in FLAMINGO on the scaling relations in Section 3. We study thermodynamic profiles, and their susceptibility to different weighting schemes and cluster selections in Section 4. Finally, we discuss our results and offer conclusions in Section 5.

Simulations overview
FLAMINGO (Full-hydro large-scale structure simulations with allsky mapping for the interpretation of next generation observations) is a large suite of hydrodynamical cosmological simulations, covering enormous cosmic volumes.The flagship run comprises a region of (2.8 Gpc) 3 , which is ideal for the statistical studies of clusters.The simulations are described in detail in Schaye et al. (2023).A unique feature is the machine learning-aided calibration of the stellar and AGN feedback to reproduce the low-redshift cluster gas fractions at  500c and the  = 0 galaxy stellar mass function (Kugel et al. 2023).
Variations on the fiducial physical model in (1 Gpc)3 volumes are made by shifting the observed gas fractions (or galaxy stellar mass function) up or down with a multiple of their uncertainty () (see Table 1) and recalibrating the model to fit those new data points.For reference, at  500c = 10 14 M ⊙ models fgas+2 and fgas−8 have gas mass fractions of ≈ 0.10 and 0.05, respectively, while the fiducial model has 0.08 in agreement with observations.We note that if we interpret the gas fraction variations as horizontal shifts in the plot of the gas fraction as a function of mass, then the mass  500c where the gas fraction equals half the cosmic baryon fraction (  gas ≈ 0.08), is about 0.2 dex lower and 0.5 dex higher for the models with the highest (fgas+2) and lowest (fgas−8) gas fractions, respectively.Besides the feedback variations, FLAMINGO also includes cosmology variations in 1 Gpc volumes.We find that the cosmology variations have a negligible impact on the cluster properties investigated here (see Schaye et al. (2023) for a comparison of cluster scaling relations), and hence we do not include them in this paper.
In this work we make use of the fiducial resolution FLAMINGO simulations ( gas = 1.07×10 9 M ⊙ ), and use the high ( gas = 1.34× 10 8 M ⊙ ), and low ( gas = 8.56 × 10 9 M ⊙ ) resolution simulations for convergence testing in Appendix C, where we show the convergence is generally excellent for the thermodynamic profiles.The largest volume and all the model variations use the fiducial resolution.An overview of the simulations used in this work is presented in Table 1.
The mass ranges over which the simulations have been calibrated depend on the resolution and the box size used for the calibration runs.The m9 simulations have been calibrated to the  = 0 galaxy stellar mass function in the stellar mass range 10 9.92 − 10 11.5 M ⊙ and low- cluster gas fractions in the mass range  500c = 10 13.5 − 10 14.36 M ⊙ .The lower (m10) and higher (m8) resolution simulations have the same upper limit for the stellar mass range but lower limit of, respectively, 10 11.17 M ⊙ and 10 8.67 M ⊙ .The cluster gas fraction are calibrated from the same lower limit, but have different upper limits, respectively, 10 14.53 M ⊙ and 10 13.73 M ⊙ (Kugel et al. 2023).
FLAMINGO uses the open source simulation code swift (Schaller et al. 2023) and the sphenix Smooth Particle Hydrodynamics (SPH) scheme (Borrow et al. 2022).Neutrinos are included as particles (with a summed mass of 0.06 eV using the f method (Elbers et al. 2021)).Initial conditions are generated with a modified version of monofonic (Hahn et al. 2021;Elbers et al. 2022), and the default cosmology is the '3x2pt + all external constraints' from the dark energy survey year 3 (Ω m = 0.306, Ω b = 0.0486,  8 = 0.807, H 0 = 68.1, s = 0.967; Abbott et al. 2022).
An example of a FLAMINGO halo is shown in Fig. 1, where we show mass-weighted projected thermodynamic quantities.

X-ray luminosities
We generate tables of collisional-and photo-ionized X-ray model spectra using Cloudy (version 17.02;Ferland et al. 2017).Contrary to the common practice of using separate packages for the radiative cooling rates used in the simulation and for X-ray spectra generated in post-processing, we self-consistently use Cloudy for both.To compute the X-ray spectra we follow the methods used by Ploeckinger & Schaye (2020;hereafter PS20) for the radiative cooling rates used in FLAMINGO.We use the same density (10 −8 cm −3 ≤  H ≤ 10 6 cm −3 ), redshift (0 ≤  ≤ 9) and metallicity (primordial to 3 ⊙ ) ranges, but add information about the X-ray luminosity.We limit the analysis to gas with temperatures 10 5 K ≤  ≤ 10 9.5 K, as this is the relevant regime for X-ray clusters.In line with PS20, we use the UV and X-ray background from Faucher-Giguère (2020) modified at  > 3 to match the effective photo-ionisation and photo-heating rates before helium reionisation.
To compute the element-by-element emissivities, we generate differential spectra isolating the contribution of a single metal (analogous to the Wiersma et al. 2009a approach).We do this for 9 metals (C, N, O, Ne, Mg, Si, Ca, S, and Fe), 7 of which are explicitly traced by FLAMINGO (Ca and S are assumed to trace Si with mass ratios of 0.094 and 0.605) 3 .We make use of the solar abundance table of Asplund et al. (2009).A final spectrum can then be generated by taking the base spectrum, which lacks all 9 metals, and putting their respective contributions back in according to their individual abundances.This allows for non-solar relative abundances to be captured.These differential spectra are generated for every grid-point in metallicity, density, temperature, and redshift.Helium is treated differently because its contribution to the free electron abundance is non-negligible, every metallicity bin corresponds to a different helium abundance, and for every such bin the differential spectra are generated.
Metals which are not taken out individually are always assumed to be at the metallicity corresponding to the helium fraction assuming a solar abundance pattern4 .
We employ a single zone cloudy model, which is appropriate for unshielded or ionised gas, and analogous to the Wiersma et al. (2009a) approach.PS20 account for self-shielding in dense, cool gas, but self-shielding is not relevant for X-ray clusters.A cosmic ray background is used with the default value from Cloudy.PS20 change the value of the cosmic ray background depending on the shielding column, but this is again not relevant for X-ray clusters (increasing the cosmic ray background by an order of magnitude has a < 0.3% effect on the result and only at the lowest temperatures).No grains are used, since we are only working in the T > 10 5 K regime.This is a safe assumption as figure 7 of PS20 shows that there is assumed to be no dust at these temperatures.

X-ray table generation
At each redshift, the spectra computed by Cloudy are integrated over the  = 0 observer band.The table stores both the photon emissivity [s −1 cm 3 ], and energy emissivity [erg s −1 cm 3 ].For the energy version of the table the last column of the Cloudy Table 1.Hydrodynamical simulation runs.The columns list the simulation name, the number of standard deviations by which the observed stellar masses are shifted before calibration Δ * , the number of standard deviations by which the observed cluster gas fractions are shifted before calibration Δ  gas , the comoving box side length , the AGN feedback implementation (thermal or jets), the number of baryonic particles  b (which equals the number of cold dark matter particles  CDM ), the number of neutrino particles   , the initial mean baryonic particle mass  b , the mean cold dark matter particle mass  CDM , the comoving gravitational softening length  com , and the maximum proper gravitational softening length  prop .output file "diffuse.dat" is used.This contains the total diffuse emisssivity 4   (erg cm −3 s −1 ) (continuum + lines), where  is the frequency and   the specific intensity.For the photon count version of the table the output from the Cloudy files "XSPEC_diffuse_reflected.FITS", "XSPEC_diffuse.FITS", "XSPEC_lines.FITS" and "XSPEC_lines_reflected.FITS" (all in units of photons s −1 cm −2 bin −1 ) is summed.To obtain a photon emissivity from this flux, the value is divided by the depth of the shell.

Identifier
Calling the emissivity spectrum obtained from Cloudy   and the midpoint energy of the spectral bins   , the following integral is performed.
where  [ − − + ] (erg cm −3 s −1 ) is the total emissivity over the energy range  − to  + .The boundary values of this energy range are fixed for an observer at  = 0, and hence they are defined as We create integrated broadbands covering energy ranges corresponding to observational campaigns: This work makes use only of the ROSAT band.For each grid point in the tables, the emissivity of these bands is integrated using the trapezoidal rule.Finally, to account for the large dynamic range in emissivities due to the scaling with the square of the density, we divide all the emissivities by the square of the hydrogen density giving units of [erg s −1 cm 3 ].We store log 10 Ĵ[ − − + ] in the tables.

X-ray interpolation
With the Cloudy-generated table we interpolate the logarithmic X-ray broadband linearly to the logarithms of the density, temperature and elemental abundances of the simulation resolution elements (particles in the case of FLAMINGO), following the five steps below.
Step 1 The density bin (idx_n) and temperature bin (idx_T) corresponding to the resolution element are found and the absolute offset from the bin centres (dx_T, dx_n).
Step 2 The abundances and ratio-to-solar are computed.The abundance of each element (helium, carbon, nitrogen, oxygen, neon, magnesium, silicon and iron) is computed using where   and  H are the element and hydrogen mass fractions taken from the simulation and   the element atomic mass in units of the proton mass.Ratio-to-solar () is the ratio of the element abundance from the simulation to the solar abundance Note that the simulations predict absolute abundances and that the results therefore do not depend on the assumed solar abundances of the elements tracked by the simulations provided the same solar abundances are applied throughout the analysis.
Step 3 The appropriate helium bin and the relative offset to that bin centre are found, using the abundance from Step 2.
Step 4 The interpolation over the helium, temperature, density and redshift axes of the table is performed for all 9 differential values of the luminosity with a single metal missing (9 entries with only a single metal present, 1 entry with no metals present) Step 5 The contributions from all metals are added based on their ratio-to-solar.With  no metals representing the interpolated emissivity with all 9 tracked metals removed, and   the interpolated emissivity where metal  is present, this is computed as log 10  final = log 10  no metals + ∑︁  log 10   ×   (7)

Halo selection
In this work we select all haloes for which the mass,  500c , is larger than 10 13 M ⊙ from the halo catalogs of the FLAMINGO simulations.These are generated using the 6D Friends-Of-Friends subhalo finder Velociraptor (Elahi et al. 2019).Physical quantities are then computed using the Spherical Overdensity and Aperture Processor (SOAP), a tool developed for the FLAMINGO project.SOAP computes (sub)halo properties using as an input only the (sub)halo centers and particle membership from the halo finder.
For the fiducial model, the L2p8_m9 and L1_m9 simulations have, respectively, 2,031,904 and 92,015 haloes with  500c > 10 13 M ⊙ , 113,168 and 5,164 haloes with  500c > 10 14  ⊙ , and 461 and 15 haloes with  500c > 10 15 M ⊙ at  = 0.The model comparison in Section 4.4 uses the L1_m9 simulations, all other parts of this work are based of the L2p8_m9 simulation.
We remove all star forming gas particles from our analysis, because their temperature reflects the imposed entropy floor and represents the average of the unresolved multiphase interstellar medium, in almost all cases they are anyway below 10 5 K. Particles heated directly by an AGN in the last 15 Myr are also removed.Because the AGN feedback model is likely not realistic on small scales, recently heated particles may have unrealistic temperatures and densities.Both the removal of star forming and recently heated gas particles have a negligible effect on the results.To compare with observations, we furthermore only consider particles which could reasonably emit non-negligible amounts of X-ray photons by imposing a temperature floor of 10 5 K. Without such a temperature floor, mass-and volumeweighted thermodynamic quantities would be less informative for our purposes.

Thermodynamic profiles
To create radial profiles, we consider all particles eligible according to the above criteria.These are assigned to 30 logarithmically spaced radial bins between 0.01 R 500c and 3 R 500c .For each radial bin, we compute either the volume-, X-ray-or mass-weighted quantity of interest, where   is the value of the quantity for particle ,   is the particle mass,   its SPH volume (  =   /  ) and  , its X-ray luminosity in the ROSAT band (0.5 − 2.0 keV).The summation runs over all particles contained within the radial bin.All our profiles are 3D profiles, we have checked the effect of projection in Appendix B and find it to be non-negligible.However, lacking a method to undo the deprojection in observational data, our work is limited to 3D space.When presenting the radial profiles, we ensure that at least 16% of haloes have at least 1 particle in each bin, to ensure a proper representation of the 1 error.This leads to different inner stopping radii for different mass haloes.

Normalisation
Density The density profile is computed from the simulation using the free electron number density of each particle.The electron number density is computed internally in FLAMINGO by interpolating the tables from Ploeckinger & Schaye (2020), which take into account the ionisation state given the temperature-, density-, metallicity-and redshift-dependent UV and X-ray background and interstellar radiation field, for each particle.
Pressure The hot gas electron pressure is defined as where  e is the free electron number density,  the temperature of a particle in our simulation, and  B the Boltzmann constant.The pressure is normalised relative to the virial equilibrium expectation  500 (see Appendix D1 for a derivation) where  is the gravitational constant, and  () the Hubble parameter.We use  B = Ω b /Ω m = 0.159 the cosmological baryon fraction in FLAMINGO,  = 0.59 the mean particle mass, and  e = 1.14 the mean particle mass per free electron.
Temperature We normalise the temperature profiles to  500 (see Appendix D2 for a derivation) with  p the proton mass.

Entropy
The entropy is defined as where both the temperature and electron number densities are taken from the particles in our simulation.This is normalised relative to  500 (see Appendix D3 for a derivation) Metallicity We define the metallicity as the abundance of iron relative to hydrogen in solar units, where  Fe and  H are the mass fractions of respectively iron and hydrogen.We use the solar abundance ratio from Asplund et al. (2009),  Fe,⊙ / H,⊙ = 10 7.5 /10 12 .

EFFECT OF MODEL VARIATIONS ON THE SCALING RELATIONS
One of the primary purposes of FLAMINGO is the study of galaxy groups and clusters.Even though the simulations have been calibrated to match observed gas fractions at R 500c for clusters at  ≈ 0.1 − 0.3 and  500c = 10 13.5 − 10 14.36 M ⊙ , this does not guarantee the 10 34 10 35 10 36 10 37 L X (erg s 1 kpc 2 ) 10 2 10 3 10 4 10 5 Pressure (K cm 3 ) 10 7 10 8 Temperature (K) 10 3 4x10 3 8x10 3 Entropy (keV cm 2 ) 1.5 1.0 0.5 0.0 [Fe/H] 10 5 10 4 10 3 10 2 n e (cm 3 )   reproduction of cluster scaling relations, as those depend on density, temperature, and metallicity profiles, as well as on the clumpiness and multiphase nature of the ICM.Furthermore, calibration does not include all masses, and is only performed at low redshift, which also makes the scaling relations a test for the predictive power of the simulations.Schaye et al. (2023) already showed that at  = 0 the different resolutions match various observed scaling relations well across more than two decades in mass.Here we study the model dependence of the same scaling relations.
Figure 2 shows the median X-ray luminosity -halo mass (left), and temperature -halo mass (right) relations for all groups and clusters in the different models in a (1 Gpc) 3 volume at the fiducial resolution.Temperatures are mass-weighted and include all particles with  > 105 K, we have compared with emission-weighted temperatures and found the differences to be within 10% for these global properties.The small horizontal arrows indicate the systematic shift applied to observational data with hydrostatic-equilibrium inferred masses, which corresponds to the value found for the hydrostatic mass bias (0.743) during the calibration in Kugel et al. (2023), which agrees with their assumed priors based on the observations of Hoekstra et al. (2015) and Eckert et al. (2016).We neglect any corrections on the observed quantities due to the hydrostatic bias as these are negligible (e.g. the X-Ray luminosity will be dominated by radii much smaller than  500c ) The observed X-ray luminosities have been shifted to the 0.5-2.0keV band using PIMMS 5 (Mukai 1993).Observational data was grouped into a limited number of bins per data set, where the error bars show the scatter between individual objects in that bin.
The different models show a consistent offset between models with different gas fractions, which is most pronounced in the luminosity -mass relation.Apart from the offset, the slope of the scaling relation varies with the model variations.The normalisation and slope, measured by fitting a single power law between 10 13.5 − 10 15 M ⊙ are tabulated in Table 3, according to with Y either the logarithmic luminosity or temperature,  the slope and  the normalisation.
As can be seen from Fig. 2 the temperature -mass relation shows little variation in the slope and normalisation above 10 13.5 M ⊙ , but the luminosity -mass relation varies by 40% in slope and by and order of magnitude in normalisation.The offset can be explained by the different gas fractions, with less gas producing less X-ray emission.The change in slope can be understood in the light of less massive objects having larger fractional changes in their gas fractions between the model variations.The paucity of gas, and hence X-ray luminosity, at lower masses, combined with the relatively small change at higher masses, yields a change in slope as seen in these relations.It should be noted that, except for the most extreme model variations, the medians of all models are within 1 of each other for both scaling relations across the entire mass range.
We have verified that using core-excised quantities (i.e.removing radii  < 0.15  500c ) or spectroscopic-like temperatures (Mazzotta et al. 2004;Vikhlinin 2006) does not significantly change the scaling relations, nor their agreement with data.In summary, the simulated cluster scaling relations are numerically converged and agree well with observations, the FLAMINGO models provide observationally motivated variations which are measurably different in these integrated quantities.

THERMODYNAMIC PROFILES
In the previous section we have shown that the scaling relations obtained from FLAMINGO clusters match observational data, and do so across a wide range of halo masses.Scaling relations are, however, by their very nature limited to comparing global properties of haloes.For many observables, global properties can be dominated by one region of the halo, such as the outer halo for mass-weighted temperatures, whereas X-ray properties tend to be sensitive to the core.There are scenarios where clusters fall on the scaling relations, but have unrealistic radial profiles.This section compares the density, temperature, metallicity, entropy, and pressure profiles from FLAMINGO with observational data.
As described in Section 2.5, we exclude star forming particles, particles that have recently been directly heated by AGN feedback, and cool gas ( < 10 5 K).We find that the exclusion of these three categories of particles has a very small impact on most quantities.Only the innermost part of the mass-and volume-weighted temperature profile is significantly affected when cool particles are included.
Using the above selection criteria, we compute the profiles for all  500c > 10 13 M ⊙ haloes.This is done using 30 equally spaced logarithmic bins between 0.01  500c and 3  500c .In each bin we compute the volume-, emission-, or mass-weighted average of each physical quantity.Volume weighting is performed using the SPH volume of each particle (eq.8).Emission-and mass-weighting are performed analogously, replacing the SPH volume by either the [0.5 − 2.0] keV X-ray luminosity (eq.9), or the mass of each particle (eq.10).Using spectroscopic-like temperatures (Mazzotta et al. 2004;Vikhlinin 2006) results in negligible changes to the temperature profile over the entire radial range compared to X-ray luminosity weighting, hence we choose not to show it.
The smallest radii at which we show the median profile is determined by the radius at which at least 16% of the haloes no longer have particles in the radial bin, this is the point where the 1 error can no longer be reliably determined.This results in mass-, and resolution-dependent cut-off radii.For higher mass objects, we can follow the curves to smaller / 500c simply because those objects have a higher density.
The observations with which we compare our simulated cluster sample are listed in Table 2, where we tabulate the sample size, mass range, median mass, and the redshift range of the observed sample.All observed normalised radii ( = / 500c ) have been multiplied to account for the hydrostatic mass bias  hs / = (0.743) 1/3 obtained during the calibration of FLAMINGO, as well as their normalised observed values.The metallicities found by Ghizzardi et al. (2021) have been corrected to the solar iron abundance in Asplund et al. (2009).All datasets have been re-scaled to the FLAMINGO baryon fraction (  B = 0.159), mean particle mass ( = 0.59) and mean particle mass per free electron ( e = 1.14).
To interpret the simulated thermodynamic profiles, and to compare them with observations, we first have to understand the impact of different mass selection effects and weighting choices.With this in mind, we will first study the mass dependence of simulated thermodynamic profiles in Section 4.1, then consider the different choices for weighting particles in Section 4.2, and finally the difference between cool-core and non-cool-core clusters in Section 4.3.These three topics will set the stage for the FLAMINGO model comparison, which we show in Section 4.4.

Mass dependence
The temperature, pressure, and entropy profiles in this work are normalised by the average value of each quantity expected within  500c assuming an isothermal sphere in hydrostatic equilibrium and the virial relations ( 500c ,  500c , and  500c ).This normalisation is naturally mass and redshift dependent as the total gravitational force exerted by a virialized structure depends on its mass, and the critical density evolves with redshift.In order to compare the profiles of different mass haloes, we use the normalised radius (/ 500c ).If all objects, of all masses, would be spherically symmetric, self-similar, virialised structures in hydrostatic equilibrium, then the normalized thermodynamic profiles would be mass independent, assuming that there are no non-gravititional contributions to the heating, cooling or gas kinetics.In reality this is most likely not the case, as the halo formation time and merger rate, and hence the departure from sphericity and equilibrium, are mass dependent, as well as the normalised radius out to which galaxy baryonic processes can have a discernible and significant impact.In addition, halo concentrations and stellar-to-halo mass ratios vary systematically with halo mass.
We find that the normalised profiles still show a strong mass dependence.Since observational campaigns tend to target, and are thus biased towards the higher-mass objects, and since mass estimation from observations is fraught with uncertainties, the mass dependence is important to keep in mind when comparing profiles from simulated objects with observations.
Figure 3 shows five  = 0 median thermodynamic profiles (temperature, density, pressure, entropy, and metallicity) in 5 different mass bins (each 0.5 dex wide), stretching from  500c = 10 13 M ⊙ to 10 15.5 M ⊙ .The 16 th -84 th percentile are shown for the 10 14 M ⊙ to 10 14.5 M ⊙ mass bin, but are comparable in both shape and magnitude for all other masses.The observations we compare with are coloured by their median mass in an identical fashion to the simulations.All profiles are clearly mass dependent, despite the applied normalizations.From one mass bin to the next, the magnitude of the profiles at fixed radius changes by ≈ 0.1 − 0.3 dex.The normalized temperature, normalized entropy and metallicity decrease with mass, whereas the density and normalized pressure increase.The peak of the temperature profile moves to smaller normalized radii when the mass increases.For the entropy profile, we observe that the flat part of the profile shifts to smaller normalized radii for higher-mass objects.
Flattening of the entropy profile compared to observations has previously been discussed in detail by Altamura et al. (2023), who found it to be common in simulations, and stronger than our simulations, but leave open the exact origin of the phenomenon.Compared to recent work, we find a similar mass dependence to MilleniumTNG (Pakmor et al. 2023), and are in broad agreement with TNG-Clusters (Lehle et al. 2023).Compared to those works, FLAMINGO clusters show a strong drop in the temperature in the core, as well as a drop in the core entropy, as opposed to flat profiles.In Section 4.3 we will show that this could be due to a different coolcore fraction.The mass dependence is also similar to what is seen in clusters from The Three Hundred project (Li et al. 2023), when taking into account that profiles evolve self-similarly as we will show in Section 4.5.
Except for the metallicity, which is ≈ 0.3 dex higher than observed, our simulated clusters follow the observed relations down to small radii, when comparing with observations which have a similar median mass.The metallicity is also high compared to previous simulations which tend to produce less iron in the core compared to observations (Vogelsberger et al. 2018;Pearce et al. 2021).
The extrapolations of the temperature and entropy profiles to smaller radii than sampled by the simulation do not seem to match the observations by Sun et al. (2009) ( 500c ≈ 10 13.5 M ⊙ , green-brown line) and McDonald et al. ( 2014)6 ( 500c ≈ 10 14.75 M ⊙ , light-blue line).However, we will show that at the smallest radii the profiles of cool-core and non-cool-core clusters diverge, with the non-cool-core objects agreeing with data.
For all masses, the simulated pressure profiles slightly undershoot the observed pressures.In addition, the sloped are slightly too shallow, underprediciting the pressure in the centers and overpredicting it in the outskirts.
The consistent discrepancy between the predicted and observed metallicity could be explained by the assumed nucleosynthetic yields being too high, or by the total stellar masses being too high in the simulation.Figure 11 of Schaye et al. (2023) compares the cluster stellar masses to observations.For  500c > 10 14 M ⊙ the total stellar masses are indeed too high, however, if only the stellar mass within 50 kpc apertures is included, then the stellar masses are too low, which implies that the comparison is sensitive to the difficult to observe extended low surface brightness stellar envelopes of galaxies.
For the different panels, where different observational data sets at the same mass are available, the difference between consecutive mass bins is of roughly the same magnitude as the difference between those data sets (with the exception of the metallicity profiles).This implies that any differences seen between the simulated profiles and observations, could potentially be explained by uncertainties in the observations.

Weighting scheme
When calculating a cluster thermodynamic profile from simulations, an important choice is how to weigh the contribution of individual particles.The arguably most intuitive method, to weigh the contribution of each particle by the fractional volume it occupies in a spherical shell (see eq. 8), is not necessarily representative of what would be  L2p8_m9) and massive clusters ( 500c = 10 14.5 − 10 15.0 M ⊙ ), using 3 different methods to weigh the particles when constructing the profiles.X-ray weighting is indicated by the solid line, mass weighting by the dashed line, and volume weighting by the dash-dotted line.measured observationally.Because thermodynamic profiles are inferred from X-ray observations, X-ray bright gas will dominate the inference.The comparison with volume-weighted profiles from simulations is only fair if all volume elements within a spherical shell have the same X-ray brightness.
To test this assumption, and to better approximate observational inferences, we compare three different weighting schemes.Volumeweighting, mass-weighting (see eq. 10), and X-ray-weighting (see eq. 9) in which particles are weighted by their density-and temperature-dependent X-ray luminosity.We expect the latter, which was used in Figure 3, to be closest to what is measured observationally.As a caveat we note that observationally the density, temperature and metallicity are measured from the radially-binned spectrum, we plan to include this in future work producing realistic synthetic X-ray observations.
In Figure 4 we compare the resulting profiles for the fiducial model (L2p8_m9) and objects with masses between  500c = 10 14.5 and 10 15.0 M ⊙ (a total of 12354 objects).The physical radius (top axis) is computed for the median mass of the sample.The differences between the different weighting schemes are generally small, and comparable to or smaller than the scatter between different observations.The scatter in the simulations is similar for all weighting schemes and comparable in magnitude to what is shown in Fig. 3.The differences due to weighting are generally much smaller than an 0.5 dex change in mass, as seen in Figure 3.However, for the temperature, density, and entropy the weighting causes larger differences in the core ( ≪  500c ) than a 0.5 dex change in mass.Compared to volume-weighting, for mass-and particularly for X-ray-weighting, the inner temperature and entropy are lower, the density and metallicity higher, and the pressure is unaffected.Differences in the temperature profile are limited to  < 0.1  500c .For the density and entropy profiles the differences increase towards both small and large radii.The pressure and metallicity profiles show the weakest dependence on the weighting scheme and an increasing difference with radius.
These tendencies are in line with expectations.Since mass and X-ray weighting favour denser gas, they yield a higher density, and this denser gas tends to have a lower temperature, lower entropy and higher metallicities.The differences increase towards small radii, where the gas tends to be more multiphase.Because the different phases tend to be in pressure equilibrium, the pressure profiles are nearly the same for the different weighting schemes.The difference at larger radii seen for the density and entropy might be due to particles bound to satellites, which will have higher densities.The differences between applying different weights are within the scatter between different observational data sets.We note, however, that when discrepancies with the observational data exist, they tend to diminish when moving from volume to mass weighting, and even more so when changing to X-ray luminosity weighting.

Cool cores
We distinguish cool core (CC) and non-cool core (NCC) clusters in our simulations by measuring the radiative cooling time within 0.048 500c .This is in line with Hudson et al. (2010), who compared 16 CC metrics and found this one to be the most distinguishing feature.We compute the cooling rate from the cooling tables used for FLAMINGO (Ploeckinger & Schaye 2020).These are interpolated to the X-ray-luminosity-weighted temperature, density, and metallicity of particles within 0.048 500c (the median NCC cluster has 10 particles within this radius), where we make the same selection to exclude recently heated, star forming and cool particles as in the rest of this work.Using the radiative cooling rate Λ, we compute the  L2p8_m9) for massive clusters ( 500c = 10 14.5 −10 15.0 M ⊙ ), comparing cool-core (dashed) and non-cool-core samples (solid).The two samples show large differences in the cluster cores.
cooling time 2 2 e Λ( e , , , ) , with  the gas metallicity.CC clusters are often defined as objects that have a central cooling time below a critical value, which tends to be set to a value between 1 and 5 Gyr.To discern strong CC clusters from NCC clusters, we define the CC sample to be objects with  cool < 1 Gyr.Conversely, the NCC sample has  cool > 5 Gyr.
We compare the (emission-weighted) cluster gas profiles for CC and NCC clusters of mass between  500c = 10 14.5 and 10 15.0 M ⊙ in Fig. 5 7 .In this mass bin, this selects 16% of the total number of clusters as NCC (1988) and 22% as CC (2785) (see also Fig. 8).We note that the sample median central cooling time of all objects is 2 Gyr, hence the total sample is more similar to CC than to NCC clusters.Since the cooling time decreases with density, we expect CC clusters to have denser cores.Cooling is also more efficient at sub-virial temperatures, and at higher metallicity, particularly below 1 keV where metal-line cooling becomes important.This gives rise to the expectation that CC clusters have cooler, higher metallicity gas in their cores compared to their NCC counterparts, something borne out in observations (e.g.De Grandi & Molendi 2001;Lovisari & Reiprich 2019).We reproduce that trend in Fig. 5, where the median relations for NCC and CC are shown.The differences between CC and NCC clusters are only large for  ≲ 0.1 500c .The offset seen in the pressure profile at small radii could be explained if CC objects have more concentrated mass profile, which is borne out by the higher density in the core.We note that for temperature, density and entropy, the difference is the core surpasses the 1 region shown as a shaded band.
CC clusters are sometimes defined as clusters with a low central ( ≲ 0.012  500c ) entropy, where the threshold tends to be in the range 30 − 50 keVcm 2 (McDonald et al. 2013).We see that our CC sample, which is selected based on the central cooling time, has a much lower central entropy compared to its NCC counterpart.The entropy at the smallest radius to which our profiles extend is below 50 keVcm 2 for the CC sample, while it is more than double that value for the NCC sample.Our classification is thus in line with this alternative definition of CC clusters.In Appendix E we show that the cooling time criterion used clearly separates the CC and NCC populations in density, temperature and entropy, which shows that this is a robust way of identifying CC clusters in FLAMINGO.

Evolution of the cool core fraction
We study the evolution of the CC fraction of clusters over a wide mass range.As before, we define a CC cluster by its central ( < 0.048  500c ) cooling time (eq.18).At  = 0 we used 1 Gyr as the critical value for an object to be classified as a CC.The CC fraction we obtain is similar to that obtained by TNG-Cluster (Lehle et al. 2023), and slightly higher than IllustrisTNG (Barnes et al. 2018).Assuming the radiative cooling is dominated by bremsstrahlung, and that the virial temperature evolves self-similarly as  () 2/3 (eq.13), the cooling time of a cluster will evolve as  cool ∝  () −5/2 .Employing a self-similarly evolving critical value, Fig. 6 show that the CC fraction is almost constant with redshift.Only at the highest masses does it decrease slightly.The dashed lines show the result when not accounting for self-similar evolution, in that case a large fraction of objects is classified as CC at higher redshift.This strong evolution of the CC fraction is analogous to what has been found with a slightly different cooling time criterion by Barnes et al. (2018) for the Illus-trisTNG simulation.However, it seems in conflict with observations which find an almost non-evolving CC fraction for a non-evolving cooling time criterion (McDonald et al. 2017;Ruppin et al. 2021).

Model comparison
Using the knowledge from the previous three sections, we compare thermodynamic profiles for the different model variations in the FLAMINGO suite.We keep in mind the mass-dependence, which, depending on the quantity and radius (see Fig. 3) can give a 20-100% difference between consecutive mass-bins; the weighting dependence, which showed a strong radial dependence and can, at specific radii, change the result by 50% (see Fig. 4); and the coolcore non-cool-core distinction, which is especially relevant for the innermost 0.1 500c of the thermodynamic profiles (see Fig. 5).
Fig. 7 shows the different model variations for two different mass bins, large groups ( 500c = 10 13.5 − 10 14 M ⊙ ) in the left panel, and clusters ( 500c = 10 14.5 − 10 15 M ⊙ ) in the right panel.For all models, the scatter is comparable to what is shown in Fig. 3.The observational data we have selected is only available for massive objects, hence we do not show it for the lower mass bin.First, we focus on the gas fraction model variations, which are calibrated to observed mass-dependent gas fractions shifted up or down by different numbers of observational  (see Table 1).These are shown in different shades of blue in Fig. 7.For  ≳ 0.1  500c , we see the expected trends: lower gas fractions correspond to lower densities, which implies a lower gas mass, and a lower gas pressure.The changes in density and pressure then determine the change in temperature which can go up or down depending on whether the density or pressure changes more.In FLAMINGO we find that the peak temperature decreases with increasing gas fraction.From the combination of a lower tem-perature and higher density with increasing gas fractions it follows that entropies are lower (see eq. 14).
We furthermore find that the gas fraction variations show the largest differences at intermediate radii between 0.1 − 0.4 500c , and converge at both larger and smaller radii (only visible for the cluster mass bin), except for the pressure profiles.The pressure profiles show an increasing difference between the models for decreasing radii, with lower gas fractions yielding up to 40% lower pressures in the core of clusters, and a factor 4 lower pressure in the core of groups.
The lower gas fraction variations have a slightly stronger AGN feedback.The gas pushed out by the AGN is seemingly predominantly deposited between 0.1 − 0.4 500c , as the difference between the models is largest there.The stronger AGN feedback deposits an increased amount of metals at these radii, explaining the divergence in metallicities in this radial range for clusters.At even larger radii, the models converge again.The correlation between temperature and gas fraction seen between 0.1 − 1 500c could also be due to the cumulative effect of the AGN.With the cosmic baryon fraction fixed between different gas fraction variations, a lower cluster gas fraction implies more gas has been expelled at early time.While clusters grow, this hot expelled gas will flow towards the cluster, naturally leading to higher temperatures in the outer regions.Though McCarthy et al. (2011) found that the infalling gas has higher entropy not because it was heated at earlier times, but simply because it is taking the place of the low entropy material that was ejected.
The effect of a change in the cluster stellar fraction is seen in the stellar mass variation model, shown by the orange line.This model has a lower stellar mass, but the same gas fraction.Though almost identical to the fiducial model in temperature, density, pressure, and entropy, the decreased production of metals, due to a lower stellar fraction, yields a consistently lower metallicity at all radii.The lower metallicity for groups in the jet_fgas−4 model can be explained by the fact that at this halo mass the stellar mass of objects in this model variation is even slightly below the stellar mass variation model ( * − ).
The two jet models show very different behaviour.Compared with the fiducial thermally driven AGN feedback, the kinetic jet feedback clearly expels more gas from the cluster core.This is seen in a depressed central density for clusters.The remaining gas in the core is heated to a significantly higher temperature compared to all thermal feedback models.The combination of a higher temperature and a lower density leads to a significantly higher entropy in the core, while the pressure remains similar.This is not true for groups, which can be explained by the gas fraction for the Jet model being similar to the fgas+2 variation in this mass range, and consequently, the Jet model is more similar to the increased gas fraction variation for groups.The jet models also show a strong deficit of metals in the core of clusters compared to the fiducial model.We interpret this as a sign that our kinetic jet feedback model is more effective at expelling enriched gas from the cluster core.Previous work based on the OWLS simulations (Schaye et al. 2010), which use thermally-driven AGN feedback, has found that the profiles are mostly affected by ejection of gas from the high- progenitors of  = 0 clusters (McCarthy et al. 2011), which raises the interesting question of exactly when and where the feedback was actually injected.Though the relative lack of metals in the core is similar for the reduced gas fraction jet model (Jet_fgas − 4), the overall metallicity is significantly lower for this model, especially at larger radii.We note that for this latter model, the stellar fraction is significantly lower (Schaye et al. 2023), possibly due the expulsion of low-entropy gas limiting radiative cooling and thus reducing the overall central stellar production, and hence the  metal production.However, we remind the reader that all variations have been calibrated to match the stellar mass function over the same range of masses.
The combination of distinctly higher temperatures, lower densities and higher entropies, is a typical signature for NCC clusters seen when comparing with CC clusters in observations (e.g.Hudson et al. 2010).In Fig. 8 we compare the CC fraction of the jet models with all other models for clusters, finding that their fraction is indeed suppressed by more than 50%.The different CC to NCC ratio in the jet sample could then possibly account for some of the difference we see in the temperature, density and entropy profiles.NCC clusters are also observed to have lower core metallicities (Lovisari & Reiprich 2019), in line with the difference between the Jet models and the thermal AGN models in Fig. 7.

Matching the cooling-time distributions
Since the model variations have different strengths of their AGN and stellar feedback, they will also have different central cooling times.We have shown in Fig. 5 that the thermodynamic profiles of CC and NCC clusters differ, especially in the core region.The core region is also where we see a strong difference between the median profile for the jet models and all other variations, which begs the question of whether these different profiles give rise to different CC fractions between the model variations, and whether correcting for this would give consistent profiles between all models.
To study this effect, we first compute the cumulative distribution function (CDF) of the central-cooling times for clusters ( cool ( < 0.048  500c ), eq.18).We plot these distributions in Fig. 8 which clearly shows that the jet models have much larger central cooling times compared to all other variations.The median profiles for the jet models will thus be much more NCC-like compared to the other variations.
To test whether this can explain the large difference seen in the cores for the median thermodynamic profiles of clusters, we create a  matched sample of clusters across all model variations for the cluster mass bin ( 500c = 10 14.5 − 10 15 M ⊙ ), selecting for every object in the fiducial model (L1_m9) an object in each variation that has a central cooling time differing by less than 10% without replacement.
If such a matching object cannot be found for every model variation, the halo is discarded.In this manner, we match ≈ 1/3 of all haloes in the simulations.This matched sample should have nearly identical distributions of central cooling times, and we have checked that their CDFs are indeed nearly identical.For this matched sample of haloes, we recompute the median thermodynamic profiles for all variations.Fig. 9 shows these profiles for the temperature and entropy, for which the biggest difference with the two jet variations were seen.We see that after matching the cooling time distributions, the median profiles for the Jet and Jet_fgas − 4 models are more similar to all other variations, but still show differences in the core.This suggests that the differences in Fig. 7 are only partly due to the different CC/NCC sample composition.

Redshift evolution at fixed mass
We study the evolution of cluster gas profiles with redshift in Fig. 10 for the L2p8_m9 simulation, selecting objects with the same mass ( 500c ) at  = 0, 0.5, 1.0, and 1.5.Because haloes accrete mass between these redshifts, the most massive objects at  = 1.5 may no longer be in the same mass bin at lower redshifts.However, by selecting the same mass range, we can compare the physical properties of haloes at a fixed mass over cosmic time.At  = 1.5 there are 17 objects in this mass bin.At all redshifts, including  = 1.5, the scatter is almost identical to the error band shown for  = 0. We find that the temperature, density, pressure, and entropy profiles show little to no evolution with redshift, when normalised to their redshift-dependent  500 quantity, or the redshift-dependent critical density.This is expected since, although the Universe has expanded since  = 1.5, leading to the gradual dilution of gas, the virial quantities  500c account for the effect of this expansion.It suggests that in FLAMINGO non-gravitational physics may impact (proto-)clusters mostly at high redshifts, as there is almost no deviation from selfsimilar evolution below  = 1.5.In Appendix A we show that scaling relations also evolve close to self-similarly.
Observations of high-redshift galaxy clusters find a deviation from self-similar expectation in the core ( < 0.1  500c ), with the cores evolving less since at least  = 1 (Ghirardini et al. 2021).When using normalised quantities, this results in higher-redshift profiles being lower in the core compared to lower-redshifts, as the normalisation over-compensates for evolution.We see this most strongly in our temperature and entropy profiles, though the effect is not significant.Outside 0.1  500c they are consistent with self-similarity, but inside that radius, the evolution is slightly slower.Unlike what is seen in observations (McDonald et al. 2017;Ruppin et al. 2021), the density and pressure profiles do not show deviations from self-similarity in the core.
The metallicity shows a clear evolution with redshift.Higher redshift clusters have a lower gas metallicity, particularly in the core.This is expected, since fewer stars have been formed by those earlier times, resulting also in fewer SNe Ia, which means that less metals have been injected into the ICM.The decrease in metallicity with redshift at large radii is consistent with what has previously been found in IllustrisTNG, C-EAGLE, and Magneticum (Vogelsberger et al. 2018;Pearce et al. 2021;Angelinelli et al. 2023).In contrast to those simulations, we find that the core metallicity also decreases with redshift.

CONCLUSIONS
In this work we explored groups and clusters in the FLAMINGO cosmological hydrodynamical simulations, investigating the effect of mass selection, particle weighting, and cool-core selection, comparing their temperature, density, pressure, entropy, and iron abundance profiles to observations, and studied the evolution of cluster scaling relations and profiles.This has allowed us to conclude that, with the exception of the metallicity of the ICM, FLAMINGO clusters are a good match to objects in the real Universe, showing similar radial profiles and scaling relations.In particular, we found that (i) Even when normalised to virial expectations, thermodynamic profiles are sensitive to the mass of the selected objects.The normalised temperature, normalised entropy and iron abundance decrease with increasing mass, whereas the density and normalised pressure increase (Fig. 3).In particular, the density profile shows a radius-independent offset of ≈ 0.2 dex per 0.5 dex in mass, which is consistent with the view that more massive objects originate from larger fluctuations in the background density field, yielding higher overall densities at all radii.Offsets in other thermodynamic profiles are of similar size but appear only at either small or large radii.
(ii) Our fiducial weighting of particles when constructing the profiles is by X-ray luminosity.We find that using mass or volume weighting instead of X-ray weighting has a large, and non-trivial influence on the thermodynamic profiles (Fig. 4).While the temperature (pressure) profiles only show a difference of 0.3 dex and 0.1 dex in the innermost (outermost) regions, the density and entropy show a difference of 0.2 dex at both large and small radii, with a decreased difference at intermediate radii.
(iii) We define cool-core clusters using the central cooling-time, and show that in FLAMINGO they have lower core temperatures and entropies, but higher densities, pressures, and metallicities (Fig. 5).This is in line with observational results, and implies that the sample composition is important to consider when comparing simulated clusters with observed clusters as the latter tend to have sample compositions biased towards more cool-core clusters.The cool-core fraction evolves strongly with redshift when using a non-evolving cool-core definition.When correcting for self-similar evolution, we find the cool-core fraction to be almost constant with redshift (Fig. 6).
(iv) The models using jet-like AGN feedback instead of our fiducial thermally-driven AGN feedback give thermodynamic profiles with qualitatively different shapes (Fig. 7).Cluster cores are significantly hotter and less dense in these models, indicating that our kinetic jet feedback is effective at shaping the cluster cores.
(v) Except for the iron abundance, FLAMINGO galaxy clusters show little evolution in their thermodynamic profiles when normalised to their respective virial quantities (Fig. 10).This indicates that massive galaxy clusters evolve self-similarly.
We have shown that galaxy clusters from the FLAMINGO simulation evolve self-similarly, in line with expectations from observations.The thermodynamic profiles are sensitive to mass, even after normalisation to virial expectations, potentially due to the importance of non-gravitational physics for galaxy groups and clusters.For the cores of clusters, the weighting scheme chosen to measure physical quantities from simulations is important to consider.Finally, we have shown that the cool-core to non-cool-core composition of galaxy cluster samples should be matched to perform reliable comparisons of thermodynamic profiles.In future work it would be useful to compare virtual observations to the data, which would allow for a fairer comparison and enable accounting for selection and projection effects.and simulated profiles.At larger radii the differences shrink, but remain non-negligible.

APPENDIX C: RESOLUTION
We test the convergence properties of our simulations for our fiducial mass bin  500c = 10 14.5 − 10 15.0 M ⊙ , and show the results in Figure C1.We compare the three simulations L1_m8, L1_m9 and L1_m10, which use the same box size but differ in mass (spatial) resolution by factors of 8 (2) (see Table 1).It should be kept in mind that although all resolutions have been calibrated to the observed low-redshift galaxy mass function and cluster gas fractions, they are not the same physical models.
For the pressure profile the resolution has a negligible effect.The temperature profile shows converged results at large radii, while at small radii the peak of the temperature shifts inwards and increases slightly in normalisation for the highest resolution.Something similar can be seen for the entropy profile, where the plateau extends to smaller radii for L1_m8.The density profile shows a slight suppression at  <  500c but only for the highest resolution.The metallicity profiles are somewhat shallower for increasing resolution.In all cases, at small radii the highest resolution model is closer to the observations than an extrapolation of the intermediate resolution model.

Figure 1 .
Figure 1.Projected maps for a massive cluster,  500c = 10 15.03 M ⊙ , of 4 500c (6.3 Mpc) on a side from the L1_m9 simulation at  = 0. Clockwise from the top left the different panels show projected X-ray surface luminosity density in the ROSAT band (0.5 − 2.0 keV), pressure, temperature, free electron density, iron abundance and gas entropy, where the last five panels show mass weighted means.

Figure 3 .
Figure 3. From top to bottom, the different panels show the median temperature, density, pressure, entropy, and metallicity profile for galaxy groups and clusters in the fiducial model (L1_m9) of the FLAMINGO simulations at  = 0.All quantities are X-ray-luminosity-weighted.The color scale encodes halo mass ( 500c ), as indicated by the color bar.Different lines are for different 0.5 dex mass bins.Data points correspond to observational data sets listed in Table 2, and have a colour corresponding to their median mass.Error bars indicate the 16th and 84th percentiles when available, or the observational 1 error.The shaded region indicates the 16th-84th percentiles for the  500c = 10 14.0 − 10 14.5 M ⊙ mass bin.

Figure 4 .
Figure 4. Similar to Fig 3, but only showing the fiducial model (L2p8_m9) and massive clusters ( 500c = 10 14.5 − 10 15.0 M ⊙ ), using 3 different methods to weigh the particles when constructing the profiles.X-ray weighting is indicated by the solid line, mass weighting by the dashed line, and volume weighting by the dash-dotted line.

Figure 6 .
Figure 6.Fraction of objects classified as CC ( cool, (<0.048R 500c ) < 1 Gyr ×  () −5/2 ) as a function of redshift.When correcting for self-similar evolution, there is no change in the CC fraction with redshift.The dashed lines show the result when using a non-evolving critical cooling time of 1 Gyr.Shaded regions indicate the 1 bootstrap error on the CC fraction.

Figure 7 .
Figure 7. Similar to Fig. 4, but showing different FLAMINGO models in the same box size and for the same resolution as L1_m9 for groups with mass  500c = 10 13.5 − 10 14 M ⊙ in the left panel, and clusters with mass  500c = 10 14.5 − 10 15 M ⊙ in the right panel.Except for the jet models and the most extreme gas fraction variation (fgas − 8), all model medians for clusters are within the 1 scatter of the fiducial model.The two jet models show distinctively different behaviour in the cores of clusters.For groups, the effect of model variations is much larger.

Figure 8 .
Figure 8. Cumulative distribution function of the central cooling times ( cool ( < 0.048R 500c ) ) of massive clusters ( 500c = 10 14.5 − 10 15 M ⊙ ) for all L1_m9 model variations.The two jet models have significantly greater central cooling times.

Figure 9 .
Figure 9. Temperature and entropy profiles for the central-cooling-time matched sample of massive clusters ( 500c = 10 14.5 − 10 15 M ⊙ ).After matching, the Jet and Jet_fgas − 4 model models are closer to all other variations, but still discrepant.

Figure 10 .
Figure 10.As figure Cluster gas profiles for cluster with mass ( 500c = 10 14.5 − 10 15 M ⊙ ) at different redshifts indicated by the colour.Observations are coloured by their mean redshift.The physical radius (top axis) for  = 0 objects is shown.Except for temperature and metallicity very little evolution is seen.

Table 2 .
The observations to which the FLAMINGO data is compared.From left to right the columns list: Reference, sample size, median halo mass ( 500c ), and the redshift of the observed objects.

Table 3 .
Slope () and normalisation () of the scaling relations for simulations calibrated to different gas fractions (see eq. 17).