# ARCHITECTS I: impact of subgrid physics on the simulated properties of the circumgalactic medium

Maxime Rey, Jérémy Blaizot, Taysun Kimm, Joakim Rosdahl, Léo

Michel-Dansac

► To cite this version:

Maxime Rey, Jérémy Blaizot, Taysun Kimm, Joakim Rosdahl, Léo Michel-Dansac. ARCHITECTS I: impact of subgrid physics on the simulated properties of the circumgalactic medium. Monthly Notices of the Royal Astronomical Society, 2025, 543 (1), pp.12-27. 10.1093/mnras/staf1372. hal-05352139

HAL Id: hal-05352139

https://hal.science/hal-05352139v1

Submitted on 6 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# ARCHITECTS I: impact of subgrid physics on the simulated properties of the circumgalactic medium

Maxime Rey<sup>1,2,★</sup> Jérémy Blaizot<sup>1,2</sup> Taysun Kimm<sup>1</sup>,<sup>1</sup> Joakim Rosdahl<sup>2</sup> and Léo Michel-Dansac<sup>3</sup>

Accepted 2025 August 14. Received 2025 August 14; in original form 2025 April 22

## ABSTRACT

Galaxy evolution is shaped by star formation and stellar feedback at scales unresolved by current high-resolution cosmological simulations. Precise subgrid models are thus necessary, and different approaches have been developed. However, they are degenerate and often primarily calibrated to reproduce stellar masses from observations. To explore these degeneracies, we perform three cosmological zoom-in radiation-hydrodynamics simulations of the same galaxy within a $5 \times 10^{11}$ M$_\odot$ dark matter halo at $z \sim 1$, each with a different subgrid model: mechanical feedback, a combination of mechanical feedback and thermal feedback, and delayed cooling. We calibrate the simulations to match in stellar mass, isolating the effect of the models on the circumgalactic medium (CGM). Our findings demonstrate that despite producing galaxies with comparable stellar masses, the three models lead to distinct feedback modes, resulting in notable variations in the CGM properties. The delayed cooling run is dominated by ejective feedback and exhibits high burstiness, whereas mechanical and the hybrid models primarily feature preventive feedback, respectively, acting at the galaxy and halo scales. Delayed cooling reduces the baryon mass to half the universal baryon fraction, while mechanical feedback retains most baryons, with the hybrid model standing in between. Delayed cooling also ejects significantly more metals into the CGM than both other models. While for delayed cooling and mechanical feedback metals are almost evenly distributed in the CGM, they are concentrated around satellites in the hybrid model. These discrepancies emphasize the need to design an appropriate subgrid model to understand how stellar feedback regulates galaxy growth.

**Key words:** methods: numerical – galaxies: evolution.

## 1 INTRODUCTION

By allowing the simultaneous modelling of non-linear processes across vast ranges of spatial and temporal scales, numerical simulations have become an essential tool to study the complex processes that govern galaxy formation and evolution. In recent years, significant progress has been made in creating simulations that describe large representative volumes of the Universe, yielding results consistent with many observations (e.g. Dubois et al. 2014; Schaye et al. 2015; Davé et al. 2019; Nelson et al. 2019a). However, these simulations rely on macroscopic subgrid models that describe the behaviour of the interstellar medium (ISM) as a whole, prompting the development of a new generation of simulations, aimed at a deeper understanding of the fundamental physics shaping galaxy evolution (Hopkins et al. 2014, 2018, 2023; Wang et al. 2015; Grand et al. 2017; Agertz et al. 2020). These simulations focus on zooming in on specific haloes within a cosmological box, preserving large-scale structures while achieving higher resolution in the galaxy. While these simulations reach much higher resolutions, they require the development of new models to describe still unresolved processes

like star formation and feedback. Accurately capturing supernovae and mitigating the overcooling problem requires a higher resolution than typically achieved in modern zoom simulations: $\approx 4$ pc (Kim & Ostriker 2015; Simpson et al. 2015) in Eulerian simulations. If these conditions are not met, most of the supernova's energy is lost before a shockwave can form (Katz 1992; Cho & Kang 2008; Ceverino & Klypin 2009; Chaikin et al. 2022). One method used to overcome this limitation is delayed cooling (Stinson et al. 2006; Governato et al. 2010; Agertz, Teyssier & Moore 2011; Teyssier et al. 2013), where cooling is temporarily disabled to prevent spurious energy loss immediately after a supernova. More recently, Kimm & Cen (2014) and Kimm et al. (2015) introduced mechanical feedback (see also Hopkins et al. 2014), where energy derived from high-resolution simulations is injected as momentum. Some approaches, such as that of Kretschmer & Teyssier (2020), combine thermal energy and mechanical feedback.

The impact of subgrid models is often tested in either high-resolution idealized galaxy simulations (Hopkins, Narayanan & Murray 2013; Rosdahl et al. 2017; Smith, Sijacki & Shen 2018; Smith 2021; Andersson et al. 2023) or in zoom-in cosmological simulations (Kimm & Cen 2014; Nuñez-Castiñeyra et al. 2021; Azartash-Namin et al. 2024; Kang et al. 2025). Although different models can yield similar galaxy properties, they can still produce distinct feedback

<footnote n="1"><sup>1</sup>Center for Galaxy Evolution Research and Department of Astronomy, Yonsei University, 50 Yonsei-ro, Seodaemun-gu, Seoul 03722, Republic of Korea</footnote>

<footnote n="2"><sup>2</sup>Centre de Recherche Astrophysique de Lyon UMR5574, ENS de Lyon, Univ. Lyon, Univ Lyon 1, CNRS, 9 av. Charles André, F-69230 Saint-Genis-Laval, France</footnote>

<footnote n="3"><sup>3</sup>Aix Marseille Univ., CNRS, CNES, Laboratoire d'Astrophysique de Marseille, F-13388 Marseille, France</footnote>

<footnote>*E-mail: maximerey.astro@gmail.com</footnote>
mechanisms such as ejective, where gas is expelled from the galaxy, preventive feedback, which prevents gas from falling on to the galaxy or galactic fountains where gas that has been ejected eventually falls back on to the galaxy (Christensen et al. 2016; Mitchell & Schaye 2022b). Discriminating between these models and refining them remains a key challenge for new-generation simulations. Most of these comparative studies focus solely on the stellar or gaseous properties of the galactic disc, often overlooking the crucial role of the circumgalactic medium (CGM). This highly multiphase multiscale environment (Tumlinson, Peebles & Werk 2017; Faucher-Giguère & Oh 2023; Peroux & Nelson 2024; Sameer et al. 2024) surrounds galaxies and this is where inflows from cosmological accretion interact with feedback-driven outflows (Crain et al. 2013; Hafen et al. 2019; Afruni, Fraternali & Pezzulli 2021; Saeedzadeh et al. 2023; Decataldo et al. 2024). Understanding how these flows are regulated is essential to unveil the processes that govern star formation beyond ISM scales. Given its critical role, many studies investigate the properties of the CGM (Mitchell et al. 2018; Afruni et al. 2021, 2023; Decataldo et al. 2024), and observations of the CGM are becoming increasingly detailed (Chen & Zahedy 2024; Fumagalli 2024), offering complementary constraints that could help differentiate between models.

Simulations studying the CGM must incorporate a cosmological context, as seen in Davies et al. (2020), Oppenheimer et al. (2020), and Medlock et al. (2025), where CGM properties were compared across the EAGLE (Crain et al. 2015; Schaye et al. 2015), SIMBA (Davé et al. 2019), and ILLUSTRISTNG (Pillepich et al. 2019; Nelson et al. 2019b) cosmological simulations. However, full cosmological simulations typically lack ISM resolution, and zoom-in simulations address this by focusing on a single halo at high resolution, though at the cost of the loss of large-number statistics. Kelly et al. (2022) compared the CGM of two simulations suites using this approach, namely the AURIGA (Grand et al. 2017) and APOSTLE simulations (Fattahi et al. 2016; Sawala et al. 2016). Nevertheless, the minimum softening length in these simulations is only about 300 pc. To the best of our knowledge, this represents the highest resolution comparison of the impact of subgrid models on the CGM to date. Moreover, in CGM comparisons from both full cosmological simulations (Davies et al. 2020; Oppenheimer et al. 2020; Medlock et al. 2025) and zoom-in simulations (Kelly et al. 2022), not only do subgrid physics change, but also the entire simulation setup, such as the code, the settings (e.g. the resolution), and initial conditions, making it difficult to isolate the specific impact of subgrid models on the CGM. In this paper, we take a similar approach but focus on the impact of star formation and feedback models on the CGM with much higher resolution. We do not aim to identify the best subgrid model, but rather to emphasize the critical role these models play in shaping the properties of the CGM, even when they produce similar stellar masses.

To this end, we introduce the ARCHITECT suite of simulations, designed to Assess the Robustness of Computational Hydrodynamics and ISM-scale Theories by Exploiting Circumgalactic Tracers. We run three high-resolution radiation-hydrodynamics simulations, based on identical initial conditions, differing only in the subgrid models employed. Each of the three simulations uses a different set of star formation and supernova feedback models often used in the literature, which we refer to as ME (mechanical feedback; Kimm et al. 2017), MT (mechanical + thermal; Kretschmer & Teyssier 2020), and DC (delayed cooling; Teyssier et al. 2013), as described in Section 2. Since CGM properties are expected to depend on galaxy mass, as indicated by both observational (Bordoloi et al. 2018; Dutta et al. 2021) and theoretical studies (Nelson et al. 2020; Anand, Nelson & Kauffmann 2021; Cook et al. 2024), we calibrate the stellar

mass of the three galaxies to isolate the impact of feedback. With similar galaxies in place, we present their properties and compare the effects of the three subgrid models on the CGM in Section 3, before concluding on our findings in Section 4.

## 2 NUMERICAL METHODS

In this section, we briefly describe the code we use to run our simulations and how we generate the initial conditions. We also list the physics included in the simulations and detail how the subgrid models work, focusing on the differences between the three approaches.

### 2.1 Numerical setup

The simulation setup and physics included closely follow those used in SPHINX (Rosdahl et al. 2018). We perform our simulations with the hydrodynamical code RAMSES (Teyssier 2002, 2010), relying on an Adaptive Mesh Refinement strategy and based on the Euler conservation equations in the presence of self-gravity and cooling. Non-equilibrium cooling is included alongside ionizing radiative transfer (Rosdahl et al. 2013; Rosdahl & Teyssier 2015).

Dark matter and stars are modelled as collisionless particles interacting gravitationally through a particle-mesh solver and cloud-in-cell interpolation (Guillet & Teyssier 2011). Their minimum masses are, respectively, $m_{\text{DM}} = 3.49 \times 10^5 \text{ M}_\odot$ and $m_* = 3.2 \times 10^3 \text{ M}_\odot$.

