\textsc{The Three Hundred} project: The \textsc{Gizmo-Simba} run

We introduce \textsc{Gizmo-Simba}, a new suite of galaxy cluster simulations within \textsc{The Three Hundred} project. \textsc{The Three Hundred} consists of zoom re-simulations of 324 clusters with $M_{200}\gtrsim 10^{14.8}M_\odot$ drawn from the MultiDark-Planck $N$-body simulation, run using several hydrodynamic and semi-analytic codes. The \textsc{Gizmo-Simba} suite adds a state-of-the-art galaxy formation model based on the highly successful {\sc Simba} simulation, mildly re-calibrated to match $z=0$ cluster stellar properties. Comparing to \textsc{The Three Hundred} zooms run with \textsc{Gadget-X}, we find intrinsic differences in the evolution of the stellar and gas mass fractions, BCG ages, and galaxy colour-magnitude diagrams, with \textsc{Gizmo-Simba} generally providing a good match to available data at $z \approx 0$. \textsc{Gizmo-Simba}'s unique black hole growth and feedback model yields agreement with the observed BH scaling relations at the intermediate-mass range and predicts a slightly different slope at high masses where few observations currently lie. \textsc{Gizmo-Simba} provides a new and novel platform to elucidate the co-evolution of galaxies, gas, and black holes within the densest cosmic environments.


ABSTRACT
We introduce G -S , a new suite of galaxy cluster simulations within T T H project. T T H consists of zoom re-simulations of 324 clusters with 200 10 14.8 drawn from the MultiDark-Planck -body simulation, run using several hydrodynamic and semi-analytic codes. The G -S suite adds a state-of-the-art galaxy formation model based on the highly successful S simulation, mildly re-calibrated to match = 0 cluster stellar properties. Comparing to T T H zooms run with G -X, we find intrinsic differences in the evolution of the stellar and gas mass fractions, BCG ages, and galaxy colour-magnitude diagrams, with G -S generally providing a good match to available data at ≈ 0. G -S 's unique black hole growth and feedback model yields agreement with the observed BH scaling relations at the intermediate-mass range and predicts a slightly different slope at high masses where few observations currently lie. G -S provides a new and novel platform to elucidate the co-evolution of galaxies, gas, and black holes within the densest cosmic environments.

INTRODUCTION
Galaxy clusters are a key class of objects for many astrophysical areas. On large scales, they are useful for constraining cosmological models via their abundance and evolution. On halo scales, they are interesting sites for environmental studies of galaxies along with the evolution of the hot intracluster medium (ICM). On galactic scales, they are ★ E-mail: weiguang.cui@ed.ac.uk important for studying the oldest galaxy stellar populations and most massive galaxies, along with the impact of supermassive black holes on galaxies and surrounding gas. For these reasons, clusters are much investigated both observationally and theoretically (see Allen et al. 2011;Kravtsov & Borgani 2012;Walker et al. 2019, for reviews).
Interpreting observations of galaxy clusters from the radio to the X-rays within a structure formation context is challenging, because clusters contain numerous components interacting over a wide range of scales. Thus models must capture both the large-scale structure within which clusters grow, while including many small-scale physical processes. Cosmologically-situated numerical simulations have played an increasingly important role in holistically understanding the physics driving clusters. Unfortunately, clusters are rare objects, so representative cosmological volumes that are able to model all the relevant small-scale physics are extremely challenging computationally. Many studies have therefore focused on using the zoom simulation technique, where individual clusters are re-simulated with full galaxy formation physics after being extracted from a large (typically dark matter-only) parent simulation. Zoom simulations must be done one object at a time, but with a sufficiently large sample they can cover the full parameter space of real clusters.
Previous cluster zoom simulations with only dark matter particles, such as Phoenix (Gao et al. 2012), Rhapsody (Wu et al. 2013), and ZOMG (Borzyszkowski et al. 2017) can elucidate the detailed internal structures of the clusters, but cannot directly model the galaxies and gas. Hydrodynamic cluster zoom simulation suites such as Dianoga (Planelles et al. 2013), MACSIS (Barnes et al. 2017a), C-EAGLE (Barnes et al. 2017b), Hydrangea (Bahé et al. 2017), MUSIC (Sembolini et al. 2013), and nIFTy (Sembolini et al. 2016a), are able to investigate detailed baryonic properties, and to compare with observations more directly. These are complemented by hydrodynamic simulations that have representative volumes, typically focusing more on the group to poor cluster regime, such as EAGLE (Schaye et al. 2015), Magneticum (Dolag et al. 2016), BAHAMAS (McCarthy et al. 2017), IllustrisTNG (Pillepich et al. 2018), FABLE (Henden et al. 2018, which also included zoom regions for galaxy clusters), and S Robson & Davé 2020). Thus there is great interest in producing state-of-the-art simulations of clusters, particularly with clusters being a target for numerous forefront observational facilities such as Euclid, the Dark Energy Survey, eROSITA, Sunyeav-Zeldovich (SZ) telescopes, and the Square Kilometre Array and its precursors.
T T H project 1 occupies a unique niche among cluster simulation suites. Other suites of cluster simulations have typically focused on a handful of objects, with zoom regions covering only the cluster and immediate surroundings. In contrast, T T H re-simulates a mass-complete sample of 324 galaxy clusters extracted from the MultiDark cosmological simulation, using the zoom region that extends out to many virial radii. The penalty for having so many clusters with large zoom regions is that the numerical resolution is necessarily lower owing to computational limitations. However, the benefit is that it covers a relatively wide and complete halo mass range, enables larger-scale cosmic web studies around clusters, and provides good statistics along with the ability to investigate rare systems. Furthermore, another interesting feature of T T H project is that all these clusters have been run with several different galaxy evolution codes. These include the cosmological hydrodynamic codes G -MUSIC (Sembolini et al. 2013) and G -X (Rasia et al. 2015), as well as three different semi-analytical models (SAMs): G (Benson 2012), SAGE (Croton et al. 2016) and SAG (Cora et al. 2018). This enables crosscomparisons between models employing different input physics, to better understand the sensitivity to the various physical processes and the robustness of the resulting predictions.
In this paper, we introduce another set of hydrodynamic runs to T T H suite, namely the G -S runs 2 . This 1 https://the300-project.org/. the300 is also used for short. 2 Note that to distinguish from the S simulation ) which is a 100 ℎ −1 Mpc cosmological hydrodynamic simulation run with the same suite uses the G cosmological hydrodynamics code in its Meshless Finite Mass (MFM) solver mode, as opposed to the G based runs which use Smoothed Particle Hydrodynamics (SPH). It further includes a suite of galaxy formation physics similar to that in the recent S simulation ) that yields an excellent match to a wide range of galaxy, black hole, and intergalactic medium properties. Its novel input physics modules such as torquelimited black hole growth, stably bipolar jet feedback, and on-the-fly dust tracking, make it a valuable addition to the existing T T H suite. This paper is organised in the following order: we first introduce the T T H project in §2. Then, we present the details of the new G -S run in §3. The general comparisons to the other models and observation results are shown in §4. We also include the gas scaling relations in the appendixes which seem less affected. At last, we conclude and discuss our results in §5.

