#  The Three Hundred project: The Gizmo-Simba runs

Weiguang Cui,1 Romeel Dave,1 Alexander Knebe,2,3,4 Elena Rasia,5,6 Meghan
Gray,7 Frazer Pearce,7 Chris Power,3 Gustavo Yepes,2,3 Dhayaa Anbajagane,8,9
Daniel Ceverino,2,3 Ana Contreras-Santos,2 Daniel de Andres,2,3 Marco De
Petris,10 Stefano Ettori,11,12 Roan Haggar,7 Qingyang Li13 Yang Wang,14,15
Xiaohu Yang,13,16 Stefano Borgani,5,6,17 Klaus Dolag,18,19 Ying Zu,13,16
Ulrike Kuchner,7 Rodrigo Cañas,2,3 Antonio Ferragam,10 Giulia Gianfagna20
1Institute for Astronomy, University of Edinburgh, Royal Observatory,
Edinburgh EH9 3HJ, United Kingdom
2Departamento de Física Teórica, Módulo 15, Facultad de Ciencias, Universidad
Autónoma de Madrid, 28049 Madrid, Spain
3Centro de Investigación Avanzada en Física Fundamental (CIAFF), Facultad de
Ciencias, Universidad Autónoma de Madrid, 28049 Madrid, Spain
4International Centre for Radio Astronomy Research, University of Western
Australia, 35 Stirling Highway, Crawley, Western Australia 6009, Australia
5IFPU - Institute for Fundamental Physics of the Universe, Via Beirut 2, 34014
Trieste, Italy
6INAF Osservatorio Astronomico di Trieste, via Tiepolo 11, I-34131, Trieste,
Italy
7School of Physics & Astronomy, University of Nottingham, Nottingham NG7 2RD,
United Kingdom
8Department of Astronomy and Astrophysics, University of Chicago, Chicago, IL
60637, USA
9Kavli Institute for Cosmological Physics, University of Chicago, Chicago, IL
60637, USA
10Dipartimento di Fisica, Sapienza Universitá di Roma, Piazzale Aldo Moro,
5-00185 Roma, Italy
11INAF, Osservatorio di Astrofisica e Scienza dello Spazio, via Piero Gobetti
93/3, 40129 Bologna, Italy
12INFN, Sezione di Bologna, viale Berti Pichat 6/2, 40127 Bologna, Italy
13Department of Astronomy, School of Physics and Astronomy and Shanghai Key
Laboratory for Particle Physics and Cosmology, Jiao Tong University, Shanghai
200240, China
14Department of Mathematics and Theories, Peng Cheng Laboratory, No.2, Xingke
1st Street, Nanshan District, Shenzhen 518000, Guangdong Province, P. R. China
15CSST Science Center for Guangdong-Hong Kong-Macau Great Bay Area, Zhuhai
519082, China
16Tsung-Dao Lee Institute and Key Laboratory for Particle Physics,
Astrophysics and Cosmology, Ministry of Education, Jiao Tong University,
Shanghai 200240, China
17University degli Studi di Trieste, Trieste, Italy
18Universitäts-Sternwarte, Fakultät für Physik, Ludwig-Maximilians-Universität
München, Scheinerstr.1, 81679 München, Germany
19Max-Planck-Institut für Astrophysik, Karl-Schwarzschild-Straße 1, 85741
Garching, Germany
20INAF, Istituto di Astrofisica e Planetologia Spaziali, via Fosso del
Cavaliere 100, 00133 Rome, Italy
E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

We introduce Gizmo-Simba, a new suite of galaxy cluster simulations within The
Three Hundred project. 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 Gizmo-Simba suite adds a state-of-the-art galaxy formation
model based on the highly successful Simba simulation, mildly re-calibrated to
match $z=0$ cluster stellar properties. Comparing to The Three Hundred zooms
run with Gadget-X, we find intrinsic differences in the evolution of the
stellar and gas mass fractions, BCG ages, and galaxy colour-magnitude
diagrams, with Gizmo-Simba generally providing a good match to available data
at $z\approx 0$. Gizmo-Simba’s unique black hole growth and feedback model
yields agreement with the observed BH scaling relations at intermediate mass
range, and predicts a slightly different slope at high masses where few
observations currently lie. Gizmo-Simba provides a new and novel platform to
elucidate the co-evolution of galaxies, gas, and black holes within the
densest cosmic environments.

###### keywords:

keyword1 – keyword2 – keyword3

††pubyear: 2020††pagerange: The Three Hundred project: The Gizmo-Simba runs–B

##  1 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 the
longest structure formation histories and environmental studies of galaxies
with the evolution of the hot intracluster medium (ICM). On galactic scales,
they are important for the study of the oldest galaxy stellar populations and
most massive galaxies with the impact of supermassive black holes on galaxies.
For these reasons, clusters have been the object of much study from both an
observational and theoretical standpoint (see Allen et al., 2011; Kravtsov &
Borgani, 2012; Walker et al., 2019, for different reviews).

Interpreting and understanding observations of galaxy clusters from the radio
to the X-rays are 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. As such, 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), 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),
FABLE (Henden et al., 2018), MUSIC (Sembolini et al., 2013), and the nIFTy
cluster comparison project (Sembolini et al., 2016), are able to investigate
detailed baryon 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), IllustrisTNG (Pillepich et al., 2018), Magneticum
(Dolag et al., 2016), BAHAMAS (McCarthy et al., 2017), and Simba (Davé et al.,
2019; 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.

The Three Hundred project111https://the300-project.org/, we also use the300
for short. occupies a unique niche among cluster simulation suites. The Three
Hundred 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. Other suites have typically focused on a handful of
clusters, and smaller zoom regions. 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.

Another feature of The Three Hundred project is that all these clusters have
been run with several different galaxy evolution codes. These include the
cosmological hydrodynamic codes Gadget-MUSIC (Sembolini et al., 2013) and
Gadget-X (Rasia et al., 2015), as well as three different semi-analytical
models (SAMs): Galacticus (Benson, 2012), SAGE (Croton et al., 2016) and SAG
(Cora et al., 2018). This enables cross-comparisons between models employing
different input physics, to better understand the sensitivity to the various
physical processes and the robustness of the resulting predictions.

To date, The Three Hundred has been used in over 30 papers investigating
galaxy clusters and their environs. These include the detailed relationship
between the central cluster and connecting filaments (Rost et al., 2021;
Kuchner et al., 2020, 2021b), the feeding of the galaxy clusters(Kuchner et
al., 2021a; Kotecha et al., 2021), cluster backsplash galaxies (Haggar et al.,
2020; Knebe et al., 2020), and the virial shock radius (Baxter et al., 2021;
Anbajagane et al., 2021). The advanced input physics in the hydrodynamic
simulations allow us to perform a detailed investigation on the cluster
properties, such as cluster profiles (Mostoghiu et al., 2019; Li et al.,
2020), substructure and its baryonic content (Arthur et al., 2019; Haggar et
al., 2021; Mostoghiu et al., 2021b, a), dynamical state and morphology
(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., 2021), the effects of
mergers on the BCG properties (Contreras-Santos et al., 2022) and the
hydrostatic equilibrium mass bias (Ansarifard et al., 2020; Li et al., 2021;
Gianfagna et al., 2022). 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 that allows constraints on the dark matter
cross-section (Vega-Ferrero et al., 2021); and even to examine chameleon
gravity (Tamosiunas et al., 2021). With many more projects in the works, it is
valuable to continue to update The Three Hundred runs with state-of-the-art
physical models to expand its range and robustness.

In this paper, we introduce another set of hydrodynamic runs to The Three
Hundred suite, namely the Gizmo-Simba runs 222Note that to distinguish from
the Simba simulation (Davé et al., 2019) which is a 100${{\,h^{-1}{\rm Mpc}}}$
cosmological hydrodynamic simulation run with the same code, this run for The
Three Hundred clusters is referred to as Gizmo-Simba run.. This suite uses the
Gizmo cosmological hydrodynamics code in its Meshless Finite Mass (MFM) solver
mode, as opposed to the Gadget-based runs which use Smoothed Particle
Hydrodynamics (SPH). It further includes a suite of galaxy formation physics
similar to that in the recent Simba simulation (Davé et al., 2019) that yields
an excellent match to a wide range of galaxy, black hole, and intergalactic
medium properties. Moreover, it includes some novel input physics modules such
as torque-limited black hole growth, stably bipolar jet feedback, and on-the-
fly dust tracking, making it a valuable addition to the existing The Three
Hundred suite. Gizmo-Simba thus provides a state-of-the-art new set of
simulations to investigate and compare to the various other models in The
Three Hundred project. Here we present the Gizmo-Simba results and its
calibration based on the stellar properties within galaxy clusters, after
which we make extensive comparisons to the Gadget-X runs which are mostly
calibrated by gas properties and represents the other state-of-the-art
hydrodynamic galaxy formation model within The Three Hundred. This work thus
sets the stage for many future investigations of cluster physics across
multiple hydrodynamic and semi-analytic models.

This paper is organised in the following order: we first introduce the The
Three Hundred project in §2. Then, we present the details of the new Gizmo-
Simba 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.

##  2 The Three Hundred project

The Three Hundred project (Cui et al., 2018b) is a set of cluster-scale zoom
simulations based on a mass-complete sample of 324 most massive galaxy
clusters ($M_{\rm vir}\gtrsim 8\times
10^{14}{{\,h^{-1}{\rm{M_{\odot}}}}}$)333The halo mass is defined as the mass
enclosed inside an overdensity of $\delta$ times the critical density of the
universe: $\delta=\sim 98$ for virial mass at $z=0$ (Bryan & Norman, 1998) and
$M_{200,500}$ is with $\delta=200,500$ respectively. Similarly, $R_{500}$ is
the radius at which the overdensity $\delta=500$ is reached. 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\,h^{-1}~{}{\rm Gpc}$ containing $3840^{3}$ DM
particles having a mass of $1.5\times 10^{9}{{\,h^{-1}{\rm{M_{\odot}}}}}$
each. Each cluster region was selected to have a comoving radius of
$15{{\,h^{-1}{\rm Mpc}}}$ (over $5\times R_{200}$) for re-simulation with
different baryonic models: Gadget-MUSIC (Sembolini et al., 2013), Gadget-X
(Rasia et al., 2015; Beck et al., 2016), and now Gizmo-Simba (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),Galacticus (Benson, 2012).

The re-simulation regions are generated with the parallel GINNUNGAGAP
code444https://github.com/ginnungagapgroup/ginnungagap: the highest resolution
Lagrangian regions share the same mass resolution as the original MDPL2
simulation with gas particles ($M_{\rm gas}=2.36\times
10^{8}{{\,h^{-1}{\rm{M_{\odot}}}}}$) split from DM particles. The outside
regions are degraded in multiple layers (with a shell thickness of $\sim
4{{\,h^{-1}{\rm Mpc}}}$) with lower mass resolution particles (mass increased
by eight times for each layer) that eventually provide the same tidal fields
yet at a much lower computational costs than in the original simulation.

##  3 The Gizmo-Simba run

###  3.1 The Simba model

The Gizmo-Simba runs of The Three Hundred are performed with the Gizmo code
(Hopkins, 2015) with the state-of-the-art galaxy formation subgrid models
following the Simba simulation (Davé et al., 2019). We refer the interested
reader to Davé et al. (2019) for full details of all of Gizmo-Simba’s
features, and here focus on its more unique aspects relevant for clusters. We
also describe our modifications to the Simba model parameters utilised for the
Gizmo-Simba runs of The Three Hundred clusters, which required re-tuning owing
to the lower numerical resolution in Gizmo-Simba relative to the original
Simba simulation.

Owing to the MFM solver implemented in Gizmo, 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 self-shielding on the fly based on the
prescription in Rahmati et al. (2013). An $\rm H_{2}$-based star formation
model is taken from its predecessor simulation Mufasa (Davé et al., 2016),
which is calibrated to match the Schmidt (1959) law. Here $\rm H_{2}$ is
estimated from the local column density and metallicity following the Krumholz
et al. (2009); Krumholz & Gnedin (2011) prescription. Besides requiring that
$H_{2}$ be present, an additional minimum density cut of $n_{\rm H}>0.1~{}{\rm
cm}^{-3}$ and a minimum metallicity of 0.05 for the $f_{\rm H_{2}}$ in its
formation (Krumholz et al., 2009) are required for active star formation;
these numbers are slightly modified from the original Simba simulation of
$0.1~{}{\rm cm}^{-3}$ and 0.01, respectively.

Star formation-driven galactic winds also shares the same decoupled two-phase
model in Mufasa, 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):