We use MUSIC (Hahn & Abel 2011, 2013) to generate cosmological initial conditions for the simulations presented here, assuming cosmological parameters $\Omega_\Lambda = 0.6825$, $\Omega_m = 0.3175$, and $\Omega_b = 0.049$ (Planck Collaboration XVI 2014). Our zoom-in simulations describe a halo of $M_h = 5.33 \times 10^{11} \text{ M}_\odot$ at $z=0$, which was selected in a $30 \text{ cMpc h}^{-1}$ box to be relatively isolated, with no massive sub-structure and no neighbour more massive than $0.2 \times M_h$ within 3 $R_{\text{vir}}$, with $R_{\text{vir}}$ the virial radius. The zoom-in region extends beyond 3 $R_{\text{vir}}$ at all redshifts down to $z=1$. The resolution of our simulation ranges from level 7 ($128^3$ particles/cells at the coarse level) to level 20 ($8192^3$ effective resolution, $\Delta x \approx 40$ pc) in the zoom region. To illustrate how our simulations compare to other simulations in terms of resolution, we show in Fig. 1 a projection of the cell width of one of our runs (ME), as well as the cell width profile of the three ARCHITECT simulations as a function of radius. We normalize the x-axis of the bottom panel by $R_{200}$, the radius within which the mean density of the medium is 200 times the critical density of the Universe. We compare our resolution to Suresh et al. (2019) (the L11.12 simulation) and Ramesh & Nelson (2024) (the GIBLE simulation), two zoom simulations that aimed at increasing the CGM resolution to resolve cold gas structures. Beyond $R_{200}$ and down to $\approx 0.5 R_{200}$, our simulations are close to those of L11.12, although ME and DC have somewhat worse resolution than MT. In the inner CGM ($\lesssim 0.5 R_{200}$), our simulation has a cell width twice lower than L11.12, and two to eight times lower than GIBLE, the highest difference being close to the galaxy, at $0.2 R_{200}$. At lower radii, while the resolution decreases sharply in GIBLE to reach 450 pc in the ISM, our simulations exhibit a sharp increase in resolution down to a roughly 10 times higher resolution than GIBLE, at 40 pc. L11.12 exhibits similar ISM resolution to our simulation, albeit a bit lower, around 100 pc.

We use spectral energy distributions from the BINARY POPULATION AND SPECTRAL SYNTHESIS code (v2.2.1, Stanway & Eldridge 2018) and trace the ionizing radiation in the simulations using three photon groups defined by the ionization energies of H, He, and He: [13.6, 24.59], [24.59, 54.42], [54.42, $\infty$). We do not account for
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Top: projection of the cells width for one of the ARCHITECT simulations (ME) at redshift $z = 1$ in parsec. We trace two circles in dashed and dotted lines, with a respective radius of 0.1 $R_{200}$ and $R_{200}$. Bottom: cell width as a function of radius in units of $R_{200}$ for the three ARCHITECT simulations (ME, MT, and DC, presented hereafter). The solid line shows the median resolution over 1 Gyr (from $z = 1.3$ to $z = 1$) and the shaded area corresponds to the 15.9 and 84.1 percentiles. We compare our simulation's resolution to Suresh et al. (2019) and Ramesh & Nelson (2024).</figcaption></figure>

radiation below the hydrogen ionization energy. The interactions simulated for these groups are photoionization, cooling and heating, and momentum transfer (Rosdahl & Teyssier 2015). Radiation flows through the grid and thus constrains the Courant condition. We hence use the reduced speed of light approximation with $c_{sim} = 0.005$ c (Gnedin & Abel 2001; Rosdahl et al. 2013) to limit the computational cost. We turn on a uniform ionization UV field at $z = 8.5$, modelling emission from unresolved surrounding galaxies, following Faucher-Giguère et al. (2009), and take into account self-shielding with an exponential damping factor above densities of $10^{-2}$ cm$^{-3}$ (Schaye 2001; Rosdahl & Blaizot 2012).

We trace the ionization states of atomic hydrogen and helium, and track metals through the single variable Z, describing the metal mass fraction. Non-equilibrium cooling is computed for hydrogen and helium via photoionization, collisional ionization, and excitation, recombination, bremsstrahlung, homogeneous Compton cooling/heating off the cosmic microwave background, and dielectronic recombination. Above $10^4$ K, metal cooling is modelled following CLOUDY tables (Ferland et al. 1998, version 6.02) assuming photoionization equilibrium with the UV background from Haardt &

Madau (1996). Below $10^4$ K, fine structure metal cooling rates from Rosen & Bregman (1995) are used down to 10 K.

In the zoom-in region, two criteria can trigger refinement in the AMR grid. The first criterion follows density and increases the resolution when a cell's mass (DM and baryonic) is greater than a selected threshold. For the cell to be refined, the condition is $m_{DM,cell} + \frac{\Omega_m}{\Omega_b} m_{b,cell} \ge 8m_{DM}$, where $m_{DM,cell}$ is the dark matter mass in the cell, $m_{b,cell}$ the baryonic mass in a cell (gas + stellar mass), and $m_{DM}$ the mass of a single dark matter particle.

The second refinement criterion aims to keep the local thermal Jeans length $\lambda_J$ resolved by at least four cells, down to the resolution limit. Thus, a cell is refined whenever

$$4\Delta x \ge \lambda_J = \sqrt{\frac{\pi c_s^2}{\rho G}}, \qquad (1)$$

with $\rho$ the gas density, $G$ the gravitational constant, and $c_s$ the sound speed.

Feedback from massive stars in the form of photoionization, heating, and radiation pressure from young stars is accounted for in our simulations via the radiation-hydrodynamics solver (Rosdahl et al. 2013; Rosdahl & Teyssier 2015). We nevertheless need two subgrid models to describe the formation of stars and the effect of explosions of Type II supernovae. In this paper, we test and compare three models for this, which we describe in Sections 2.2–2.4.

## 2.2 ME simulation

ME uses a Schmidt-based thermo-turbulent subgrid model for star formation (Kimm et al. 2017), and mechanical supernova feedback by Kimm et al. (2015). This model is used in many simulations such as SPHINX (Rosdahl et al. 2018, 2022), SPHINX-MHD (Katz et al. 2021), NewHorizon (Dubois et al. 2014), OBELISK (Trebitsch et al. 2021), or the NEWHORIZON2 model (Yi et al. 2024). We first summarize how the multi-free-fall model for star formation is applied in ME. Then, we describe the radiative feedback used in all our simulations alongside the modelling of supernovae feedback in ME.

### 2.2.1 Star formation

Star formation can be broadly described in three steps: (1) determining whether a cell is star-forming, (2) calculating the local SFR, and (3) implementing the SFR over a time-step dt using finite-mass star particles.

In ME, we allow star formation if a cell has a gas hydrogen number density higher than $10 \text{ cm}^{-3}$, is a local maximum in density compared to its six neighbouring cells, is in a converging flow ($\nabla \cdot (\rho \vec{v}) \le 0$), and has a size $\Delta x \ge \lambda_{J,turb}/4$, where $\lambda_{J,turb}$ is the turbulent Jeans length:

$$\lambda_{J,turb}(\rho) = \frac{\pi \sigma_{\text{gas}}^2 \pm \sqrt{36\pi c_s^2 G \Delta x^2 \rho + \pi^2 \sigma_{\text{gas}}^4}}{6G \Delta x \rho}. \qquad (2)$$

The three-dimensional gas velocity dispersion $\sigma_{\text{gas}}$ in equation (2) is computed locally as the norm of the gradient of the velocity field from the neighbouring cells.

Once a cell is flagged for star formation, it will form stars at a rate given by

$$\dot{\rho}_{\star} = \rho \frac{\epsilon_{\text{ff}}}{t_{\text{ff}}}, \qquad (3)$$

with $t_{ff} = [3\pi/(32G\rho)]^{1/2}$ the free-fall time. The star formation efficiency per free-fall time ($\epsilon_{ff}$) is modelled following the multi-
free-fall approach laid out by Hennebelle & Chabrier (2011) and Federrath & Klessen (2012). The main assumption is that star formation occurs in a supersonic turbulent medium so that the unresolved density distribution in the cell follows a lognormal function. Then, this distribution can be integrated over a critical density $s_{crit}$ to define the star formation efficiency

$$ \epsilon_{ff} = \epsilon_{loc} \exp \left( \frac{3}{8} \sigma_s^2 \right) \left[ 1 + \operatorname{erf} \left( \frac{\sigma_s^2 - s_{crit}}{\sqrt{2\sigma_s^2}} \right) \right]. \quad (4) $$

Here, erf is the Gauss error function and $\sigma_s$ is the standard deviation of the density distribution. We write $\sigma_s^2 = \ln(1 + b^2 M^2)$, where $M$ is the Mach number, and $b = 0.4$ the turbulent forcing parameter that describes the ratio of solenoidal to compressive modes. The Mach number $M = \sigma_{gas}/c_s$ depends on the sound speed $c_s$ and the velocity dispersion $\sigma_{gas}$. The model from Padoan & Nordlund (2011), which we use in ME, uses

$$ s_{crit} = \ln[0.62\alpha_{vir}M^2], \quad (5) $$

where

$$ \alpha_{vir} = \frac{5}{\pi \rho G} \frac{\sigma_{gas}^2 + c_s^2}{\Delta x^2} \quad (6) $$

is the virial parameter describing the gravitational stability of the cloud. The local star formation efficiency $\epsilon_{loc} = 0.1425$ is a constant used to account for unresolved physics such as protostellar jets removing gas from the star-forming region and the estimated error on the time-scale for the gas to become unstable (Krumholz & McKee 2005). When the Mach number gets large, turbulent compression gets more efficient and the global star formation efficiency $\epsilon_{ff}$ increases. Conversely, as the virial parameter gets larger, the cloud gets more gravitationally unstable and $\epsilon_{ff}$ decreases (see fig. 1 in Federrath & Klessen 2012).

Finally, the mass of stars formed is $N$ times $m_*$, where $N$ is an integer number drawn from a Poisson distribution of mean $SFR \times \Delta t$ (Rasera & Teyssier 2006). The stellar particle created represents a simple stellar population defined through its mass, metallicity, age, position, and velocity. We set the minimum mass of the stellar particles in our simulation to $\approx 3200 M_\odot$.

### 2.2.2 Feedback

We assume that all Type-II supernovae explode 10 Myr after a stellar particle is spawned with a rate of two supernovae per $100 M_\odot$ of stellar mass formed.

We artificially boost the supernova rate by a factor of 2 to better reproduce the stellar mass to halo mass relation and match the stellar mass of the other models. Boosting feedback by, e.g. increasing the supernova frequency (Rosdahl et al. 2018), decreasing the minimum mass of stars undergoing supernovae (Dubois et al. 2021), or increasing the supernovae efficiency (Li, Gnedin & Gnedin 2018; Yi et al. 2024) is common practice in supernovae models as galaxies otherwise tend to be too massive. This calibration can be seen as compensation for unaccounted physics such as cosmic rays pressure, proto-stellar jets, or to model the unresolved clustering of supernovae, which can increase the momentum transfer to the ISM (Keller et al. 2014; Kim & Ostriker 2015).

The method used in Kimm et al. (2015) determines in which phase the supernova explosion is when reaching the neighbour cells, depending on the local gas conditions, and injects momentum and energy accordingly.

In the snowplough limit where much gas has been swept up ($\sim 100$ times the ejecta), the model injects (Kimm et al. 2017, and references therein)

$$ p_{\text{rad}} = \left[ 2.5 \times 10^5 N_{\text{SN}}^{16/17} n_{\text{H}}^{-2/17} Z'^{-0.14} \right] e^{-\Delta x/R_{\text{st}}} + \left[ 5 \times 10^5 N_{\text{SN}}^{16/17} Z'^{-0.14} \right] (1 - e^{-\Delta x/R_{\text{st}}}), \quad (7) $$

where $N_{SN}$ is the number of supernovae at a given time-step, $n_H$ the density of the neighbouring cell, $Z' = \max(Z, 0.01)$, with $Z$ the metallicity in solar units, $\Delta x$ the size of the host cell, and $R_{st}$ the Strömgren radius. If the Strömgren radius is not sufficiently resolved, the impact of photoionization is likely to be underestimated. In this case, the second term of $p_{rad}$, fitted from Geen et al. (2015), will dominate the equation. Conversely, if the Strömgren radius is well-resolved, the first term based on Thornton et al. (1998) will prevail.

In the other limit, where little gas has been swept up, the Sedov-Taylor phase (energy conserving) is well resolved and momentum is injected as