T T H PROJECT
T T H project (Cui et al. 2018b) is a set of clusterscale zoom simulations based on a mass-complete sample of 324 most massive galaxy clusters ( vir 8 × 10 14 ℎ −1 M ) 3 drawn from the MultiDark simulation (MDPL2, Klypin et al. 2016). MDPL2 assumes cosmological parameters from Planck (Planck Collaboration et al. 2016), and has a periodic cube of comoving length 1 ℎ −1 Gpc containing 3840 3 DM particles having a mass of 1.5 × 10 9 ℎ −1 M each. Each cluster region was selected to have a comoving radius of 15 ℎ −1 Mpc (over 5 × 200 ) for re-simulation with different baryonic models: G -MUSIC (Sembolini et al. 2013), G -X (Rasia et al. 2015; Beck et al. 2016), and now G -S (this work). Additionally, galaxy catalogues in the same cluster regions are extracted from three different SAMs that were run on MDPL2 (Knebe et al. 2018): SAG (Cora et al. 2018), SAGE (Croton et al. 2016), and G (Benson 2012). The re-simulation regions are generated with the parallel GINNUNGAGAP code 4 : the highest resolution Lagrangian regions share the same mass resolution as the original MDPL2 simulation with gas particles ( gas = 2.36 × 10 8 ℎ −1 M ) split from DM particles. The outside regions are degraded in multiple layers (with a shell thickness of ∼ 4 ℎ −1 Mpc) with lower mass resolution particles (mass increased by eight times for each layer) that eventually provide the same tidal fields at a much lower computational costs than in the original simulation.
The aforementioned unique features of T T H project has enabled many studies of various aspects of galaxy clusters to be carried out. To date, T T H has been used in over 30 papers investigating galaxy clusters and their environs. These include studying the detailed relationship between the central cluster and connecting filaments (Rost et al. 2021;Kuchner et al. 2020Kuchner et al. , 2021, the feeding of the galaxy clusters (Kuchner et al. 2022;Kotecha et al. 2022), cluster backsplash galaxies (Haggar et al. 2020;Knebe et al. 2020), and the virial shock radius (Baxter et al. 2021; code, this run for T T H clusters is referred to as G -S run. 3 The halo mass is defined as the mass enclosed inside an overdensity of times the critical density of the universe: =∼ 98 for virial mass at = 0 (Bryan & Norman 1998) and 200,500 is with = 200, 500 respectively. Similarly, 500 is the radius at which the overdensity = 500 is reached. 4 https://github.com/ginnungagapgroup/ginnungagap T T H : G -S 3 Anbajagane et al. 2022a). The advanced input physics in the hydrodynamic simulations further allow detailed investigations on cluster properties, such as cluster profiles (Mostoghiu et al. 2019;Li et al. 2020), substructure and baryon content (Arthur et al. 2019;Haggar et al. 2021;Mostoghiu et al. 2021b,a), dynamical state and morphologies (Capalbo et al. 2021;De Luca et al. 2021;Capalbo et al. 2022;Zhang et al. 2021), ICM (non-)thermalization (Sayers et al. 2021;Sereno et al. 2021), the fundamental plane (Díaz-García et al. 2022), the effects of mergers on the BCG properties (Contreras-Santos et al. 2022), and various methods for estimating galaxy cluster masses, namely dynamics (Ansarifard et al. 2020;Li et al. 2021Li et al. , 2022, hydrostatic equilibrium (Gianfagna et al. 2022), and machine learning (de Andres et al. 2022b,a). Lastly, comparing to the void/field region runs in this project allows us to study the effect of environment (Wang et al. 2018); a self-interacting dark matter run was done that allows constraints on the dark matter cross-section ; and even chameleon gravity was examined (Tamosiunas et al. 2022). With many more projects in the works, it is valuable to continue to update T T H runs with state-of-the-art physical models in order to expand its range and robustness.

The S model
The G -S runs of T T H are performed with the G code (Hopkins 2015) with the state-of-the-art galaxy formation subgrid models following the S simulation . We refer the interested reader to Davé et al. (2019) for full details of all of S 's features, and here focus on its more unique aspects relevant for clusters. We also describe our modifications to the S model parameters utilised for the G -S runs of T T H clusters, which required re-tuning owing to the lower numerical resolution in G -S relative to the original S simulation. Owing to the Meshless Finite Mass (MFM) solver implemented in G , gas particles are evolved following an accurate description of shocks and shear flows, without the need for any artificial viscosity. This feature improves the description of shocks and flows with high Mach number, which provides a realistic simulation of outflows and jets. It also provides improved handling of contact discontinuities relative to SPH. See Hopkins (2015) for a full discussion of the differences of MFM with respect to other hydrodynamics methods.
Radiative cooling and photon-heating/ionization processes of gas are implemented using the Grackle-3.1 library (Smith et al. 2017), which also accounts for metal cooling with non-equilibrium primordial chemistry treatment. SIMBA adopts a spatially-uniform (Haardt & Madau 2012) ultraviolet background model, accounting for selfshielding on the fly based on the prescription in Rahmati et al. (2013). An H 2 -based star formation model is taken from its predecessor simulation M (Davé et al. 2016), which is calibrated to match the Schmidt (1959) law. Here H 2 is estimated from the local column density and metallicity following the Krumholz et al. (2009);Krumholz & Gnedin (2011) prescription. Besides requiring that 2 be present, an additional minimum density cut of H > 0.1 cm −3 and a minimum metallicity of 0.05 for the H 2 in its formation (Krumholz et al. 2009), compared to the original S simulation of 0.1 cm −3 and 0.01, respectively, are required for active star formation.
Star formation-driven galactic winds also shares the same decoupled two-phase model in M , but the mass loading factor scaling with stellar mass is based on the Feedback in Realistic Environments (FIRE) zoom simulations of Anglés-Alcázar et al. (2017b): Here, 0 = 2 × 10 9 M ; this is slightly different than the original S simulation, as we will motivate later. The ejection velocity is based on scalings from Muratov et al. (2015) as in M : where Δ (0.25 vir ) is the velocity corresponding to the potential difference between the launch point and one-quarter of the virial radius. Again, this has been changed relative to S , who used a normalisation of 1.7; the normalisation used here in G -S is the original one proposed by Muratov et al. (2015). Identically to S , galaxies are identified with the on-the-fly approximate friends-of-friends (FOF) finder for star, dense gas and BH particles in Davé et al. (2016), allowing galaxy properties such as * to be computed on-the-fly, with circ obtained from a scaling based on the observed baryonic Tully-Fisher relation (McGaugh 2012).
The chemical enrichment model tracks eleven elements (H, He, C, N, O, Ne, Mg, Si, S, Ca, Fe), with metals from supernovae type Ia (Iwamoto et al. 1999) and type II (Nomoto et al. 2006), and Asymptotic Giant Branch (AGB) stars (Oppenheimer & Davé 2006). Furthermore, S also includes metal-loaded winds, i.e. metals in the wind particle are enhanced, and correspondingly subtracted from nearby gas in a kernel-weighted manner. All this is identical to the original S model; see Davé et al. (2019) for details. S seeds black hole (BH) particles based on the host galaxy stellar mass, * > BH × seed . If the galaxy meets the aforementioned condition and does not already contain a black hole particle, then the star particle closest to the centre of mass of the galaxy is converted into a black hole particle. For this G -S run, we employ seed = 10 5 ℎ −1 M and BH = 3 × 10 5 , which sets the galaxy stellar mass threshold for seeding BH is * ≈ 10 10.5 ℎ −1 M . This is 10 times higher than original S model, owing to the lower resolution. By assuming the dynamical friction is efficient enough to maintain black holes near the host galaxy's centre (within 4 times the size of the BH kernel, 0 , considered for the accretion model), black hole particles are re-positioned to the location of the potential minimum within the FOF host group at each time-step. Furthermore, any two black holes located within 0 are allowed to merge instantaneously if their relative velocity is lower than three times their mutual escape velocity.
The BH accretion follows a dual model: The cold accretion mode is described with a torque-limited accretion model for the cold gas ( ≤ 10 5 ), driven by disk gravitational instabilities arising from galactic scales down to the accretion disk around the central BH (Hopkins & Quataert 2011; see also Anglés-Alcázar et al. 2013. Hot gas ( > 10 5 ) is accreted based on the Bondi rate (Bondi 1952). We reduce the Bondi accretion rate to the rate for a BH = 10 9 ℎ −1 M BH, no matter how big it gets; in S , this was set to 10 10 . The black hole accretion kernel has a distance enclosing 256 baryonic particles or 0 = 6 ℎ −1 kpc (comoving), whichever is smaller, within which the gas quantities are calculated; the latter was set to 2 ℎ −1 kpc in S . The total accretion rate for a given black hole is then the sum of Torque and Bondi , times an additional constant 1 − , with the radiative efficiency assumed to be = 0.1 ). This total accretion rate is used to determine the AGN feedback modes detailed in the following paragraph.
There are three different AGN feedback modes: A kinetic subgrid model for both the 'radiative mode' and 'jet mode', and a mostly kinetic X-ray feedback mode accounting for radiation pressure from X-rays off the accretion disc broadly following Choi et al. (2012). The 'radiative mode' feedback is turned on when the BH is accreting at a high Eddington ratio ( Edd ≡ BH / Edd > 0.2). The radiative velocity for wind particles, which is based on ionized gas linewidth observations of X-ray detected AGN from SDSS (see Fig. 8 in Perna et al. 2017), scales as: , Radiative = 500 + 500 3 log 10 BH M + 6 km s −1 .
When the BH's in low Eddington accretion mode, Edd < 0.2, the wind begins to transition into a jet mode, with the velocity scaled with Edd as follows: , Jet = , Radiative + 15000 log 10 0.2 Edd km s −1 .
Note that the wind in both modes is ejected in the form of purely bipolar outflows, based on the angular momentum of gas and stars within 0 . The wind velocity in the jet mode is capped at 15000 km/s (as opposed to 7000 km/s in original S ) when the Eddington rate drops below Edd ≤ 0.02. There is another condition to trigger jet mode -the minimum BH mass has to be greater than 10 7.5 M . X-ray feedback only operates when the 'jet mode' AGN feedback is in action, and it further has to meet another two conditions: * > 10 9 M and gas / baryon < 0.2. The * condition is raised with respect to the original S model, owing to the lower resolution. Another novel feature of the S simulation, the on-the-fly dust production and destruction model, is also employed in G -S , unchanged. This model is quite successful in reproducing galaxy dust properties over cosmic time .

Re-calibration for T T H
We now summarise the modifications from the original S simulation, and describe the datasets used for the re-calibration. The set of parameter changes is tabulated in Table 1. To reiterate, this is required because T T H simulation has about 10 times worse mass resolution than the S simulation. It has been found that the S model is reasonably well converged towards higher resolutions than in the original S simulation, but it has so-called weak convergence (Schaye et al. 2015) towards poorer resolution, and must be re-calibrated to achieve an equivalently good match to observations.
We use three observational relations focusing on the stellar component in galaxy clusters to calibrate the model parameters for our G -S runs, all at ≈ 0: (C1) total stellar mass fraction within 500 ; (C2) BCG stellar mass -halo mass relation; and (C3) the satellite galaxy stellar mass function (SSMF) in galaxy clusters. The denotations C1, C2, and C3 are further noted in the three subsection titles in §4 where we show comparisons to some calibration datasets. We note that the calibration was done mostly by trial-anderror using intuition to guide the variations, so should not be regarded as a unique parameter set that achieves agreement.
The calibration was not done over the entire sample, but rather only using a single cluster region at = 0 whose largest object has 500 ≈ 5 × 10 14 , chosen to be typical of T T H sample 5 . Additionally, all uncontaminated halos (i.e. those that do not contain any low-resolution particles) within that region are considered (which includes about 10 more halos down to ≈ 10 13 ), although typically these do not add much information since the observational constraints tend to get weaker towards lower masses. Once sufficiently calibrated based on this single region, all the parameters in Table 1 were frozen and run for all the remaining galaxy clusters. As such, while the agreements to the datasets used for C1, C2, and C3 are to be considered as tuned and not an intrinsic success of the model, it is also the case that in principle other cluster regions could have shown large variations with respect to the one region used for tuning, particularly at different masses. Hence the agreements with C1, C2, and C3 over the entire mass range (and to higher redshifts, where applicable) may still be regarded as a modest success.
We now describe the motivation for the individual changes, namely what went wrong with the cluster run when we adopted the original S parameters. At this low resolution, we found that the star formation is insufficient within satellite galaxies, resulting in a much lower satellite stellar mass function (SSMF) compared to observations (C3). The changes for H , Metallicity floor, 0 in SN mass loading factor and galaxy stellar mass limit for X-ray feedback serve to boost the SF in these satellite galaxies. However, as a consequence of these changes, the total stellar mass and BCG stellar mass ended up significantly higher compared to observations (C1 and C2). Thus we strengthened the AGN jet feedback, which has been shown as the key to quench massive galaxies , by increasing the maximum jet speed; this reduces the stellar mass of BCG (and hence total). Finally, we choose to use a fixed comoving softening length of 5 ℎ −1 kpc for consistency with the other T T H runs, as opposed to a variable softening length with a minimum of 0.5 ℎ −1 kpc as in original S . The BH accretion kernel maximum radius was commensurately increased in G -S owing to the coarser spatial resolution, and the BH seeding stellar mass was increased by an order of magnitude to reflect the order of magnitude poorer particle mass resolution.

Detailed differences between G -X and G -S
As noted in the Introduction, we will focus on the differences between the results from G -X and G -S runs, given that they are the two hydrodynamic models that match a reasonably wide suite of observations. Though both simulation codes are well presented in a series of literature papers, it is useful to relist their key features again in Table 2 for comparison. For more detailed models of G -X, we refer to Cui et al. (2018b) and references therein. Although we give descriptions of the S model in this section 3, interested readers are referred to Davé et al. (2019) and Davé et al. (2016) for further information. Besides these model differences as shown in Table 2, we emphasise that G -X is quite successful in reproducing the observed gas properties and relations, while the original S simulation was primarily tuned to reproduce galaxy stellar properties. Here, we follow S 's process by calibrating G -S according to the observed stellar properties, with no regard to gas properties. Comparing the two runs in different properties will help us to better constrain galaxy formation models. the mean or median values to rerun the whole cluster data set. However, this requires an infeasible amount of computation time. Furthermore, the observation results have a large scatter, and the sub-grid models are not expected to perfectly match all calibration data. Thus, we only require the calibrated parameters yield results in rough agreement with observations. 0.001 0.05 0 in SN mass loading factor 5.2 × 10 9 M 2 × 10 9 M Galaxy stellar mass limit for seeding BH 10 9.5 ℎ −1 M 10 10.5 ℎ −1 M BH mass for Bondi accretion rate cap 10 10 ℎ −1 M 10 9 ℎ −1 M BH accretion kernel radius 2 ℎ −1 kpc 6 ℎ −1 kpc Cap wind velocity limit in jet mode 7000 km/s 15000 km/s Galaxy stellar mass limit for X-ray feedback 0 10 9 M Gravitational softening length 0.5 ℎ −1 kpc minimum 5 ℎ −1 kpc fixed Table 1. Summary of parameter changes with respect to the original S simulation.  (Springel & Hernquist 2003) Two-phase winds with mass loading factor depending on * baryon models -BH BH seeding condition > 8 × 10 11 ℎ −1 M & * > 1.6 × 10 10 ℎ −1 M Galaxy stellar mass > 3 × 10 10 ℎ −1 M BH seed mass 5 × 10 6 ℎ −1 M 10 5 ℎ −1 M BH accretion Bondi accretion torque-limited and Bondi accretion models BH feedback Thermal feedback Kinetic feedback + X-ray feedback Number of stars that per gas particle can spawn. Note that the detailed implementations are different, see Tornatore et al. (2007) and Davé et al. (2016) for more information. Note here * is the total stellar mass within the FoF halo. And another two conditions have to be met as well: * > 0.05 & > 0.10 * The feedback energy is coming from both mechanical and radiative modes. Different outflow velocities are used to mimic the jet and radiative mode AGNs in observation. Thermal or thermal+kinetic feedback is used for the X-ray heating depending on whether the surrounding gas is non-ISM or ISM, respectively.

The AHF halo and C galaxy catalogues
Two object catalogues are generated from the suite of the300 galaxy cluster runs: the AHF 6 (Knollmann & Knebe 2009) halo catalogues and the C 7 galaxy/FoF halo catalogues. AHF provides halo and subhalo (and thus galaxy) catalogue generated using a spherical overdensity (SO) algorithm, while C also provide halo catalogue generated using a 3-D FoF algorithm along with a matched galaxy catalog using a 6-D FoF. C further provides a large range of pre-computed physical and photometric properties (with and without dust extinction) for each object. 6 http://popia.ft.uam.es/AHF 7 https://github.com/dnarayanan/caesar In this paper, we use as many (uncontaminated) objects as possible (if not specified) to do the investigation because T T H runs have a much larger radius, thus many smaller mass halos besides the central cluster in each region. If not specified, we always use the halo mass defined as 500 with the quantities are always calculated within 500 . Therefore, halo properties from the AHF catalogue are used to compare with the global clusters properties from observations, while galaxy properties from C which has a 6D (in both spatial and velocity field) galaxy finder, is used. We further match the galaxies from C to the halos from AHF, by simply taking all the C galaxies within the AHF halo radius ( 500 ). The BCG is selected as the most massive C galaxy that lies close to the AHF distinct halo centre. As the 6D galaxy finder makes no distinction between central and satellite galaxies, it doesn't matter whether the galaxy is coming from a FoF halo or a SO halo.
To track cluster growth histories, we use the cluster main progenitors determined by the MERGERTREE package integrated into the AHF program. The main progenitors are selected based on the matched dark matter particle IDs.

Baryon fractions
We first focus on the total gas and stellar components within the clusters (mostly within 500 ) in this subsection. We will also detail their evolution and show the differences between G -X and G -S runs.

C1: The gas and stellar fractions
The first quantity we examine in G -S clusters is the abundance of stars and gas. This global quantity provide the overall measure of the cluster's baryonic content that can be directly compared to observational results (see Oppenheimer et al. 2021, for a recent review). As have been revealed by both observation and theoretical works (e.g. Behroozi et al. 2013;Yang et al. 2013), more gas is consumed and converted into stars (i.e. star formation efficiency is higher) in less massive halos through the group and cluster regime. Therefore, low-mass halos also tend to have somewhat higher stellar mass fractions. To achieve this general trend, AGN feedback is invoked in all current models (Somerville & Davé 2015) in order to solve the cooling flow problem (e.g. Fabian 1994) in which too many stars are formed, especially in the BCG, due to gas cooling being very efficient in the centres of galaxy clusters (Kravtsov & Borgani 2012). In models, AGN feedback from the BH of the central galaxy is invoked to counteract cooling and/or expel gas, thereby quenching the galaxy. Using this, hydrodynamic simulations of galaxy clusters can roughly reproduce the correct stellar mass fraction as a function of halo mass.
In the case of G -S , we have used the stellar fraction to constrain our baryon parameters (C1). But since the calibration was only done for a single object, it is still interesting to examine this relation over all the G -S cluster regions. The gas fractions were not used for calibration, so they represent an independent prediction.
In Figure 1 we present comparisons of both stellar and gas fractions within 500 between our simulated clusters at = 0, versus recent simulations (FABLE, Henden et al. 2018 and C-EAGLE, Barnes et al. 2017b) and observational data at 0.1 Laganá et al. 2011;Sanderson et al. 2013;Gonzalez et al. 2013;Kravtsov et al. 2018;Chiu et al. 2018). The gas fraction is calculated using all gas particles, but as indicated in Li et al. (2020), the cold gas only contributes a very small faction to the total gas mass, hence it is reasonable to compare to the results from observations which mainly use hot gas. Statistical or best-fit results are presented for Lovisari et al. instead of X-ray) that also examined clusters at ∼ 0.1 to investigate similar fractions, are not included in this comparison due to various reasons, such as no 500 or * being available, or that the sample is dominated by 0 clusters. However, it has been suggested that there is almost no redshift evolution in both gas and stellar fractions within ∼ 1 (see more discussion below), so we include the best fitting results from the recent work of Akino et al. (2021) out to modest redshifts using weak-lensing masses (Umetsu et al. 2020) for the 136 XXL clusters in the HSC-SSP survey, shown as the grey shaded band.
As seen by comparing the gas fractions in the left panel of Figure 1, among the simulations G -X, FABLE and C-EAGLE are in very good agreement with each other for both the median and scatters. In contrast, G -S shows a steeper slope, steepening further for 500 < ∼ 10 14 . A similarly steep trend was found within the group regime in the original S simulation (Robson & Davé 2020). The difference could owe to the fact that the first three simulations employ a thermally-based AGN feedback scheme, while G -S employs a kinetic scheme; we leave a detailed exploration into the origin of such differences for future work.
Turning to the gas fraction observations, both individual clusters and sample fits show a wide range of gas fractions particularly towards lower masses. All the simulation predictions are contained within the observational scatter, but the scatter in any given model is much smaller than in the data. The fitting function from Lovisari et al. (2015) is the flattest among observations, and agrees well with G -X results. Meanwhile, Dietrich et al. (2019) and Akino et al. (2021) are in better agreement with G -S at 500 > 10 14 M ; measurements to even lower masses are as yet highly uncertain. The downturn for 500 < 10 14 * in G -S owes to the high jet velocity which blows gas particles well outside of these low-mass halos' virial radius (Sorini et al. 2021). While these observations (and to a lesser extent the simulations) assume slightly different cosmological parameters, most are broadly consistent with a Planck cosmology (Planck Collaboration et al. 2016). Given the complexity in deriving these observed gas masses and 500 , it is not obvious how to make a correction for cosmology, so we take the data as-is and simply note that the predicted and observed ranges are likely to be much larger than differences due to cosmology.
The right-hand panel of Figure 1 shows the stellar mass fractions versus 500 . Here, G -X and G -S show very similar results, and IllustrisTNG (Pillepich et al. 2018) shows a similar slope but a slightly higher amplitude that is in better agreement with the observations from Akino et al. (2021). In contrast, FABLE shows higher stellar fractions at group scales ( 500 < 10 14 M ), and thus a steeper slope that matches better with the individual observations shown (also the results in Andreon 2010). We note that the low resolution of T T H will tend to suppress stellar fractions since they do not resolve as far down the mass function as Illus-trisTNG; if the shape of the galaxy stellar mass function is not a strong function of 500 , this can explain the constant ∼ 20% offset between these models. Overall, we note that where the observations are most robust at 500 > ∼ 10 14.5 M , all models are within the range of the observations. The scatter in the simulations is lower than in observations, which may reflect observational uncertainties in addition to intrinsic scatter. Finally, we note that Anbajagane et al. (2020) compared several cosmological simulations -BAHAMAS, Magneticum and TNG -and found that TNG has the lowest total stellar mass within 200 . Thus, we suspect that the stellar mass fractions from BAHAMAS and Magneticum would also lie above the TNG line in the right panel of Figure 1.
In conclusion, for massive (∼ 10 15 ) clusters all simulations are generally in agreement with each other for both gas fraction (differences within ∼5 per cent) and stellar fractions (difference within ∼ 1 per cent), and in reasonable agreement with observations. For G -S , the latter was mostly achieved via tuning parameters as described in §3. The models tend to differ more significantly to- Note that all the fitting results only cover the region of the observed data points. These two plots show that G -X is very similar to the FABLE and the C-EAGLE simulations in both fractions; there is little difference between G -X and G -S in the stellar fraction, but the gas fraction from G -S shows a much steeper slope. Both fractions from the observational data present a large scatter.
wards lower mass halos, with G -S producing particularly low gas fractions in groups, and lower stellar fractions compared to observations of individual 10 14 systems. The origin of the differences between G -X and G -S are being investigated by studying their density profiles in Li et al. (in preparation). Observationally, deeper and more precise estimation of the gas fraction from next-generation surveys such as NIKA2 (Adam et al. 2018), CMB-S4 (Abazajian et al. 2016) using the Sunyaev-Zeldovich effect 8 and ATHENA (Nandra et al. 2013) and Lynx (The Lynx Team 2018) in the X-rays, will be required to test the input galaxy formation physics. In the meantime, we can look into other properties to distinguish between these simulations.