$\eta(M_{*})\approx\begin{cases}9(\frac{M_{*}}{M_{0}})^{-0.317}&\text{if
$M_{*}<M_{0}$}\\\ 9(\frac{M_{*}}{M_{0}})^{-0.761}&\text{if $M_{*}\geq
M_{0}$}\end{cases}.$ (1)

Here, $M_{0}=2\times 10^{9}\,\rm{M_{\odot}}$; this is slightly different than
the original Simba simulation, as we will motivate later. The ejection
velocity is based on scalings from Muratov et al. (2015) as in Mufasa:

$v_{w}=0.854\left(\frac{v_{\rm circ}}{200\ {\rm kms^{-1}}}\right)^{0.12}v_{\rm
circ}+\Delta v(0.25R_{\rm vir}),$ (2)

where $\Delta v(0.25R_{\rm 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 Simba, who used a
normalisation of 1.7; the normalisation used here in Gizmo-Simba is the
original one proposed by Muratov et al. (2015). Identically to Simba, 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 $M_{*}$ are also computed on-the-fly, with $v_{\rm circ}$
obtained from a scaling based on the baryonic Tully–Fisher relation (McGaugh,
2012).

The chemical enrichment model tracks eleven elements (H, He, C, N, O, Ne, Mg,
Si, S, Ca, Fe) 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, Simba also consider metal-loaded winds, i.e. metal
of the wind particle is subtracted from nearby gas in a kernel-weighted
manner. All this is identical to the original Simba model; see Davé et al.
(2019).

Simba seeds black hole (BH) particles based on the host galaxy stellar mass,
$M_{*}>\gamma_{\rm BH}\times M_{\rm 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 Gizmo-Simba run, we employ
$M_{\rm seed}=10^{5}{{\,h^{-1}{\rm{M_{\odot}}}}}$ and $\gamma_{\rm BH}=3\times
10^{5}$, which sets the galaxy stellar mass threshold for seeding BH is
$M_{*}\approx{}10^{10.5}{{\,h^{-1}{\rm{M_{\odot}}}}}$. This is 10 times higher
than original Simba model. 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, $R_{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 $R_{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 ($T\leq 10^{5}K$),
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; Anglés-Alcázar et al. 2015, 2017a); The hot-
accretion mode is based on the Bondi accretion from the hot ($T>10^{5}K$) gas
component (Bondi, 1952). We reduce the Bondi accretion rate to the rate for a
$M_{\rm BH}=10^{9}{{\,h^{-1}{\rm{M_{\odot}}}}}$ BH, no matter how big it gets;
in Simba, this was set to $10^{10}M_{\odot}$. The black hole accretion kernel
has a distance enclosing 256 particles or $R_{0}=6{{\,h^{-1}{\rm kpc}}}$
(comoving), whichever is smaller, within which the gas quantities are
calculated; the latter was $2{{\,h^{-1}{\rm kpc}}}$ in Simba. The total
accretion rate for a given black hole is then the sum of $\dot{M}_{\rm
Torque}$ and $\dot{M}_{\rm Bondi}$ times an additional constant $1-\eta$ with
the radiative efficiency $\eta=0.1$ (Davé et al., 2019). 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 energy
feedback mode accounting for 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 ($f_{\rm Edd}\equiv\dot{M}_{\rm
BH}/\dot{M}_{\rm Edd}>0.2$). The radiative velocity for wind particles scales
as:

$v_{w,\ \rm Radiative}=500+\frac{500}{3}\left(\log\frac{M_{\rm
BH}}{\,\rm{M_{\odot}}}+6\right)\,\rm km\,s^{-1}.$ (3)

When the BH’s in low Eddington accretion mode, $f_{\rm Edd}<0.2$, the wind
begins to transition into a jet mode, with the velocity scaled with $f_{\rm
Edd}$ as follows:

$v_{w,\ \rm Jet}=v_{w,\ \rm Radiative}+15000\log\left(\frac{0.2}{f_{\rm
Edd}}\right)\,{\rm km\,s^{-1}}.$ (4)

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 $R_{0}$. The
wind velocity in the jet mode is capped at $15000\,\rm km/s$ (as opposed to
$7000$ km/s in original Simba) when the Eddington rate drops below $f_{\rm
Edd}\leq 0.02$. There is another condition to trigger jet mode – the minimum
BH mass has to be greater than $10^{7.5}\,\rm{M_{\odot}}$. X-ray feedback only
operates when the ‘jet mode’ AGN feedback is in action, and it further has to
meet another two conditions: $M_{*}>10^{9}\,\rm{M_{\odot}}$ and $M_{\rm
gas}/M_{\rm baryon}<0.2$. The $M_{*}$ condition is raised with respect to the
original Simba model, owing to the lower resolution.

Another new feature of the Simba simulation, the on-the-fly dust production
and destruction model, is also employed in Gizmo-Simba, unchanged. This model
is quite successful in reproducing galaxy dust properties over cosmic time (Li
et al., 2019).

Parameter | Original Simba simulation |  Gizmo-Simba run for the300 clusters
---|---|---
$n_{\rm H}$ threshold for SF | 0.13 | 0.1
Metallicity floor for $f_{\rm H2}$ | 0.001 | 0.05
$M_{0}$ in SN mass loading factor | $5.2\times 10^{9}\,\rm{M_{\odot}}$ | $2\times 10^{9}\,\rm{M_{\odot}}$
Galaxy stellar mass limit for seeding BH | $10^{9.5}{{\,h^{-1}{\rm{M_{\odot}}}}}$ | $10^{10.5}{{\,h^{-1}{\rm{M_{\odot}}}}}$
BH mass for Bondi accretion rate cap | $10^{10}{{\,h^{-1}{\rm{M_{\odot}}}}}$ | $10^{9}{{\,h^{-1}{\rm{M_{\odot}}}}}$
BH accretion kernel radius | $2{{\,h^{-1}{\rm kpc}}}$ | $6{{\,h^{-1}{\rm kpc}}}$
Cap wind velocity limit in jet mode | $7000\,\rm km/s$ | $15000\,\rm km/s$
Galaxy stellar mass limit for X-ray feedback | 0 | $10^{9}\,\rm{M_{\odot}}$
Gravitational softening length |  $0.5{{\,h^{-1}{\rm kpc}}}$ minimum |  $5{{\,h^{-1}{\rm kpc}}}$ fixed
Table 1: Summary of parameter changes with respect to the original Simba
simulation.

###  3.2 Re-calibration for The Three Hundred

We now summarise the modifications from the original Simba 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 The
Three Hundred simulation has about 10 times worse mass resolution than the
Simba simulation. It has been found that the Simba model is reasonably well
converged towards higher resolutions than in the original Simba 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 Gizmo-Simba runs,
all at $z\approx 0$: (C1) total stellar mass fraction within $R_{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-and-error 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 whose largest object has $M_{500}\approx 5\times
10^{14}M_{\odot}$, chosen to be typical of The Three Hundred sample555Ideally,
we could calibrate the parameters for each cluster region, then used 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.. 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 $\approx 10^{13}M_{\odot}$), 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 Simba parameters. At
this low resolution, we found that the star formation is insufficient for
satellite galaxies, resulting in a much lower satellite stellar mass function
(SSMF) compared to observations (C3). The changes for $n_{\rm H}$, Metallicity
floor, $M_{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 (Cui et al., 2021), 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{{\,h^{-1}{\rm kpc}}}$ for consistency
with the other The Three Hundred runs, as opposed to a variable softening
length with a minimum of $0.5{{\,h^{-1}{\rm kpc}}}$ as in original Simba. The
BH accretion kernel maximum radius was commensurately increased in Gizmo-Simba
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.

###  3.3 Detailed differences between Gadget-X and Gizmo-Simba

Table 2: The detailed model differences between Gadget-X and Gizmo-Simba. modules | Gadget-X | Gizmo-Simba
---|---|---
| Gravity and hydro solvers |
Simulation code | Gadget-3P | GIZMO (Gadget-3 based)
Gravity solver | TreePM | TreePM
Hydro solver | SPH | MFM
Kernel for hydro solver | Wendland C4 | Cubic spline
| baryon models – gas |
Gas cooling | Metal cooling table (Wiersma et al., 2009) | Grackle-3.1 library (Smith et al., 2017)
UV/X-ray background radiation | Haardt & Madau (2001) | Haardt & Madau (2012)
Self-shielding | No | Rahmati et al. (2013)
| baryon models – star |
Star formation model | Tornatore et al. (2007) | Davé et al. (2016)
Star formation rate | Gas density and temperature based |  $\mathrm{H}_{2}$-based
Star generationsa | Multiple | Single
IMF | Chabrier (2003) | Chabrier (2003)
Chemical enrichment model | 11 elements from SN-II, SN-Ia, and AGB stars | The sameb
Stellar feedback | kinetic feedback (Springel & Hernquist, 2003) | Two-phase winds with mass loading factor depending on $M_{*}$
| baryon models – BH |
BH seeding condition | $M_{FoF}>8\times 10^{11}{{\,h^{-1}{\rm{M_{\odot}}}}}\ \&\ M_{*}>1.6\times 10^{10}{{\,h^{-1}{\rm{M_{\odot}}}}}$ | Galaxy stellar mass $>3\times 10^{10}{{\,h^{-1}{\rm{M_{\odot}}}}}^{c}$
BH seed mass | $5\times 10^{6}{{\,h^{-1}{\rm{M_{\odot}}}}}$ | $10^{5}{{\,h^{-1}{\rm{M_{\odot}}}}}$
BH accretion | Bondi accretion | torque-limited and Bondi accretion models
BH feedback | Thermal feedbackd | Kinetic feedbacke \+ X-ray feedbackf
a Number of stars that per gas particle can spawn.
b Note that the detailed implementations are different, see Tornatore et al.
(2007) and Davé et al. (2016) for more information.
c Note here $M_{*}$ is the total stellar mass within the FoF halo. And another
two conditions have to be met as well: $M_{*}>0.05M_{DM}\ \&\
M_{gas}>0.10M_{*}$
d The feedback energy is coming from both mechanical and radiative modes.
e Different outflow velocities are used to mimic the jet and radiative mode
AGNs in observation.
f Thermal or thermal+kinetic feedback is used for the X-ray heating depending
on whether the surrounding gas is non-ISM or ISM, respectively.

As noted in the Introduction, we will focus on the differences between the
results from Gadget-X and Gizmo-Simba 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 Gadget-X, we refer to Cui et al. (2018b) and
references therein. Although we give descriptions of the Simba 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 Gadget-X is quite successful in reproducing the
observed gas properties and relations, while the original Simba simulation was
primarily tuned to reproduce galaxy stellar properties. Here, we follow
Simba’s process by calibrating Gizmo-Simba 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.

###  3.4 The AHF halo and Caesar galaxy catalogues

Two object catalogues are generated from the suite of the300 galaxy cluster
runs: the AHF666http://popia.ft.uam.es/AHF (Knollmann & Knebe, 2009) halo
catalogues and the Caesar 777https://github.com/dnarayanan/caesar galaxy/FoF
halo catalogues. AHF provides halo and subhalo (and thus galaxy) catalogue
generated using a spherical overdensity (SO) algorithm, while Caesar also
provide halo catalogue generated using a 3-D FoF algorithm along with a
matched galaxy catalog using a 6-D FoF. Caesar further provides a large range
of pre-computed physical and photometric properties (with and without dust
extinction) for each object.

In this paper, we use as many (uncontaminated) objects as possible (if not
specified) to do the investigation because The Three Hundred 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
$M_{500}$ with the quantities are always calculated within $R_{500}$.
Therefore, halo properties from the AHF catalogue are used to compare with the
global clusters properties from observations, while galaxy properties from
Caesar which has a 6D (in both spatial and velocity field) galaxy finder, is
used. We further match the galaxies from Caesar to the halos from AHF, by
simply taking all the Caesar galaxies within the AHF halo radius ($R_{500}$).
The BCG is selected as the most massive Caesar 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.

##  4 results

###  4.1 Baryon fractions

We first focus on the total gas and stellar components within the clusters
(mostly within $R_{500}$) in this subsection. We will also detail their
evolution and show the differences between Gadget-X and Gizmo-Simba runs.

####  4.1.1 C1: The gas and stellar fractions

Figure 1: The baryon fractions within $R_{500}$: gas fractions on the left-
hand side panel and stellar fractions on the right-hand side panel at $z=0$.
Observational and different simulation results can be crossly viewed from
legends in both figures. The symbols and/or line styles for the same reference
result are the same in both figures. Therefore, we only show them once in
either legend. The statistical results from The Three Hundred project:
Gadget-X and Gizmo-Simba, are presented with both symbols and lines. Note that
the error bars for the two runs are marking the $16^{th}-84^{th}$ percentiles.
For the statistical/fitting results, we only include the errorbar (light
shaded areas in both figures) for the most recent result – Akino et al.
(2021). The vertical dashed lines indicate the mass completeness for the
clusters from The Three Hundred project (see Cui, 2022, for details) and the
horizontal line in the left-hand-side figure is the cosmological baryon
fraction from the Planck cosmology (Planck Collaboration et al., 2016). Note
that all the fitting results only cover the region of the observed data
points. These two plots show that Gadget-X is very similar to the FABLE and
the C-EAGLE simulations in both fractions; there is little difference between
Gadget-X and Gizmo-Simba in the stellar fraction, but the gas fraction from
Gizmo-Simba shows a much steeper slope. Both fractions from the observational
data present a large scatter.

The first quantity we examine in Gizmo-Simba 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, the 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 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 Gizmo-Simba, 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 Gizmo-
Simba 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
$R_{500}$ between our simulated clusters at $z=0$, versus recent simulations
(FABLE, Henden et al. 2018 and C-EAGLE, Barnes et al. 2017b) and observational
data at $z\lesssim 0.1$ (Zhang et al., 2011; 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. (2015); Eckert et al. (2016); Dietrich et al. (2019) since
these do not provide individual cluster data. Note that other works, such as
Lin & Mohr (2004); Gonzalez et al. (2007); Andersson et al. (2011); Budzynski
et al. (2014) and also the recent works Lim et al. (2020); Chen et al. (2022)
(estimating gas fractions with SZ signal instead of X-ray) that also examined
clusters at $z\sim 0.1$ to investigate similar fractions, are not included in
this comparison due to different reasons, such as no $M_{500}$ or $M_{*}$
being available, or that the sample is dominated by $z\gg 0$ clusters.
However, it has been suggested that there is almost no redshift evolution in
both gas and stellar fractions within $z\sim 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 Gadget-X, FABLE and C-EAGLE are in very good agreement with
each other for both the median and scatters. In contrast, Gizmo-Simba shows a
steeper slope, steepening further for $M_{500}\la 10^{14}M_{\odot}$. The
difference could owe to the fact that the first three simulations employ a
thermally-based AGN feedback scheme, while Gizmo-Simba employs a kinetic
scheme (Robson & Davé, 2020); 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 Gadget-X results. Meanwhile, Dietrich et al. (2019) and Akino
et al. (2021) are in better agreement with Gizmo-Simba at
$M_{500}>10^{14}\,\rm{M_{\odot}}$; measurements to even lower masses are as
yet highly uncertain. The downturn for $M_{500}<10^{14}{{\,M_{*}}}$ in Gizmo-
Simba 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 $M_{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
$M_{500}$. Here, Gadget-X and Gizmo-Simba 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 ($M_{500}<10^{14}\,\rm{M_{\odot}}$), and thus a slope that matches
better with the individual observations shown. We note that the low resolution
of The Three Hundred will tend to suppress stellar fractions since they do not
resolve as far down the mass function as IllustrisTNG; if the shape of the
galaxy stellar mass function is not a strong function of $M_{500}$, this can
explain the constant $\sim 20\%$ offset between these models. Overall, we note
that where the observations are most robust at $M_{500}\ga
10^{14.5}\,\rm{M_{\odot}}$, 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 $R_{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 ($\sim 10^{15}M_{\odot}$) clusters all simulations
are generally in agreement with each other for both gas fraction (differences
within $\sim$5 per cent) and stellar fractions (difference within $\sim 1$ per
cent), and in reasonable agreement with observations. For Gizmo-Simba, the
latter was mostly achieved via tuning parameters as described in §3. The
models tend to differ more significantly towards lower mass halos, with Gizmo-
Simba producing particularly low gas fractions in groups, and lower stellar
fractions compared to observations of individual $10^{14}M_{\odot}$ systems.
The origin of the differences between Gadget-X and Gizmo-Simba 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 effect888Interested readers are
referred to Yang et al. (2022) on how the next-generation SZ observations can
be used to distinguish different baryon models. 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.

####  4.1.2 The evolution of the baryon fractions at the same halo mass

Figure 2: The gas fraction (left-hand-side panel) and stellar fraction (right-
hand-side panel) binned by $M_{500}$ at different redshifts. To clearly view
the fraction changes along redshift, we use linear scale for the y-axes. The
same symbols and colourful line as in Figure 1 for the two runs are adopted,
while different line styles are for different redshifts, see the legend on the
left-hand-side panel for details. We don’t include any observation data in
this busy plot because there lacks statistical results and weak or no redshift
evolution is claimed. It seems that Gizmo-Simba shows much clear redshift
evolution for both gas (much clearer at lower halo mass end) and stellar
fractions compared to Gadget-X, of which shows a weak redshift evolution in
their stellar mass fractions.

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. 2018
from observation side, or Planelles et al. 2013; Truong et al. 2018 from
simulations). This is something of a paradoc in hierarchical structure
formation models: how can the baryonic scaling relations vary steeply with
mass and still remain roughly constant in time as structures grow? In
$\Lambda$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 highlights how baryon fractions change
for a sample selected at a given mass, while the latter explicitly shows the
real 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 Gadget-X and Gizmo-Simba, 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 Gadget-X
and Gizmo-Simba owes to Gadget-X tending to have a much higher gas density in
the halo centre than Gizmo-Simba. This points to galactic feedback processes
being the primary driver of the model differences. This is corroborated by the
most significant drop occurring between $z=2.5\to 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
Gizmo-Simba that in Gadget-X down to $z=0$ in $M_{500}\la 10^{14.5}M_{\odot}$
halos.

Although for clarity we don’t show higher-redshift observation data in Figure
2, Gizmo-Simba has the closest gas fractions compared to Chiu et al. (2016a),
around 10 per cent at $z\approx 0.9$ and $M_{500}\approx 6\times
10^{14}\,\rm{M_{\odot}}$, while the gas fraction from Chiu et al. (2018)
(around 12 per cent at $z\approx 0.6$ and $M_{500}\approx 4.8\times
10^{14}\,\rm{M_{\odot}}$) lies between Gadget-X and Gizmo-Simba. 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. Gizmo-Simba predicts some evolution
over this redshift range, but at a $\la 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 galaxy group masses
$5-10\times 10^{13}\,\rm{M_{\odot}}$.

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 $M_{500}\la 10^{14.5}M_{\odot}$, and by $z=0$ they are fairly
similar. However, the evolution is quite different, with again Gizmo-Simba
showing much more evolution than Gadget-X. Neither simulation evolves much at
the massive end; we reiterate that Gizmo-Simba was calibrated to match $z=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 $z\sim
3$ (see the following subsection and §4.2.4 for details). This results in mild
but still visible redshift evolution – about 0.5 per cent from $z=1.5\to 0$.

In smaller clusters and groups, jet feedback in Gizmo-Simba is once again
implicated by the fact that the most significant change in stellar fraction
happens between $z=2.5$ and $z=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 $z=1\to 0$ the evolution is relatively modest. However,
it is still significantly more than in Gadget-X, which predicts essentially no
evolution in stellar mass fractions since $z\sim 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-$z$ stellar fraction
determinations on this plot, we comment on some results below in relation to
our predictions. Overall, the mild evolution is in agreement with Henden et
al. (2020) which also report a marginally significant change for
$M_{500}\gtrsim 3\times 10^{14}\,\rm{M_{\odot}}$ up to $z\sim 1.2$.
Nevertheless, they reported an anti redshift evolution at lower halo masses,
much stronger at group scales, opposite to Gizmo-Simba predictions. Although
the stellar fraction of 0.023 at $8\times 10^{13}\,\rm{M_{\odot}}$ with the
median redshift around 0.5 from Chiu et al. (2016b) is in good agreement (see
also Chen et al. 2021 for a similar result at lower redshift and $M_{200}$)
with Gizmo-Simba the stellar fractions from Chiu et al. (2018) – 0.0083 at
$M_{500}\approx 4.8\times 10^{14}\,\rm{M_{\odot}}$ and $z\approx 0.6$, and
from Chiu et al. (2016a) – 0.011 at $M_{500}\approx 6\times
10^{14}\,\rm{M_{\odot}}$ and $z\approx 0.9$, are much lower than both Gadget-X
and Gizmo-Simba. Controversially, Decker et al. (2021) report higher
$f_{*}\approx 0.025$ for 12 clusters at $z=0.95-1.43$, which seems in good
agreement with Gizmo-Simba. Besides these results, most observations claim
there is no clear redshift evolution of the cluster stellar fraction, (see Lin
et al., 2012, 2017; Chiu 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 $M_{500}$ means we don’t take the halo pseudo evolution into
account in either observation and simulations, so this should not have an
effect on our conclusions. Another point to keep in mind is that the high
redshift halos in zoom simulations are necessarily the progenitors of the
$z=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 mild decrease in both gas and stellar
fractions with time at lower halo mass ($M_{500}\lesssim
10^{14.5}{{\,h^{-1}{\rm{M_{\odot}}}}}$) in both Gadget-X and Gizmo-Simba.
However, the rate of decrease is much stronger in Gizmo-Simba, which shows
much higher stellar fractions at high-$z$ and much lower gas fractions at
low-$z$. This differences can be expected from Gizmo-Simba’s kinetic AGN
feedback which imparts much stronger feedback than Gadget-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.

####  4.1.3 The evolution of the cluster baryon fractions by tracking

Figure 3: The redshift evolution of gas fraction (left-hand-side panel),
stellar fraction (middle panel) and total/baryon masses (upper and lower
family of lines in the right-hand-side panel, respectively) within $R_{500}$
for three halo mass bins at $z=0$. The halos at $z=0$ are separated into the
same three halo mass bins as indicated in the legends. Only the median values
of the halos and their progenitors in the same mass bins are shown. Errorbars
are not shown because of the very dense data points and the curves are very
close to each other. Through this tracking which also goes the highest
redshift in the simulation snapshots, we can view the detailed differences
between Gadget-X and Gizmo-Simba.

Simulations provide a view that cannot be matched 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 $z\sim 1.5-2.5$ in Gizmo-Simba. To see how this is reflected in individual
halos, we can place clusters into three mass bins at $z=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 Gizmo-Simba (green curves) and Gadget-X (blue).

The gas fraction evolution is shown in the left panel of Figure 3. Gadget-X
and Gizmo-Simba 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 ($z\approx 5-6$), Gadget-X and Gizmo-
Simba reach a gas fraction near the cosmic baryon fraction value, with Gizmo-
Simba’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 $z\sim 2$, gas fractions starts to decline very quickly (more strongly
in Gizmo-Simba), until they reach a minimum at $z\sim 2-3$ (slightly later in
Gizmo-Simba). 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 Gizmo-Simba and Gadget-X are caused by the different star
formation and stellar feedback models. The net result is that the gas fraction
in Gizmo-Simba is about 2% lower than Gadget-X, persisting until today,
indicating that the sensitivity to feedback for halo gas contents is
established early on. After $z\sim 2$, gas fractions start to increase
gradually, levelling off after $z\sim 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 $z\sim 4$, the stellar fractions start
to drop, and level off after $z\sim 1$. Gizmo-Simba has higher stellar
fractions at early times versus Gadget-X; at $z=4$, Gizmo-Simba has produce
50% more stars than Gadget-X in these halos, but by $z=0$ they are the same.
This shows that Gizmo-Simba has earlier stellar formation times in groups and
clusters compared to Gadget-X. Gizmo-Simba grows early galaxies faster than
Gadget-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 Gizmo-Simba limit the
available fuel for star formation, resulting in less rapid stellar growth, and
the stellar fraction returns to meet Gadget-X by $z=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-$M_{*}$ 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
Gizmo-Simba, 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 $z\sim 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 Gizmo-Simba having lower gas fractions.

The differences in the gas fractions (and associated slight differences in the
stellar fractions) emerging at $z\sim 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 this redshift when Gizmo-Simba’s jet
mode AGN feedback turns on in these massive halos, which is what drives galaxy
quenching (Davé et al., 2019). 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. Gizmo-Simba and
Gadget-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 Gadget-X and Gizmo-Simba are
established at fairly early epochs. The gas fractions, remain are
significantly different at all redshifts with Gizmo-Simba having lower gas
fractions by $\sim 2\%$ of the cosmic baryon fraction. Gizmo-Simba 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,
Gizmo-Simba has an earlier stellar growth phase relative to Gadget-X, with 50%
more stars formed than in Gadget-X at $z\sim 3-5$. This is caused by Gizmo-
Simba 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 Gizmo-Simba curtail stellar
growth more than in Gadget-X, so that by $z\la 1$ the stellar fractions in the
two models are similar. The overall differences in the halo baryon masses are
small but noticeable by $z=0$. While both models produce a similar global
stellar fraction by $z=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.

###  4.2 BCG Properties

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

####  4.2.1 C2: The $M_{BCG,*}$ \- $M_{halo}$ relation

Figure 4: The BCG mass -$M_{500}$ relation from the Gadget-X and Gizmo-Simba
runs of The Three Hundred. From left to right panels, we show the BCG mass
estimated within different 3D apertures: 30 kpc, 50 kpc and 0.1$*R_{500}$. The
observation results from Kravtsov et al. (2018) are shown by black filled
stars. The red, blue and green lines with error bars are the median values
with $16^{th}-84^{th}$ percentiles from Gadget-X and Gizmo-Simba,
respectively. Generally, Gadget-X and Gizmo-Simba are inline with the Kravtsov
et al. (2018).

As we have discussed in section 3, the BCG – halo mass relation is one of the
three constraints for the Gizmo-Simba 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–$M_{500}$ relations from all the
uncontaminated halos in our 324 regions at $z=0$. Here, the BCG stellar masses
are estimated at three different 3D apertures: 30 kpc (left), 50 kpc (middle)
and $0.1\times R_{500}$ (right). Gizmo-Simba results are in green, Gadget-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 $M_{500}$. The
general trends are very similar regardless of the chosen aperture. However, in
detail there are some differences. Generally, Gizmo-Simba produces somewhat
lower BCG masses, particularly at $M_{500}\lesssim 10^{14}M_{\odot}$. Also,
Gadget-X’s BCG masses show a clearly larger dispersion at a fixed halo mass
bin. At all apertures, the BCG masses from Gizmo-Simba are in somewhat better
agreement with the Kravtsov et al. (2018) data. It seems that both models
predicts a little steeper slope for the relation $M_{BCG,0.1\times
R_{500}}-M_{500}$ than Kravtsov et al. (2018), which could be because we are
using 3D apertures, and the smallest halos should suffer the most projection
effects.

As shown in Henden et al. (2020), the clusters from FABLE simulation have a
systematically higher BCG stellar mass compared to observations (by
$\sim$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) indicated
that 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).

####  4.2.2 BH-BCG-halo relations

Figure 5: The BH – BCG – halo relation comparisons between the Gizmo-Simba run
and observations. From left to right panels, we show the BH mass as a function
of BCG stellar mass $M_{*}$, BCG stellar velocity dispersion $\sigma_{*}$ and
halo mass $M_{500}$, respectively. We refer to the legends in each panel for
the detailed observation result with the Gizmo-Simba run is also shown in
green dotted-dashed lines. Again, the errorbars show $16^{th}-84^{th}$
percentiles. Shaded regions are the errorbars coming from Gaspari et al.
(2019). Note that Bandara et al. (2009) – BCS thereafter, used halo mass
$M_{200}$ in their fitting result, we simply convert the $M_{200}$ to
$M_{500}$ by assuming a fixed concentration of 3. Also note that all the
fitting results only roughly cover the region of the observed data points.

All BCGs are believed to contain large supermassive black holes. As shown by
previous studies (e.g. Martizzi et al., 2012; Martizzi et al., 2014; Frigo 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 Gizmo-Simba’s BH mass $M_{\bullet}$ vs. stellar mass
$M_{*}$ (left), stellar velocity dispersion $\sigma_{*}$ (middle), and
$M_{500}$ (left) versus observations. Note here $M_{*}$ and $\sigma_{*}$ are
obtained from the Caesar 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. Gadget-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;
Savorgnan et al., 2016; Bogdán et al., 2018; Kormendy & Ho, 2013; Savorgnan &
Graham, 2016). While we don’t specifically select early-type galaxies, we
expect that our BCGs are all early-types since they are massive, quenched and
bulge dominated.

The $M_{\bullet}-M_{*}$ relation is shown on the left panel of Figure 5.
Gizmo-Simba shows good agreement with Gaspari et al. (2019) (and Kormendy & Ho
2013, not shown) at $M_{*}\approx 10^{11.2}-10^{11.7}\,\rm{M_{\odot}}$.
However, the BH masses for galaxies for $M_{*}\ga 10^{11.7}\,\rm{M_{\odot}}$
galaxies are above these observations. In contrast, Sahu et al. (2019a) found
a steeper result that is in very good agreement with Gizmo-Simba. It is beyond
the scope of this work to delve into the differences between these
observations results, so we simply note that Gizmo-Simba predicts that more
massive galaxies will follow the same trend from Sahu et al. (2019a), all the
way up to $M_{*}\approx 10^{12.5}\,\rm{M_{\odot}}$.

Looking at the middle panel of Figure 5, there is more deviation between the
Gizmo-Simba run and the observation results for the $M_{\bullet}-\sigma_{*}$
relation. At $\log\sigma_{*}\lesssim 2.6$, Gizmo-Simba 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 $M_{\bullet}-\sigma$ relation from the original
Simba 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
$\ga\log\sigma_{*}\gtrsim 2.6$, the $M_{\bullet}-\sigma_{*}$ relation from
Gizmo-Simba 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. Hence examining
the BH mass–halo mass connection in Gizmo-Simba is interesting.

Nevertheless, we argue here that the $M_{\bullet}-M_{halo}$ relation exists
for these BCGs in galaxy clusters which are mostly ellipticals.

The result for $M_{\bullet}-M_{500}$ from Gizmo-Simba is shown in the right
panel of Figure 5. We compare to Bandara et al. (2009) (BCS2009) and Gaspari
et al. (2019). Gizmo-Simba is in very good agreement with the observed ones at
group scales. However, there is a hint that at galaxy cluster scales, Gizmo-
Simba 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
commensurate growth in the black hole.

####  4.2.3 BCG colours

Figure 6: The g - r colour - magnitude diagram for the BCGs in the mass-
complete 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^{th},50^{th}$ and $84^{th}$.

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 Gizmo-Simba
and Gadget-X using the Pyloser code implemented in Caesar, which uses the age
and metallicity of each star within a galaxy to generate a spectrum based on a
simple stellar population model from
FSPS999https://github.com/cconroy20/fsps(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. Pyloser 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 $g-r$ colour as a
function of its magnitude in SDSS $r$-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 Bernardi et al., 2007) and the max-BCGs (cyan
points Koester et al., 2007), and from (Yang et al., 2018, coloured contours)
with BCGs from their group catalogue. Note that the $g-r$ colour is corrected
in rest frame for all the observation data based on Chilingarian & Zolotukhin
(2012).

The $g-r$ colour from the observed galaxies typically $\approx 0.8$, slightly
higher for Yang et al. (2018). This is well in line with the BCG colours from
Gizmo-Simba, particularly for the Roche et al. (2010) sample which extends to
overlapping magnitudes. In contrast, Gadget-X BCGs have a much bluer BCG
colour, typically $g-r\approx 0.5-0.7$. This indicates that BCG colours, and
hence their star formation histories, are more realistic in Gizmo-Simba than
in Gadget-X.

That said, for both Gizmo-Simba and Gadget-X, the BCG brightnesses extend up
farther than the data. The C4-BCG and max-BCG can reach the magnitude of
$r\approx-24.5$ comparable to the brightest galaxies in Gizmo-Simba, but
Gadget-X shows BCGs extending to brighter than $r\approx-26$. For Gizmo-Simba,
however, the mild differences can probably be attributed to sample selection,
in that The Three Hundred focuses on only the most massive clusters.
Particularly for Gadget-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.

####  4.2.4 Halo and BCG formation time

Figure 7: The correlation between the BCG age and halo formation redshift for
the mass-complete cluster sample. As indicated in the legend, different
symbols are for the BCGs from different versions of simulations and SAMs. We
further detail their BCG age distributions by showing the histograms in
redshift bins with the same colour to the right side panel.

The different BCG colours in Gizmo-Simba and Gadget-X suggest that the BCG
formation times are different in these models. The $z=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 the stars in the BCG at $z=0$ are
formed? (2) when the progenitor galaxies of the BCG are 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 as the star
formation time. The second one tightly connects with the hierarchical
structure formation theory, which can be hard to probe in observation, we
refer as the assembly formation time. We expect this latter one may be more
related to halo assembly, which will be quite similar within Gizmo-Simba and
Gadget-X.

For the first question, the BCG 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 $\sim$13 Gyr (see also Loubser et al., 2009), which
is in agreement with the SAM 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 $\sim$ 5 Gyr, which can be caused by
more recent in-situ star formation ($\sim$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 observation to roughly constrain their assembly time. Collins et al. (2009)
found several BCGs at z $\sim 1-1.5$ which have stellar masses comparable to
the most-massive galaxies in the Universe, suggesting a very early assembly
time of the BCG. More theoretical studies with 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
falls only into place typically by z $\sim$ 1.5. Similarly, the BCGs in FABLE
which are in good agreement with several observational inferences at
$z\lesssim 1$ (Henden et al., 2020), have moderate stellar mass growth, about
a factor of 1.5 from $z=1$ to 0.3 and effectively halts at $z\lesssim 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 ( z $>$ 3) rather than at low redshifts (z < 1). Thus, highly
overdense protoclusters assemble their stellar mass into brightest cluster
galaxies within $\sim$1 Gyr of evolution - producing massive blue elliptical
galaxies at high redshifts (z $>$ 1.5). Due to our limited study in this
paper, we only focus on the first question and try to see 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 $z=0$ in Gizmo-Simba (green points) and Gadget-X (blue). The star
formation time is the mean age of the stellar population, while the halo
formation time here is the commonly used $z_{\rm 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 $z\sim 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. It is has been sometimes
called "anti-hierarchical", 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 Gizmo-Simba
and Gadget-X. Clearly, Gizmo-Simba has older BCGs, with the formation redshift
peaking at $z\sim 3$, and little star formation after $z\sim 2$. The mean BCG
ages from Gadget-X is younger, as young as $\sim 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 scatter in mean stellar ages.

These results are consistent with our earlier findings that the stellar
assembly within halos is shifted towards earlier epochs in Gizmo-Simba vs.
Gadget-X. This owes to both the jet feedback, as well as Gizmo-Simba’s X-ray
feedback that is responsible for quenching star formation completely Cui et
al. (2021). In the original Simba, the SFR function at $z\sim 2-4$ extends to
very high SFRs, $>10^{3}\,\rm{M_{\odot}}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 difficult time matching these
constraints. For instance, Gadget-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 $z\sim 2-4$.

###  4.3 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., 2022;
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., 2019; Lotz et al., 2021, for example). Therefore it is important to check
how well our simulations reproduce observations of the satellite population in
groups and clusters.

####  4.3.1 C3: Satellite stellar mass function

Figure 8: Satellite stellar mass functions. Only the satellite galaxies from
the mass-complete clusters are used in this satellite galaxy stellar mass
function. We only show the median values in each stellar mass bin of the
selected clusters. While the observation results within two particular halo
mass bins are shown in black and gray thick lines with errorbars. This plot is
different to Fig. 7 in Cui et al. (2018b) by removing the other models results
except Gadget-X and adding the new Gizmo-Simba result with errorbars, which is
in very good agreement with the observation results.

We first investigate whether The Three Hundred simulations produce the correct
amount of satellites at different stellar mass bins in Figure 8. The satellite
stellar mass functions (SSMFs) for Gizmo-Simba (green) and Gadget-X (blue) are
shown, along with observations of the SSMF in clusters from Yang et al. (2018)
within two mass bins spanning much of The Three Hundred sample. Recall that
the satellite stellar mass function is used to constrain the feedback
parameters in Gizmo-Simba (cf. constraint C3), so the agreement with Gizmo-
Simba should not be considered a true success of the model, although the
tuning was only done on a single cluster. Detailed comparisons and discussions
for the other models can be found in Cui et al. (2018b).

The median satellite stellar mass function from Gizmo-Simba 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. Gizmo-
Simba and Gadget-X agree at $M_{*}\ga 10^{10}{{\,h^{-1}{\rm{M_{\odot}}}}}$,
below which Gadget-X falls below the observations. There are some interesting
features, such as Gadget-X showing a strong peak at $M_{*}\ga
10^{10.3}{{\,h^{-1}{\rm{M_{\odot}}}}}$ and Gizmo-Simba showing a steep decline
at slightly higher masses; we suspect that these arise due to the onset of
strong AGN feedback at these masses (Davé et al., 2019). This is because the
BH is seeded and AGN feedback is in action when the galaxy stellar mass reach
$M_{*}\sim 10^{10.5}{{\,h^{-1}{\rm{M_{\odot}}}}}$, 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 The Three Hundred. 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 ($\sim$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 The Three Hundred.

####  4.3.2 Satellite galaxy colour distribution

Figure 9: 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 $M_{200}\gtrsim 6\times 10^{14}{{\,h^{-1}{\rm{M_{\odot}}}}}$. Again
the three contour lines/maps mark the $16^{th},50^{th}$ and $84^{th}$
percentiles. Note here that we apply the same stellar mass cut $M_{*}\geq
10^{10}\,\rm{M_{\odot}}$ to both simulation and observed dataset for a
consistent comparison.

Similar to the colour-magnitude diagram for the BCGs, the satellite colour-
magnitude diagram gives insights into their stellar growth histories. With
Gizmo-Simba 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 $g-r$ colours vs. $r$ magnitudes are shown from
Gizmo-Simba (green contours) and Gadget-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 frame101010Note 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 ($M_{*}\geq 10^{9}{{\,h^{-1}{\rm{M_{\odot}}}}}$)
applied. By comparing to Fig. 8 in Cui et al. (2018b), Gadget-X tends to have
a bluer colour, which suggests that less massive satellite galaxies should be
dominated by red colour..

Satellite galaxy colours from Gizmo-Simba show very good agreement with Yang
et al. (2018), with again typically $g-r\approx 0.8$ just as with the BCGs.
Hence the stellar populations in both the central and satellite galaxies in
Gizmo-Simba tend to be quite old. In contrast, Gadget-X results, meanwhile,
show a somewhat bluer $g-r$ 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 Gizmo-Simba’s distribution of colours has a small
scatter compared to Yang et al. (2018). In particularly, Gizmo-Simba lacks any
bluer satellites with $g-r\la 0.7$. While these are rare in the observations,
they do occur. This could be caused by two reasons: (1) the strong AGN
feedback quenched all the 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. In the future we plan
to run higher resolution The Three Hundred runs to check whether such
resolution effects can impact these comparisons.

##  5 Conclusions and discussions

In this paper, we introduce a new suite of simulated galaxy clusters in The
Three Hundred project called Gizmo-Simba, run using a feedback model based on
the successful Simba simulation (Davé et al., 2019). We re-tune its feedback
model parameters (mostly due to the coarse resolution of The Three Hundred
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 Gadget-X runs which is also successful in reproducing many
galaxy cluster properties and relations. By comparing the two hydrodynamic in
The Three Hundred we aim to understand the origin of observational properties
of galaxy clusters under different galaxy formation models. We note that a
comparison between Gadget-X and other galaxy formation models in The Three
Hundred was done in Cui et al. (2018a). Our main results are as follows:

  * •

Similar stellar mass fractions $f_{*}(<R_{500})$ are found at $z=0$ between
Gadget-X and Gizmo-Simba, which are also in agreement with observational
results. However, their gas fractions at $z=0$ are different ($\sim 5$ per
cent in absolute mass fraction) at $M_{500}\sim 10^{14}\,\rm{M_{\odot}}$, 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 Gizmo-Simba’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 $z=0$ in both models,
Gizmo-Simba has a much stronger evolution of this fraction versus redshift
than Gadget-X, especially at high redshift $z\gtrsim 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 Gizmo-Simba 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 Gadget-X by $z=0$.

  * •

Using BCG stellar masses within fixed apertures, both Gadget-X and Gizmo-Simba
are in agreement with observations from Kravtsov et al. (2018). However, the
BCGs (identified by the 6D Caesar galaxy finder) from Gadget-X have much
brighter and bluer colour compared to Gizmo-Simba; observations agree better
with Gizmo-Simba. This is corroborated by the mean BCG ages in Gadget-X being
systematically lower than in Gizmo-Simba, with a long tail of younger BCGs.

  * •

To test Gizmo-Simba’s black hole accretion and feedback models, we compare
three BH scaling relations from Gizmo-Simba versus $M_{*}$, $\sigma_{*}$, and
$M_{500}$. Gizmo-Simba matches observations generally well, although it
predicts higher black hole masses in the most massive galaxies, a shallower
slope for the $M_{\bullet}-\sigma_{*}$ at the low-$\sigma_{*}$ end, and a
shallower $M_{\bullet}-M_{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, Gizmo-Simba shows better agreement to observation than Gadget-X at
the lower mass end, which is only partly a result of tuning. The satellite
galaxies in Gizmo-Simba agree very well with the $g-r$ vs. $r$ colour-
magnitude diagram of the SDSS galaxies, while Gadget-X satellites are somewhat
bluer. However, Gizmo-Simba shows less scatter than observations, with no
satellites bluer than $g-r\la 0.7$ whereas observations indicate some.

The gas scaling relations are compared in the appendixes. As a result of
Gizmo-Simba’s strong feedback, its temperature–mass ($T-M$) relation tends to
be a little bit higher than Gadget-X. However, in combination with its
slightly lower gas fractions, the resulting integrated Sunyaev-Zeldovich
decrement vs. mass ($Y-M$) relation is quite similar to that in Gadget-X.

The Gizmo-Simba runs in The Three Hundred 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 would like to thank The Three Hundred collaboration for all the help on
this project and this paper.

As part of The Three Hundred project, this work has received financial support
from the European Union’s Horizon 2020 Research and Innovation programme under
the Marie Sklodowskaw-Curie grant agreement number 734374, the LACEGAL project

The simulations used in this paper have been performed in the MareNostrum
Supercomputer at the Barcelona Supercomputing Center, thanks to CPU time
granted by the Red Española de Supercomputacíon. The CosmoSim database used in
this paper is a service by the Leibniz-Institute for Astrophysics Potsdam
(AIP). The MultiDark database was developed in cooperation with the Spanish
MultiDark Consolider Project CSD2009-00064.

WC and RD are supported by the Science and Technology Facilities Council
(STFC) AGP Grant ST/V000594/1. WC and YZ acknowledge the science research
grants from the China Manned Space Project with NO. CMS-CSST-2021-A01 and CMS-
CSST-2021-B01. AK, GY, DC, RC and AC are supported by the Ministerio de
Ciencia, Innovación y Universidades (MICIU/FEDER) under research grant
PGC2018-094975-C21. AK further thanks Lana Del Rey for video games. DA is
supported by the National Science Foundation Graduate Research Fellowship
under Grant No. DGE 1746045. MDP acknowledges support from Sapienza Università
di Roma thanks to Progetti di Ricerca Medi 2019, RM11916B7540DD8D. SE
acknowledges financial contribution from the contracts ASI-INAF Athena
2019-27-HH.0, “Attività di Studio per la comunità scientifica di Astrofisica
delle Alte Energie e Fisica Astroparticellare” (Accordo Attuativo ASI-INAF n.
2017-14-H.0), INAF mainstream project 1.05.01.86.10, and from the European
Union’s Horizon 2020 Programme under the AHEAD2020 project (grant agreement n.
871158). YW is supported by NSFC grant No.11733010 and NSFC grant No.11803095.
XY is supported by the national science foundation of China (Nos. 11833005,
11890692) and by the China Manned Space Project with No. CMS-CSST-2021-A02. DC
is a Ramon-Cajal researcher. KD acknowledges support by the Deutsche
Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s
Excellence Strategy - EXC-2094 - 390783311 and support through the COMPLEX
project from the European Research Council (ERC) under the European Union’s
Horizon 2020 research and innovation program grant agreement ERC-2019-AdG
882679. RH acknowledges support from the STFC through a studentship. UK
acknowledges support from the STFC through grant number RA27PN. YZ
acknowledges the support by the National Science Foundation of China
(11873038, 12173024).

## Data Availability

The results shown in this work use data from The Three Hundred galaxy clusters
sample. These data are available on request following the guidelines of The
Three Hundred 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, 2013),
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).

## References

  * Adam et al. (2018) Adam R., et al., 2018, A&A, 609, A115
  * Akino et al. (2021) Akino D., et al., 2021, arXiv e-prints, p. arXiv:2111.10080
  * Allen et al. (2011) Allen S. W., Evrard A. E., Mantz A. B., 2011, ARA&A, 49, 409
  * Anbajagane et al. (2020) Anbajagane D., Evrard A. E., Farahi A., Barnes D. J., Dolag K., McCarthy I. G., Nelson D., Pillepich A., 2020, MNRAS, 495, 686
  * Anbajagane et al. (2021) Anbajagane D., et al., 2021, arXiv e-prints, p. arXiv:2111.04778
  * Anbajagane et al. (2022) Anbajagane D., et al., 2022, MNRAS, 510, 2980
  * Andersson et al. (2011) Andersson K., et al., 2011, ApJ, 738
  * Anglés-Alcázar et al. (2013) Anglés-Alcázar D., Özel F., Davé R., 2013, ApJ, 770, 5
  * Anglés-Alcázar et al. (2015) Anglés-Alcázar D., Özel F., Davé R., Katz N., Kollmeier J. A., Oppenheimer B. D., 2015, ApJ, 800, 127
  * Anglés-Alcázar et al. (2017a) Anglés-Alcázar D., Davé R., Faucher-Giguère C.-A., Özel F., Hopkins P. F., 2017a, MNRAS, 464, 2840
  * Anglés-Alcázar et al. (2017b) Anglés-Alcázar D., Faucher-Giguère C.-A., Kereš D., Hopkins P. F., Quataert E., Murray N., 2017b, MNRAS, 470, 4698
  * Ansarifard et al. (2020) Ansarifard S., et al., 2020, A&A, 634, A113
  * Arthur et al. (2019) Arthur J., et al., 2019, MNRAS, 484, 3968
  * Astropy Collaboration et al. (2013) Astropy Collaboration et al., 2013, A&A, 558, A33
  * Astropy Collaboration et al. (2018) Astropy Collaboration et al., 2018, AJ, 156, 123
  * Bahé et al. (2017) Bahé Y. M., et al., 2017, MNRAS, 470, 4186
  * Baldi et al. (2018) Baldi A. S., De Petris M., Sembolini F., Yepes G., Cui W., Lamagna L., 2018, MNRAS, 479, 4028
  * Bandara et al. (2009) Bandara K., Crampton D., Simard L., 2009, ApJ, 704, 1135
  * Barnes et al. (2017a) Barnes D. J., Kay S. T., Henson M. A., McCarthy I. G., Schaye J., Jenkins A., 2017a, MNRAS, 465, 213
  * Barnes et al. (2017b) Barnes D. J., et al., 2017b, MNRAS, 471, 1088
  * Bassini et al. (2019) Bassini L., et al., 2019, A&A, 630, A144
  * Bassini et al. (2020) Bassini L., et al., 2020, A&A, 642, A37
  * Baxter et al. (2021) Baxter E. J., Adhikari S., Vega-Ferrero J., Cui W., Chang C., Jain B., Knebe A., 2021, MNRAS, 508, 1777
  * Beck et al. (2016) Beck A. M., et al., 2016, MNRAS, 455, 2110
  * Behroozi et al. (2013) Behroozi P. S., Wechsler R. H., Conroy C., 2013, ApJ, 770, 57
  * Behroozi et al. (2015) Behroozi P., et al., 2015, MNRAS, 454, 3020
  * Benson (2012) Benson A. J., 2012, New Astron., 17, 175
  * Bernardi et al. (2007) Bernardi M., Hyde J. B., Sheth R. K., Miller C. J., Nichol R. C., 2007, AJ, 133, 1741
  * Bernardi et al. (2011) Bernardi M., Roche N., Shankar F., Sheth R. K., 2011, MNRAS, 412, 684
  * Bogdán & Goulding (2015) Bogdán Á., Goulding A. D., 2015, ApJ, 800, 124
  * Bogdán et al. (2018) Bogdán Á., Lovisari L., Volonteri M., Dubois Y., 2018, ApJ, 852, 131
  * Böhringer et al. (2012) Böhringer H., Dolag K., Chon G., 2012, A&A, 539, A120
  * Bondi (1952) Bondi H., 1952, MNRAS, 112, 195
  * Booth & Schaye (2010) Booth C. M., Schaye J., 2010, MNRAS, 405, L1
  * Borzyszkowski et al. (2017) Borzyszkowski M., Porciani C., Romano-Díaz E., Garaldi E., 2017, MNRAS, 469, 594
  * Bryan & Norman (1998) Bryan G. L., Norman M. L., 1998, ApJ, 495, 80
  * Budzynski et al. (2014) Budzynski J. M., Koposov S. E., McCarthy I. G., Belokurov V., 2014, MNRAS, 437, 1362
  * Bulbul et al. (2019) Bulbul E., et al., 2019, ApJ, 871, 50
  * Cañas et al. (2020) Cañas R., Lagos C. d. P., Elahi P. J., Power C., Welker C., Dubois Y., Pichon C., 2020, MNRAS, 494, 4314
  * Capalbo et al. (2021) Capalbo V., De Petris M., De Luca F., Cui W., Yepes G., Knebe A., Rasia E., 2021, MNRAS, 503, 6155
  * Capalbo et al. (2022) Capalbo V., et al., 2022, in European Physical Journal Web of Conferences. p. 00008 (arXiv:2111.01411), doi:10.1051/epjconf/202225700008
  * Chabrier (2003) Chabrier G., 2003, PASP, 115, 763
  * Chen et al. (2021) Chen X., Zu Y., Shao Z., Shan H., 2021, arXiv e-prints, p. arXiv:2112.03934
  * Chen et al. (2022) Chen Z., Zhang P., Yang X., 2022, arXiv e-prints, p. arXiv:2201.12591
  * Chilingarian & Zolotukhin (2012) Chilingarian I. V., Zolotukhin I. Y., 2012, MNRAS, 419, 1727
  * Chiu et al. (2016a) Chiu I., et al., 2016a, MNRAS, 455, 258
  * Chiu et al. (2016b) Chiu I., et al., 2016b, MNRAS, 458, 379
  * Chiu et al. (2018) Chiu I., et al., 2018, MNRAS, 478, 3072
  * Choi et al. (2012) Choi E., Ostriker J. P., Naab T., Johansson P. H., 2012, ApJ, 754, 125
  * Collins et al. (2009) Collins C. A., et al., 2009, Nature, 458, 603
  * Conroy & Gunn (2010) Conroy C., Gunn J. E., 2010, ApJ, 712, 833
  * Conroy et al. (2009) Conroy C., Gunn J. E., White M., 2009, ApJ, 699, 486
  * Contreras-Santos et al. (2022) Contreras-Santos A., et al., 2022, MNRAS,
  * Cora et al. (2018) Cora S. A., et al., 2018, preprint,  (arXiv:1801.03883)
  * Croton et al. (2016) Croton D. J., et al., 2016, ApJS, 222, 22
  * Cui (2022) Cui W., 2022, in European Physical Journal Web of Conferences. p. 00011 (arXiv:2111.01889), doi:10.1051/epjconf/202225700011
  * Cui et al. (2014) Cui W., et al., 2014, MNRAS, 437, 816
  * Cui et al. (2017) Cui W., Power C., Borgani S., Knebe A., Lewis G. F., Murante G., Poole G. B., 2017, MNRAS, 464, 2502
  * Cui et al. (2018a) Cui W., Knebe A., Yepes G., Yang X., Borgani S., Kang X., Power C., Staveley-Smith L., 2018a, MNRAS, 473, 68
  * Cui et al. (2018b) Cui W., et al., 2018b, MNRAS, 480, 2898
  * Cui et al. (2021) Cui W., Davé R., Peacock J. A., Anglés-Alcázar D., Yang X., 2021, Nature Astronomy,
  * Davé et al. (2016) Davé R., Thompson R., Hopkins P. F., 2016, MNRAS, 462, 3265
  * Davé et al. (2019) Davé R., Anglés-Alcázar D., Narayanan D., Li Q., Rafieferantsoa M. H., Appleby S., 2019, MNRAS, 486, 2827
  * De Luca et al. (2021) De Luca F., De Petris M., Yepes G., Cui W., Knebe A., Rasia E., 2021, MNRAS, 504, 5383
  * De Lucia & Blaizot (2007) De Lucia G., Blaizot J., 2007, MNRAS, 375, 2
  * Decker et al. (2021) Decker B., et al., 2021, arXiv e-prints, p. arXiv:2112.12239
  * Dekel & Birnboim (2006) Dekel A., Birnboim Y., 2006, MNRAS, 368, 2
  * Díaz-García et al. (2021) Díaz-García L. A., Umetsu K., Rasia E., Cui W., Meneghetti M., 2021, arXiv e-prints, p. arXiv:2108.13498
  * Dietrich et al. (2019) Dietrich J. P., et al., 2019, MNRAS, 483, 2871
  * Dolag et al. (2010) Dolag K., Murante G., Borgani S., 2010, MNRAS, 405, 1544
  * Dolag et al. (2016) Dolag K., Komatsu E., Sunyaev R., 2016, MNRAS, 463, 1797
  * Dullo et al. (2021) Dullo B. T., Gil de Paz A., Knapen J. H., 2021, ApJ, 908, 134
  * Eckert et al. (2016) Eckert D., et al., 2016, A&A, 592, A12
  * Eckert et al. (2021) Eckert D., Gaspari M., Gastaldello F., Le Brun A. M. C., O’Sullivan E., 2021, Universe, 7, 142
  * Edwards et al. (2020) Edwards L. O. V., et al., 2020, MNRAS, 491, 2617
  * Fabian (1994) Fabian A. C., 1994, ARA&A, 32, 277
  * Ferragamo et al. (2022) Ferragamo A., De Petris M., Yepes G., Macías-Pérez J., Cui W., Jiménez-Muñoz A., 2022, in European Physical Journal Web of Conferences. p. 00018 (arXiv:2111.01724), doi:10.1051/epjconf/202225700018
  * Frigo et al. (2019) Frigo M., Naab T., Hirschmann M., Choi E., Somerville R. S., Krajnovic D., Davé R., Cappellari M., 2019, MNRAS, 489, 2702
  * Gao et al. (2012) Gao L., Navarro J. F., Frenk C. S., Jenkins A., Springel V., White S. D. M., 2012, MNRAS, 425, 2169
  * Gaspari et al. (2019) Gaspari M., et al., 2019, ApJ, 884, 169
  * Gianfagna et al. (2022) Gianfagna G., Rasia E., Cui W., De Petris M., Yepes G., 2022, in European Physical Journal Web of Conferences. p. 00020 (arXiv:2111.01903), doi:10.1051/epjconf/202225700020
  * Gonzalez et al. (2007) Gonzalez A. H., Zaritsky D., Zabludoff A. I., 2007, ApJ, 666, 147
  * Gonzalez et al. (2013) Gonzalez A. H., Sivanandam S., Zabludoff A. I., Zaritsky D., 2013, ApJ, 778, 14
  * Haardt & Madau (2001) Haardt F., Madau P., 2001, in Neumann D. M., Tran J. T. V., eds, Clusters of Galaxies and the High Redshift Universe Observed in X-rays. (arXiv:astro-ph/0106018)
  * Haardt & Madau (2012) Haardt F., Madau P., 2012, ApJ, 746, 125
  * Haggar et al. (2020) Haggar R., Gray M. E., Pearce F. R., Knebe A., Cui W., Mostoghiu R., Yepes G., 2020, MNRAS, 492, 6074
  * Haggar et al. (2021) Haggar R., Pearce F. R., Gray M. E., Knebe A., Yepes G., 2021, MNRAS, 502, 1191
  * Henden et al. (2018) Henden N. A., Puchwein E., Shen S., Sijacki D., 2018, MNRAS, 479, 5385
  * Henden et al. (2020) Henden N. A., Puchwein E., Sijacki D., 2020, MNRAS, 498, 2114
  * Hopkins (2015) Hopkins P. F., 2015, MNRAS, 450, 53
  * Hopkins & Quataert (2011) Hopkins P. F., Quataert E., 2011, MNRAS, 415, 1027
  * Hunter (2007) Hunter J. D., 2007, Computing In Science & Engineering, 9, 90
  * Iwamoto et al. (1999) Iwamoto K., Brachwitz F., Nomoto K., Kishimoto N., Umeda H., Hix W. R., Thielemann F.-K., 1999, ApJS, 125, 439
  * Kereš et al. (2005) Kereš D., Katz N., Weinberg D. H., Davé R., 2005, MNRAS, 363, 2
  * Klypin et al. (2016) Klypin A., Yepes G., Gottlöber S., Prada F., Heß S., 2016, MNRAS, 457, 4340
  * Knebe et al. (2018) Knebe A., et al., 2018, MNRAS, 474, 5206
  * Knebe et al. (2020) Knebe A., et al., 2020, MNRAS, 495, 3002
  * Knollmann & Knebe (2009) Knollmann S. R., Knebe A., 2009, ApJS, 182, 608
  * Koester et al. (2007) Koester B. P., et al., 2007, ApJ, 660, 239
  * Kormendy & Bender (2011) Kormendy J., Bender R., 2011, Nature, 469, 377
  * Kormendy & Ho (2013) Kormendy J., Ho L. C., 2013, ARA&A, 51, 511
  * Kotecha et al. (2021) Kotecha S., et al., 2021, arXiv e-prints, p. arXiv:2110.13419
  * Kravtsov & Borgani (2012) Kravtsov A. V., Borgani S., 2012, ARA&A, 50, 353
  * Kravtsov et al. (2018) Kravtsov A. V., Vikhlinin A. A., Meshcheryakov A. V., 2018, Astronomy Letters, 44, 8
  * Krumholz & Gnedin (2011) Krumholz M. R., Gnedin N. Y., 2011, ApJ, 729, 36
  * Krumholz et al. (2009) Krumholz M. R., McKee C. F., Tumlinson J., 2009, ApJ, 699, 850
  * Kuchner et al. (2020) Kuchner U., et al., 2020, MNRAS, 494, 5473
  * Kuchner et al. (2021a) Kuchner U., et al., 2021a, arXiv e-prints, p. arXiv:2111.11467
  * Kuchner et al. (2021b) Kuchner U., et al., 2021b, MNRAS, 503, 2065
  * Laganá et al. (2011) Laganá T. F., Zhang Y.-Y., Reiprich T. H., Schneider P., 2011, ApJ, 743, 13
  * Li et al. (2019) Li Q., Narayanan D., Davé R., 2019, MNRAS, 490, 1425
  * Li et al. (2020) Li Q., et al., 2020, MNRAS, 495, 2930
  * Li et al. (2021) Li Q., Han J., Wang W., Cui W., Li Z., Yang X., 2021, MNRAS, 505, 3907
  * Lim et al. (2020) Lim S. H., Mo H. J., Wang H., Yang X., 2020, ApJ, 889, 48
  * Lin & Mohr (2004) Lin Y.-T., Mohr J. J., 2004, ApJ, 617, 879
  * Lin et al. (2012) Lin Y.-T., Stanford S. A., Eisenhardt P. R. M., Vikhlinin A., Maughan B. J., Kravtsov A., 2012, ApJ, 745, L3
  * Lin et al. (2017) Lin Y.-T., et al., 2017, ApJ, 851, 139
  * Lotz et al. (2019) Lotz M., Remus R.-S., Dolag K., Biviano A., Burkert A., 2019, MNRAS, 488, 5370
  * Lotz et al. (2021) Lotz M., Dolag K., Remus R.-S., Burkert A., 2021, MNRAS, 506, 4516
  * Loubser et al. (2009) Loubser S. I., Sánchez-Blázquez P., Sansom A. E., Soechting I. K., 2009, MNRAS, 398, 133
  * Lovell et al. (2021) Lovell C. C., Geach J. E., Davé R., Narayanan D., Li Q., 2021, MNRAS, 502, 772
  * Lovisari et al. (2015) Lovisari L., Reiprich T. H., Schellenberger G., 2015, A&A, 573, A118
  * Lovisari et al. (2021) Lovisari L., Ettori S., Gaspari M., Giles P. A., 2021, Universe, 7, 139
  * Lustig et al. (2022) Lustig P., et al., 2022, arXiv e-prints, p. arXiv:2201.09068
  * Marasco et al. (2021) Marasco A., Cresci G., Posti L., Fraternali F., Mannucci F., Marconi A., Belfiore F., Fall S. M., 2021, MNRAS, 507, 4274
  * Martizzi et al. (2012) Martizzi D., Teyssier R., Moore B., 2012, MNRAS, 420, 2859
  * Martizzi et al. (2014) Martizzi D., Jimmy Teyssier R., Moore B., 2014, MNRAS, 443, 1500
  * McCarthy et al. (2017) McCarthy I. G., Schaye J., Bird S., Le Brun A. M. C., 2017, MNRAS, 465, 2936
  * McConnell & Ma (2013) McConnell N. J., Ma C.-P., 2013, ApJ, 764, 184
  * McGaugh (2012) McGaugh S. S., 2012, AJ, 143, 40
  * Millman & Aivazis (2011) Millman K. J., Aivazis M., 2011, Computing in Science and Engineering, 13, 9
  * Mostoghiu et al. (2019) Mostoghiu R., Knebe A., Cui W., Pearce F. R., Yepes G., Power C., Dave R., Arth A., 2019, MNRAS, 483, 3390
  * Mostoghiu et al. (2021a) Mostoghiu R., et al., 2021a, MNRAS, 501, 5029
  * Mostoghiu et al. (2021b) Mostoghiu R., et al., 2021b, A&A, 652, A10
  * Munari et al. (2013) Munari E., Biviano A., Borgani S., Murante G., Fabjan D., 2013, MNRAS, 430, 2638
  * Murante et al. (2007) Murante G., Giovalli M., Gerhard O., Arnaboldi M., Borgani S., Dolag K., 2007, MNRAS, 377, 2
  * Muratov et al. (2015) Muratov A. L., Kereš D., Faucher-Giguère C.-A., Hopkins P. F., Quataert E., Murray N., 2015, MNRAS, 454, 2691
  * Nagarajan et al. (2018) Nagarajan A., et al., 2018, preprint,  (arXiv:1804.03671)
  * Nandra et al. (2013) Nandra K., et al., 2013, arXiv e-prints, p. arXiv:1306.2307
  * Nomoto et al. (2006) Nomoto K., Tominaga N., Umeda H., Kobayashi C., Maeda K., 2006, Nuclear Phys. A, 777, 424
  * Ogiya et al. (2022) Ogiya G., van den Bosch F. C., Burkert A., 2022, MNRAS, 510, 2724
  * Oliphant (2007) Oliphant T. E., 2007, Computing in Science and Engg., 9, 10
  * Oppenheimer & Davé (2006) Oppenheimer B. D., Davé R., 2006, MNRAS, 373, 1265
  * Oppenheimer et al. (2021) Oppenheimer B. D., Babul A., Bahé Y., Butsky I. S., McCarthy I. G., 2021, Universe, 7, 209
  * Pérez & Granger (2007) Pérez F., Granger B. E., 2007, Computing in Science and Engineering, 9, 21
  * Pillepich et al. (2018) Pillepich A., et al., 2018, MNRAS, 475, 648
  * Planck Collaboration et al. (2014) Planck Collaboration et al., 2014, A&A, 571, A20
  * Planck Collaboration et al. (2016) Planck Collaboration et al., 2016, A&A, 594, A13
  * Planelles et al. (2013) Planelles S., Borgani S., Dolag K., Ettori S., Fabjan D., Murante G., Tornatore L., 2013, MNRAS, 431, 1487
  * Puchwein et al. (2010) Puchwein E., Springel V., Sijacki D., Dolag K., 2010, MNRAS, 406, 936
  * Ragone-Figueroa et al. (2013) Ragone-Figueroa C., Granato G. L., Murante G., Borgani S., Cui W., 2013, MNRAS, 436, 1750
  * Ragone-Figueroa et al. (2018) Ragone-Figueroa C., Granato G. L., Ferraro M. E., Murante G., Biffi V., Borgani S., Planelles S., Rasia E., 2018, preprint,  (arXiv:1803.08049)
  * Rahmati et al. (2013) Rahmati A., Pawlik A. H., Raičević M., Schaye J., 2013, MNRAS, 430, 2427
  * Rasia et al. (2015) Rasia E., et al., 2015, ApJ, 813, L17
  * Rennehan et al. (2020) Rennehan D., Babul A., Hayward C. C., Bottrell C., Hani M. H., Chapman S. C., 2020, MNRAS, 493, 4607
  * Robson & Davé (2020) Robson D., Davé R., 2020, MNRAS, 498, 3061
  * Robson & Davé (2021) Robson D., Davé R., 2021, arXiv e-prints, p. arXiv:2107.01206
  * Roche et al. (2010) Roche N., Bernardi M., Hyde J., 2010, MNRAS, 407, 1231
  * Rost et al. (2021) Rost A., et al., 2021, MNRAS, 502, 714
  * Rudick et al. (2011) Rudick C. S., Mihos J. C., McBride C. K., 2011, ApJ, 732, 48
  * Sabra et al. (2015) Sabra B. M., Saliba C., Abi Akl M., Chahine G., 2015, ApJ, 803, 5
  * Sahu et al. (2019a) Sahu N., Graham A. W., Davis B. L., 2019a, ApJ, 876, 155
  * Sahu et al. (2019b) Sahu N., Graham A. W., Davis B. L., 2019b, ApJ, 887, 10
  * Sanderson et al. (2013) Sanderson A. J. R., O’Sullivan E., Ponman T. J., Gonzalez A. H., Sivanandam S., Zabludoff A. I., Zaritsky D., 2013, MNRAS, 429, 3288
  * Savorgnan & Graham (2016) Savorgnan G. A. D., Graham A. W., 2016, ApJS, 222, 10
  * Savorgnan et al. (2016) Savorgnan G. A. D., Graham A. W., Marconi A., Sani E., 2016, ApJ, 817, 21
  * Sayers et al. (2021) Sayers J., Sereno M., Ettori S., Rasia E., Cui W., Golwala S., Umetsu K., Yepes G., 2021, MNRAS, 505, 4338
  * Schaye et al. (2015) Schaye J., et al., 2015, MNRAS, 446, 521
  * Schmidt (1959) Schmidt M., 1959, ApJ, 129, 243
  * Sembolini et al. (2013) Sembolini F., Yepes G., De Petris M., Gottlöber S., Lamagna L., Comis B., 2013, MNRAS, 429, 323
  * Sembolini et al. (2016) Sembolini F., et al., 2016, MNRAS, 457, 4063
  * Sereno et al. (2021) Sereno M., Lovisari L., Cui W., Schellenberger G., 2021, MNRAS, 507, 5214
  * Smith et al. (2017) Smith B. D., et al., 2017, MNRAS, 466, 2217
  * Somerville & Davé (2015) Somerville R. S., Davé R., 2015, ARA&A, 53, 51
  * Sorini et al. (2021) Sorini D., Dave R., Cui W., Appleby S., 2021, arXiv e-prints, p. arXiv:2111.13708
  * Springel & Hernquist (2003) Springel V., Hernquist L., 2003, MNRAS, 339, 289
  * Sunyaev & Zeldovich (1970) Sunyaev R. A., Zeldovich Y. B., 1970, Ap&SS, 7, 3
  * Tamosiunas et al. (2021) Tamosiunas A., Briddon C., Burrage C., Cui W., Moss A., 2021, arXiv e-prints, p. arXiv:2108.10364
  * The Lynx Team (2018) The Lynx Team 2018, arXiv e-prints, p. arXiv:1809.09642
  * Thomas et al. (2019) Thomas N., Davé R., Anglés-Alcázar D., Jarvis M., 2019, MNRAS, 487, 5764
  * Tornatore et al. (2007) Tornatore L., Borgani S., Dolag K., Matteucci F., 2007, MNRAS, 382, 1050
  * Truong et al. (2018) Truong N., et al., 2018, MNRAS, 474, 4089
  * Umetsu et al. (2020) Umetsu K., et al., 2020, ApJ, 890, 148
  * Vega-Ferrero et al. (2021) Vega-Ferrero J., Dana J. M., Diego J. M., Yepes G., Cui W., Meneghetti M., 2021, MNRAS, 500, 247
  * Vikhlinin et al. (2006) Vikhlinin A., Kravtsov A., Forman W., Jones C., Markevitch M., Murray S. S., Van Speybroeck L., 2006, ApJ, 640, 691
  * Vikhlinin et al. (2009) Vikhlinin A., et al., 2009, ApJ, 692, 1033
  * Vulcani et al. (2011) Vulcani B., et al., 2011, MNRAS, 412, 246
  * Walker et al. (2019) Walker S., et al., 2019, Space Sci. Rev., 215, 7
  * Wang & White (2012) Wang W., White S. D. M., 2012, MNRAS, 424, 2574
  * Wang et al. (2018) Wang Y., et al., 2018, ApJ, 868, 130
  * Wiersma et al. (2009) Wiersma R. P. C., Schaye J., Smith B. D., 2009, MNRAS, 393, 99
  * Wu et al. (2013) Wu H.-Y., Hahn O., Wechsler R. H., Mao Y.-Y., Behroozi P. S., 2013, ApJ, 763, 70
  * Yang et al. (2009) Yang X., Mo H. J., van den Bosch F. C., 2009, ApJ, 695, 900
  * Yang et al. (2013) Yang X., Mo H. J., van den Bosch F. C., Bonaca A., Li S., Lu Y., Lu Y., Lu Z., 2013, ApJ, 770, 115
  * Yang et al. (2018) Yang X., et al., 2018, ApJ, 860, 30
  * Yang et al. (2022) Yang T., Cai Y.-C., Cui W., Davé R., Peacock J. A., Sorini D., 2022, arXiv e-prints, p. arXiv:2202.11430
  * Zhang et al. (2011) Zhang Y.-Y., Laganá T. F., Pierini D., Puchwein E., Schneider P., Reiprich T. H., 2011, A&A, 535, A78
  * Zhang et al. (2021) Zhang Z., Wang H., Luo W., Zhang J., Mo H. J., Jing Y., Yang X., Li H., 2021, arXiv e-prints, p. arXiv:2112.04777
  * van den Bosch (2016) van den Bosch R. C. E., 2016, ApJ, 831, 134
  * van der Walt et al. (2011) van der Walt S., Colbert S. C., Varoquaux G., 2011, CoRR, abs/1102.1523

##  Appendix A The gas T-M relation

In this and the following appendixes, we briefly check the gas properties from
Gizmo-Simba by comparing the gas T-M (Appendix A) and Y-M (Appendix B) scaling
relations (see Lovisari et al., 2021, for a recent review). The T-M relation
will test the gas temperature in Gizmo-Simba is affected by its strong
feedback or not; while the second relation is a simply check to see if the
differences in gas fraction and temperature can be directly viewed from the SZ
observations.

Figure 10: The gas mass-weighed temperature $T_{500}$ \- halo mass $M_{500}$
relation. As indicated in the legend at top-left corner, data points from
Gadget-X and Gizmo-Simba 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 Gadget-X and Gizmo-
Simba are presented by blue dotted and lime dashed lines, respectively. The
thick solid black line shows the self-similar relation $\log T_{500}\propto
2/3\log M_{500}$ is derived from self-similarity cluster approximations (see
Böhringer et al., 2012, for example).

Following Cui et al. (2018b), we calculate the mass-weighted gas temperature
within $R_{500}$ and add the scaling relation from Gizmo-Simba in Figure 10.
At the most cluster mass range, Gizmo-Simba is in line with Gadget-X and the
observational results, albeit a slightly higher temperature (see the zoomed-in
panel for details). This suggests that the gas heating inside these most
massive clusters is reasonable. However, the median data points from Gizmo-
Simba at the group halo mass range, $M_{500}<10^{14}\,\rm{M_{\odot}}$, seem to
be a little bit higher than the other results albeit a very large scatter.
This can be related to the strong AGN feedback which is also efficiently
expelling the gas out of $R_{500}$ as shown in Figure 1. By applying a
weighted fitting to $\log T_{500}=A+B*\log\frac{M_{500}}{6\times
10^{14}\,\rm{M_{\odot}}}$, in which the weight is given by the completeness of
the halo number, Gizmo-Simba ($B=0.54$) seems deviating from the self-similar
slope much further than Gadget-X ($B=0.61$). This is mainly because the high
temperature in Gizmo-Simba at lower halo mass. Note that the median data
points deviating from the fitting line at the lower halo mass 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 Gizmo-Simba is adopted (see Cui (2022) for the changes).
This indicates the fitting results are very sensitive to the weights.

##  Appendix B The gas Y-M relation

Figure 11: The updated Y-M relation from Gizmo-Simba. As indicated in the
legend, two fitting results (Planck Collaboration et al., 2014; Nagarajan et
al., 2018) from observations are shown by brown solid and black dotted lines,
respectively. While the weighted fitting results with the parameters from the
two hydro-simulations are shown by blue dashed and green dotted-dashed lines
as before.

Besides the gas fraction and temperature, it would be also interesting to see
how the overall simulated signal, the integrated SZ $Y_{500}$ value in this
case, is comparable to the one which can be directly obtained from
observation. Thus, we apply the PyMSZ package111111Publicly available at
https://github.com/weiguangcui/pymsz, see Baldi et al. (2018) for the details
of generating the kinetic SZ signal with this package. to Gizmo-Simba
simulated halos to generate mock $y-$maps along z-direction. Then, $Y_{500}$
the total $y$ values within $R_{500}$. Note that the ling-of-sight depth is
set to $2\times R_{500}$ to be consistent with Gadget-X, and we refer to Yang
et al. (2022) for the projection effect.

As shown in Figure 11, the median data points from Gizmo-Simba are comparable
to the ones from Gadget-X. This is not surprising as the SZ$-y$ signal is
proportional to gas density times temperature (Sunyaev & Zeldovich, 1970, and
references therein). The lower gas fraction in Gizmo-Simba 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 Gizmo-Simba at the low
halo mass range are also shown in this $Y_{500}-M_{500}$ relation. We also
apply a weighted fitting of this $Y_{500}-M_{500}$ relation to $\log
d_{A}^{2}Y_{500}=A+B*\log\frac{M_{500}}{6\times 10^{14}\,\rm{M_{\odot}}}$ for
the two hydro runs, where $d_{A}$ is the angular diameter distance. It is not
surprise, as it is consistent with these data points, to see such a good
agreement on the fitting parameters between Gizmo-Simba and Gadget-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 self-similar value. Adopting
the weighted fitting makes the slope dominated by these massive halos and
results in a value close to 5/3. While at the low mass end, the slope from the
data points is much steeper in both Gadget-X and Gizmo-Simba. If we do an
unweighted fitting to these data points, we get a stepper slope which is
comparable to the result from Planck. Actually, the different slopes for
$Y_{500}-M_{500}$ relation at different mass range have been found in several
works using different hydrodynamic simulations, for example, Yang et al.
(2022). Therefore, it would be better to use the fitting function with halo-
mass dependent slopes to overcome this problem.