$$ p_{\text{ad,ME}} = B_{\text{ph}} \sqrt{2\chi m_{\text{ej}} f_e E_{\text{SN}}}, \quad (8) $$

where $f_e$ is a coefficient ensuring a smooth transition between the two propagation phases, $B_{ph}$ is a boost factor accounting for unresolved photoionization (Geen et al. 2015), $\chi$ is the ratio of the mass swept up by the supernova to the mass ejected in each cell, $m_{ej}$ is the total mass of the ejecta, and $E_{SN}$ is the energy released by a single supernova.

Along with energy and momentum, SN explosions release metals, and we model the mass fraction of metal in the ejecta with a yield of 7.5 per cent, as $Z_{ej} = Z_{star} + 0.075(1 - Z_{star})$.

## 2.3 MT simulation

We now detail the subgrid models used in the MT simulation, which come from Kretschmer & Teyssier (2020). This model is used for example in the MIGA simulations (Kretschmer, Dekel & Teyssier 2022), the EMP-Pathfinder suite (Reina-Campos et al. 2022), and the SPICE suite of simulations (Bhagwat et al. 2024).

### 2.3.1 Star formation

Here, the model for star formation is also an implementation of the multi-free-fall model, with a number of subtle differences compared to the model from Kimm et al. (2017). One major difference is the definition of $\sigma_{gas}$, which corresponds to the one-dimensional turbulence. Kretschmer & Teyssier (2020) compute turbulence from the velocity and density fields, relying on numerical diffusion to provide an implicit subgrid model (Schmidt, Niemeyer & Hillebrandt 2006; Semenov, Kravtsov & Gnedin 2016).

Furthermore, this models sets slightly different parameters values in equation (4) (see Table 1), and relies on a critical density given by (Krumholz & McKee 2005)

$$ s_{\text{crit}} = \ln \left[ \alpha_{\text{vir}} \left( 1 + \frac{2M^4}{1 + M^2} \right) \right]. \quad (9) $$

While the role of $\alpha_{vir}$ in $s_{crit}$ is the same in both MT and in ME, the Mach number has a distinct impact in the two models. In particular, for large $M$, $s_{crit}$ in MT approaches the form used in ME, as $s_{crit} \rightarrow \ln[2\alpha_{vir}M^2]$, albeit higher by 1.17. Consequently, MT limits star formation to denser media than ME, leading to a seemingly smaller star formation efficiency. However, this is counterbalanced by a higher local star formation efficiency $\epsilon_{loc} = 0.5$, as shown in Table 1.
<table><caption>**Table 1.** Summary of the main differences between the subgrid models used in ME, MT, and DC. The three first rows describe the differences in the star formation subgrid models, i.e. the critical density $s_{crit}$, the local star formation efficiency $\epsilon_{loc}$, and the method to estimate turbulence. The following rows are the differences in the feedback subgrid models, namely the time sampling of the supernovae explosions and its boundaries, the average supernova rate (the number of supernovae per 100 M$_{\odot}$), the method for the thermal energy injection, the equation for the momentum injected in the adiabatic phase of the supernova, the transition and the radiative phase, where this momentum is injected, and the metal yield of the ejecta.</caption><thead><tr><th></th><th>ME</th><th>MT</th><th>DC</th></tr></thead><tbody><tr><td>s<sub>crit</sub></td><td>ln [0.62α<sub>vir</sub>M<sup>2</sup>]</td><td>ln [α<sub>vir</sub> (1 + <sup>2M<sup>4</sup></sup>/<sub>1+M<sup>2</sup></sub>)]</td><td>ln [0.62α<sub>vir</sub>M<sup>2</sup>]</td></tr><tr><td>ε<sub>loc</sub></td><td>0.1425</td><td>0.5</td><td>0.1425</td></tr><tr><td>Turbulence σ<sub>1D</sub></td><td>Local estimation</td><td>Subgrid turbulence model</td><td>Local estimation</td></tr><tr><td>Supernova sampling</td><td>Single explosion</td><td>Uniform time sampling</td><td>Single explosion</td></tr><tr><td>Time of supernova</td><td>10 Myr</td><td>0 to 20 Myr</td><td>10 Myr</td></tr><tr><td>Average supernova rate</td><td>2 SN/100 M<sub>☉</sub></td><td>8 SN/100 M<sub>☉</sub></td><td>2 SN/100 M<sub>☉</sub></td></tr><tr><td>Thermal energy injection</td><td>None</td><td>Random cell within the 8 parent cell</td><td>In the central cell + cooling turned off</td></tr><tr><td>Momentum in the adiabatic phase</td><td>B<sub>ph</sub>√2χm<sub>ej</sub>f<sub>e</sub>E<sub>SN</sub></td><td>None</td><td>None</td></tr><tr><td>Transition</td><td>χ ≡ dm<sub>W</sub>/dm<sub>ej</sub></td><td>R<sub>cool</sub> = 43.6n<sub>H</sub><sup>-0.42</sup>Z<sup>-0.05</sup> pc</td><td>None</td></tr><tr><td>Momentum in the radiative phase [M<sub>☉</sub> km s<sup>-1</sup>]</td><td>2.5 × 10<sup>5</sup>N<sub>SN</sub><sup>16/17</sup>n<sub>H</sub><sup>-2/17</sup>Z<sup>-0.14</sup>e<sup>-Δx/R<sub>st</sub></sup> + 5 × 10<sup>5</sup>N<sub>SN</sub><sup>16/17</sup>Z<sup>-0.14</sup>(1 - e<sup>-Δx/R<sub>st</sub></sup>)</td><td>2.66 × 10<sup>5</sup>η<sub>R</sub>N<sub>SN</sub>Z<sup>-0.114</sup>n<sub>H</sub><sup>-0.19</sup></td><td>None</td></tr><tr><td>Momentum injection sites</td><td>Neighbouring cells</td><td>Neighbour cells + source term</td><td>None</td></tr><tr><td>Metal yield</td><td>7.5 per cent</td><td>10 per cent</td><td>7.5 per cent</td></tr></tbody></table>

### 2.3.2 Feedback

The original model of Kretschmer & Teyssier (2020) includes a subgrid description for radiative feedback. As we already include RT, we turn this component of the subgrid models off to avoid double-counting. Individual supernova explosions are modelled as taking place for each stellar particle at ages ranging from 3 to 20 Myr, with a uniform distribution in time. The location at which the supernovae explode is chosen randomly over the eight cells comprised in the oct where the stellar particle is located. At each supernova explosion, $10^{51}$ erg is injected as thermal energy into the host cell. Kretschmer & Teyssier (2020) also compute a local cooling scale (Martizzi, Faucher-Giguère & Quataert 2015)

$$R_{\text{cool}} = 43.6 n_{\text{H}}^{-0.42} Z^{-0.05} \text{ pc}, \qquad (10)$$

which they compare to the cell size to decide when unresolved cooling losses are important. When $R_{\text{cool}} > 4\Delta x$, MT also inject the momentum expected at the radiative phase as

$$p_{\text{rad,MT}} = 2.66 \times 10^5 \eta_R N_{\text{SN}} Z^{-0.114} n_{\text{H}}^{-0.19} M_{\odot} \text{ km s}^{-1}. \qquad (11)$$

$\eta_R$ is a factor dependent on how resolved the cooling radius is. As done for the previous model, we enhance the effect of supernova feedback, this time by boosting the supernova rate by a factor of 4 compared to the fiducial value from Kretschmer & Teyssier (2020), resulting in eight supernovae explosions per 100 M$_{\odot}$. A further difference with ME is that this momentum is not simply injected as momentum in the surrounding cells but also treated as a source term at the solver level, following Agertz et al. (2013). The ejected mass is obtained by taking the product of the number of supernovae and their average mass and has a metal yield of 10 per cent.

## 2.4 DC simulation

### 2.4.1 Star formation

We use the same star formation recipe as for ME here.

### 2.4.2 Feedback

We refer to this last model as DC, as it follows the delayed cooling feedback model from Teyssier et al. (2013). In this model, the

supernova energy is injected in the form of thermal energy, and the overcooling problem is solved by inhibiting cooling locally. This prescription is used for example in the MIRAGE simulations (Perret et al. 2014) as well as in many comparative studies of feedback models (Dale 2015; Gabor et al. 2016; Rosdahl et al. 2017; Nuñez-Castiñeyra et al. 2021), and similar implementations are commonly used in major simulation suites with other codes, such as NIHAO (Wang et al. 2015).

As in the ME simulation, star particles undergo all their supernovae explosions at once, 10 Myr after they were created. The total return mass and supernovae energy are deposited in the host cell with the same yield as in ME. In order to prevent excessive cooling, cooling is turned off for a given amount of time after the supernova energy is injected. The duration depends on an additional variable $\epsilon_{DC}$, which increases every time a supernova occurs. This quantity is advected with the gas and is only used as a passive scalar that is exponentially damped at each time-step with a characteristic dissipation time-scale of 10 Myr (Teyssier et al. 2013).

We summarize the main differences between ME, MT, and DC in Table 1. The primary distinction for star formation lies in the varying methods used to estimate the velocity dispersion $\sigma_{1D}$. For feedback, major changes are the energies and supernova rates used in each model, the time sampling of the supernova (single explosion in DC and ME sampled in time in MT), and how the supernova is modelled: mechanical feedback (ME), thermal feedback combined with mechanical feedback (MT), or delayed cooling (DC). ME also includes compensation for unresolved H II regions. The observed differences observed in the resulting galaxies might thus come from different implementations, but could also be due to variations in the parameter choices. As this paper aims to highlight differences between the models as a whole, we do not try to find the exact mechanism causing these differences.

## 3 GLOBAL PROPERTIES

Our objective is to create galaxies with similar stellar masses with the three different models and test whether these models produce different CGM properties. In this section, we first focus on the stellar properties of the calibrated galaxies, before comparing the global gas properties at halo scale.
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Star formation histories (top) for mechanical feedback (ME, blue), thermal feedback combined with mechanical feedback (MT, orange), and delayed cooling (DC, green), and burstiness (bottom). The star formation history is computed using a 100 Myr moving average for clarity, and the burstiness is quantified through a 1 Myr moving average to highlight short-term fluctuations. Following equation (12), we compute the burstiness from 2 Gyr to the end of the simulation (markers) over moving periods of 30 and 500 Myr (left and right side of the violins). The violins show the PDF of the burstiness, and the horizontal bar corresponds to the median.</figcaption></figure>

## 3.1 Star formation histories

In the top panel of Fig. 2, we show the star formation rate as a function of time. We compute it from the age and initial mass of all stellar particles within the galaxy at the last snapshot, and smooth it over 100 Myr. Here and for the rest of the paper, we define all stellar properties within 0.1 $R_{200}$. The three simulations have a comparable level of star formation at late times ($z < 3$) but their star formation histories show distinct features. The three simulations start forming stars simultaneously, with an initial star formation rate sharply increasing. The three models exhibit contrasting behaviour, with MT reaching $\approx 3.4 M_{\odot} \text{yr}^{-1}$, and both ME and DC reaching $\approx 1.6 M_{\odot} \text{yr}^{-1}$. The three simulations also reveal different levels of burstiness, with DC showing the widest variations, in the form of deep dips where the star-forming activity drops by one order of magnitude on a short-but-significant time-scale.

Analysing the same models (albeit using feedback sampled in time with ME) with much greater time resolution around a star formation burst in a similar simulation, Rey (2022) found that stellar particles form through distinct pathways with different models. While stars form in highly supersonic environments due to turbulent compression in MT, they form in mildly turbulent environments with much lower virial parameters in ME and DC (see fig. 6.3 in Rey 2022).