The evolution of the baryon fractions at the same halo mass
It has been suggested that cluster baryon fractions (both gas and stellar) depend weakly on redshift, but strongly on cluster mass (see Chiu et al. 2018from observation side, or Planelles et al. 2013Truong et al. 2018 from simulations). At first glance this seems paradoxical within hierarchical structure formation models: how can the baryonic scaling relations vary steeply with mass and still remain 8 Interested readers are referred to  on how the nextgeneration SZ observations can be used to distinguish different baryon models. roughly constant in time as structures grow? In ΛCDM, big halos are mostly formed later by merging with smaller halos. If we ignore the baryon processes and halo accretion at low redshift, given that the smaller (group-sized) halos have low gas fractions and high stellar fractions, the later-formed massive (cluster) halos should have even lower gas fractions and higher stellar fractions, opposite to what is observed. In this and the following section, we will study the redshift evolution of the galaxy clusters in two ways, by binning at the same halo mass at all redshifts and by tracking individual halos. The former, discussed in this section, highlights how baryon fractions change for a sample selected at a given mass, while the latter, discussed in the next section explicitly shows the true evolution of baryon fractions as halos grow hierarchically.
For the gas fraction evolution shown in the left panel of Figure 2, at all redshifts the clusters show increasing gas fractions with mass, approaching but not reaching the full expected baryonic budget (horizontal dashed line). There is modest evolution, with the highest gas fractions at high redshifts. Higher gas fractions are expected at early epochs when cooling is rapid and feedback processes are dominated by star formation whose energetics are typically not sufficient to unbind gas from protoclusters.
Comparing between G -X and G -S , the latter has significantly more evolution at lower masses, and overall shows lower gas fractions. In Li et al. (2022, in prep.) we identify that the gas fraction difference between G -X and G -S owes to G -X tending to have a much higher gas density in the halo To clearly view the fraction changes along redshift, we use linear scale for the y-axes. The star symbol with dotted lines are for G -X and square with solid lines are for G -S . Different colourful lines are used to highlight the redshift evolution, see the legend on the left-hand-side panel for details. We don't include any observation data (especially at high redshift) in this busy plot because there lacks statistical results and weak or no redshift evolution is claimed. It seems that G -S shows much clear redshift evolution for both gas (much clearer at lower halo mass end) and stellar fractions compared to G -X, of which shows a weak redshift evolution in their stellar mass fractions. centre than G -S . This points to galactic feedback processes being the primary driver of the model differences. This is corroborated by the most significant drop occurring between = 2.5 → 1.5, which is when Simba's jet mode AGN feedback which becomes important prominent in this mass range (Robson & Davé 2021). The gas fraction continues to drop more rapidly in G -S that in G -X down to = 0 in 500 < ∼ 10 14.5 halos. Although for clarity we don't show higher-redshift observation data in Figure 2, G -S has the closest gas fractions compared to Chiu et al. (2016a), around 10 per cent at ≈ 0.9 and 500 ≈ 6 × 10 14 M , while the gas fraction from Chiu et al. (2018) (around 12 per cent at ≈ 0.6 and 500 ≈ 4.8 × 10 14 M ) lies between G -X and G -S . However, Chiu et al. (2016a) suggests no statistically significant redshift trend at fixed mass (see also Chiu et al. 2018;Bulbul et al. 2019) when accounting for a 15 per cent systematic mass uncertainty. G -S predicts some evolution over this redshift range, but at a < ∼ 1% level which is much smaller than the uncertainties in Chiu et al. (2018)). This is also suggested by Henden et al. (2020) and the redshift evolution is even stronger at poor group masses of 0.5 − 1 × 10 14 M .
The stellar fraction-halo mass relation at various redshifts is shown in the right figure of Figure 2. Both simulations predict a dropping stellar fraction with mass for 500 < ∼ 10 14.5 , and by = 0 they are fairly similar. However, the evolution is quite different, with again G -S showing much more evolution than G -X. Neither simulation evolves much at the massive end; we reiterate that G -S was calibrated to match = 0 observations in this mass range.
The rapid drop in stellar fraction indicates that galaxies in these systems tend to form their stars very early on, so that over time the halo mass grows but the stellar content does not keep up. This can happen because star formation is quenched early on, and also if the hierarchical growth is predominantly adding smaller systems that have lower stellar fractions (well below group scales).
For massive cluster halos, the stellar mass should mainly come from the accreting small halos, as the central galaxy is typically quenched by ∼ 3 (see the following subsection and §4.2.4 for details). This results in mild but still visible redshift evolutionabout 0.5 per cent from = 1.5 → 0.
In smaller clusters and groups, jet feedback in G -S is once again implicated by the fact that the most significant change in stellar fraction happens between = 2.5 and = 1.5. We have argued that jet feedback is responsible for dropping the gas fraction, which then reduces the fuel for star formation and thus causes galaxy quenching that results a decrease in the stellar fraction. From = 1 → 0 the evolution is relatively modest. However, it is still significantly more than in G -X, which predicts essentially no evolution in stellar mass fractions since ∼ 1.5. This suggests that the stellar fractions within protocluster environments at high redshifts provides a significant discriminant between models.
Although we don't include observations of higher-stellar fraction determinations on this plot, we comment on some other results below in relation to our predictions. On the simulation side, the mild evolution is in agreement with Henden et al. (2020) who also report a marginally significant change for 500 3 × 10 14 M up to ∼ 1.  2021) report higher * ≈ 0.025 for 12 clusters at = 0.95 − 1.43, which seems in good agreement with G -S . Besides these results, most observations claim there is no clear redshift evolution of the cluster stellar fraction, (see Lin et al. 2012Lin et al. , 2017Chiu et al. 2018, for example). This likely owes to small sample sizes, the difficulty in measuring stellar fractions in distant systems, and the larger uncertainty in mass estimation. Therefore, such a very small fraction change predicted in simulations can be difficult to detect. Note that using 500 means we don't account for halo pseudo-evolution, but this is the case for both observations and simulations so this should not have an effect on our comparisons. Another point to keep in mind is that the high redshift halos in zoom simulations are necessarily the progenitors of the = 0 halos, unlike in observations. Therefore, the redshift evolution can appear more significant than comparing to random samples from a cosmological volume.
In conclusion, we find that there is a decrease in both gas and stellar fractions with time at lower halo mass ( 500 10 14.5 ℎ −1 M ) in both G -X and G -S . However, the rate of decrease is much stronger in G -S , which shows much higher stellar fractions at high-and much lower gas fractions at low-. This differences can be expected from G -S 's kinetic AGN feedback which imparts much stronger feedback than G -X's thermal AGN feedback. For high mass halos, the AGN feedback becomes relatively unimportant, and both models predict similar gas and stellar fractions at all redshifts (see Eckert et al. 2021, for a recent review on the AGN feedback on galaxy groups). Thus the greatest discrimination between models occurs at high redshifts for poor clusters and groups, which motivates future X-ray and SZ surveys that can probe this regime.

