ELUCID. VIII. Simulating the Coma Galaxy Cluster to Calibrate Model and Understand Feedback

We conducted an investigation of the Coma cluster of galaxies by running a series of constrained hydrodynamic simulations with GIZMO-SIMBA and GADGET-3 based on initial conditions reconstructed from the SDSS survey volume in the ELUCID project. We compared simulation predictions and observations for galaxies, intracluster medium (ICM) and intergalactic medium (IGM) in and around the Coma cluster to constrain galaxy formation physics. Our results demonstrate that this type of constrained investigation allows us to probe in more detail the implemented physical processes, because the comparison between simulations and observations is free of cosmic variance and hence can be conducted in a “one-to-one” manner. We found that an increase in the earlier star formation rate and the supernova feedback of the original GIZMO-SIMBA model is needed to match observational data on stellar, interstellar medium, and ICM metallicity. The simulations without active galactic nucleus (AGN) feedback can well reproduce the observational ICM electron density, temperature, and entropy profiles, ICM substructures, and the IGM temperature–density relation, while the ones with AGN feedback usually fail. However, one requires something like AGN feedback to reproduce a sufficiently large population of quiescent galaxies, particularly in low-density regions. The constrained simulations of the Coma cluster thus provide a test bed to understand processes that drive galaxy formation and evolution.


Introduction
Observed galaxies are diverse in mass, color, morphology, size, star formation, nuclear activity, gas content, metallicity, and environment.These properties are usually related to each other, indicating that a complex network of processes drives the formation and evolution of galaxies.Various approaches have been used to understand these processes, and some key aspects have been identified, including the growth of the cosmic web and dark matter halos, gas cooling and condensation, ram pressure and tidal stripping, galaxy merger/interaction and instability, star formation and evolution, and feedback from stars and active galactic nuclei (AGNs; see, e.g., Mo et al. 2010).Many of the processes such as feedback, although originating from small scales, are expected to have significant impacts on the gas on larger scales, such as the intergalactic medium (IGM), the intracluster medium (ICM) in clusters of galaxies, as well as the circumgalactic medium (CGM; Cen & Ostriker 2006;Fabian 2012;Heckman & Best 2014;Somerville & Davé 2015;Tumlinson et al. 2017;Nelson et al. 2019a;Cui et al. 2021;Eckert et al. 2021;Boselli et al. 2022;Yang et al. 2024).Such a coupling between large and small scales makes it difficult to model the underlying processes accurately on the one hand, but provides the link between theoretical models and astronomical observations on the other.
Cosmological hydrodynamic simulations can, in principle, implement many of the relevant processes.However, owing to limits on the numerical resolution and computational power, some of the key processes can only be implemented as subgrid prescriptions with parameters that need to be calibrated by observational data and/or higher-resolution simulations.These simulations have been quite successful in reproducing many properties of the gas and stellar components of the universe at different cosmic epochs, such as the stellar mass/luminosity functions, the star formation rates (SFRs) and gas content of galaxies, the stellar mass-supermassive black hole (SMBH) mass relation, and stellar mass-halo mass relation (e.g., Dubois et al. 2014;Hirschmann et al. 2014;Hopkins et al. 2014;Vogelsberger et al. 2014;Schaye et al. 2015;Wang et al. 2015;Pillepich et al. 2018;Davé et al. 2019;Cui et al. 2022).
However, these apparently "successful" models usually make different assumptions for the same subgrid physics and adopt different prescriptions to treat the same physical processes, indicating that the details of these processes are still poorly understood and modeled.
Clearly, to differentiate different hypotheses, realistic models and reliable calibrations of model parameters are needed (e.g., Davé et al. 2019).Traditionally, free parameters describing the subgrid physics are calibrated using small, high-resolution simulations.For example, the galaxy model for the Illustris simulations (Vogelsberger et al. 2014), which the IllustrisTNG simulations (Weinberger et al. 2017;Pillepich et al. 2018;Springel et al. 2018;Nelson et al. 2019b) were also based on, was calibrated using simulations with a side length of 25 h −1 Mpc (Vogelsberger et al. 2013), and the EAGLE galaxy formation model was calibrated and tested using simulations with side lengths of 50 and 25 h −1 Mpc, respectively (Crain et al. 2015).On a scale of 25 h −1 Mpc, the typical density fluctuation at z = 0 is already greater than 70% (Somerville et al. 2004;Chen et al. 2019), indicating that cosmic variance is a significant issue affecting the accuracy of the calibration.Cosmic variance also affects observational statistics used for the calibration of models.For example, Chen et al. (2019) found that the stellar mass function (SMF) of Sloan Digital Sky Survey (SDSS) galaxies is significantly underestimated for faint galaxies that are observed in relatively small volumes.
On the other hand, the potential of observational data has not yet been fully exploited in the model-data comparison.Statistical quantities of galaxies commonly used for model calibrations are the results of various mechanisms on various scales, and they alone may not be able to provide full information about the underlying processes.Indeed, to disentangle different processes, it is necessary to compare observations and simulations in different mass components in different environments separately.For example, the dominant process driving the evolution of the gas on large scales (1 Mpc) is gravitational and can now be modeled accurately by numerical simulations.However, processes operating on smaller scales, in particular those related to feedback, which are difficult to simulate accurately, can also affect gas on large scales.Thus, the gas properties on large scales provide an important avenue to understand the underlying feedback processes.Unfortunately, few studies in the literature have directly used observations of the gas distribution on large scales to constrain galaxy formation models.Second, conventional simulations, which assume random phases of Fourier modes in a relatively small simulation box, may not be able to match the ecosystems for galaxies in the observed volume, owing to cosmic variance as mentioned above.In contrast to traditional simulations, constrained simulations (CSs hereafter) are designed to reproduce the real large-scale structures found in a given volume of the universe, which allows simulated galaxies and structures to reside in the same environments and ecosystems as their observed counterparts, making the comparison between simulation results and observational data less affected by potential cosmic variance and hence more accurate.For example, one may select a particularly interesting volume, where rich multiband observational data are available, to carry out CSs, and use comparisons between simulations and observations to calibrate models of galaxy formation.Since such "one-to-one" comparisons are, in principle, free of the cosmic variance, the results are thus much more powerful in probing the underlying physical processes.In addition, such CSs also reconstruct the formation histories of the objects in question, which allows us to investigate the effects of history on the "remnants" we observe today.
The initial condition (IC) of a CS is usually reconstructed from galaxy redshift surveys and/or peculiar velocity surveys, and many reconstruction methods have been developed (e.g., Hoffman & Ribak 1991;Nusser & Dekel 1992;Frisch et al. 2002;Klypin et al. 2003;Jasche & Wandelt 2013;Kitaura 2013;Wang et al. 2013Wang et al. , 2014;;Sorce et al. 2016;Modi et al. 2018;Bos et al. 2019;Horowitz et al. 2019;Kitaura et al. 2021;Li et al. 2022b).Wang et al. (2014) developed a method that implements particle-mesh dynamics to predict the final density field from a given initial density field, and uses a Hamiltonian Monte Carlo Markov Chain to obtain the posterior model parameters (the set of Fourier modes representing the initial density field).Their tests based on N-body simulations show that the method can accurately recover the formation history of large-scale structures, even at scales of ∼2 h −1 Mpc.This method has been successfully applied to the SDSS volume, as shown in Wang et al. (2016).
The Coma galaxy cluster (CM) is one of the most interesting structures in the local Universe suitable for our "one-to-one" study.This cluster has attracted a lot of attention because it is one of the most massive, nearby galaxy clusters with a mass of M 200c ≈ 6.2 × 10 14 h −1 M e (Okabe et al. 2014) and a redshift of 0.0241 (Yang et al. 2007).Large amounts of high-quality observational data have been accumulated for the Coma cluster, from the radio band all the way to the X-ray band (Brown & Rudnick 2011;Matsushita 2011;Petropoulou et al. 2012;Akamatsu et al. 2013;Ogrean & Brüggen 2013;Planck Collaboration et al. 2013;Simionescu et al. 2013;Mirakhor & Walker 2020;Churazov et al. 2021).These observations can be used to understand the details of different mass components in and around the Coma cluster.For example, a radio relic around the virial radius, which is likely caused by a merger shock, has been detected by radio observations (Brown & Rudnick 2011); a strong bow-like shock in the inner region of the cluster has been revealed by X-ray observations (Churazov et al. 2021); and optical observations show that the Coma cluster is connected to several massive filaments and contains two bright central galaxies (BCGs).All these provide important information about the formation and evolution of this particular cluster and about the underlying processes that drive the evolution of galaxies and gas components in general.
CSs have already been used to study the Coma cluster.Malavasi et al. (2023) used CSs to study the connection between the Coma cluster and the surrounding filaments, such as their spatial configuration, kinematic state, and evolution.Dolag et al. (2016) studied the thermal Sunyaev-Zeldovich (SZ) effect in their simulated Coma cluster and found that the predicted Compton y-profile is slightly lower than that observed.Based on the ELUCID reconstruction (Wang et al. 2014(Wang et al. , 2016)), Li et al. (2022a) found that their predicted Compton y and X-ray brightness profiles are in good agreement with the observed profiles within about half of the virial radius.In particular, their CS can reproduce bow-like shocks in the same locations as those revealed in observations, and these shocks are associated with recent mergers (see also Zhang et al. 2019).The CS used by Li et al. (2022a) does not include AGN feedback, commonly believed to have a significant impact on both galaxies and the ICM, and hence, the problem needs to be revisited.
In this paper, we use ICs obtained from the ELUCID project (Wang et al. 2014(Wang et al. , 2016) ) to carry out hydrodynamic simulations in a volume around the Coma cluster.These simulations implement different galaxy formation models, so that we can see directly how galaxy formation physics is reflected in the predicted galaxy population and gas components in different environments.This paper is organized as follows.In Section 2, we present observational data for the Coma cluster, including the galaxy catalog and the reduced data for the interstellar medium (ISM) and ICM.We also introduce an observational result for the IGM obtained in the SDSS region.In Section 3, we describe the ICs for the CSs, the two simulation codes used, GIZMO (GZ)-SIMBA and GADGET-3 (G3), and the method to identify halos and galaxies.We compare the predictions of the fiducial models of the GZ-SIMBA and G3 simulations with observational data in Section 4. In Section 5, we use the observational data for the Coma cluster to calibrate model parameters, and discuss the origin of the discrepancy between observations and simulations in connection to the role of AGN feedback.Finally, in Section 6, we summarize our main results and make some further discussions.