To quantify how the burstiness varies among the three simulations, we compute a burstiness parameter inspired by Goh & Barabasi

(2006), as

$$B_{\Delta t} = \frac{\sigma_{\Delta t} / \mu_{\Delta t} - 1}{\sigma_{\Delta t} / \mu_{\Delta t} + 1}, \quad (12)$$

with $\mu_{\Delta t}$ and $\sigma_{\Delta t}$ the mean and standard deviation of the star formation rate over a time $\Delta t$. When the standard deviation is much smaller than the mean star formation rate ($\sigma_{\Delta t} \ll \mu_{\Delta t}$), the burstiness $B$ approaches $-1$, indicating no burstiness. Conversely, when the standard deviation is large, $B$ tends to 1, representing maximal burstiness. We compute the mean and standard deviation of the star formation rate as

$$\mu_{\Delta t} = \frac{1}{\Delta t} \sum_{i \in [t \pm \Delta t/2]} \text{SFR}_i^1, \quad (13)$$

$$\sigma_{\Delta t}^2 = \frac{1}{N} \sum_{i \in [t \pm \Delta t/2]} (\text{SFR}_i^1 - \mu_{\Delta t})^2, \quad (14)$$

where $N$ is the number of bins over $[t \pm \Delta t/2]$, and $\text{SFR}_i^1$ is the star formation rate smoothed with a 1 Myr sliding average. In all our computations, we restrict the star formation rates to $t > 2$ Gyr, when the galaxy in ME becomes more massive than 20 per cent of its final stellar mass.

In the bottom panel of Fig. 2, we show the resulting estimate of burstiness for $\Delta t = 30, 500$ Myr and over the whole simulation, with $t > 2$ Gyr. Over the entire simulation (shown with crosses), we find a burstiness of $B \approx -0.14$ and $B \approx 0.02$ for ME and MT, and $B \approx 0.12$ for DC, confirming the overall higher burstiness observed in the star formation rate with DC.

To get a sense of how the burstiness differs, we additionally compare the burstiness over $\Delta t = 500$ Myr, roughly encapsulating the large bursts observed in DC, and a shorter time-scale of $\Delta t = 30$ Myr. The left (right) side of each violin represents the probability density function (PDF) of the burstiness over moving windows of 30 Myr (500 Myr). Focusing first on the larger bursts, we find that ME and MT exhibit a similar main peak of burstiness at $B_{500} \approx -0.3$. MT exhibits another region of even lower burstiness, corresponding to times when the SFR is almost flat, such as over 3.2–3.8 Gyr. DC stands far above ME and MT with a positive median burstiness.

Focusing on the smaller time-scale ($\Delta t = 30$ Myr), we find that the distribution is quite different and much lower for all simulations. Although the median burstiness in ME decreases, ME's burstiness distribution is less impacted by the change of time-scale than the distributions of other models. MT drops to much lower values. This possibly hints at supernovae disrupting star-forming clouds and hindering star formation faster in ME than in MT. This could be because energy injection is uniformly sampled up to 20 Myr in MT, whereas in ME, energy is injected instantaneously at 10 Myr. It is also likely that the sampled injection times in MT slow down star formation in the cloud so that it can form stars for a longer period before reaching a critical state and exploding. For ME, there is no supernova feedback for the first 10 Myr, and star formation can thus proceed very rapidly, producing after 10 Myr a large and synchronous number of supernovae events that disrupt the cloud. While the median burstiness in DC also decreases, it remains far above other models.

We thus find very distinct behaviour in the three simulations. DC is significantly more bursty than both other models at all time-scales considered, a direct consequence of more disruptive feedback leading to frequent suppression of star formation. We also find contrasting behaviour between ME and MT. ME exhibits a slightly lower burstiness on 30 Myr time-scales compared to 500 Myr time-scales, whereas in MT the difference is much more pronounced. This highlights the distinct star formation regulation in the three models, and the role of
<figure><img src="image_4.png" /></figure>

**Figure 3.** Stellar mass-to-halo mass relation. Each marker represents a different simulation, as indicated in the legend. The colour code shows the redshift, which goes down to $z=1$ (blue). The solid blue line corresponds to the empirical model from Behroozi et al. (2019), with the blue shaded area showing the 16th–84th percentiles range. The grey shaded area indicates the broad contour of predictions of Moster et al. (2021) through reinforced learning¹. Both predictions are at redshift $z=1$.

both star formation and feedback subgrid models, although we see with DC that feedback plays a critical role.

In Fig. 3, we show the stellar mass-to-halo mass relation for our three simulations. The halo mass is measured from a twin dark-matter-only simulation to match the approach of Behroozi et al. (2019) and Moster et al. (2021) since it can vary with baryonic physics (Velliscig et al. 2014).

Due to stronger feedback, the stellar mass remains roughly twice lower in DC than in the other two models down to the end of the simulations at $z=1$. Near redshift $z=2.5$, numerous mergers in rapid succession lead to a quick increase in the stellar mass of the galaxy in all simulations and results in the step that can be seen at a halo mass of $M_{\text{halo,200}} \approx 2-3 \times 10^{11} M_{\odot}$.

Our simulations are above the stellar mass-to-halo mass mass relation inferred from observations at $z=1$. As galaxies get more massive, their gravitational potential becomes too strong for the supernovae feedback models to regulate star formation, and AGNs (which are not simulated here) may take over as the driver of feedback to solve this problem (Silk & Rees 1998; Somerville & Davé 2015; Naab & Ostriker 2017). Alternatively, different star formation models leading to bursty star formation may produce stronger outflows that inhibit galaxy growth (Kang et al. 2025). The high stellar mass might be explained by the history of the simulated galaxy, as we only follow a single set of initial conditions, or by the stochasticity of galaxy formation (Genel et al. 2019; Keller et al. 2019). Lastly, some uncertainties still weigh on these relations, as we can see with the large error bars from Moster et al. (2021). Sources of such uncertainties include underestimating the observed stellar mass (Munshi et al. 2013).

We saw that our three simulations produce galaxies with similar stellar masses, slightly above the stellar mass-to-halo mass relation inferred from observations. The ME and MT runs yield very similar masses, despite slightly different, moderately bursty star formation histories. The DC model has very strong feedback, which shuts down star formation over a significant fraction of the time, resulting in a galaxy about half the mass of the other runs. These properties being understood, we now move on to analyse the properties of the gas in these three simulations.

<figure><img src="image_5.png" /></figure>

**Figure 4.** Phase diagrams of the gas contained within $R_{200}$ for each simulation. The result consists in 100 snapshots stacked over 1 Gyr from $z=1.3$ to $z=1$. From top to bottom are ME, MT, and DC with a colour code corresponding to the hydrogen mass normalized probability density function. The solid line contour encompasses 90 per cent of the mass within $R_{200}$, and the outline in the orange dotted line is that of ME. We show the mass-weighted temperature and density probability distribution functions on the left and lower sides of each panel.

## 3.2 The multiphase nature of ISM and CGM gas

In Fig. 4, we compare the density–temperature distributions of the gas in the dark matter haloes of ME (top), MT (middle), and DC (bottom). Each panel shows the stacked phase diagram of 100 outputs from redshift $z=1.3$ down to redshift $z=1$.

The gas in all simulations is split into three main phases. The first phase is the hot ($10^5–10^6$ K) and low-density ($10^5–10^3$ cm$^{-3}$) shock-heated gas. At higher densities, cooling becomes efficient, and the gas quickly drops to $10^4$ K to develop the second phase of the density–temperature diagram in equilibrium between cooling and photoionization heating from the UV background. This gas is spread over a wide range of densities. At higher densities, cooling becomes rapid again and is responsible for the third phase: dense ($30–300$ cm$^{-3}$) and cold gas ($10^2–10^3$ K) in which star formation can occur. For the rest of the paper, we categorize the gas into three phases based on temperature: cold ($T < 10^{3.5}$ K), warm ($10^{3.5}$ K <
$T < 10^{4.5}$ K), and hot ($T > 10^{4.5}$ K), corresponding to the three distinct peaks observed in the phase diagrams.

We find a striking difference between DC and the two other models, as gas in DC populates the intermediate- to high-density ($n_H \gtrsim 10^{-2}$ cm$^{-3}$) high-temperature ($T \gtrsim 10^5$ K) part of the diagram where cooling is supposed to be highly efficient. This is a direct consequence of the delayed cooling model, and the thermal state of these cells should be taken with caution. We have checked that the mass (volume) of gas in these cells represents less than 0.1 per cent of the total gas mass (volume) within $0.1 R_{200}$ at $z = 1$ and much less in the CGM. This hot dense gas is concentrated around stars and does not affect our results in the next sections. However, it might play a role in star formation itself by keeping gas from cooling and forming a new generation of stars.

Focusing now on the density distribution (the solid black lines at the bottom of each panel), we find that DC exhibits a notable difference relative to the other simulations. For both ME and MT, the density distribution has two peaks, roughly containing an equal mass of gas. These peaks respectively correspond to gas in the CGM (largely hot and very diffuse) and gas in the ISM (predominantly cold and relatively dense). This starkly contrasts with the results of DC, which has most of its gas in the low-density high-temperature CGM. In this simulation, the ISM contains much less gas than in the other two models, which is another manifestation of the very efficient (and ejective) feedback obtained with delayed cooling.

In the ISM, a large fraction of gas is expected to be in the warm phase (Ferrière 2001; Troland 2004). Recent high-resolution simulations focused on idealized galaxies (e.g. Katz et al. 2022) or the solar neighbourhood (e.g. Kim & Ostriker 2017) have observed comparable distribution. In contrast, the mass fraction of the warm phase within the ISM over $z = 1.3 - 1$ in ME and MT is $\approx 5$ per cent. The cold phase dominates the ISM gas with a mass fraction of $\approx 94$ per cent. Conversely, the cold and warm gas phases over that period, respectively, represent 49 per cent and 31 per cent of the ISM mass in DC. Although the mass distribution across the different phases remains stable over that period in ME and MT, it largely fluctuates in DC with the cold and warm gas mass fractions oscillating between 18–66 per cent and 20–54 per cent. The hot gas mass fraction in DC can reach up to 60 per cent, whereas it stays below 2 per cent in ME and MT throughout the entire period. This behaviour is anticipated and is a direct result of delayed cooling, which rapidly heats large regions of gas that do not cool down quickly. Despite the similar stellar masses of the galaxies in our three simulations, the varying subgrid models for feedback and star formation lead to distinct differences in the gas content of the galaxy and its surrounding halo.

### 3.3 Outflow properties

Having understood how the properties of the ISM and CGM differ in each simulation (i.e. stellar properties, temperature), we investigate the main culprit in driving these properties: galactic outflows. At a given radius $r$, we compute the instantaneous outflow rates using a method similar to Dalla Vecchia & Schaye (2008) and Nelson et al. (2019b), as

$$ \dot{M}(r) = 4\pi r^2 \frac{\sum_i^{N_{\text{shell}}} (m_i v_{r,i})}{\sum_i^{N_{\text{shell}}} \Delta x_i^3}, $$

with $4\pi r^2$ the surface of the spherical shell at a distance $r$ from the galaxy centre, and $v_{r,i} = v_i \cdot r_i / |r_i|$ the radial velocity of the gas, where the velocity is obtained by subtracting the mean mo-