The evolution of the cluster baryon fractions by tracking
Simulations provide a view that cannot be seen in observations, in that they can track the evolution of individual halos over time. In comparing baryon fractions at a fixed halo mass over time, we saw a significant change at ∼ 1.5 − 2.5 in G -S . To see how this is reflected in individual halos, we can place clusters into three mass bins at = 0, and then track the individual systems back to the earliest epochs, showing how a specific set of clusters evolves. We track each cluster's main progenitor using the AHF MERGERTREE catalogue. We show their median baryon fraction values as a function of time (and redshift) in Figure 3, for our three chosen mass bins as indicated in the legend, for G -S (green curves) and G -X (blue). The gas fraction evolution is shown in the left panel of Figure 3. G -X and G -S show qualitatively similar evolution, but quantitatively there are differences in the values as well as shifts in the maxima and minima of the evolution. At very high redshift ( ≈ 5−6), G -X and G -S reach a gas fraction near the cosmic baryon fraction value, with G -S 's occurring slightly earlier. The lower values prior to this likely owe to the halos not being properly resolved at these very early epochs. From then until ∼ 2, gas fractions starts to decline very quickly (more strongly in G -S ), until they reach a minimum at ∼ 2 − 3 (slightly later in G -S ). At these redshifts, the feedback models are relatively simple, since there is little AGN feedback and wind recycling is not yet common. Hence differences between G -S and G -X are caused by the different star formation and stellar feedback models. The net result is that the gas fraction in G -S is only slightly lower (a few per cent) than G -X, persisting until today, indicating that the sensitivity to feedback for halo gas contents is established early on. After ∼ 2, gas fractions start to increase gradually, levelling off after ∼ 1 at approximately their present-day value.
The middle panel shows the stellar fraction evolution. Here, there is sharp early rise as rapid cooling in the dense early universe is able to drive copious gas to the halo centres. After ∼ 4, the stellar fractions start to drop, and level off after ∼ 1. G -S has higher stellar fractions at early times versus G -X; at = 4, G -S has produce ∼ 50 per cent more stars than G -X in these halos, but by = 0 they are the same. This shows that G -S has earlier stellar formation times in groups and clusters compared to G -X. G -S grows early galaxies faster than G -X. Comparing these first two panels, it is clear that the sharp early drop in gas content accompanies the sharp increase in stellar content. This implicates rapid gas consumption owing to star formation as the main driver of these early trends. At later epochs, the lower gas contents in G -S limit the available fuel for star formation, resulting in less rapid stellar growth, and the stellar fraction returns to meet G -X by = 0. Comparing the three mass bins, until 2 Gyrs there is little difference among the mass bins in either run, in either gas or stellar fractions. This shows that at early times, galaxies within these regions grow self-similarly along a linear SFR- * relation, so there is no significant trend with halo mass. But as time goes on, the galaxies in more massive halos break away from this as they quench owing to AGN feedback. These differences are most apparent in G -S , where the low-mass halos have lower gas fractions because AGN feedback can more easily remove baryons. Interestingly, the stellar fractions show much less trends with halo mass. This is because once these halos quench at ∼ 2 − 3, the stellar content doesn't grow much, while the halo mass continues to grow.
The growth of the total and baryon halo mass is illustrated in the right panel of Figure 3. As expected, the total halo masses show no differences in evolution between the two runs, as this is driven by hierarchical structure formation. There is about 2 orders of magnitude growth of the total halo mass in the first 2 Gyrs, while the rest of the growth in halo mass (also about 2 orders of magnitude) takes the rest 10 Gyrs. The total baryon mass shows more differences, owing to G -S having lower gas fractions. The differences in the gas fractions (and associated slight differences in the stellar fractions) emerging at ∼ 2−3 can be understood from an interplay between halo growth and star formation quenching. As shown in Cui et al. (2021); Robson & Davé (2021), it is at these redshifts when G -S 's jet mode AGN feedback turns on in these massive halos, which is what drives galaxy quenching . With star formation curtailed and the halo having a stable virial shock (Kereš et al. 2005;Dekel & Birnboim 2006), accreted gas can be held up and accumulated in the halo, causing the gas fractions to increase and the stellar fractions to drop. G -S and G -X have their AGN feedback in operation at a similar time, which is what is required to produce today's quenched galaxy population. Therefore, both models present a qualitatively similar picture, although the quantitative details of their evolutions differ.
In summary, by tracking individual halos within mass bins, we see that the differences in stellar and gas fractions between G -X and G -S are established at fairly early epochs. The gas fractions, remain are significantly different at all redshifts with G -S having low gas fractions at about 98 per cent of the cosmic baryon fraction. G -S also shows a stronger dependence with halo mass owing to the ability of its AGN jet feedback model to evacuate  gas particularly from smaller halos. Meanwhile, G -S has an earlier stellar growth phase relative to G -X, with ∼ 50% more stars formed than in G -X at ∼ 3 − 5. This is caused by G -S consuming more gas into stars. After this, the stellar fractions drop in both models, as AGN feedback quenches galaxy stellar growth while the halo mass continues to grow. The lower gas fractions in G -S curtail stellar growth more than in G -X, so that by < ∼ 1 the stellar fractions in the two models are similar. The overall differences in the halo baryon masses are small but noticeable by = 0.
While both models produce a similar global stellar fraction by = 0, they arrive at this by different histories. Hence next we examine BCG evolution in more detail to understand how such differences manifest in their stellar properties.