Observational Data
The galaxy sample we used is taken from the New York University Value-Added Galaxy Catalog (Blanton et al. 2005) based on SDSS DR7, which is also used for the ELUCID reconstruction (see below).The stellar masses of galaxies are obtained using the relation between the stellar mass-to-light ratio (M/L) and the color based on model magnitudes (see Bell et al. 2003;Yang et al. 2007).The initial mass function (IMF) used to calculate the stellar mass is the Kroupa IMF (Kroupa 2001), which differs from the Chabrier IMF (Chabrier 2003) used in the simulations.However, both Kroupa and Chabrier IMFs have a stellar M/L that is about 0.25 dex lower than the Salpeter IMF (Salpeter 1955).Thus, no correction is needed when comparing simulated and observed galaxies, as discussed in Borch et al. (2006).
The specific star formation rate (sSFR), defined as the ratio between the SFR of a galaxy and its stellar mass, is taken from Brinchmann et al. (2004), derived either directly from emission lines or indirectly from the 4000 A ̊break.Star-forming and quiescent galaxies are separated by a value of logsSFR = − 11, the same as that used for the simulations.The stellar metallicity, Z * , of galaxies is taken from the eBOSS Firefly Value-Added Catalog (Comparat et al. 2017) DR16.12It is measured assuming a Chabrier IMF and the MILES library (Sánchez-Blázquez et al. 2006).In this paper, we take the solar metallicity Z e to be the conventional value of 0.02, following previous studies (e.g., Bruzual & Charlot 2003;Le Borgne et al. 2003;Ma et al. 2016). Tremonti et al. (2004) derived the gas-phase oxygen abundance 12+log(O/H) using optical nebular emission lines for star-forming galaxies in the SDSS.We use their data to compare with the ISM oxygen abundance measured in simulated galaxies.We take the gas-phase oxygen abundance as the metallicity of the ISM (Z ISM ) hereafter.
We corrected the redshift distortion of these galaxies using the method shown in Wang et al. (2016).This method, which employs linear theory to correct for the Kaiser effect (e.g., Wang et al. 2009Wang et al. , 2012) ) and uses the group catalog of Yang et al. (2007) to correct for the finger-of-God effect, has been shown to be capable of correcting for the anisotropy in the twodimensional correlation function (e.g., Shi et al. 2016).Based on the corrected galaxy positions, we select two galaxy samples according to their distances from the Coma cluster.The first sample of galaxies resides in a high-density region (HDR), located within 3R 200c (see the definition of R 200c in Section 3.3) of the Coma cluster center, and the second sample contains galaxies in a low-density region (LDR) between 3R 200c and 7R 200c from the cluster center.There are 135 (371) and 166 (94) star-forming (quiescent) galaxies with the HDR and LDR, respectively.These counts are reduced to 43 and 107 for star-forming galaxies with available data for the oxygen abundance.
In the following, we briefly describe the observational data of the ICM used in this paper.Matsushita (2011) analyzed 28 galaxy clusters, observed with XMM-Newton, including the Coma cluster, our subject of study in this paper.They integrated the spectra in each of the six concentric annular regions centered on the X-ray peak, and fitted the composite spectra with a multitemperature model in the 0.5−10 keV band.They then determined the Fe abundance from the flux ratios of Fe lines to the continuum within an energy range of 3.5 −6 keV.Finally, they obtained the ICM Fe-abundance and temperature profiles for the Coma cluster with projected radius up to 0.5R 180c .We correct their assumed solar abundance value of 0.0133 to our assumed value of 0.02.Simionescu et al. (2013) presented X-ray data from Suzaku for the Coma cluster.The data consist of 24 pointings covering the E, NW, and SW directions contiguously out to a radius of 2°, combined with more pointings toward the NE and W from the cluster center.They modeled the hot gas as a thermal plasma with a single temperature in collisional ionization equilibrium within each shell.Furthermore, they performed a fit to the spectra in each annulus to measure the metallicity.They adopted a solar metallicity of 0.017, which we correct to 0.02.Both temperature and metallicity profiles obtained along azimuths not aligned with the infalling southwestern subcluster (i.e., E+NE+NW+W) are given in their paper.With the assumption of spherical symmetry, they derived deprojected radial profiles of the electron density and specific entropy along the more relaxed NW+W and E+NE directions.In this paper, we only show their averaged profiles along the NW+W and E +NE directions.
Mirakhor & Walker (2020) presented a new extended XMM-Newton mosaic of the Coma cluster that extends beyond its virial radius with almost complete azimuthal coverage.After combining it with the SZ observation from the Planck, they obtained the thermodynamic properties of the ICM in an azimuthally averaged profile as well as in 36 angular sectors.They folded an APEC model through the XMM-Newton response in XSPEC and used the 0.7-1.2keV count rate to derive the APEC normalization used to calculate the electron density.The derived electron density profile is deprojected by assuming spherical symmetry.They then derived the pressure profile P(R) by fitting the projected y profile to a universal pressure formula.The corresponding threedimensional temperature and entropy profiles are calculated using , respectively.We also use observational data for the IGM.Lim et al. (2018b) derived the IGM temperature as a function of local mass density in the volume covered by SDSS DR7.They assumed a broken power-law relation between the electron pressure and the reconstructed mass density field (Wang et al. 2009(Wang et al. , 2016) ) to predict the SZ Compton y parameter map over the SDSS sky.They then performed a Monte Carlo Markov Chain to constrain the free parameters in the broken power law to yield the best match to the observed y map.The temperature-density relation is then derived by assuming a mean cosmic baryon fraction.Although the relation is obtained in the whole SDSS volume rather than the region around the Coma cluster, it is still interesting to compare with our simulation results, as we will show below.

ELUCID Reconstruction of Initial Conditions and Zoom-in Realizations
The ELUCID project aims to reconstruct the ICs of the local Universe.The reconstruction method consists of four steps: (i) identifying galaxy groups; (ii) correcting for redshift distortion; (iii) reconstructing the present-day density field; and (iv) recovering the ICs.A large constrained N-body simulation in the volume of the SDSS DR7 (Abazajian et al. 2009) up to z = 0.12 has been carried out, as shown in Wang et al. (2016).
In what follows, we refer to this CS as the original CS (hereafter OCS).We refer the reader to Wang et al. (2016) and references therein for details of the method and the OCS.In this paper, we focus mainly on one particularly interesting region, centered on the CM at a redshift of 0.0241.We adopt a zoomin technique, which simulates structures of interest with a higher resolution than other regions in the simulation volume (Katz & White 1993).
Details of the method to generate the zoom-in ICs are presented in Li et al. (2022a).Here, we present a brief description.First, we select a high-resolution region (HIR) containing the CM from the z = 0 snapshot of the OCS.The HIR is chosen to be a spherical volume centered on the simulated Coma.We follow all the dark matter particles in the HIR back to the initial time and select an initial HIR (iHIR) that contains all the HIR particles.The iHIR has the same comoving center and geometry as the HIR.To prevent lower-resolution particles from entering the HIR, we set up a buffer region with a size equal to 5% of the iHIR.The cosmic density field outside the iHIR is sampled with particles with three different levels of lower resolution.Gas particles are only placed in the iHIR.We split each high-resolution particle into two, one dark matter particle and one gas particle with a mass ratio equal to (Ω m,0 − Ω b,0 )/Ω b,0 , and separate the two particles in a random direction by a distance equal to half of the mean particle separation.
We adopt the same cosmological parameters (WMAP5; Dunkley et al. 2009) as the OCS: Ω Λ,0 = 0.742, Ω m,0 = 0.258, Ω b,0 = 0.044, h = H 0 /(100 km s −1 Mpc −1 ) = 0.72, σ 8 = 0.80.The mass resolution of the HIR is 8 times as high as the OCS simulation.The corresponding masses of dark matter and gas particles in the iHIR are 3.20 × 10 7 h −1 M e and 6.58 × 10 6 h −1 M e , respectively.The radii of the spherical volume of these simulations are set to be either 20 or 30 h −1 Mpc, much larger than the virial radius of the Coma cluster (see below).This enables us to study galaxies in both high-and low-density environments around the Coma cluster.In addition, we also have CSs for a void (VD) in the SDSS region at z ≈ 0.05 to check the cosmic variance when necessary.The HIR for the VD simulation is a spherical volume with a radius of 40 h −1 Mpc.See Li et al. (2022a) for a detailed definition of the VD.
Information about the simulations is presented in Table 1.The name of each simulation consists of three parts.The first part indicates the code used, either GZ or G3.The second part denotes the galaxy formation model, which is discussed in the following subsections.The third part specifies the simulated region, either the CM or the VD.For example, GZ-SIMBA (SB)-CM is a CS of the Coma cluster run with GZ and SIMBA physics (Davé et al. 2019; see Section 3.2.1),and G3-H-VD is a CS of the VD with G3 and the galaxy formation model presented in Huang et al. (2020; see Section 3.2.2).The first two parts of the name can also be used to indicate the model used.For example, the GZ-SBnA model is the GZ-SB model with the AGN feedback switched off, and the G3-H model is the model presented in Huang et al. (2020).These simulations are run with the support of the Jiutian simulation project.