tion of the dark matter halo. The two summations are done over $N_{\text{shell}}$ cells within a shell defined by $r_s \in [r \pm \Delta x_i / 2]$. $\Delta x_i$, $m_i$, $r_i$, and $v_i$, respectively, describe the cells' size, mass, position, and velocity in Cartesian coordinates. To minimize errors due to Cartesian grid structures, we compute the outflows using the factor $\dot{M}(r) = 4\pi r^2 / \sum_i^{N_{\text{shell}}} \Delta x_i^3$ instead of simply dividing the total momentum by the thickness of the shell with $1/\Delta x$. Additionally, we exclude all cells within $0.2 R_{\text{vir,sub}}$ of any subhaloes containing stars, with $R_{\text{vir,sub}}$ defined as the virial radius for subhaloes computed by the ADAPTAHOP halo finder (Tweed et al. 2009). Lastly, to avoid counting turbulent motions, we restrict data to cells with radial velocities $v_{r,i} > 0.5 v_{\text{circ}}$ for outflows ($v_{r,i} < -0.5 v_{\text{circ}}$ for inflows), with $v_{\text{circ}} = \sqrt{\frac{GM_{\text{halo}}}{R_{200}}}$, and $M_{\text{halo}}$ the total mass within the halo. Over $z = 1 - 1.3, 128 \text{ km s}^{-1} < v_{\text{circ}} < 139 \text{ km s}^{-1}$.

In Fig. 5, we show the outflow rates of the three simulations as a function of radius.

The three models produce total outflow rates within 1 dex of each other, with DC higher than ME, which is higher than MT. For all three models, the CGM mass flux increases with radius, a sign of outflows entraining CGM gas, as explained in Mitchell et al. (2020). For MT, ME, and DC, the outflow rates, respectively, range from 2, 3.5, and 8 $M_{\odot} \text{yr}^{-1}$ at $0.1 R_{200}$ to 7, 14, and 24 $M_{\odot} \text{yr}^{-1}$ at $R_{200}$, revealing a factor of 3–4 increase in flux, the highest increase being seen in ME and the smallest in DC. By the time the outflows reach the halo edge, at least 2/3 of the gas is entrained from the CGM.

We now focus on the different phases of the outflows and find that in all three simulations, the total outflow rates are widely dominated by hot gas. The entrained gas is hence mainly hot gas from the CGM.

In warm gas, we find contrasting behaviour between the three models. In the CGM, DC exhibits continuous and overall higher outflow rates than ME and MT, consistent with time. This warm outflowing gas is relatively conserved up to $\approx 25$ kpc, from whence its median value start decreasing. Conversely, ME and MT vary significantly over time (see the shaded area), with the median in ME dropping sharply in the CGM and MT showing a single peak of outflowing gas at $\approx 0.5 R_{200}$. This highlights the different behaviour of warm outflows in the three simulations. In DC, warm outflows are continuously driven in time, while in ME and MT, warm outflows are episodic, oscillating between states in which the galaxy is sending plumes of hot gas propagating outwards and having low warm outflow rates.

Lastly, there is no outflowing cold gas in any of our simulations.

Several recent studies suggest that hot outflows entrain a population of cold clouds from the ISM (Tan & Fielding 2024), which then grow (Gronke & Oh 2023) at a rate set by mixing (Tan, Oh & Gronke 2021) through cooling-induced pulsation (Gronke & Oh 2020) if they manage to survive. Our result possibly shows how cold gas fails to be entrained as it is heated up while being ejected out of the galaxy and thermally mixed with hotter gas. The current main explanation in the literature for this inability to entrain cold gas into the CGM is a lack of resolution. The resolution in our CGM is at kiloparsec scale, while Gronke & Oh (2018) and Gronke et al. (2022) suggest a need for resolution better than 100 pc, which corresponds roughly to the estimated size of CGM cold clouds (McCourt et al. 2018). Alternatively, non-thermal support provided by cosmic rays (Salem, Bryan & Corlies 2016; Ji et al. 2020) or magnetic fields (McCourt et al. 2015; Ramesh, Nelson & Pillepich 2023; Ramesh & Nelson 2024) could explain the lack of cold gas outflows and survival. However, cosmic ray properties are still poorly constrained (Ruszkowski, Yang & Zweibel 2017; Hopkins et al. 2020, 2021), and their role is greatly affected by the parameters used (Butsky et al.
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Spherical instantaneous outflow rates as a function of radius in the CGM of ME, MT, and DC. We show total outflow rates (top) and split it between hot ($T> 10^{4.5}$ K), warm ($10^{3.5} < T < 10^{4.5}$ K), and cold ($T< 10^{3.5}$ K) outflow phases (from top to bottom). The solid lines correspond to the median of each shell over 1 Gyr from $z = 1.3$ to $z = 1$, and the shaded area corresponds to the 15.9 and 84.1 percentiles. We split the figure between the ISM and the CGM by placing a vertical grey line at 0.1 $R_{200}$.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 6.** Spherical instantaneous inflow rates as a function of radius in the CGM of ME, MT, and DC. We show total inflow rates (top) and split it between hot ($T> 10^{4.5}$ K), warm ($10^{3.5} < T < 10^{4.5}$ K), and cold ($T< 10^{3.5}$ K) inflow phases (from top to bottom). The solid lines correspond to the median of each shell over 1 Gyr from $z = 1.3$ to $z = 1$, and the shaded area corresponds to the 15.9 and 84.1 percentiles. We split the figure between the ISM and the CGM by placing a vertical grey line at 0.1 $R_{200}$.</figcaption></figure>

2022; DeFelippi et al. 2024). The impact of magnetic fields on the CGM also remains debated (Li et al. 2020; Das & Gronke 2024).

In Fig. 6, we show the inflow rates of the three simulations. At first glance, we find once again total inflow rates within 1 dex of each other among the three simulations, with DC once again exhibiting higher fluxes than ME, which in turn are higher than MT. However, there are differences between the three models. The total inflow rate in MT decreases by a factor slightly higher than two from $R_{200}$ to 0.1 $R_{200}$, going from $7.5 M_{\odot} yr^{-1}$ to $3.5 M_{\odot} yr^{-1}$. In contrast, DC shows the smallest relative decrease in total inflow rates with

decreasing radius, going from $27 M_{\odot} yr^{-1}$ to $18.5 M_{\odot} yr^{-1}$, although the absolute decrease is the highest. This highlights how infalling gas in MT is more efficiently prevented from falling on to the galaxy, possibly stopped by outflows. By looking at the inflow temperature decomposition, we find that the hot inflowing phase fully drives this behaviour. This is expected as stellar feedback and gravitational shock heat up the halo gas (Birnboim & Dekel 2003; Dekel & Birnboim 2006; Ocvirk, Pichon & Teyssier 2008).

From $R_{200}$ to ≈ 5 kpc, the inflow rate in DC appears relatively constant. In reality, the inflow rates decrease down to ≈ 25 kpc, from
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Fraction of baryons in different media of the simulation for ME (left), MT (centre), and DC (right). The total baryonic mass considered is given by $f_b M_{200,\text{min}}$, with $f_b \approx 0.154$ and $M_{200,\text{min}}$ the dark matter mass within $R_{200,\text{min}}$. We split it into the stellar component (blue), the ISM (orange), the CGM (red), and the difference with the total (white). The difference between the baryonic mass of the simulation and that expected from the dark matter halo mass corresponds to the *missing baryons*. The result consists of 100 snapshots stacked over 1 Gyr from $z=1.3$ to $z=1$. For each simulation, we use the smallest $R_{200}$ among the three simulations $R_{200,\text{min}}$ to define the different regions at each time-step. The ISM contains all the gas within $0.1 R_{200,\text{min}}$, and the CGM corresponds to gas between $0.1 R_{200,\text{min}}$ and $R_{200,\text{min}}$.</figcaption></figure>

which they increase again to reach the same value as measured at $R_{200}$. It means that at large radii the gas is indeed being prevented from falling, but below $\approx 25$ kpc, there is added gas falling on the galaxy, likely in the form of recycled gas. By looking at the metallicity of the inflowing gas (not shown), we confirm that gas is being recycled as its metallicity increases with decreasing radius. By looking at the temperature of the inflows, we find that in DC, inflow rates are dominated by hot gas in the outer CGM and warm gas in the inner CGM. It thus seems that hot gas in DC is also being prevented from falling into the galaxy while warm gas is being recycled into the galaxy. This is consistent with the warm outflows observed at $\approx 25$ kpc, where warm gas is either slowed down or transitions to a hotter phase. In ME, the inflow rate increases from $R_{200}$ to $\approx 0.5 R_{200}$ possibly showing gas recycled at large radii, which is then conserved as falls into the galaxy. In MT, the warm inflowing rate decreases with decreasing radius down to $\approx 25$ kpc, where a sudden five-fold increase occurs. The warm inflow rate is then conserved at lower radii and exhibits a much lower time variability, possibly hinting at more consistent recycling flows. However, as the warm inflow rates of ME and MT are only a fraction of their respective total inflow rate, we cannot conclude what causes this behaviour. It might indeed show recycling, but it could also correspond to a fraction of the hot inflows being converted into warm gas, or more coherent inflows.

Lastly, the inflow rates of cold gas are negligible in all three simulations, showing that the cold gas that fuels star formation in the galaxy is either cold gas infalling at velocities below $0.5v_{\text{circ}}$, cold gas accreted through satellite mergers, or warm recycled gas cooling down once in the galaxy. By looking at inflows without removing the satellite contribution and without a velocity cut, we confirm that cold gas falls on to the galaxy from satellites, but there is no slow-moving infalling cold gas. We also compare both inflow and outflow rates of all gas phases without any velocity cut. As the velocity cut reduces the flow rates similarly in all three simulations, our conclusions on their comparative behaviour hold. However, while the outflow rate in DC slightly decreases (as all models do in the other gas phases), it extends to somewhat larger radii in ME. In MT, the outflow rates increase and become similar to ME, rapidly decreasing and extending to $\approx 25$ kpc.

This is consistent with the picture drawn from Fig. 6 in which warm outflowing gas slows down as it extends and is eventually recycled into the galaxy.

In summary, both outflows and inflows are highly dominated by hot gas, except for inflowing gas in DC, which is dominated by warm gas in the inner CGM. For all phases, DC exhibits the highest flux in the CGM, while it is lowest in MT. In all three simulations, hot outflows entrain hot gas from the CGM, while inflows are being prevented from infalling on to the galaxy. Warm outflowing gas in ME and MT exhibit very sporadic outflows. In DC, warm outflows are stable through time and decrease as a function of radius as they get recycled and fall back on to the galaxy.

## 3.4 The baryon distribution

In Fig. 7, we now investigate the location of the baryons *expected* within the halo, given by $f_b M_{200,\text{min}}$, with $f_b \approx 0.154$ the universal baryon fraction, and $M_{200,\text{min}}$ the dark matter mass within $R_{200,\text{min}}$. We define $R_{200,\text{min}}$ as the smallest $R_{200}$ among the three simulations. We stack outputs over the redshift range $z = 1-1.3$, and split $f_b M_{200,\text{min}}$ between stars, the ISM, the CGM, and the missing baryons (i.e. the difference between the expected amount of baryons and the sum of stars, ISM, and CGM masses).

This figure gives a visual impression of marked differences between the three simulations. Feedback processes reduce the baryon content to $\approx 56$ per cent of the universal fraction in the DC run, $\approx 79$ per cent in MT and 94 per cent in the ME simulation. ME and MT have comparable distribution of gas fractions in the halo components, with most of the baryons in stellar form, and a similar fraction of CGM and ISM gas. Conversely, DC has a much reduced stellar and ISM fraction and a CGM gas fraction higher than the other two simulations.