BCG Properties
We now look into the key stellar component in the galaxy clusters: the brightest cluster galaxy. We will study several relations for BCGs, from their black holes to their halo properties. We further look into the BCG colour-magnitude and age-halo formation time relations. Note that various definitions of BCGs are used in this section at times, to better compare with observational results.

C2: The
, *ℎ relation As we have discussed in section 3, the BCG -halo mass relation is one of the three constraints for the G -S parameter calibration. One well-known issue for the BCGs in galaxy clusters is that they intermingle with the intra-cluster light (ICL) in their outskirts. It is not trivial to separate a BCG's stars from its surrounding ICL in hydrodynamic simulations (see Murante et al. 2007;Dolag et al. 2010;Cañas et al. 2020, for examples). A new advanced method has been proposed in Cañas et al. 2022 (in prep.), which is robust and independent of simulation resolution. However, the BCG identified with these more physical motivated methods is not directly comparable with the BCG as identified in observations (see Cui et al. 2014;Rudick et al. 2011, for more discussions). Therefore, it is not easy to make apples-to-apples comparisons. Here, we adopt a simple definition of the enclosed mass within different 3-D apertures, to be able to compare with the observation results from Kravtsov et al. (2018). Note here that the observed masses are within projected 2D apertures, so the simulations are expected to have a slightly lower BCG mass compared to the observed one.
In Figure 4, we show the BCG stellar mass-500 relations from all the uncontaminated halos in our 324 regions at = 0. Here, the BCG stellar masses are estimated at three different 3D apertures: 30 kpc (left), 50 kpc (middle) and 0.1 × 500 (right). G -S results are in green, G -X in blue, and observations from Kravtsov et al. (2018) are shown as black stars. Errorbars enclose 16-84% of the simulated values.
Overall, there is broad agreement with the Kravtsov et al. (2018) data for both models, with a mild increase in BCG mass as a function of 500 . The general trends are very similar regardless of the chosen aperture. However, in detail there are some differences. Generally, G -S produces somewhat lower BCG masses, particularly at 500 10 14 . Also, G -X's BCG masses show a clearly larger dispersion at a fixed halo mass bin. At all apertures, the BCG masses from G -S are in somewhat better agreement with the Kravtsov et al. (2018) data. To quantify aperture effects for 3D versus 2D, we have checked the projection effect on the BCG masses within different apertures and find that the BCG mass increases are mostly less than 20 per cent. This increase will not change the agreement with the BCG masses from Kravtsov et al. (2018). It seems that both models predicts a slightly steeper slope for the ,0.1× 500 − 500 relation than Kravtsov et al. (2018), which could be caused by different amount of ICL contribution which will be studied in detail in Cañas et al. 2022 (in prep.).
As shown in Henden et al. (2020), the clusters from FABLE simulation have a systematically higher BCG stellar mass compared to observations (by ∼0.2-0.3 dex, which is still a significant improvement from Puchwein et al. 2010). The IllustrisTNG and C-EAGLE clusters (orange dashed line and orange circles in Fig. 5 of Henden et al. 2020, see also Pillepich et al. 2018) also possess significantly more massive BCGs than observed despite their more realistic galaxy sizes. Furthermore, Anbajagane et al. (2020)  Kravtsov+18 GIZMO-SIMBA GADGET-X the BCG mass (within 100 kpc) from Magneticum is even higher than IllustrisTNG, while BAHAMAS seems to have a slightly lower stellar mass compared IllustrisTNG. Ragone-Figueroa et al. (2018) show that their BCG masses are in very good agreement with observations, which could be the improvement on the BH's cold and hot modes gas accretion compared to their previous result (Ragone-Figueroa et al. 2013). In comparison to these works, the BCG masses in G -S seem somewhat more realistic.