Notes.
a The name of the large-scale structure that we simulate.The Coma galaxy cluster is located at α J = 194.8,δ J = 27.9, and z = 0.0241, and the mean location of the void is α J = 222.6,δ J = 45.5, and z = 0.0519.b The HIR radius is in units of h −1 Mpc.c This is the fiducial GIZMO-SIMBA model (Davé et al. 2019).
d Both star formation and supernova feedback are enhanced.
e The jet velocity is similar to GZ-SB-CM, but we allow the jet at a lower velocity to heat the surrounding gas to the virial temperature to shift the jet mode-on time to higher redshifts.
f The jet velocity is almost doubled to increase the AGN feedback strength.
g The star formation and stellar feedback model is presented in Huang et al. (2020).

GIZMO-SIMBA
SIMBA (Davé et al. 2019) is a successor to the earlier MUFASA simulations (Davé et al. 2016); both were developed on top of a branched version of the gravity plus hydrodynamics code GZ (Hopkins 2015) in its meshless finite mass mode.SIMBA's default setup, which was used to run the 100 h −1 Mpc boxsize cosmological simulation in Davé et al. (2019), is strictly adopted here and named GZ-SB.Here, we provide brief descriptions of this fiducial model together with our modified models.They are also listed in Table 1 for reference.
This fiducial model includes both photoionization heating and radiative cooling, using the GRACKLE-3.2library (Smith et al. 2017) in its nonequilibrium mode.The neutral hydrogen fraction, with self-shielding following the Rahmati et al. (2013) prescription, is computed self-consistently within GRACKLE, with the metagalactic ionizing flux attenuated according to the gas density.Furthermore, the ionizing background is assumed to be uniform and to follow the Haardt & Madau (2012) Schmidt (1959) law is used to model the star formation: , where the H 2 density is computed from the metallicity and local column density using the subgrid prescription of Krumholz & Gnedin (2011).
Although the newer SIMBA-C (Hough et al. 2023) adopts the advanced Kobayashi chemical model to track 34 different elements, we stick to the original chemical enrichment model in SIMBA, which tracks 11 elements (H, He, C, N, O, Ne, Mg, Si, S, Ca, Fe) from Type II supernovae (SN II), Type Ia supernovae (SN Ia), and asymptotic giant branch (AGB) stars, as described in Oppenheimer & Davé (2006).Yield tables for SN II from Nomoto et al. (2006), for SN Ia from Iwamoto et al. (1999) with 1.4M e of metals per supernova (SN), and for AGB from Oppenheimer & Davé (2006), are adopted in the simulation (assuming a helium fraction of 36%; and a Chabrier 2003 IMF).The dust model, as described in Li et al. (2019), broadly follows that given in McKinnon et al. (2016) with the metal mass passively advected along with the gas.The growth of dust particle mass is through condensation, and the accretion of gasphase metals is via two-body collisions.The dust inside a gas particle can be destroyed by SN shocks, AGN feedback, and collisions with thermally excited gas elements, and its mass and metals are then returned to the gaseous phase.
The stellar feedback is prescribed by two-phase, decoupled, and metal-enriched winds, with 30% of the wind particles ejected hot and the rest ejected at T ≈ 10 3 K.The mass loading factor scales with the stellar mass of the galaxy based on the mass outflow rates from the Feedback In Realistic Environments (FIRE; Hopkins et al. 2014) simulations (Anglés-Alcázar et al. 2017b), where the stellar mass is provided by an on-thefly friends-of-friends galaxy finder.The wind velocity scaling is also based on results from the FIRE simulations (Muratov et al. 2015), with modifications adopted in the SIMBA simulation.Ejected wind particles are temporarily hydrodynamically decoupled to avoid numerical inaccuracies.The cooling and other baryon processes are also paused during the decoupled time.Wind particles are enriched with metals from their surroundings at the launch time according to Nomoto et al. (2006).The star-formation-driven wind mechanisms are described in more detail in Appleby et al. (2021).
SIMBA employs a novel method for SMBH growth and feedback.SMBHs are seeded with M seed ∼ 10 4 M e into galaxies with M *  10 9.5 M e in the SIMBA simulation.However, the iHIR has a slightly higher resolution than SIMBA.Therefore, we seed an SMBH particle when the host galaxy stellar mass is larger than 10 9 M e by scaling the SMBH seeding mass to ∼3 × 10 3 M e .The SMBH mass growth is through a two-mode accretion prescription in the simulation: For cold gas (T < 10 5 K) surrounding the SMBH, torque-limited accretion is implemented following the prescription of Anglés-Alcázar et al. (2017a), which is based onHopkins & Quataert (2011).For hot gas (T > 10 5 K), classical Bondi accretion is adopted.
Kinetic AGN feedback, similar to stellar feedback, is implemented in two modes: "radiative" and "jet" at high and low Eddington ratios separated by f Edd = 0.2.This is designed to mimic the dichotomy in real AGN (e.g., Heckman & Best 2014).The radiative and jet feedback outflows are ejected in a ± direction along the angular momentum vector of the inner disk with zero opening angle and different velocities.The mass outflow rate varies so that the momentum output of the SMBH is 20L/c, where L is the bolometric AGN luminosity.In the radiative mode, winds are ejected at the ISM temperature with an SMBH-mass-dependent velocity calibrated using X-ray detected AGN in SDSS (Perna et al. 2017): 2, the transition to the jet mode begins for galaxies with M BH > 10 7.5 M e , with winds particles ejected at the virial temperature of the halo.Jet outflows receive a velocity boost that increases with decreasing f Edd , up to a maximum of 7000 km s −1 above the radiative mode velocity at f Edd = 0.02: In addition, SIMBA includes an X-ray mode of the AGN feedback, following Choi et al. (2012), to mimic X-ray heating from the accretion disk.This X-ray feedback model is only turned on when the host galaxies have f gas < 0.2, and the AGN feedback is in the full-velocity jet mode.As shown in Cui et al. (2021), jet-mode feedback is the key to quenching galaxies while the X-ray feedback enhances color bimodality of the galaxy population.
In addition to the model GZ-SB that adopts the original SIMBA setup, we have other two models, GZ-SBrw and GZ-SBrs, based on the SIMBA simulation.In Section 5.1, we describe the details of the two models.In general, both star formation and SN feedback are enhanced in the two models, while the former tries to shift the onset of the jet mode to higher redshift, and the latter tries to simply increase the jet velocity.Moreover, we ran a simulation with the fiducial setup but without the SMBH part.This simulation, referred to as GZ-SBnA, is used to compare with the G3 runs that also neglect the SMBH growth and AGN feedback.