In Fig. 8, we look at the radial density distribution of gas in the three simulations. We weigh the y-axis by the radius squared to compare the relative mass fractions. In ME and MT, we find densities of $\approx 5 \text{ cm}^{-3}$ in the ISM (at 1 kpc), which then decrease sharply by four orders of magnitude, reaching densities of $\approx 6-10 \times 10^{-4} \text{ cm}^{-4}$ at $0.1 R_{200}$.
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Radial density profile of all the gas in the halo (top), as well as the hot, warm, and cold phases (from top to bottom) in ME (blue), MT (orange), and DC (green). The solid lines correspond to the median of each shell over 1 Gyr from $z = 1.3$ to $z = 1$, and the shaded area corresponds to the 15.9 and 84.1 percentiles. We also show with dashed lines the radial density profile of the last snapshot, at $z = 1$. Lastly, the figure is split between the ISM, the CGM and the IGM with two grey lines placed at $0.1 R_{200}$ and $R_{200}$.</figcaption></figure>

As expected from Fig. 7, the ISM gas densities are lower in DC, with no steep decrease at the galaxy edge, decreasing by slightly less than three dex over the same range, starting from ≈ 1 cm⁻³. The transition between the ISM and the CGM gas profile is more continuous in DC than in ME and MT and is caused by DC’s feedback removing gas efficiently from the ISM and ejecting it into the CGM and beyond, as observed in Fig. 7. In the CGM, the total gas density of the three models decreases similarly, roughly following an inverse-square law with radius.

Focusing on cold gas (lower panel of Fig. 8), we first find that DC exhibits much lower densities than both ME and MT, which both have comparable density profiles. However, all three simulations have a negligible amount of cold gas in the CGM as the amount of cold gas drops at $0.1 R_{200}$. The percentile distribution shown by the shaded area highlights short-lived cold gas in the form of satellites falling on to the main galaxy. We highlight their presence through the dashed lines, showing the density profile of a single snapshot. This is consistent with the absence of outflowing and inflowing cold gas observed in the CGM. As mentioned earlier, this possibly highlights how all models fail to entrain cold gas out into the CGM.

For warm gas, DC exhibits the highest densities up to ≈ 30 kpc. In the ISM, we find that the warm gas content in DC is quite significant and comparable to the cold gas content. In the CGM, the warm gas radial density also roughly evolves at a rate similar to the inverse square of the radius for ME and MT, but does so more steeply in DC. We have seen in Figs 5 and 6 that DC shows signs of warm gas recycling in the CGM with outflows and inflows both decreasing as a function of radius. This leads to a larger amount of hot gas in the inner CGM and a steeper decrease in the warm gas radial density profile with increasing radius.

Lastly, DC shows higher ISM densities in the hot phase than the other models. We have seen in Fig. 4 that gas in DC unrealistically populates the high-density ($n_H \gtrsim 10^{-2}$ cm⁻³), high-temperature ($T \gtrsim 10^5$ K) part of the halo phase diagram. This gas is close to star formation sites in which density is still high and where cooling has been turned off when supernova explosions are initiated. Thus, although this gas represents a small fraction of the total gas mass, it is expectedly dominant in radial density plots focusing on the hot gas phase. In the CGM, all three models show similar behaviour, with a decrease that scales inversely with radius. Outflows in all three models are dominated by hot gas, carrying it further than warm and cold gas. As a result, hot gas dominates the halo and displays the shallowest radial density decline among the three gas phases. However, note that while feedback can temporarily affect the pressure equilibrium of the halo, the radial density profile of the hot gas is constrained by the requirement that the total pressure gradient (bulk flows + turbulence + thermal) needs to balance gravity.

## 3.5 The metal budget

Having looked at the density properties of the halo, we now turn to metallicity. We show in Fig. 9 the stellar mass–metallicity relation over the redshift range $z = 1–1.3$. The stellar metallicity is mass-weighted, and computed within $0.1 R_{200}$. We find that ME and DC are both compatible with observations. As the galaxies have similar stellar masses, they should have produced a comparable amount of metals. However, this is not the case for the MT model, which uses a higher metal yield than ME and DC (0.1, compared to 0.075). This leads to a stellar metallicity much higher in MT than in ME, above observational constraints. There is however a factor of 4 difference between the stellar metallicity in ME and DC, which is not explained by a different yield factor. This difference is due to the feedback
<figure><img src="image_10.png" /><figcaption>**Figure 9.** Mass-weighted stellar metallicity relation for ME, MT, and DC, shown with markers. The lines associated with the markers show the metallicity evolution from $z = 1.3$ to $z = 1$. We compare our results to Gallazzi et al. (2014) at $0.65 \le z \le 0.75$ and Carnall et al. (2022) at $1 < z < 1.3$.</figcaption></figure>

models themselves, showing us that they regulate metals differently. DC is very efficient in removing metals from the galaxy and the halo, while ME and MT both retain most of the metals locked in stars.

For similar halo masses as ours, the EAGLE and the SIMBA simulations have reported that the mass of metals within stars in haloes over the total mass of metals in the same halo are respectively ≈ 60 per cent and ≈ 70 per cent (Appleby et al. 2021; Mitchell & Schaye 2022a). If we do the same calculation, we find that our simulations bracket theirs (≈ 80 per cent for ME and MT, ≈ 50 per cent for DC). However, it is more instructive to look at the mass of metals in the different regions of the halo over the mass of metals *released* by stars in the halo than only the metals within the halo.

We show this in the top panel of Fig. 10, with the location of all the metals produced by the galaxy. More specifically, we count the mass of metals locked up in stars $M_{Z,star}$, the mass of metals in the ISM $M_{Z,ISM}$ (i.e. with $r < 0.1 R_{200,\text{min}}$) and within the CGM $M_{Z,CGM}$ (i.e. with $0.1 R_{200,\text{min}} < r < R_{200,\text{min}}$). We normalize these quantities to the total mass of metals produced by the stars in the

simulated galaxy $M_{Z,\text{prod}}$, and compute the mass of metals that have been ejected out of the DM halo as $M_{Z,\text{out}} = M_{Z,\text{prod}} - M_{Z,\text{star}} - M_{Z,\text{ISM}} - M_{Z,\text{CGM}}$. We refer to them as *missing metals*. On average, the mass of metals produced by the galaxy over the redshift range $z = 1.3-1$ is $M_{Z,\text{prod,ME}} = 6.7 \times 10^8 M_{\odot}$ for ME, $M_{Z,\text{prod,MT}} = 7.6 \times 10^8 M_{\odot}$ for MT and $M_{Z,\text{prod,DC}} = 3.2 \times 10^8 M_{\odot}$ for DC.

Strikingly, the galaxies in ME and MT retain ≈ 90 per cent of the metals they produce within their stars or ISM, while the DC galaxy only keeps ≈ 25 per cent of its metal production within $0.1 R_{200,\text{min}}$. The very large fraction of ejected metals in the DC simulation makes it clear that the difference in the fraction of metals in stars or the ISM is due to DC very efficiently expelling metals from their production sites to large distances, beyond the virial radius, from where they do not return. This is possibly due to a collective enhancement of the supernovae (the supernovae in a stellar particle all occur at the same time in DC) and to the model producing stronger and hotter outflows (Rosdahl et al. 2017). ME and MT exhibit comparable behaviour regarding metal distribution. With ME’s feedback, metals that leave the ISM also mostly leave the DM halo for good, such that there is only a small fraction of metals present in the CGM. In MT, there is only a small fraction of metals in the CGM, but the feedback also ejects very little matter beyond the virial radius. Despite these differences, both ME and MT exhibit a similar total amount of metals in the CGM ($2.3 \text{ per cent} \times M_{Z,\text{prod,ME}} = 1.54 \times 10^7 M_{\odot}$ and $2.1 \text{ per cent} \times M_{Z,\text{prod,MT}} = 1.59 \times 10^7 M_{\odot}$). The CGM in DC exhibits both a high gas content and a high metallicity, resulting in a higher total amount of metal ($12.5 \text{ per cent} \times M_{Z,\text{prod,DC}} = 3.95 \times 10^7 M_{\odot}$). Despite a disparate fraction of baryons in the CGM, these results are within the scatter found in the EAGLE simulations. Our results are however inconsistent with IllustrisTNG, in which the CGM gas mass fraction is predominantly higher than 30 per cent for similar halo masses (Davies et al. 2020).

Observers often characterize the CGM through radial column densities, and all ions commonly used for CGM measurements depend on metallicity. Understanding the distribution of metals is thus crucial in understanding the CGM. Furthermore, we found in Fig. 10 that while the distributions are different, the metal mass also varies a lot between the three models due to different feedback modes, with DC having much more metals in the CGM than ME and MT. In

<figure><img src="image_11.png" /><figcaption>**Figure 10.** Fraction of metals in different media of the simulation for ME (left), MT (centre), and DC (right). The total mass of metals considered here is the mass of metals created through the simulation by the stars within the halo at each time-step. We split it into the stellar component (blue), the ISM (orange), the CGM (red), and the difference with the total (white). This difference corresponds to the fraction of metals ejected out of the halo. The ISM, CGM, and $R_{200,\text{min}}$ are defined as in Fig. 7. The result consists of 100 snapshots stacked over 1 Gyr from $z = 1.3$ to $z = 1$.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**Figure 11.** Mass of metals as a function of radius for each simulation. We bin the data radially and plot the median over 1 Gyr (from $z = 1.3$ to $z = 1$) in solid lines. The shaded areas correspond to the 15.9 and 84.1 percentiles. We also show with dashed lines the radial mass of metals in the last snapshot, at $z = 1$. Lastly, the figure is split between the ISM, the CGM and the IGM with two grey lines placed at 0.1 $R_{200}$ and $R_{200}$.</figcaption></figure>

Fig. 11, we finally look at the radial density distribution of metals in the CGM. We find that the metal mass in DC is much lower than in ME and MT in the galaxy, yet stands much higher than ME and MT in the CGM. Similarly to the total radial gas density observed in Fig. 8, we observe that the radial distribution of metals is more continuous in DC, while we notice a large drop at the galaxy edge in ME and MT. The absence of such a decline in DC shows how the galaxy in DC is unable to retain its metals in the ISM against the powerful feedback produced, unlike in ME and MT. Although ME and MT have a similar fraction of metals in the CGM (respectively, 2.3 per cent and 2.1 per cent, cf Fig. 10), ME exhibits higher metal densities in the CGM than MT. This is due to both the spherical binning and the temporal evolution of the CGM. The solid lines correspond to the binned median over all time-steps within $1.3 > z > 1$. We also show in dashed lines the metal distribution of the last output and find that metals are indeed especially concentrated around satellites in MT, while ME has an overall higher metal density in the CGM, yet lower metal densities around its satellites. As the satellites fall into the central galaxy, the metal distribution is smoothed in MT and brought to similar values as in ME. In the last snapshot, all three simulations are on the lower part of the shaded areas due to recent lower star formation rates and a thus lower supply of metals. The lines at higher redshift, e.g. at $z = 1.3$ (the first snapshot of our stack) are, conversely, higher than the shaded area in all three simulations.

## 3.6 Different feedback modes

Feedback can change the galaxy and its CGM properties in multiple ways, depending on how outflows affect the gas they interact with (Faucher-Giguère & Oh 2023). We mainly distinguish two feedback modes in the form of preventive and ejective feedback (Christensen et al. 2016; Mitchell & Schaye 2022b). While ejective feedback removes fuel for star formation by expelling gas out of the galaxy, preventive feedback impedes gas infall before it can reach the galaxy.

With an extreme case of preventive feedback, one could expect to find a metal-rich high-pressure CGM filled with ejected gas, keeping pristine gas at larger distances and shutting down accretion. With extreme ejective feedback, one could expect the opposite behaviour, with all the enriched gas ejected out of the halo while pristine gas can freely infall. In reality, it is likely that both regimes coexist in a complex manner, and the final state depends on how the gas