BH-BCG-halo relations
All BCGs are believed to contain large supermassive black holes (BH). As shown by previous studies (e.g. Martizzi et al. 2012Martizzi et al. , 2014Frigo et al. 2019), AGN feedback is a key to have simulated galaxy clusters reproduce observed ones. Since AGN feedback is tied to the BH accretion rate which also sets the final BH mass, comparisons of the BH-galaxy-halo relations will help us to test and constrain the AGN feedback in models. In Figure 5 we compare G -S 's BH mass • vs. stellar mass * (left), stellar velocity dispersion * (middle), and 500 (left) versus observations. Note here * and * are obtained from the C catalogue instead of in a fixed aperture as in the previous section. In addition, only central galaxies are considered in Figure 5; i.e. no satellite galaxies are included. G -X is excluded from this comparison. For the observation results, we show data for early-type massive galaxies from Gaspari et al. (2019) (with BH properties from van den Bosch 2016), Sahu et al. (2019a), Sahu et al. (2019b), and Bandara et al. (2009), which are similar to or supersede various other data (McConnell & Ma 2013;Bogdán et al. 2018;Kormendy & Ho 2013;. While we don't specifically select early-type galaxies, we confirmed that our BCGs all have high bulge-to-total ratios and very low rotational support, which is expected since they are massive, quenched and appear bulge-dominated. The • − * relation is shown on the left panel of Figure 5. G -S shows good agreement with Gaspari et al. (2019) (and Kormendy & Ho 2013, not shown) at * ≈ 10 11.2 − 10 11.7 M . However, the BH masses for galaxies for * > ∼ 10 11.7 M galaxies are above these observations. In contrast, Sahu et al. (2019a) found a steeper result that is in very good agreement with G -S . It is beyond the scope of this work to delve into the differences between these observations results, so we simply note that G -S predicts that more massive galaxies will follow the same trend from Sahu et al. (2019a), all the way up to * ≈ 10 12.5 M .
Looking at the middle panel of Figure 5, there is more deviation between the G -S run and the observation results for the • − * relation. At log 10 * 2.6, G -S has a flatter slope compared to the observed results. Although earlier observations indicate a shallower slope (e.g. Sabra et al. 2015), the most recent observations prefer a steeper relation; for example Dullo et al. (2021) predicts a even higher slope than Sahu et al. (2019b) for early-type galaxies. Note that Thomas et al. (2019) showed that the • − relation from the original S simulation is in very good agreement with Kormendy & Ho (2013). This suggests that the re-tuned AGN feedback and lower resolution has more effects on internal dynamics than on the galaxy stellar mass -it seems to reduce the galaxy velocity dispersion at lower stellar or halo masses. Meanwhile at log 10 * 2.6, the • − * relation from G -S is basically following the extensions of Kormendy & Ho (2013); Gaspari et al. (2019).
There has been some controversy regarding the connection between BH mass and host halo mass, with some claiming a tight relationship (Booth & Schaye 2010;Bogdán & Goulding 2015;Bogdán et al. 2018;Bassini et al. 2019;Marasco et al. 2021), but others claiming it is incidental (e.g. Kormendy & Bender 2011) and are only related to classical bulges. For instance, Zhang et al. (2021) suggests that the BH and AGN feedback is not relevant for the formation history of massive star forming disks which show very high gas-to-star conversion rates, suggesting small BHs in their large halos. Nevertheless, we argue here that the • − halo relation exists for these BCGs in galaxy clusters which are mostly ellipticals. Hence examining the BH mass-halo mass connection in G -S is interesting.
The result for • − 500 from G -S is shown in the right panel of Figure 5. We compare to Bandara et al. (2009 and Gaspari et al. (2019). G -S is in very good agreement with the observed ones at group scales. However, there is a hint that at galaxy cluster scales, G -S predicts that the slope that becomes slightly shallower; observations do not extend into this regime. This is understandable as the BCGs in more massive clusters get quenched earlier (see more discussions in subsubsection 4.2.4), so like with the stellar mass, the halo continues to grow without  Figure 6. The g -r colour -magnitude diagram for the BCGs in the masscomplete cluster sample compared to the observation results in cluster mass range. The filled colourful contours are the results from SDSS (Yang et al. 2018). All the other simulation and observation data can be deciphered from the legend. The contours for the SDSS, simulation and SAM results are at the same percentiles: 16 ℎ , 50 ℎ and 84 ℎ . commensurate growth in the black hole. Interestingly, the relation is quite tight with 500 , even though neither BH accretion nor feedback in G -S is directly governed by any halo property.

BCG colours
BCG galaxies tend to be red in colour, indicating predominantly older stellar populations. The colour thus provides a constraint on the formation history of the BCG's stars. We determine the photometry of the galaxies in G -S and G -X using the P code implemented in C , which uses the age and metallicity of each star within a galaxy to generate a spectrum based on a simple stellar population model from FSPS 9 (Conroy et al. 2009;Conroy & Gunn 2010). We adopt a Chabrier (2003) initial mass function for the stars, matching the one used in the simulation. P also computes the dust extinction to each star based on the line-of-sight dust column, but for BCGs the dust model has a negligible effect since they have little cold ISM gas or dust.
In Figure 6, we compare the distributions of the BCG's − colour as a function of its magnitude in SDSS -band. The observational results are taken from Bernardi et al. (2011, best-fit shown as lime dashed line with dotted lines for the rms scatter); from Roche et al. (2010) for the results of both the C4-BCG (magenta points from Bernardi et al. 2007) and the max-BCGs (cyan points from Koester et al. 2007) samples, and from (Yang et al. 2018, coloured contours) with BCGs from their group catalogue. Note that the − colour is corrected in rest frame for all the observation data based on Chilingarian & Zolotukhin (2012).
The − colour from the observed galaxies typically ≈ 0.8, slightly higher for Yang et al. (2018). This is well in line with the BCG colours from G -S , particularly for the Roche et al. (2010) sample which extends to overlapping magnitudes. In contrast, G -X BCGs have a much bluer BCG colour, typically − ≈ 0.5−0.7. This indicates that BCG colours, and hence their star formation histories, are more realistic in G -S than in G -X. That said, for both G -S and G -X, the BCG brightnesses extend up farther than the data. The C4-BCG and max-BCG can reach the magnitude of ≈ −24.5 comparable to the brightest galaxies in G -S , but G -X shows BCGs extending to brighter than ≈ −26. For G -S , however, the mild differences can probably be attributed to sample selection, in that T T H focuses on only the most massive clusters. Particularly for G -X, the brighter BCGs go hand-in-hand with the bluer colour, since younger stellar populations are more luminous; this is unlikely to be reconciled by the aforementioned effect.

T T H
: G -S 13

Halo and BCG formation time
The different BCG colours in G -S and G -X suggest that the BCG stellar formation times are different in these models. The = 0 stellar populations arise from a combination of the star formation histories in the BCG's progenitors, plus mergers bringing in stellar mass. Hence in these simulations we can ask two distinct questions: (1) When were the stars in the BCG at = 0 formed? and (2) When were the progenitor galaxies of the BCG assembled? The first one tightly connects with the detailed star formation history of the BCG, which can be revealed from its mean stellar age; we refer to this as the star formation time. The second one is tightly connected to hierarchical structure formation, which is difficult to probe observationally; we refer to this as the assembly formation time. We expect this latter one may be more related to halo assembly, which will be quite similar within G -S and G -X by construction.
For the first question, the BCG stellar age has been estimated to be generally very old in observations. For example, Edwards et al. (2020) found the luminosity-weighted BCG core age of ∼13 Gyr (see also Loubser et al. 2009), which is in agreement with the semi-analytic model prediction from e.g. De Lucia & Blaizot (2007). However, note that two BCGs (Abell 671 and Abell 602) in the sample of Edwards et al. (2020) have quite young ages of ∼ 5 Gyr, which can be caused by more recent in-situ star formation (∼25% of the growth is in-situ; Ragone-Figueroa et al. 2018). By studying early-quenched galaxies in Magneticum and TNG simulations, Lustig et al. (2022) reported about a factor of 2 older for the average ages of simulated quiescent galaxies than observed ones.
For the second question, one can examine the most massive BCG at high redshift in observations to roughly constrain their assembly time. Collins et al. (2009) found several BCGs at z ∼ 1 − 1.5 that have stellar masses comparable to the most massive galaxies in the Universe, suggesting a very early assembly time of their BCGs. Similar suggestion is also found in Andreon (2013) which proposed with a even higher redshift ∼ 2.5 for the BCG mass build up time. More theoretical studies using hydro simulations find similar conclusions. For example, Ragone-Figueroa et al. (2018) found the assembly of the BCG occurs over an extended time-span and half of the BCGs' stellar mass only falls into place typically by ∼ 1.5. Similarly, the BCGs in FABLE, which are in good agreement with several observational inferences at 1 (Henden et al. 2020), have moderate stellar mass growth, about a factor of 1.5 from = 1 to 0.3, and effectively halts at 0.3 as suggested by observation. Furthermore, Rennehan et al. (2020) found that the stellar assembly time of a sub-set of brightest cluster galaxies occurs at high redshifts ( > 3) rather than at low redshifts ( < 1). Thus, highly overdense protoclusters assemble their stellar mass into brightest cluster galaxies within ∼1 Gyr of evolution, producing massive blue elliptical galaxies at high redshifts ( > 1.5). Here we only focus on the first question to determine if there is any connection between the BCG star formation time and the halo formation time.
In Figure 7, we show the halo formation time versus the BCG star formation time at = 0 in G -S (green points) and G -X (blue). The star formation time is the mean mass-weighted age of the stellar population, while the halo formation time here is quantified by the commonly used half , i.e. the redshift by which the halo had accreted half of its final total mass.
Interestingly, there is no obvious correlation between the BCG star formation time and the halo formation redshift, in either model. The halo formation redshifts are fairly recent, generally after ∼ 1. In contrast, the BCG mean stellar age shows the stars are fairly old. This is a manifestation of stellar population downsizing, where the stars forming in massive halos are old even though the most massive halos assemble late. This is has sometimes been called "antihierarchical" behaviour, but it is exactly as expected in hierarchical structure formation models where the largest density perturbations collapse first and hence can form the oldest stars, even though the halos assemble late.
In contrast, the BCG age distribution is quite different between G -S and G -X. Clearly, G -S has older BCGs, with the formation redshift peaking at ∼ 3, and little star formation after ∼ 2. The mean BCG ages from G -X is younger, as young as ∼ 5 − 6 Gyrs ago, indicating late star formation which can be caused by either the in-situ star formation shown in Ragone-Figueroa et al. (2018) or rejuvenation at low redshift. This also results in a much larger tail towards young mean stellar ages in G -X. These results are consistent with our earlier findings that the stellar assembly within halos is shifted towards earlier epochs in G -S vs. G -X. This owes to both the jet feedback, as well as G -S 's X-ray feedback that is responsible for quenching star formation completely ). In the original S , the SFR function at ∼ 2 − 4 extends to very high SFRs, > 10 3 M yr −1 , which turns out to be critical for reproducing the sub-millimetre galaxy (SMG) population at these epochs (Lovell et al. 2021), which are likely to be the progenitors of BCGs today; other simulations tend to have a more a difficult time matching these constraints. It seems that G -X may also under-predict the SFR in protoclusters at high redshift, especially the starburst population, compared to observations (see Bassini et al. 2020, for more details). It appears that such rapid early stellar growth within massive halos is crucial both for matching BCG colours today, as well as SMGs at ∼ 2 − 4.

Satellite galaxies
Satellite galaxies provide another important constraint on the stellar properties of galaxy clusters. They can connect with the cluster total mass through cluster richness (for example Anbajagane et al. 2020) or velocity dispersion (for example Munari et al. 2013;Anbajagane et al. 2022b;Ferragamo et al. 2022), with the cluster's dynamical state through its distribution and mass fraction (for example Cui et al. 2017), and with cluster shapes and orientations. Successfully reproducing the satellite galaxy properties in galaxy cluster is very important for studying environment quenching processes such as ram-pressure stripping and harassment (see Lotz et al. 2019Lotz et al. , 2021. Therefore it is important to check how well our simulations reproduce observations of the satellite population in groups and clusters.

C3: Satellite stellar mass function
We first investigate whether T T H simulations produce the correct amount of satellites at different stellar mass bins in Figure 8. The satellite stellar mass functions (SSMFs) for G -S (green) and G -X (blue) are shown, along with observations of the SSMF in clusters from Yang et al. (2018) within two mass bins spanning much of T T H sample. Recall that the satellite stellar mass function is used to constrain the feedback parameters in G -S (cf. constraint C3), so the agreement with G -S should not be considered a true success of the model, although the tuning was only done on a single cluster. Detailed comparisons and discussions among the other T T H models can be found in Cui et al. (2018b). The median satellite stellar mass function (SSMF) from G -S show very good agreement with the observation results, which confirms that the tuning down for one region was successful even when considering all 324 regions. G -S and G -X agree at * > ∼ 10 10 ℎ −1 M , below which G -X falls below the observations. There are some interesting features, such as G -X showing a strong peak at * > ∼ 10 10.3 ℎ −1 M and G -S showing a steep decline at slightly higher masses; we suspect that these arise due to the onset of strong AGN feedback at these masses . This is because the BH is seeded and AGN feedback is in action when the galaxy stellar mass reach * ∼ 10 10.5 ℎ −1 M , which then halts galaxy growth and piles up galaxies just below this threshold. A more gradual onset of AGN feedback with halo mass may produce a smoother SSMF; we leave this for a future work.
Using the Hydrangea simulations, Bahé et al. (2017) reported very good agreement when compared with satellite stellar mass functions from various observations (Yang et al. 2009;Vulcani et al. 2011;Wang & White 2012), down to lower masses than we can probe in T T H . Although FABLE did not show its satellite stellar mass function, Henden et al. (2020) found that the fraction of the total stellar mass in satellite galaxies is significantly smaller in FABLE clusters (∼40 per cent) than observed. This could be caused by the stars in their satellite galaxies being stripped at an accelerated rate which would also result in a larger ICL mass fraction. It is generally the case that it is easier to strip stars in lower resolution simulations, owing to less deep potential wells, hence this effect may also be relevant for T T H . Planned higher-resolution runs will be able to quantify this in the future.

Satellite galaxy colour distribution
Similar to the colour-magnitude diagram for the BCGs, the satellite colour-magnitude diagram gives insights into their stellar growth histories. With G -S in particular using the SSMF as a constrain, it is important to see whether it also produces a reasonable distribution of satellite colours as a complementary check on the model.
As shown in Figure 9, the − colours vs. magnitudes are shown from G -S (green contours) and G -X (blue), with the observations from (Yang et al. 2018, coloured shading). As with the BCGs, their galaxy colour has been corrected into the rest frame 10 . 10 Note that this plot is slightly different to the middle panel of Fig. 8 in Cui et al. (2018b) due to a different stellar mass cut ( * ≥ 10 9 ℎ −1 M ) applied. By comparing to Fig. 8 in Cui et al. (2018b), G -X tends to have a bluer colour, which suggests that less massive satellite galaxies should be dominated by red colour. Satellite galaxy colours from G -S show very good agreement with Yang et al. (2018), with again typically − ≈ 0.8 just as with the BCGs. Hence the stellar populations in both the central and satellite galaxies in G -S tend to be quite old. In contrast, G -X results, meanwhile, show a somewhat bluer − colour, mimicking the trend seen in the BCG colours. The satellites are somewhat redder than the BCGs in this model.
It is notable however, the G -S 's distribution of colours has a small scatter compared to Yang et al. (2018). In particularly, G -S lacks any bluer satellites with − < ∼ 0.7. While these are rare in the observations, they do occur. This could be caused by two reasons: (1) the strong AGN feedback quenches all the satellite galaxies instead of only most of them; (2) the low resolution means that gas stripping is too efficient, which quenches satellite galaxies progenitors prior to entering into the cluster or too quickly within the cluster. In the future we plan to run higher resolution T T H runs to check whether such resolution effects can impact these comparisons.

CONCLUSIONS AND DISCUSSIONS
In this paper, we introduce a new suite of simulated galaxy clusters in T T H project called G -S , run using a feedback model based on the successful S simulation . We re-tune its feedback model parameters (mostly due to the coarse resolution of T T H clusters) by calibrating them with the observed stellar properties using one random cluster, specifically the stellar fraction, the BCG mass-halo mass relation and the satellite stellar mass function. Using this calibration we run all 324 clusters, and make detailed comparisons to both observational results and the G -X runs which is also successful in reproducing many galaxy cluster properties and relations. By comparing the two hydrodynamic in T T H we aim to under- . Satellite galaxy colour-magnitude diagram. This plot is similar to Figure 6, but for all the satellite galaxies in the mass-complete cluster sample. Contour maps are the observation results from Yang et al. (2018) with clusters 200 6 × 10 14 ℎ −1 M . Again the three contour lines/maps mark the 16 ℎ , 50 ℎ and 84 ℎ percentiles. Note here that we apply the same stellar mass cut * ≥ 10 10 M to both simulation and observed dataset for a consistent comparison. stand the origin of observational properties of galaxy clusters under different galaxy formation models. We note that a comparison between G -X and other galaxy formation models in T T H was done in Cui et al. (2018a). Our main results are as follows: • Similar stellar mass fractions * (< 500 ) are found at = 0 between G -X and G -S , which are also in agreement with observational results. However, their gas fractions at = 0 are different (∼ 5 per cent in absolute mass fraction) at 500 ∼ 10 14 M , with agreement for the most massive haloes and larger differences for low mass haloes. This difference is seen at all redshifts, either at a fixed halo mass or by tracking progenitors. This could be rooted mainly in the feedback scheme, of which G -S 's kinetic feedback is very strong and efficiently blows gas out of haloes. More accurate gas fractions from observation are required to distinguish between models.
• Though the stellar mass fraction is in agreement at = 0 in both models, G -S has a much stronger evolution of this fraction versus redshift than G -X, especially at high redshift 1.5. This could be caused by two reasons: (1) the star formation model is more efficient; (2) the quenching of galaxies is very quick (at mostly high redshift) and thorough. The high star formation rate in G -S produces a higher stellar mass fraction at high redshift, while at later epochs the complete shutdown of star formation makes its stellar mass fraction similar to G -X by = 0. • Using BCG stellar masses within fixed apertures, both G -X and G -S are in agreement with observations from Kravtsov et al. (2018). However, the BCGs (identified by the 6D C galaxy finder) from G -X have much brighter and bluer colour compared to G -S ; observations agree better with G -S . This is corroborated by the mean BCG ages in G -X being systematically lower than in G -S , with a long tail of younger BCGs.
• To test G -S 's black hole accretion and feedback models, we compare three BH scaling relations from G -S versus * , * , and 500 . G -S matches observations generally well, although it predicts higher black hole masses in the most massive galaxies, a shallower slope for the • − * at the low- * end, and a shallower • − 500 relation at the more massive end. A more careful comparison versus these data mimicking aperture and selection effects is required to assess the significance of these discrepancies.
• For satellite galaxies, we focus on satellite stellar mass function and colour-magnitude diagram. Although both agree with data for more massive galaxies, G -S shows better agreement to observation than G -X at the lower mass end, which is only partly a result of tuning. The satellite galaxies in G -S agree very well with the − vs. colour-magnitude diagram of the SDSS galaxies, while G -X satellites are somewhat bluer. However, G -S shows less scatter than observations, with no satellites bluer than − < ∼ 0.7 whereas observations indicate some.
The gas scaling relations are compared in the appendixes. As a result of G -S 's strong feedback, its temperature-mass ( − ) relation tends to be a little bit higher than G -X. However, in combination with its slightly lower gas fractions, the resulting integrated Sunyaev-Zeldovich decrement vs. mass ( − ) relation is quite similar to that in G -X. The G -S runs in T T H provide a state of the art galaxy formation model with which to examine the evolution of galaxies, gas, and black holes within the largest virialized structures in the universe. Many follow-up works using this suite are already underway. For instance, for the BCG and satellite galaxies, we don't provide details on their evolution in this paper because (1) to disentangle BCG with ICL needs a careful study, which we will present in a companion paper (Cañas et al. 2022); (2) we need to do a careful tracking to fully look into the the evolution history of satellite galaxies in hydro-simulations because as discussed in Behroozi et al. (2015), halo and galaxy tracking is complicated within the cluster environment to properly account for mergers and tidal stripping (see Ogiya et al. 2022, and references therein, for the formation of the ultra-diffused galaxies without dark matter). These sorts of studies can shed further light on how galaxies evolve within dense environments.

ACKNOWLEDGEMENTS
We thank the referee Prof Peter Thomas for his valuable comments and suggestions on this paper. We would like to thank T T H collaboration for all the help on this project and this paper.

DATA AVAILABILITY
The results shown in this work use data from T T H galaxy clusters sample. These data are available on request following the guidelines of T T H collaboration, at https: //www.the300-project.org. The data specifically shown in this paper will be shared upon request to the corresponding author. The plots and analyses done in this paper have extensively used these python packages: Ipython with its Jupyter notebook (Pérez & Granger 2007), astropy (Astropy Collaboration et al. 2018, NumPy (van der Walt et al. 2011) and SciPy (Oliphant 2007;Millman & Aivazis 2011). All the figures in this paper are plotted using the python matplotlib package (Hunter 2007 GADGET-X: A=0.67, B=0.61 GIZMO-SIMBA: A=0.74, B=0.54 Vikhlinin et al. 2006Vikhlinin et al. 2009Lovisari et al. 2015 GADGET-X GIZMO-SIMBA Figure A1. The gas mass-weighed temperature 500 -halo mass 500 relation. As indicated in the legend at top-left corner, data points from G -X and G -S are shown with blue and green symbols with error bars (16th -84th percentile) respectively. The solid and dotted black lines show the observational results from Vikhlinin et al. (2006) and Vikhlinin et al. (2009), respectively. The maroon dashed line shows the fitting result from (Lovisari et al. 2015). Our weighted fitting results from G -X and G -S are presented by blue dotted and lime dashed lines, respectively. The thick solid black line shows the self-similar relation log 10 500 ∝ 2/3 log 10 500 is derived from self-similarity cluster approximations (see Böhringer et al. 2012, for example). mass relation ( − ; Appendix B) (see Lovisari et al. 2021, for a recent review). The − relation will test whether the gas temperature in G -S is affected by its strong feedback or not; while the − relation is checks whether the differences in gas fraction and temperature can be directly viewed from SZ observations. Following Cui et al. (2018b), we calculate the mass-weighted gas temperature within 500 and add the scaling relation from G -S in Figure A1. In the cluster mass range, G -S is in line with G -X and the observational results, albeit a slightly higher temperature (see the zoomed-in inset panel for details). This suggests that the gas heating inside these most massive clusters is reasonable. However, the median data points from G -S at the group halo mass range, 500 < 10 14 M , seem to be slightly higher than the other results, albeit with a very large scatter. This is likely related to the strong AGN feedback which is also efficiently expelling the gas out of 500 as shown in Figure 1.
By applying a weighted fitting to log 10 500 = + * log 10 500 6×10 14 M , in which the weight is given by the completeness of the halo counts, G -S ( = 0.54) seems to be deviating from the self-similar slope much more than G -X ( = 0.61). This is mainly because of the high temperature in G -S at lower halo masses. Note that the median data points deviating from the fitting line at lower halo masses is mainly caused by the weighted fitting method. Furthermore, the slightly different fitting results between this work and Cui et al. (2018b) is because the new mass completeness based on G -S is adopted (see Cui (2022) for the changes). This indicates the fitting results are very sensitive to the weights. We note that observational results are quite uncertain within this mass regime, so it is too early to determine which model agrees better with data.

APPENDIX B: THE GAS Y-M RELATION
It is also interesting to compare the overall simulated SZ signal, the integrated SZ 500 value in this case, versus that directly obtained from observation. For this, we apply the P MSZ package 11 to G -S simulated halos to generate mock −maps along z-direction. Then, 500 is taken to be the sum of the values within 500 . Note that the ling-of-sight depth is set to 2 × 500 to be consistent with published results from G -X; we refer to  for more a detailed examination of projection effects.
As shown in Figure B1, the median data points from G -S are comparable to the ones from G -X. This is not surprising as the SZ− signal is proportional to gas density times temperature (Sunyaev & Zeldovich 1970, and references therein). The lower gas fraction in G -S seems to be compensated by its high temperature to give such an agreement. Note that the large scatters in gas mass fraction and temperature in G -S at the low halo mass range are also shown in this 500 − 500 relation. We also apply a weighted fitting of this 500 − 500 relation to log 10 2 500 = + * log 10 500 6×10 14 M for the two hydro runs, where is the angular diameter distance. Given the similarity of the predictions, it is not surprising to see such a good agreement in the fitting parameters between G -S and G -X. Regarding the slope difference to the Planck result (Planck Collaboration et al. 2014), we note that the fitting line actually deviates from the median data points at both massive and low halo mass range: the slope from the median data points is more close to or even flatter than the selfsimilar slope of 5/3. Adopting the weighted fitting makes the slope dominated by these massive halos, and results in a value close to self-similar. Meanwhile at the low mass end, the slope from the data 11 Publicly available at https://github.com/weiguangcui/pymsz, see Baldi et al. (2018) for the details of generating the kinetic SZ signal with this package. points is much steeper in both G -X and G -S . If we do an unweighted fitting to these data points, we get a steeper slope which is comparable to the result from Planck. The different slopes for 500 − 500 relation across various mass ranges have been found in several works using different hydrodynamic simulations (e.g. . Therefore, using a fitting function with a single slope is not a good choice for representing the full mass range.

APPENDIX C: RESOLUTION EFFECTS
As described in the main text, we re-tuned the baryonic feedback parameters for T T H galaxy clusters relative to S due to their poor resolution. This was because we found that we cannot match key observations using the default S setup. It would be interesting to see whether the default S setup works for the clusters which have a similar resolution as its original run. This requires doing a higher-resolution run of a T T H cluster, with a similar resolution to that in S . We show the results of such a run with 8× higher mass resolution in this section using one of our cluster regions.
In Fig. C1, we show the comparisons for both gas and star fractions between the high-resolution run with our default run from the same region. It seems that these fractions are generally in agreement between the two runs; there are no obvious systematic differences, albeit with small numbers. There may be a hint that gas fractions are slightly higher at low masses when run at higher resolution. Furthermore, we would caution that there maybe more significant differences in the detailed galaxy cluster properties, such as profiles (as suggested in previous nIFTy comparison projects Cui et al. 2016;Sembolini et al. 2016b;Elahi et al. 2016). To statistically investigate and present these differences, we need a large sample of the high-resolution cluster runs, which is currently in progress. This paper has been typeset from a T E X/L A T E X file prepared by the author.  Fig. 1, baryon fraction within 500 . This plot focuses on comparing these fractions in the same zoomed-in region but one run with the default SIMBA setup for a high-resolution IC, marked as 'HR-  in open magenta stars, and one run in this paper with about 8 time lower resolution, marked as 'LR-CL-192' with black cross. The overall statistical result from G -S , green solid line with error bars, in Fig. 1 is also shown here for reference. There is not a large systematic effect owing to resolution obvious for this single object, but a larger sample is needed to quantify the differences.