GADGET-3
We also include simulations evolved using the G3 code and the galaxy formation model presented in Huang et al. (2019Huang et al. ( , 2020)).We use names starting with "G3-H" to denote these simulations, as shown in Table 1.The G3-H-CM and G3-H-VD simulations are exactly those presented in Li et al. (2022a).Unlike GZ-SB and its other versions, the G3-H model adopts a smoothed particle hydrodynamics (SPH) technique to solve fluid equations and a galaxy formation model without AGN feedback.Therefore, comparing G3-H, GZ-SB with observations can shed light on the importance of AGN feedback.We refer readers to Huang et al. (2019Huang et al. ( , 2020) ) and Li et al. (2022a) for details of the code and the simulations.Here, we only give a brief description.
The G3 code (Huang et al. 2019(Huang et al. , 2020)), an updated version of GADGET-2 (Springel 2005), includes several numerical improvements in the SPH technique, such as using the pressure-entropy formulation (Hopkins 2013) to integrate fluid equations, using a quintic spline kernel over 128 neighboring particles to measure fluid quantities, and adopting the Cullen & Dehnen (2010) viscosity algorithm and artificial conduction to capture shocks.These updates lead to considerable improvements in the instabilities at fluid interfaces in subsonic flows (Sembolini et al. 2016a(Sembolini et al. , 2016b;;Huang et al. 2019).
The radiative cooling of the G3-H model includes cooling from hydrogen, helium, and metal lines for altogether 11 elements (Wiersma et al. 2009).The cooling rate is recalculated according to a uniform ionizing background given by Haardt & Madau (2012).Following Springel & Hernquist (2003), we use a subgrid approach to model the multiphase ISM in dense regions with n H > 0.13 cm −3 and a star formation recipe that matches the Kennicutt-Schmidt law.The G3-H model traces the enrichment of four metal species, C, O, Si, and Fe, produced from SN II, SN Ia, and AGB stars (Oppenheimer & Davé 2008).
We adopt a kinetic subgrid model for the stellar feedback (see also Huang et al. 2019Huang et al. , 2020).An SPH particle in starforming regions can be ejected as a wind particle with a probability proportional to the local SFR.We adopt the same set of wind parameters as the fiducial simulation from Huang et al. (2020), which matches a broad range of observations (e.g., Davé et al. 2013).This results in momentum-driven wind scalings for large σ gal and SN-energy-driven wind scalings for small σ gal , where σ gal is the velocity dispersion of the galaxy.These scalings are very similar to those found in very highresolution zoom simulations (Hopkins et al. 2012(Hopkins et al. , 2014;;Muratov et al. 2015).We model the launch of galactic winds from star-forming galaxies with the mass loading factor, η ≡ ejection rate/SFR.So defined, we assume that, gal h s ~bfor small σ gal , and gal h s ~afor large σ gal .The initial wind speed, v w , scales linearly with σ gal .The wind scalings are assumed at the wind launch from the star-forming regions of the galaxy while the very high-resolution zoom simulations report their wind scalings at 0.25r vir (Muratov et al. 2015).Hence, we have had to slightly increase our wind launch velocities to reproduce their behavior at 0.25r vir .We also cap the wind speed so that the energy in the winds does not exceed that available in SNe.

Galaxies and Halos around the Coma Cluster
Galaxies are identified using the Spline Kernel Interpolative Denmax algorithm (Kereš et al. 2005).This technique groups particles on the density gradient lines that converge on the same local density maxima.To do this, a smoothed density field is first constructed using all particles.Then, test particles at the position of star-forming gas and stars are moved along the positive gradient direction of the density field until a local density maxima is reached.All particles with test particles at the same local density maxima are considered as potential members of one galaxy.An unbinding procedure is then applied to the group of particles to determine the final member particles of the galaxy.The center of the galaxy is defined as the location of its density maxima.The star and/or star-forming gas particles in each galaxy can be used to calculate the stellar mass M * , stellar metallicity Z * , sSFR, and ISM oxygen abundance 12 + log(O/H).We divided the simulated galaxies into star-forming and quiescent galaxies using logsSFR = −11.
We use an FoF algorithm to identify halos.High-resolution dark matter particles are linked to each other within a link length b equal to 0.2 times the mean separation of dark matter particles.Gas, star, and black hole particles are also linked to their closest dark matter particle if their distances from the particle are less than b.The halo center is defined as the position of the dark matter particle that possesses the minimum potential among all the FoF particles.The halo virial radius, R 200c , is the radius of a sphere within which the average density is 200 times the critical density of the universe at that redshift, and the halo mass, M 200c , refers to the total mass within R 200c .
The Coma clusters in these simulations have almost the same M 200c and assembly history.At z = 0, the Coma cluster has M 200c ≈ 7.52 × 10 14 h −1 M e and R 200c ≈ 1.48 h −1 Mpc, corresponding to 70. 9 ¢ viewed from the Earth.These are consistent with those of the real Coma cluster, M 200c ≈ 6.2 × 10 14 h −1 M e (Okabe et al. 2014) and R 70 200c » ¢ (Simionescu et al. 2013).As an example, Figure 1 shows the evolution of dark matter density, gas density, gas temperature, gas metallicity, and stellar density around the Coma cluster in the GZ-SBrw-CM simulation.During the last 2.4 Gigayears, the Coma cluster experienced several massive merger events.These mergers generate multiple shocks that heat the ICM and spread far beyond the virial radius, as can be seen in the temperature map.Interestingly, our simulated Coma cluster contains two BCGs residing in the center of the cluster, similar to the observed Coma cluster.
Figure 2 shows the spatial distributions (in J2000.0 coordinates) of the simulated and observed galaxies, as seen from the Earth.Galaxies with * M h M log 9.5 2 ( )  > are complete at this redshift (see below), and thus, only these galaxies are shown.As one can see, the simulations reproduce the observed large-scale structures, such as the Coma cluster in the center, two smaller clusters at the SW and NE of the Coma cluster, and the filaments connecting the clusters.However, smaller structures are not accurately reproduced, as expected.One of the most significant discrepancies between the simulations and observation is the quiescent galaxy population in the LDR.It appears that simulations produce too few quiescent galaxies, as we will discuss in the following section.

Fiducial Simulations versus Observations
In this section, we compare the predictions of the two fiducial simulations, GZ-SB-CM and G3-H-CM, with observational data.We focus on galaxies in Section 4.1, and on the ICM and IGM in Section 4.2.).In general, galaxy properties in the two regions differ significantly.Galaxies in HDR tend to be more massive, more frequently quenched, and metal richer than those in LDR, suggesting that environmental processes have a significant impact on galaxy properties (Wang Let us first look at the SMFs.Since all galaxies are located within a small volume, the SMFs can be obtained by directly counting the number of galaxies without any correction.As shown in Figure 3, the observed star-forming galaxies are complete to 10 9.0 h −2 M e while quiescent galaxies are complete to 10 9.5 h −2 M e .So the total population is complete to 10 9.5 h −2 M e .Note that the stellar mass of simulated galaxies with 200 star particles roughly corresponds to

Properties of Galaxies
-, lower than the observational limit.In HDR, the predicted SMFs of the two simulations for the total galaxy population agree well with the observational data at , while the GZ-SB-CM SMF is lower than the observation.
Inspecting the quiescent and star-forming galaxies separately, we can see a larger difference between observational and simulation results.The G3-H model produces more starforming galaxies and fewer quiescent galaxies in both HDR and LDR than in the observation.This could be because the G3-H model does not include AGN feedback, but includes strong stellar and SN feedback to prevent the growth of massive galaxies.Interestingly, the discrepancy, in the massive galaxies (10 10 h −2 M e ), for star-forming galaxies is more prominent in HDR than in LDR, while the discrepancy for quiescent galaxies is the opposite.In particular, there are almost no quiescent galaxies in LDR, which is very different from the observation.This indicates that most of the quiescent galaxies in the G3-H-CM should be produced primarily by environmental effects, such as ram pressure and tidal stripping, while environmental effects are too weak to quench these massive galaxies in LDRs.This large discrepancy seen in LDR suggests that additional quenching mechanisms, such as AGN feedback, are probably necessary to quench massive galaxies.Moreover, the "over" quenching of low-mass galaxies (10 9 h −2 M e ) may owe to the observational and simulation limits.
The GZ-SB model produces results that are much closer to the observational data than the G3-H model.In HDR, the predicted SMFs for both star-forming and quiescent galaxies match the observational results well.In LDR, the GZ-SB model works well for star-forming galaxies, but significantly underestimates the number of quiescent galaxies, with a difference of about 0.6 dex.This discrepancy is only visible when galaxies are separated according to their environments, as the number of quiescent galaxies in LDR is much lower than that found in HDR.This suggests that modifications of the subgrid physics are needed to reproduce the observed properties of the Coma cluster and that CSs can provide important constraints on galaxy formation models, as we will demonstrate below using more observations.
Next, we examine the stellar mass-metallicity relation, shown in Figure 4. Metallicity measurements of observed galaxies have quite a large uncertainty, as shown in the top row.Hence, to compare the observation with our simulation results, we assign the same observed metallicity uncertainties to our simulated galaxies.Specifically, for each simulated galaxy of a given mass, we generate 100,000 mock galaxies with their metallicity following a Gaussian distribution.This distribution has the same mean metallicity as the simulated galaxy and a dispersion equal to the mean uncertainty of observed galaxies with the same stellar mass as the simulated galaxy.The solid line shows the median metallicity of our mock galaxies as a function of stellar mass, and the dashed lines indicate the 16th-84th percentiles of the distribution.Results for the observed galaxies are repeated in panels of the simulated results for comparison.
The median M * -Z * relation for the G3-H model in LDR is in good agreement with the observational data, and the predicted scatter is only slightly smaller than observed.However, the predicted M * -Z * relation in HDR is similar to that in LDR and is systematically lower than the observational data.This appears in conflict with the SMF result, which shows a strong environmental dependence (Figure 3).In addition, quiescent galaxies and star-forming galaxies in HDR follow the same trend, which is not consistent with observations (Peng et al. 2015).One potential explanation is that the strong stellar wind feedback in the G3-H model makes both the ISM and CGM of simulated galaxies more extended and hence more susceptible to ram pressure stripping, leading to a rapid removal of the ISM and little evolution in the stellar metallicity.This effect is more pronounced for low-mass galaxies and could owe to their shallower gravitational potential wells, which would give a higher quenching rate for low-mass galaxies than for massive galaxies in HDR (see Figure 3).
The GZ-SB model predicts M * -Z * relations that are significantly different from the observational data, particularly at the massive end.This discrepancy in the same mass range is also visible in the original SIMBA paper (Figure 10 in Davé et al. 2019).In LDR, the intrinsic scatter of the simulated relation is very small, similar to that of the G3-H model.In HDR, however, the intrinsic scatter is much larger, which is different from the G3-H model.This could owe to the fact that, in a high-density environment, the CGM of a galaxy is more easily stripped while its ISM is not, allowing the stellar metallicity to increase gradually over time owing to metals produced by SNe.This suggests that the stripping effect in the low-mass GZ-SB galaxies is weaker than that in low-mass G3-H galaxies.
Figure 5 displays the oxygen abundance of the ISM as a function of stellar mass for both HDR and LDR.We only present results for star-forming galaxies, as the ISM metallicity of real galaxies is measured using emission lines.The observed slope is steeper for low-mass galaxies than for massive galaxies, which is consistent with results found by Tremonti et al. (2004).In addition, the ISM metallicity for low-mass The vertical gray dashed lines represent the estimated complete mass limit of the observations.galaxies is slightly higher in HDR than in LDR.The ISM metallicity of simulated galaxies is estimated using starforming gas particles with a nonzero SFR, which might not be exactly the same as the observational definition.The solid lines indicate median relations, and the dashed lines represent the scatter (16th-84th percentiles) of the relation.Note that the scatter of the simulated M * -Z ISM relation also involves measurement uncertainties of the ISM metallicity, similar to the M * -Z * relation described above.
The M * -Z ISM relation for G3-H galaxies in LDR is in agreement with the observational data in terms of both the overall trend and the scatter.In particular, the G3-H model predicts a flatter relation at the massive end, consistent with the observations.In HDR, the predicted metallicity is slightly lower than that observed for star-forming galaxies with masses * M h M 9 log 10 2 ( ) , likely caused by the strong stripping effect discussed above.The simulated relation extends to more massive galaxies than the observed one, as the simulation predicts many massive star-forming galaxies at Similarly to the M * -Z * relation, the predicted M * -Z ISM relation does not show any significant dependence on the environment.Unlike the G3-H model, the GZ-SB predicts an ISM with metallicity that is significantly lower than those observed in both LDR and HDR, with the deficit being particularly large for low-mass galaxies (about 0.5 dex).This deficit likely has the same origin as the deficit in the stellar metallicity.Additionally, the predicted M * -Z ISM  relation is steeper at the high-mass end than at the low-mass end, contrary to observational results.

Intracluster and Intergalactic Media
The ICM and IGM can fuel star formation in galaxies, and the feedback from galaxies can also leave imprints on these media.Therefore, it is essential to compare simulation predictions and observations of the gas they represent.Figure 6 shows the observed electron density (n e ), electron temperature, metallicity, and entropy profiles of the ICM of the Coma cluster, taken from Matsushita (2011), Simionescu et al. (2013), and Mirakhor & Walker (2020).The electron density decreases rapidly with increasing distance from the cluster center, while the temperature decreases slowly, with a typical temperature of about 5 keV within the virial radius.The ICM entropy increases rapidly toward the cluster outskirts, and the ICM metallicity is roughly constant with radius, around 0.3Z e .Note that different observations give quite similar results, suggesting that the observational data are reliable.
The simulated profiles are depicted as lines in Figure 6.The values are the average in three-dimensional spherical shells.As discussed in Section 2, the situation in observational results is more complex.The electron density and entropy profiles in Figure 6 are deprojected, allowing direct comparisons with the simulation results.The temperature profile from Mirakhor & Walker (2020) is also deprojected, while others are not.Despite this, the profiles are similar, indicating that the projection effect is not significant.Note that the temperature profile of Simionescu et al. (2013) is measured in a narrow sky region, which may explain the large fluctuation seen in the data.The observed metallicity profiles are not deprojected.Since the metallicity gradient is small, the projection effect may be ignored.As one can see, the electron density profile predicted by the G3-H-CM simulation roughly matches the observational data, although it is slightly lower (higher) than the observation at R < ( > )0.3R 200c .The GZ-SB model produces a density profile that is similar to the observational data at R > 0.3R 200c , but has a much lower and flatter profile at R < 0.3R 200c .Recently, Robson & Davé (2020) found that massive clusters in the SIMBA flagship simulation show flat central gas cores similar to what we find here.This implies that flat cores are very common in SIMBA clusters, which could be due to its strong AGN feedback (see a similar result in Li et al. 2023b, with the 300 cluster sample with their hypothesis).The situation in the temperature profile is different.The GZ-SB-CM simulation produces a temperature profile that is in agreement with the observational data from 9 keV at 0.1R 200c to 3 keV at the virial radius.The only difference appears in the innermost region (<0.1R200c ), where the simulated temperature reaches 12 keV, higher than the observed value.The G3-H-CM simulation, on the other hand, predicts a much lower temperature within that half of the virial radius, with the difference becoming smaller near the halo boundary.This is likely because the cooling of the ICM is not effectively suppressed in G3-H-CM.
Employing the same observational method of Simionescu et al. (2013), Mirakhor & Walker (2020), we calculate the ICM entropy using the definition , where P e (R) and n e (R) are the electron pressure and number density profiles, respectively.We assume an adiabatic gas and adopt γ = 5/3, so that P e (R) = (γ − 1)u e (R), where u e (R) is the mean internal energy density of electrons at R. Both simulations generate steep entropy profiles in the cluster outskirts (R > 0.4R 200c ), in agreement with both the observations and previous results (Pratt et al. 2010).In the inner region (0.1R 200c < R < 0.4R 200c ), the GZ-SB-CM shows a much higher and more flat entropy profile than the G3-H-CM.The observational data lies between the two simulated profiles and appears to be closer to G3-H-CM.The discrepancy between the observation and GZ-SB-CM increases with decreasing radius.In the innermost region (R < 0.1R 200c ), the entropy profile predicted by G3-H-CM drops quickly, while that predicted by GZ-SB-CM remains constant.This is likely owing to the presence or absence of AGN feedback in the two simulations.
As shown in Figures 4 and 5, the stellar and interstellar metallicities predicted by the GZ-SB model are lower than the observed values.Where are the missing metals?Are they expelled into the ICM/IGM by feedback processes, or do the galaxy formation models produce too few metals?The ICM metallicity is crucial for distinguishing between these two possibilities.As shown in Figure 6, the GZ-SB-CM produces a nearly flat metallicity profile within the virial radius, and the average metallicity is around 0.16Z e , which is about half of the observed value.The discrepancy between the observation and the G3-H-CM simulation is even greater, with the predicted mean metallicity only about 0.13Z e .These results demonstrate that the star formation processes used in these models produce too few metals.In the next section, we will improve the model predictions with new simulations.Li et al. (2022a) found a bow-like shock in the G3-H-CM Coma cluster, which was also observed in the same location by eROSITA (Churazov et al. 2021).Figure 7 displays the shock feature in the 0.4−2 keV band observed by eROSITA, taken from Figure 6 in Churazov et al. (2021).To compare the simulations with the observation, we show X-ray brightness maps for the simulated Coma in the same band.The method for calculating the X-ray maps is described in Li et al. (2022a).To emphasize the shock feature, Churazov et al. (2021) applied a flat-fielding procedure by normalizing the original image with a mean model profile.To mimic the observation, a background X-ray emission was added, and a similar procedure was adopted in the simulated maps.It is evident that the bow-like shock front in G3-H-CM is very prominent, while the feature is weaker in GZ-SB-CM.The formation of a shock front is sensitive to gas properties before the collision of gas flows (Zinger et al. 2018), thus providing an opportunity to explore the ICM properties.In addition, the substructure seen in the SW of the Coma cluster in the eROSITA observations is also present in the X-ray map of the G3-H-CM, but almost absent in that of the GZ-SB-CM.
The lower right panel of Figure 6 shows that the entropy at R > 4R 200c is significantly different between the two fiducial simulations.This implies that the IGM property may also be used to distinguish different models.Combining the reconstructed density field (Wang et al. 2009(Wang et al. , 2016) ) and the Planck SZ map, Lim et al. (2018b) obtained the IGM temperature (T IGM ) as a function of mass density (ρ m ), where the temperature and density are both estimated on grids of a size of 1 h −1 Mpc.The result is shown in Figure 8, with the shaded region representing the uncertainty of the relation.The observational T IGM increases with increasing ρ m from )~at 100 m r .We note that the result of Lim et al. (2018b) was derived using data from the entire SDSS survey region, rather than just the Coma region, and so, the comparison with simulation results is with this caveat.
To compare with the observation, we divide the HIR into grids of side length 1 h −1 Mpc.For each grid, we calculate

¯=
, where P e ¯and n e ¯are the average electron pressure and electron density of the grid, respectively, and k B is the Boltzmann constant.Here, n e ¯is calculated directly from ρ m and the cosmic baryon fraction, assuming zero metallicity and full ionization (see Lim et al. 2018b, Section 4.2).This mimics the observation as the SZ y parameter is the integral of electron pressure.In this analysis, we exclude contributions of starforming gas particles and wind particles.As before, the simulation results in Figure 8  , the GZ-SB-CM predicts a T IGM that is about 3 orders of magnitude higher than the observations.In contrast, the predicted T IGM increases monotonically with ρ m in G3-H-CM, with both the amplitude and slope of the relation in good agreement with the observational data.

Using the Coma Cluster as a Test Bed to Diagnose
Feedback Processes In the preceding section, we see a complex situation in model predictions.Neither simulation can accurately reproduce all the observational results, and each of the simulations has its advantages and disadvantages.These simulations are tailored to fit observational data in a large volume, but almost no data on ICM and IGM were used to constrain model parameters.Our analyses above revealed two main discrepancies: the first is the underprediction of metallicities in stars, the ISM, and the ICM (the metal budget problem), and the second is the discrepancies in the physical properties of the ICM and IGM.In this section, we attempt to understand the underlying mechanisms that cause the discrepancies and to make improvements in the related modeling.

Metal Budget Problem and Model Calibration
We first examine the metal budget problem for the GZ-SB model.The metallicities of stars, ISM, and ICM can all be affected by a variety of processes.To determine the source of the metal budget problem, it is essential to distinguish the effects produced by different processes.As seen in Figures 4  and 5, this problem is present in both HDR and LDR, and the discrepancy between the observations and the model is similar in both regions.This suggests that internal processes are responsible for the problem, since environmental effects are expected to be much weaker in LDR than in HDR.Furthermore, the problem is significant in both low-mass and high-mass galaxies, indicating that the relevant process should not depend strongly on stellar mass.This eliminates AGN activity and feedback as the cause, as AGN feedback is expected to be important only for massive galaxies.The predicted amounts of metals in all baryonic components, i.e., stars, ISM, and ICM, are all lower than the observed values, suggesting that we should increase the total metal yield rather than redistributing it among the different baryonic components.
Thus, one potential solution to the metal budget problem is to increase the SFR.In the two modified models, GZ-SBrw and GZ-SBrs, we adopt a slightly higher star formation efficiency per freefall time, SFR/(M gas /t freefall ) = 0.03, instead of 0.02 used in the fiducial SIMBA model.Additionally, a slightly higher metallicity floor is adopted to raise the H 2 fraction, to increase the SFR at high redshift to meet the recent JWST findings of very high redshift galaxies.To counterbalance the mass growth caused by the increased SFR, the SN strength is also enhanced through the coupling factor and the wind efficiency, which should also increase the metals in the ISM and ICM.The underprediction of quiescent galaxies in LDR suggests that the AGN feedback should also be increased, which is consistent with the results from The 300 Clusters Project (Cui et al. 2018(Cui et al. , 2022)).To this end, we adopt a much higher jet velocity of 15,000 km s −1 in the strong jet model-GZ-SBrs, instead of 7000 km s −1 in the fiducial run.In the weak jet model-GZ-SBrw, we employ a different strategy that changes the AGN feedback efficiency by shifting the jet onset time to an earlier redshift, which we accomplish by heating the wind particles to the virial temperature with a jet velocity similar to the fiducial model.GZ-SBrw and GZ-SBrs have the same star formation and SN feedback model but different AGN feedback models (Table 1).In addition, we include the no AGN  (2021).The purple dashed curve marks the observational shock front with a small shift.The SW structure corresponds to the NGC 4839 subcluster.Note that there is a small offset between the simulated and real subclusters.The simulated maps adopt the same color scale as the observational map (private communication for the observational scale).The gray dashed circles represent the scale of R 200c in each simulation.
through quasar absorption lines.Figure 1 shows the evolution of the ICM metallicity of the GZ-SBrw-CM simulation.Little change is seen in the mean metallicity since z = 0.2 when the Coma cluster was being formed.This might explain the flat slope of the ICM metallicity, as the metals had already been released to the ICM before the Coma cluster was formed.
We also examine the impact of the modification on the SMF.The SMF discrepancies between the GZ-SB-CM and the observations mainly appear in LDR.In the new models, the number density of quiescent galaxies is significantly enhanced in LDR.Therefore, the deficit of quiescent galaxies becomes much smaller.In contrast, the number density in HDR is almost unchanged.Comparing the two new models, we see that enhancing the jet velocity indeed increases the quiescent population, but the effect is small.We compare the mass of star-forming gas (M SF ) in star-forming galaxies before and after the modification.In LDR, star-forming galaxies in the GZ-SBrw/rs model have lower M SF /M * than those in the GZ-SB model.The higher SFR and SN feedback may consume and expel more cold gas than in the fiducial model.This may make these galaxies more vulnerable to AGN feedback, so as to increase the number of quenched galaxies with In HDR, environmental effects may play the dominant role in quenching galaxies, so the increase of quiescent galaxies by the change in SFR and SN feedback is insignificant.
We also examine the physical properties of the ICM and IGM predicted by the new models.As seen in Figure 6, the electron density and entropy profiles of the two new models are similar to those of the original GZ-SB model.The GZ-SBrw model appears to be a better match to the observed temperature profile in the inner region (R < 0.1R 200c ) than both the GZ-SB and GZ-SBrs models, likely because of its weaker AGN feedback.The bow-like shock in the two simulations is still very weak, as shown in Figure 7.In addition, the IGM temperature at low-density environments is much higher than the observed one, similar to the GZ-SB model (Figure 8).In general, the predicted ICM and IGM properties, except for gas metallicity, are comparable with the original simulations.Thus, the modifications made do not seem to resolve the discrepancies found in the ICM and IGM (see the next subsection for further discussions).
Finally, we consider the metal budget problem for the G3-H model.The predicted M * -Z * relation and M * -Z ISM relation are in good agreement with the observations in LDR, but the discrepancy in HDR is significant smaller than that for the GZ-SB model.This is distinct from the GZ-SB model problems, and suggests that the modification strategy for the GZ-SB model would not work for the G3-H model.We believe that the issue is related to environmental effects.The similarity in the M * -Z * relation between star-forming and quiescent galaxies suggests that environmental processes can quickly strip the ISM from satellites and prevent further star formation and metal production.A similar argument may also apply to the metal deficit in the ICM.One potential cause of the strong environmental effect is that strong stellar wind feedback makes the ISM vulnerable to tidal/ram pressure stripping (Bahé & McCarthy 2015).However, decreasing the feedback strength is not a solution as it will create too many massive galaxies.An additional internal quenching mechanism seems to be necessary.

Sources of the Discrepancy between Simulations and Observations in ICM and IGM Properties
There are a number of reasons why the simulated ICM/IGM properties are different from those observed, as suggested by Figures 6, 7, and 8.The first is that the reconstruction may not be precise enough to recover substructures that depend sensitively on nonlinear processes.The second is that hydrodynamic solvers may not be able to handle complex and violent processes accurately.The third is that the subgrid physics, such as star formation, SN wind feedback, and black hole accretion and AGN feedback, may not be appropriately modeled in the simulations, owing to either a lack of resolution, failures in the hydrodynamic scheme, or inadequacies of the models, or any combination of the above.Moreover, some mechanisms that are not included here, such as cosmic rays and magnetic fields, may also affect the ICM/IGM.To better understand the effect of AGN feedback, we involve the GZ-SBnA-CM simulation (Table 1), which turns off black hole accretion and AGN feedback.As one will see, CSs based on ELUCID provide an unprecedented opportunity to assess all these effects by comparing predictions of different models with the observational data.
We first focus on the ICM within the Coma cluster.Figure 9 illustrates the development of the density, temperature, and entropy profiles at z < 0.1 for the GZ-SB, GZ-SBrw, GZ-SBnA, and G3-H models.This shows how these profiles today are established in the recent past.A similar evolution history can be observed in all four simulations.At z ∼ 0.1, the temperature profile is nearly flat.The temperature in the inner region then increases rapidly, producing a declining temperature profile at z = 0.There are only two major heating mechanisms for the ICM in the simulations: one is the shock heating caused by accretion and mergers, and the other is heating by AGN feedback.Since the two simulations without AGN feedback give comparable results, only shock heating seems to be able to drive the evolution in the temperature profile.As can be seen in both the dark matter and star maps shown in Figure 1, there is indeed a massive substructure that is moving toward the Coma center at z ∼ 0.1.
It is interesting to examine how the ICM reacts to the merger event in the different simulations.Comparing GZ-SBnA-CM with GZ-SB-CM and GZ-SBrw-CM, it is evident that GZ-SBnA-CM is similar to the other two in the outer region of the cluster (0.4R 200c < R < 1.5R 200c ) but significantly different in the inner region.Therefore, in the following analysis, we will focus on the inner region.At z ∼ 0.1, GZ-SBnA-CM has a much higher density, a slightly lower temperature, and a much lower entropy than the other two simulations.These differences likely owe to the presence or lack of AGN feedback.As the merger occurs, the density initially increases with time and then decreases, the temperature initially rises quickly with time and then remains more or less constant for the simulations with AGN feedback but quickly drops for the GZ-SBnA-CM, and the entropy increases quickly and forms a flat core in the central region in all three simulations.However, the density and entropy profiles in the inner region in GZ-SB-CM and GZ-SBrw-CM are always flat, very different from those in GZ-SBnA-CM, in which the profile slopes vary significantly with time.This suggests that the AGN heating in SIMBA stabilizes the inner density at a relatively low level, and the entropy at a relatively high level.This is also seen in other state-of-art simulations (Altamura et al. 2023;Lehle et al. 2023).
Next, let us compare G3-H-CM with GZ-SBnA-CM.At z ∼ 0.1, the G3-H-CM simulation has density, temperature, and entropy profiles that are similar to GZ-SBnA-CM: the density/ entropy profiles are steep while the temperature profile is flat.After the merger, however, differences between the two become visible.The density and entropy of G3-H-CM remain more-or-less unchanged and stay close to the observed values.The temperature of G3-H-CM first increases quickly to reach the observed profile, and then quickly drops back to the value at z ∼ 0.1.In contrast, the temperature in GZ-SBnA-CM drops much more slowly.The different responses to the merger in G3-H-CM and GZ-SBnA-CM might be caused by the fact that G3 and GZ codes handle hydrodynamics differently and that the cooling in G3 is more efficient (see also Huang et al. 2019).
The primary source of uncertainties in the reconstruction on halo scales is that highly nonlinear events, such as mergers, may not be well reproduced.In principle, one might investigate the effects of such uncertainties by simulating different realizations sampling the posterior distribution of reconstruction (e.g., Ata et al. 2022).However, running many hydrodynamic simulations, especially at the resolution in this study, to achieve this is very costly.The basic property of ELUCID (Wang et al. 2014(Wang et al. , 2016) ) is that CSs based on the reconstructed ICs can reproduce the large-scale structure at z ∼ 0. Our N-body simulations showed that the large-scale structure around, and the mass of the Coma cluster have evolved little since z ∼ 0.1, suggesting that each snapshot at z < 0.1 may be taken to represent those in the present-day universe.The changes produced by merger events at z ∼ 0.1 might then be considered as a rough measure of the fluctuations generated by reconstruction uncertainties.
In general, the density and temperature profiles predicted by the four simulations are in line with the observational data within the "uncertainties" given by the variations among the snapshots at z ∼ 0.1.In the inner region, the median density profiles of GZ-SB-CM and GZ-SBrw-CM are lower than those observed, whereas the median profiles of the two non-AGN simulations match the observational data better.The predicted median temperature of the G3-H-CM is slightly lower than the observed value, while the predictions of the other three simulations match the observation.The most significant difference is seen in the entropy profiles.The simulations that incorporate AGN feedback fail to reproduce the observed entropy profiles in the inner region, while the two simulations without AGN are in good agreement with the observational data.Despite the differences between GZ and G3, both G3-H-CM and GZ-SBnA-CM demonstrate that the CSs of the Coma cluster without including the AGN feedback can reproduce the observed ICM properties better than simulations that include the AGN feedback.Interestingly, the density, temperature, and entropy profiles at a slightly higher redshift, e.g., z ∼ 0.05 in GZ-SBnA-CM and G3-H-CM, match the observational results better than at z = 0. Thus, if the reconstruction were fine-tuned so that the predicted merger occurred slightly later, the simulation results at z = 0 would match the observational data better.
Figure 7 shows the X-ray map around the Coma cluster, and as discussed earlier, the map contains two prominent structures: a subcluster in the SW with a small offset relative to the observed one and a shock front in the central region., consistent with the halo mass of NGC 4839.As shown in Figure 7, the X-ray emission from the subcluster is quite prominent in the observational data, and in the two simulations without AGN feedback, but is quite weak in the three simulations with AGNs.The X-ray luminosities of the NGC 4938 subcluster in G3-H-CM and GZ-SBnA-CM are comparable and about 10-20 times higher than those in the other three simulations with AGN feedback.Sasaki et al. (2016) estimated the gas mass and total mass within the central 98 kpc of the NGC 4839 subcluster using X-ray data and weak lensing, respectively.They obtained a gas fraction of 6.7%.The predicted gas fraction within 98 kpc is about 7.5% for the G3-H-CM simulation and 3.7% for the GZ-SBnA-CM simulation, but only about 0.5% for the three simulations with AGN feedback.Apparently, the AGN feedback implemented in our GZ-SB series is too effective in removing the CGM in the central region of the subcluster.Our analysis suggests that the CGM in group-sized halos can be significantly affected by AGN feedback, and that the X-ray emission of the CGM is a sensitive probe of this feedback.
In the simulations, the inner shock feature seen in the X-ray maps is triggered by a collision between the main body of the Coma cluster and a subhalo.We found that the subhalo gas was largely expelled by the AGN feedback before the collision in GZ-SB-CM and GZ-SBrw/rs-CM, so that the collision cannot generate a strong shock.The shock in GZ-SBnA-CM appears weaker than that in G3-H-CM, probably because a larger portion of the gas in the subhalo is converted into stars in GZ-SBnA-CM (see SMF in Figure 3) owing to a lack of AGN feedback in GZ-SBnA-CM and weaker SN wind feedback than in the G3-H-CM simulation.The bow-like structure observed in the eROSITA X-ray data does not correspond to features generated by a directional AGN jet.Typically, AGN jets create X-ray cavities by driving the ICM away (Fabian 2012), while the observed structure is in fact brighter than the X-ray background produced by the ICM.The morphology of the structure closely resembles the shock front generated by merger events in our CSs (e.g., Li et al. 2023a).This suggests that the amount of subhalo gas in the real Universe should not have been reduced severely by feedback before the collision, which generates the shock; because otherwise, the collision would not be able to generate large amounts of X-ray.Unfortunately, a quantitative comparison with the observational data is not feasible due to the unavailability of the data.In reality, AGN activity is expected to occur sporadically and randomly, a behavior quite different from that expected from the roughly stable structure seen in the simulations.A more stochastic AGN feedback model might result in a lower coupling with nearby gas elements, thereby exerting a lesser influence on the ICM.
Finally, we attempt to understand the source of the discrepancy in IGM properties (Figure 8).We can see that GZ-SBrw also predicts a much higher T IGM than the observations at low density, similar to GZ-SB.We do not present results of the GZ-SBrs simulation, as the results for the two modified models are similar.This is expected, as the AGN feedback models in the GZ-SB model and its revisions are similar.We also display the result for GZ-SBnA-CM in Figure 8.As one can see, this result is very similar to that for G3-H-CM and the observational data.This suggests that the two codes, G3 and GZ, can generate consistent IGM properties, as neither of them has AGN feedback.In addition, GZ-SBnA-CM and G3-H-CM implement very different stellar and SN feedback models, suggesting that stellar and SN feedback has little effect on the physical properties of the large-scale IGM.
The observational temperature-density relation of the IGM is derived from the entire SDSS region (Lim et al. 2018a), while the simulated relations shown here are obtained in a small volume around the Coma cluster.It is thus possible that the relation depends on large-scale structures, and the similarity between the non-AGN models and observation is just a coincidence.As a check, we also present results for the G3-H-VD simulation (Table 1).This is a simulation for a large VD in the SDSS region (see Li et al. 2022a, for more details).Interestingly, there is no significant difference in the median relationship between G3-H-VD, G3-H-CM, and GZ-SBnA-CM, although the fraction of high-temperature IGM is slightly larger in G3-H-CM and GZ-SBnA-CM than in G3-H-VD.This suggests that, when AGN feedback is switched off, the T IGM -ρ m relation of the IGM is not strongly affected by large-scale structures, indicating that the simulated relations presented above may be valid for the global relation.It is also worth noting that a more careful study, which mimics the Planck beam, nulling filter, CMB, and infrared background subtraction, is required to properly quantify the uncertainty in the observational data.
Our findings demonstrate that, when AGNs are deactivated, the simulated T IGM -ρ m relation is not affected significantly by the large-scale structure, the stellar/SN feedback model, and the method used to solve the fluid equations.More importantly, the simulated T IGM -ρ m relations without AGN feedback are in good agreement with the observations, suggesting that, over a broad range of density, the AGN feedback may not have a significant impact on the IGM temperature.The large discrepancy between simulations containing AGN feedback and the observational data suggests that the implemented AGN feedback in these simulations is too strong in heating the IGM.

Summary and Discussion
We conducted a series of zoom-in, CSs to investigate the formation and evolution of the CM and its surroundings.We employed two different galaxy formation models, GZ-SIMBA (GZ-SB) and G3-H, which differ in that GZ-SB includes the growth of SMBHs and AGN feedback while G3-H does not.Our CSs, which accurately reproduce the large-scale structures of the local Universe, enable us to make a more detailed, "oneto-one" comparison between models and observations.To this end, we divided the simulated HIR into two regions according to their distance to the Coma cluster center and compared observations and models in the two regions separately.We also included observations of large-scale gas, such as the ICM, CGM, and IGM, to constrain and calibrate the galaxy formation model.The Coma cluster has abundant observational data for its ICM, from the radio to the X-ray band, which provides valuable information about the evolutionary history of the Coma ICM, but has rarely been used to constrain the subgrid physics in galaxy formation models.Our main results are summarized below.
Generally, the GZ-SB model is able to reproduce the observed SMF in the Coma region, except that it significantly underestimates the population of quiescent galaxies in the LDR.This discrepancy cannot be detected using the conventional approach of considering the LDRs and HDRs together.Additionally, the GZ-SB model underestimates the stellar and ISM metallicities across the entire stellar mass range in both the LDRs and HDRs.This suggests that some internal processes are not modeled correctly, as environmental effects may be negligible in the LDR.The metallicity of the ICM is also underestimated, indicating that the total metal yield of the fiducial SIMBA model is too low.To address this problem, one may weaken the AGN feedback to increase the SFR so as to produce more metals.However, this does not work for lowmass galaxies, in which AGN feedback is not important.It seems that an increase in the SFR and in SN feedback strength is needed to fix the problem.
Our new model assuming higher SFR and SN feedback strength matches the observed metallicities in stars, the ISM, and the ICM better than the fiducial model.The predicted M * -Z * relation and its scatter are also in agreement with the observational results in both LDRs and HDRs.In the HDR, owing to strangulation, quiescent galaxies have a higher M * -Z * relation than star-forming galaxies, which is consistent with the observations.In LDRs, we find that the intrinsic scatter of the M * -Z * relation is very small, indicating that these are good places to study the origin of the M * -Z * relation.The new simulations also reproduce the ICM metallicity and its radial dependence.The amplitude of the predicted M * -Z ISM relation is consistent with the observations, but the slope is very different.Interestingly, in the new simulation, the predicted SMF for quiescent galaxies in the LDR matches the observations much better because the high SFR and SN feedback reduce the amount of cold gas and make galaxies more susceptible to AGN feedback.
The G3-H model is able to accurately reproduce the total SMF in both LDRs and HDRs.However, when star-forming and quiescent galaxies are separated, large discrepancies become apparent.It predicts too many star-forming galaxies and too few quiescent galaxies.In particular, there are almost no quiescent galaxies in the LDR, which implies that environmental effects are a major factor in quenching galaxies in the G3-H model.On the other hand, the predicted M * -Z * relation and M * -Z ISM relation for star-forming galaxies in the LDR are in good agreement with the observational data in terms of amplitude, slope, and scatter.In the HDR, it appears that the environmental effect is too strong, likely owing to the strong feedback from stars and SNe, so that star-forming and quiescent galaxies share the same M * -Z * relation, in contrast to the observational result.Additionally, the G3-H model underpredicts the ICM metallicity.All of these findings suggest that the G3-H model is successful in capturing the formation history of star-forming galaxies but fails to properly quench galaxies.An additional internal quenching mechanism, such as AGN feedback, seems to be required for the G3-H model to work.
We analyzed the physical properties of the ICM of the Coma cluster, such as electron density, electron temperature, and entropy.The predictions of the simulations at z = 0 are broadly consistent with the observational results.The GZ-SB model and its modifications can reproduce the temperature profile, while the G3-H model can recover the density and entropy profiles.We found that these properties are significantly affected by a recent merger event.Our further investigations reveal that the simulation results can be improved if the merger occurs slightly later, and so, the discrepancy may be produced by the uncertainty in the reconstruction to precisely recover the highly nonlinear merger event.Allowing for this uncertainty, the two simulations without AGN feedback, G3-H-CM and GZ-SBnA-CM, can accurately reproduce the three profiles.Simulations with AGN feedback, GZ-SB-CM, GZ-SBrw/rs-CM, can accurately reproduce the temperature and density profiles, but overestimate the entropy profile in the inner region.This is because the heating by AGN feedback stabilizes the entropy in the inner region at a high level and produces a central core.
The eROSITA X-ray observation revealed a strong bow-like shock feature in the inner region and a massive subcluster in the SW of the Coma cluster.Our two non-AGN simulations, GZ-SBnA-CM and G3-H-CM, produce similar features at similar locations, while all simulations with AGN feedback basically fail to recover these features.In our simulations, the shock is produced by the collision of the halo gas between the main Coma cluster and a subhalo.When AGN feedback is implemented, the halo gas in the subhalo is severely reduced by AGN feedback prior to the collision so that a strong shock is not produced.A similar situation is also found in the SW subcluster, which has a halo mass of M 200c = 10 13.6 h −1 M e .In GZ-SB-CM and GZ-SBrw/rs-CM, strong AGN feedback expels its gas so that its X-ray emission becomes very weak, very different from what is seen in the non-AGN simulations and in the observation.
We also investigated the properties of the IGM.Observational studies revealed that the IGM temperature (T IGM ) increases with increasing cosmic density (ρ m ).The two non-AGN models predict a T IGM -ρ m relation that is consistent with the observations within the observational uncertainty.Furthermore, our tests demonstrate that the simulated T IGM -ρ m relation is insensitive to the adopted star formation model, hydrodynamic solver, and cosmic variance.On the other hand, when AGN feedback is included, the IGM in LDRs can be heated to a very high temperature, which is very different from what is observed.
Our study indicates that, to accurately recover the SMF of both star-forming and quiescent galaxies, an internal quenching mechanism, such as AGN feedback, must be implemented.However, we find that models without AGN feedback can accurately reproduce physical properties of the observed ICM and IGM, while those with strong AGN feedback usually fail.This implies that quenching mechanisms may only operate on relatively small scales and do not significantly alter the gas properties on cluster and larger scales.Moreover, models without AGNs can also reproduce the metallicity-stellar mass relations for star-forming galaxies in the LDR, indicating that AGN feedback should not significantly affect the properties of star-forming galaxies before quenching them.
Our study demonstrates that the "one-to-one" approach provided by our CSs is a powerful tool to constrain galaxy formation models.Investigating LDRs and HDRs separately can effectively differentiate between environmental and internal effects.Additionally, the "one-to-one" study with the combined multiple-band observational data can significantly increase the constraining power of observational data.For instance, considering the star, ISM, and ICM metallicity together gives us a more comprehensive picture of star formation activity.These are particularly important for the calibration of model parameters.The physical properties of the ICM and IGM are not significantly affected by star formation models or hydrodynamic solvers and thus are ideal probes of AGN feedback.Our study shows that, with a reliable CS, the properties and structures of the ICM within one single cluster, such as the Coma cluster, can put a significant constraint on the AGN feedback.Thus, CSs, such as the Coma cluster simulations presented here, provide a unique and effective platform to test, verify, and calibrate galaxy formation models.

Figures 3
Figures 3, 4, and 5 show the SMF, stellar mass-stellar metallicity relation (M * -Z * relation), and stellar mass-ISM metallicity relation (M * -Z ISM relation) of observed galaxies in two different regions, HDR (within 3R 200c of the Coma cluster) and LDR (between 3R 200c and 7R 200c).In general, galaxy properties in the two regions differ significantly.Galaxies in HDR tend to be more massive, more frequently quenched, and metal richer than those in LDR, suggesting that environmental processes have a significant impact on galaxy properties (Wang

Figure 1 .
Figure1.The evolution of the Coma cluster in the simulation GZ-SBrw-CM from z = 1 to 0. The columns from left to right show dark matter density, gas density, gas temperature, gas metallicity, and stellar density, respectively.The circles in the z = 0 panels indicate the virial radius of the cluster at the present time.

Figure 2 .
Figure 2. Two-dimensional distributions of galaxies in the J2000.0coordinate system.The horizontal and vertical axes are, respectively, the R.A. and decl. in degrees.The six panels show results from observations and five simulations as labeled.Only galaxies with * M h M log 9.5 2 ( )  > within a radius of 20 h −1 Mpc around the Coma cluster are displayed, colored by their specific star formation rate.

Figure 3 .
Figure 3. Stellar mass functions for galaxies in the high-density region (HDR; R/R 200c < 3, upper panels) and the low-density region (LDR; 3 R/R 200c < 7, lower panels).The first column shows results for all galaxies, while the second and third columns show results for star-forming and quiescent galaxies, respectively.Different lines show results for different simulations as indicated in the lower right panel, while observational results are shown by symbols with Possion error bars.The vertical gray dashed lines represent the estimated complete mass limit of the observations.

Figure 4 .
Figure 4.The M * -Z * relation of galaxies in HDR (left panels) and in LDR (right panels).Red and blue dots represent quiescent and star-forming galaxies respectively.The first row shows observational data, and the rest of rows show the results of simulations as indicated at upper right corner of each panel.The solid lines are the median relations, while the dashed lines show the 16th-84th percentiles of the distribution.The scatter in the simulation relations represents measurement uncertainties (see the text for details).For comparison, the observational median relation and scatter (brown color) are repeated in the other panels.In the third row, we also show the median relation and scatter predicted by GZ-SBrs-CM.

Figure 5 .
Figure 5. Similar to Figure 4, but for the M * -12 log O H ( ) + relation of starforming galaxies, which is referred to as the M * -Z ISM relation in this paper.

Figure 6 .
Figure 6.ICM profiles for the Coma cluster.The four panels show electron number density (top left), temperature (top right), metallicity (bottom left), and entropy (bottom right) profiles.Symbols show observational results taken from the literature as indicated in each panel, while lines plot the simulations' results.See Section 2 for details of the observational data.
are presented so that the median T IGM -ρ m relations are shown by solid lines, and their 16th and 84th percentiles are shown by dashed lines.As one can see, the IGM temperature predicted by the GZ-SB-CM first decreases and then increases with the density.The GZ-SB-CM result is consistent with the observational data at 10

Figure 7 .
Figure 7. Scaled X-ray maps in the observations and the simulations as indicated (see the text for details).The observational data are taken from Churazov et al.(2021).The purple dashed curve marks the observational shock front with a small shift.The SW structure corresponds to the NGC 4839 subcluster.Note that there is a small offset between the simulated and real subclusters.The simulated maps adopt the same color scale as the observational map (private communication for the observational scale).The gray dashed circles represent the scale of R 200c in each simulation.

Figure 9 .
Figure 9. Evolution of the Coma ICM profiles from z = 0.1 (blue) to 0 (red) as indicated by the color-coded lines.Different columns show results for different simulations, as indicated in the first row.The first row shows electron density, the second temperature, and the third entropy.Data points are observational results taken from the literature (the same as Figure 6).