phases will mix, and on the scales over which this feedback occurs. For example, if preventive feedback were effective at the galactic scale, one could expect to find a significant amount of accumulated pristine gas in the CGM. On the other hand, if feedback prevented gas infall further than the halo scale, the CGM could effectively be devoid of gas. With these arguments in mind, it is possible to discuss feedback modes through which these models affect the properties of the galaxies and their CGM.

Simulations ME and DC have a higher fraction of missing metals than that of missing baryons in their haloes. This suggests that the *ejection* of matter outside of the haloes by galactic winds is the driving process that regulates the baryon fraction in these simulated haloes. On the contrary, MT shows a relatively large missing baryon fraction ($\approx 20$ per cent) but seems to retain most of its metals within $R_{200}$ ($\approx 94$ per cent). This suggests that outflows are not efficient at pushing gas outside of the halo in this simulation, but somehow manage to slow down accretion of pristine material (see Figs 5 and 6). There is an additional complexity due to the fraction of metals locked up in stars. ME and MT both retain a large fraction of metals in stars. This suggests that metals released by stars are retained long enough in the ISM to form new stars because supernovae feedback is not strong enough to lower densities in the central region, allowing metals to mix in the ISM. On the contrary, DC keeps only a small fraction of its metals in its stars due to strong feedback that ejects a large amount of metals from the galaxy. This reminds us that accurately modelling outflows in simulations is difficult not only because of the numerically demanding nature of physical processes that may happen in the CGM (Faucher-Giguère & Oh 2023), but also because details of metal mixing at the launch sites of galactic winds determine the metallicity profile of outflows (Mac Low & Ferrara 1999).

As ME and DC share the same star formation models, the difference in metal ejection is solely driven by the feedback model. Delayed cooling is known to produce stronger outflows (Rosdahl et al. 2017), but it also forbids immediate star formation by keeping the surrounding media at high effective pressure. This leads to burstier star formation and an overall more powerful feedback than in both other models. This feedback is hence more efficient in driving far-reaching outflows, which directly correlates with a higher efficiency in driving metals out of the halo. While the star formation and feedback models in ME and MT are relatively similar, they none the less produce distinct feedback modes. By resimulating a star formation peak in similar simulations, Rey (2022) showed that they form stars in contrasting ways, either driven by turbulence or gravitation. Ensuing supernovae will thus explode in disparate media and interact with the ISM gas differently. Also, while ME only injects feedback energy as momentum in both resolved and unresolved regimes, MT also injects thermal energy in all cases, assuming it will dissipate due to the lack of resolution if momentum is also injected. With a different threshold to distinguish the two regimes and even different equations for the momentum injected, these additional distinct features in both models can collectively drive contrasting feedback modes. The most appropriate approach to conclusively determine how the feedback modes are driven would be to perform simulations of the three models in an idealized environment; however, this analysis lies beyond the scope of this study.

## 4 CONCLUSIONS

In this work, we run cosmological zoom-in simulations of the same galaxy with three different subgrid models for star formation and feedback using the radiation hydrodynamical code RAMSES-RT. We
refer to our three models as ME (Kimm et al. 2017), MT (Kretschmer & Teyssier 2020), and a variation of ME with supernova feedback modelled through delayed cooling (Teyssier et al. 2013), DC. We have seen that three simulations with different feedback models can be calibrated to produce galaxies with similar stellar masses. Despite these similar stellar properties, we find that feedback from supernovae operates very differently in each model, and leads to varying CGM properties. Our results can be summarized as follows.

(i) **Three different star formation and stellar feedback subgrid models can be calibrated in stellar mass.** The three models ME, MT, and DC, often used in the literature, are employed on the same set of initial conditions and are calibrated to produce a similar stellar mass within a factor of ~2.

(ii) **The burstiness is sensitive to subgrid models.** Due to a more powerful feedback, DC exhibits much higher burstiness than both ME and MT at all time-scales considered. ME exhibits similar burstiness on both short and long time-scales, while MT is burstier than ME and shows a higher burstiness when computing it on larger time-scales, possibly due to the implementation of subgrid turbulence in MT and the different properties of the star-forming media.

(iii) **The halo gas of the three simulations shows similar properties, but disparate distributions.** The phase diagram of the three simulations is comparable with three phases: the cold dense ISM phase, the hot and low-density ejected or shock-heated phase, and the warm phase plateauing at $T \sim 10^4$ K. However, while ME and MT have a similar amount of gas in the ISM and in the CGM, DC has much more gas in the CGM, and only a small fraction of cold gas in the ISM as of its powerful feedback strips most of the gas from the ISM and the model itself heats a fraction of the remaining gas by construction.

(iv) **Inflows and outflows are dominated by hot gas in all simulations, but warm gas shows contrasting behaviours.** In all three simulations, both inflows and outflows are dominated by hot gas, which entrains gas from the CGM and reduces the inflow rate of hot gas. However, in the inner CGM, the total inflow rate is dominated by warm gas in DC. The flows in DC exhibit signs of warm gas recycling, with both inflows and outflows decreasing as a function of radius. In ME and MT, warm outflows are sporadic in time but have otherwise comparable rates to DC.

(v) **Feedback leads to drastically different baryon distributions in the halo.** Through feedback much more efficient in removing gas from the halo, the baryonic content of the halo in DC is reduced $\approx 56$ per cent of the universal fraction, and to $\approx 79$ per cent and $\approx 94$ per cent in MT and ME. The distribution of this gas is also different among the three simulations, DC's ISM and stellar mass each constitute less than half the fractions found with ME and MT. The CGM gas fractions are also distinct, with DC exhibiting the highest fraction and MT the lowest.

(vi) **Feedback leads to drastically different metal distributions in the halo.** The metal mass in the ISM is much lower in DC than in ME and MT, while it is much higher in the CGM. DC feedback manages to eject 65 per cent of the metal mass out of the halo while ME and MT only eject 5–10 per cent. None the less, the galaxy's metallicity is compatible with observations in both DC and ME. Also, while ME and MT exhibit a similar total metal mass in the CGM, the radial mass of metals is higher in ME as metals in MT are concentrated around satellites.

(vii) **Diverse feedback modes dominate the three simulations.** DC is largely dominated by ejective feedback, leading to a metal- and gas-rich CGM and a high fraction of metals (64 per cent) and baryons (44 per cent) ejected or kept out of the halo. Conversely,

ME and MT are dominated by preventive feedback, keeping most metals in stars (72–75 per cent) and leading to a metal-poor CGM (2.3–2.1 per cent). MT keeps pristine gas beyond the halo, reaching $\approx 80$ per cent of the universal baryon fraction inside the halo. ME keeps pristine gas beyond the galaxy and thus exhibits a CGM filled with more gas, and baryon fraction close to the universal value.

We have seen in this paper that feedback modes can regulate galaxy growth through distinct means, and produce galaxies with similar stellar masses. None the less, we also found that even when producing comparable galaxies, the density, metallicity, and temperature of the CGM can greatly differ. Understanding the CGM is thus crucial to trace back and understand what dictates galaxy growth and the role of different feedback processes. As most simulations rely on distinct feedback models, it is timely and necessary to find new ways to distinguish subgrid models through new constraints extending beyond the properties of the galaxies. One such constraint exists in the form of quasar absorption lines. In this series of papers, we will investigate how our state-of-the-art simulations compare against these constraints, and examine whether CGM absorption lines can help distinguish subgrid models in the near future.

## ACKNOWLEDGEMENTS

We thank the anonymous referee for their suggestions. We gratefully acknowledge support from the CBPsmn (PSMN, Pôle Scientifique de Modélisation Numérique) of the ENS de Lyon for the computing resources. The platform operates the SIDUS solution (Quemener & Corvellec 2013) developed by Emmanuel Quemener. MR and TK were supported by the National Research Foundation of Korea (NRF) grant funded by the Korean government (RS-2022-NR070872 and RS-2025-00516961). TK is also supported by the Yonsei Fellowship, funded by Lee Youn Jae. The supercomputing time for part of the numerical tests was kindly provided by KISTI (KSC-2024-CRE-0200), and large data transfer was supported by KREONET, which is managed and operated by KISTI. This research has made use of the Astrophysics Data System, funded by NASA under Cooperative Agreement 80NSSC21M00561. We also acknowledge the use of PYTHON (Van Rossum & Drake Jr 1995), MATPLOTLIB (Hunter 2007), NUMPY (Harris et al. 2020), and ASTROPY (Astropy Collaboration 2013, 2018) for this work.

## DATA AVAILABILITY

The data generated and used in this article will be shared upon reasonable request to the corresponding author.

## REFERENCES

* Afruni A., Fraternali F., Pezzulli G., 2021, MNRAS, 501, 5575

* Afruni A. et al., 2023, A&A, 680, A112

* Agertz O., Teyssier R., Moore B., 2011, MNRAS, 410, 1391

* Agertz O., Kravtsov A. V., Leitner S. N., Gnedin N. Y., 2013, ApJ, 770, 25

* Agertz O. et al., 2020, MNRAS, 491, 1656

* Anand A., Nelson D., Kauffmann G., 2021, MNRAS, 504, 65

* Andersson E. P., Agertz O., Renaud F., Teyssier R., 2023, MNRAS, 521, 2196

* Appleby S., Davé R., Sorini D., Storey-Fisher K., Smith B., 2021, MNRAS, 507, 2383

* Astropy Collaboration, 2013, A&A, 558, A33

* Astropy Collaboration, 2018, AJ, 156, 123

* Azartash-Namin B. et al., 2024, ApJ, 970, 40

* Behroozi P., Wechsler R. H., Hearin A. P., Conroy C., 2019, MNRAS, 488, 3143
Bhagwat A., Costa T., Ciardi B., Pakmor R., Garaldi E., 2024, MNRAS, 531, 3406

Birnboim Y., Dekel A., 2003, MNRAS, 345, 349

Bordoloi R., Prochaska J. X., Tumlinson J., Werk J. K., Tripp T. M., Burchett J. N., 2018, ApJ, 864, 132

Butsky I. S. et al., 2022, ApJ, 935, 69

Carnall A. C. et al., 2022, ApJ, 929, 131

Ceverino D., Klypin A., 2009, ApJ, 695, 292

Chaikin E., Schaye J., Schaller M., Bahé Y. M., Nobels F. S. J., Ploeckinger S., 2022, MNRAS, 514, 249

Chen H.-W., Zahedy F. S., 2024, The Circumgalactic Medium

Cho H., Kang H., 2008, New Astron., 13, 163

Christensen C. R., Davé R., Governato F., Pontzen A., Brooks A., Munshi F., Quinn T., Wadsley J., 2016, ApJ, 824, 57

Cook A. W. S., van de Voort F., Pakmor R., Grand R. J. J., 2024, The Halo Mass Dependence of Physical and Observable Properties in the Circumgalactic Medium

Crain R. A., McCarthy I. G., Schaye J., Theuns T., Frenk C. S., 2013, MNRAS, 432, 3005

Crain R. A. et al., 2015, MNRAS, 450, 1937

Dale J. E., 2015, New Astron. Rev., 68, 1

Dalla Vecchia C., Schaye J., 2008, MNRAS, 387, 1431

Das H. K., Gronke M., 2024, MNRAS, 527, 991

Davé R., Anglés-Alcázar D., Narayanan D., Li Q., Rafieferantsoa M. H., Appleby S., 2019, MNRAS, 486, 2827

Davies J. J., Crain R. A., Oppenheimer B. D., Schaye J., 2020, MNRAS, 491, 4462

DeFelippis D., Bournaud F., Bouché N., Tollet E., Farcy M., Rey M., Rosdahl J., Blaizot J., 2024, MNRAS, 530, 52

Decataldo D., Shen S., Mayer L., Baumschlager B., Madau P., 2024, A&A, 685, A8

Dekel A., Birnboim Y., 2006, MNRAS, 368, 2

Dubois Y. et al., 2014, MNRAS, 444, 1453

Dubois Y. et al., 2021, A&A, 651, A109

Dutta R. et al., 2021, MNRAS, 508, 4573

Fattahi A. et al., 2016, MNRAS, 457, 844

Faucher-Giguère C.-A., Oh S. P., 2023, ARA&A, 61, 131

Faucher-Giguère C.-A., Lidz A., Zaldarriaga M., Hernquist L., 2009, ApJ, 703, 1416

Federrath C., Klessen R. S., 2012, ApJ, 761, 156

Ferland G. J., Korista K. T., Verner D. A., Ferguson J. W., Kingdon J. B., Verner E. M., 1998, PASP, 110, 761

Ferrière K. M., 2001, Rev. Mod. Phys., 73, 1031

Fumagalli M., 2024, The Multiphase Circumgalactic Medium and Its Relation to Galaxies: An Observational Perspective

Gabor J. M., Capelo P. R., Volonteri M., Bournaud F., Bellovary J., Governato F., Quinn T., 2016, A&A, 592, A62

Gallazzi A., Bell E. F., Zibetti S., Brinchmann J., Kelson D. D., 2014, ApJ, 788, 72

Geen S., Rosdahl J., Blaizot J., Devriendt J., Slyz A., 2015, MNRAS, 448, 3248

Genel S. et al., 2019, ApJ, 871, 21

Gnedin N. Y., Abel T., 2001, New Astron., 6, 437

Goh K.-I., Barabasi A.-L., 2006, EPL (Europhysics Letters)81, 48002

Governato F. et al., 2010, Nature, 463, 203

Grand R. J. J. et al., 2017, MNRAS, 467, 179

Gronke M., Oh S. P., 2018, MNRAS, 480, L111

Gronke M., Oh S. P., 2020, MNRAS, 492, 1970

Gronke M., Oh S. P., 2023, MNRAS, 524, 498

Gronke M., Oh S. P., Ji S., Norman C., 2022, MNRAS, 511, 859

Guillet T., Teyssier R., 2011, J. Comput. Phys., 230, 4756

Haardt F., Madau P., 1996, ApJ, 461, 20

Hafen Z. et al., 2019, MNRAS, 488, 1248

Hahn O., Abel T., 2011, MNRAS, 415, 2101

Hahn O., Abel T., 2013, Astrophysics Source Code Library, record ascl:1311.011

Harris C. R. et al., 2020, Nature, 585, 357

Hennebelle P., Chabrier G., 2011, ApJ, 743, L29

Hopkins P. F., Narayanan D., Murray N., 2013, MNRAS, 432, 2647

Hopkins P. F., Kereš D., Oñorbe J., Faucher-Giguère C.-A., Quataert E., Murray N., Bullock J. S., 2014, MNRAS, 445, 581

Hopkins P. F. et al., 2018, MNRAS, 480, 800

Hopkins P. F. et al., 2020, MNRAS, 492, 3465

Hopkins P. F., Chan T. K., Squire J., Quataert E., Ji S., Kereš D., Faucher-Giguère C.-A., 2021, MNRAS, 501, 3663

Hopkins P. F. et al., 2023, MNRAS, 519, 3154

Hunter J. D., 2007, Comput. Sci. Eng., 9, 90

Ji S. et al., 2020, MNRAS, 496, 4221

Kang C., Kimm T., Han D., Katz H., Devriendt J., Slyz A., Teyssier R., 2025, A&A, 693, A149

Katz N., 1992, ApJ, 391, 502

Katz H. et al., 2021, MNRAS, 507, 1254

Katz H. et al., 2022, PRISM: A Non-Equilibrium, Multiphase Interstellar Medium Model for Radiation Hydrodynamics Simulations of Galaxies

Keller B. W., Wadsley J., Benincasa S. M., Couchman H. M. P., 2014, MNRAS, 442, 3013

Keller B. W., Wadsley J. W., Wang L., Kruijssen J. M. D., 2019, MNRAS, 482, 2244

Kelly A. J., Jenkins A., Deason A., Fattahi A., Grand R. J. J., Pakmor R., Springel V., Frenk C. S., 2022, MNRAS, 514, 3113

Kim C.-G., Ostriker E. C., 2015, ApJ, 802, 99

Kim C.-G., Ostriker E. C., 2017, ApJ, 846, 133

Kimm T., Cen R., 2014, ApJ, 788, 121

Kimm T., Cen R., Devriendt J., Dubois Y., Slyz A., 2015, MNRAS, 451, 2900

Kimm T., Katz H., Haehnelt M., Rosdahl J., Devriendt J., Slyz A., 2017, MNRAS, 466, 4826

Kretschmer M., Teyssier R., 2020, MNRAS, 492, 1385

Kretschmer M., Dekel A., Teyssier R., 2022, MNRAS, 510, 3266

Krumholz M. R., McKee C. F., 2005, ApJ, 630, 250

Li H., Gnedin O. Y., Gnedin N. Y., 2018, ApJ, 861, 107

Li Z., Hopkins P. F., Squire J., Hummels C., 2020, MNRAS, 492, 1841

Mac Low M.-M., Ferrara A., 1999, ApJ, 513, 142

Martizzi D., Faucher-Giguère C.-A., Quataert E., 2015, MNRAS, 450, 504

McCourt M., O'Leary R. M., Madigan A.-M., Quataert E., 2015, MNRAS, 449, 2

McCourt M., Oh S. P., O'Leary R., Madigan A.-M., 2018, MNRAS, 473, 5407

Medlock I. et al., 2025, ApJ, 980, 61

Mitchell P. D., Schaye J., 2022a, MNRAS, 511, 2600

Mitchell P. D., Schaye J., 2022b, MNRAS, 511, 2948

Mitchell P. D., Blaizot J., Devriendt J., Kimm T., Michel-Dansac L., Rosdahl J., Slyz A., 2018, MNRAS, 474, 4279

Mitchell P. D., Schaye J., Bower R. G., Crain R. A., 2020, MNRAS, 494, 3971

Moster B. P., Naab T., Lindström M., O'Leary J. A., 2021, MNRAS, 507, 2115

Munshi F. et al., 2013, ApJ, 766, 56

Naab T., Ostriker J. P., 2017, ARA&A, 55, 59

Nelson D. et al., 2019a, Comput. Astrophys. Cosmol., 6, 2

Nelson D. et al., 2019b, MNRAS, 490, 3234

Nelson D. et al., 2020, MNRAS, 498, 2391

Nuñez-Castiñeyra A., Nezri E., Devriendt J., Teyssier R., 2021, MNRAS, 501, 62

Ocvirk P., Pichon C., Teyssier R., 2008, MNRAS, 390, 1326

Oppenheimer B. D. et al., 2020, ApJ, 893, L24

Padoan P., Nordlund Å., 2011, ApJ, 730, 40

Peroux C., Nelson D., 2024, The Multi-Scale Multi-Phase Circumgalactic Medium: Observed and Simulated

Perret V., Renaud F., Epinat B., Amram P., Bournaud F., Contini T., Teyssier R., Lambert J. C., 2014, A&A, 562, A1

Pillepich A. et al., 2019, MNRAS, 490, 3196

Planck Collaboration XVI, 2014, A&A, 571, A16

Quemener E., Corvellec M., 2013, Linux J., 2013, 3:3

Ramesh R., Nelson D., 2024, MNRAS, 528, 3320

Ramesh R., Nelson D., Pillepich A., 2023, MNRAS, 522, 1535

Rasera Y., Teyssier R., 2006, A&A, 445, 1
Reina-Campos M., Keller B. W., Kruijssen J. M. D., Gensior J., Trujillo-Gomez S., Jeffreson S. M. R., Pfeffer J. L., Sills A., 2022, MNRAS, 517, 3144

Rey M., 2022, PhD thesis, École normale supérieure de Lyon

Rosdahl J., Blaizot J., 2012, MNRAS, 423, 344

Rosdahl J., Teyssier R., 2015, MNRAS, 449, 4380

Rosdahl J., Blaizot J., Aubert D., Stranex T., Teyssier R., 2013, MNRAS, 436, 2188

Rosdahl J., Schaye J., Dubois Y., Kimm T., Teyssier R., 2017, MNRAS, 466, 11

Rosdahl J. et al., 2018, MNRAS, 479, 994

Rosdahl J. et al., 2022, MNRAS, 515, 2386

Rosen A., Bregman J. N., 1995, ApJ, 440, 634

Ruszkowski M., Yang H. Y. K., Zweibel E., 2017, ApJ, 834, 208

Saeedzadeh V. et al., 2023, MNRAS, 525, 5677

Salem M., Bryan G. L., Corlies L., 2016, MNRAS, 456, 582

Sameer et al., 2024, MNRAS, 530, 3827

Sawala T. et al., 2016, MNRAS, 457, 1931

Schaye J., 2001, ApJ, 562, L95

Schaye J. et al., 2015, MNRAS, 446, 521

Schmidt W., Niemeyer J. C., Hillebrandt W., 2006, A&A, 450, 265

Semenov V. A., Kravtsov A. V., Gnedin N. Y., 2016, ApJ, 826, 200

Silk J., Rees M. J., 1998, A&A, 331, L1

Simpson C. M., Bryan G. L., Hummels C., Ostriker J. P., 2015, ApJ, 809, 69

Smith M. C., 2021, MNRAS, 502, 5417

Smith M. C., Sijacki D., Shen S., 2018, MNRAS, 478, 302

Somerville R. S., Davé R., 2015, ARA&A, 53, 51

Stanway E. R., Eldridge J. J., 2018, MNRAS, 479, 75

Stinson G., Seth A., Katz N., Wadsley J., Governato F., Quinn T., 2006, MNRAS, 373, 1074

Suresh J., Nelson D., Genel S., Rubin K. H. R., Hernquist L., 2019, MNRAS, 483, 4040

Tan B., Fielding D. B., 2024, MNRAS, 527, 9683

Tan B., Oh S. P., Gronke M., 2021, MNRAS, 502, 3179

Teyssier R., 2002, A&A, 385, 337

Teyssier R., 2010, Astrophysics Source Code Library, record ascl:1011.007

Teyssier R., Pontzen A., Dubois Y., Read J. I., 2013, MNRAS, 429, 3068

Thornton K., Gaudlitz M., Janka H. T., Steinmetz M., 1998, ApJ, 500, 95

Trebitsch M. et al., 2021, A&A, 653, A154

Troland T. H., 2004, in The Magnetized Interstellar Medium, Eds: B. Uyaniker W. Reich and R. Wielebinski, Copernicus GmbH, Katlenburg-Lindau. p. 105

Tumlinson J., Peebles M. S., Werk J. K., 2017, ARA&A, 55, 389

Tweed D., Devriendt J., Blaizot J., Colombi S., Slyz A., 2009, A&A, 506, 647

Van Rossum G., Drake F. L., Jr, 1995, Python Reference Manual. Centrum voor Wiskunde en Informatica Amsterdam, Amsterdam, The Netherlands.

Velliscig M., van Daalen M. P., Schaye J., McCarthy I. G., Cacciato M., Le Brun A. M. C., Dalla Vecchia C., 2014, MNRAS, 442, 2641

Wang L., Dutton A. A., Stinson G. S., Macciò A. V., Penzo C., Kang X., Keller B. W., Wadsley J., 2015, MNRAS, 454, 83

Yi S. K. et al., 2024, ApJS, 271, 1

This paper has been typeset from a *TEX/LATEX* file prepared by the author.