efficient cooling and gravitational collapse leads to the initiation of nuclear fusion and the consequent birth of new stars (e.g. Kennicutt 1983 ;McKee & Ostriker 2007 ;Kennicutt & Evans 2012 ). A fraction of the gas that makes it into galaxy nuclei can also be accreted onto supermassive black holes (SMBHs), which are known to inhabit the dense central regions of many galaxies (see re vie ws by Ferrarese & Ford 2005 ;Kormendy & Ho 2013 ;Heckman & Best 2014 ).

Both star formation and accretion onto SMBH drive large-scale outflows, which play a crucial role in regulating the baryon cycle (Veilleux, Cecil & Bland-Hawthorn 2005 ;Veilleux et al. 2020 ). Stellar feedback, through supernovae explosions, radiation, and stellar winds, injects energy, momentum, and enriched material into the ISM and potentially beyond (Heckman et al. 2015 ). Active galactic nucleus (AGN) feedback, powered by accretion onto SMBHs at galactic centres, also releases substantial amounts of energy and momentum, that can drive powerful outflows, impacting the surrounding gas and quenching star formation in some cases (e.g. Fabian 2012 ;Heckman & Best 2014 ;King & Pounds 2015 ).

On the scale of the ISM, various observational techniques provide relatively well-constrained measurements of the stellar mass, galaxy MNRAS 532, 3417-3440 (2024) cold gas (H īI and H 2 ) content, and star formation rates (SFRs). Beyond the ISM, direct observational constraints on the CGM and IGM gas are associated with increasing uncertainty due to its diffuse nature (e.g. Werk et al. 2014 ;Tumlinson et al. 2017 ). Outside the Milky Way, gas accretion has been detected in a number of cases where absorption signatures are identified in galaxy spectra (e.g. Rubin et al. 2012 ;Martin et al. 2012 ;Stone et al. 2016 ; for a re vie w see Rubin 2017 ); and in a number of studies where it is possible to use a fortuitously aligned background quasar to probe circumg alactic g as flows (e.g. Lehner et al. 2013 ;Bouch é et al. 2016 ;Martin et al. 2019 ). Even with these disco v eries, the rate at which the gas accretes and its associated properties have yet to be wellquantified on a statistical basis; and are associated with considerable uncertainties.

There is a significant body of literature presenting supporting evidence for the presence of feedback-driven gas outflows from galaxies (e.g. Heckman et al. 2000Heckman et al. , 2015 ; ;Martin 2005 ;Veilleux et al. 2005 ;Feruglio et al. 2010 ;Newman et al. 2012 ;Rubin et al. 2014 ;Schroetter et al. 2016 ;Chisholm et al. 2016 ;Rupke et al. 2019 ; for a recent re vie w see Veilleux et al. 2020 ). Ho we ver, in each case, the determination of associated mass flux is affected by several systematic uncertainties. In particular, a given outflow tracer only provides information about gas on a subset of spatial scales (nominally confined to well within a given halo virial radius), and to specific gas phases. Furthermore, these measurements are biased towards g as-rich g alaxies exhibiting particularly strong outflows -either those at high redshift, or local galaxies experiencing a starburst.

Modern cosmological hydrodynamical simulations constitute a powerful tool to explore the role of different processes within the baryon cycle, and galaxy evolution more broadly (e.g. Schaye et al. 2010Schaye et al. , 2015 ; ;Dubois et al. 2014 ;Vogelsberger et al. 2014 ;Christensen et al. 2016 ;Pillepich et al. 2018 ;Hopkins et al. 2018 ;Dav é et al. 2019 ). Ho we ver, in order to accurately simulate the formation of galaxy populations within the context of large-scale structure, it is crucial to model essential small-scale processes taking place within individual galaxies. These processes, such as star formation, black hole formation and accretion, and stellar and AGN feedback, operate at scales below the resolution capabilities of the simulations. Therefore, they are commonly incorporated using 'subgrid' techniques, which ef fecti vely account for their influence on galaxy evolution in a physically motivated, but necessarily coarse-grained manner (for re vie ws, see Somerville & Dav é 2015 ;Vogelsberger et al. 2020 ). A wide range of subgrid techniques have demonstrated success in reproducing realistic galaxy populations in cosmological simulations, particularly in capturing the observed shape of the stellar mass function (e.g. Thorne et al. 2021 ;Schaye et al. 2023 , see also Section 3 ).

Despite these promising results, recent studies have suggested that the behaviour of the gaseous phase and the broader-scale baryon cycle differ significantly between different simulations and feedback models, and can be degenerate in producing the same stellar mass assembly in galaxies o v er cosmic time (Naab & Ostriker 2017 ;Mitchell et al. 2018Mitchell et al. , 2020 ; ;Pandya et al. 2020 ;Dav é et al. 2020 ;Villaescusa-Navarro et al. 2021 ;Ni et al. 2023 ;Tillman et al. 2023a ). F or e xample, Dav é et al. ( 2020 ) compare the cold gas content of galaxies in three publicly available modern cosmological hydrodynamical simulations -EAGLE, TNG100, and SIMBA (for a technical o v erview of the subgrid prescriptions in each simulation, see Section 2 ). They find that these three simulations produce relatively similar predictions for the molecular gas content of galaxies (H 2 mass function or H 2 MF) o v er cosmic time, within a range of 0 . 3 dex. This phase is strongly correlated with star formation, which also agrees fairly well between models as also established by the aforementioned stellar mass function. Predictions for the slightly more tenuous atomic phase, ho we ver, sho w greater di vergence.

Similarly, the emergent gas mass outflow rates and mass loading factors ( η ≡ Ṁ out ( r) / Ṁ ) at ISM scales that arise from different subgrid implementations of stellar-driven winds appear to differ significantly even across simulations that are all calibrated to reproduce the observed z = 0 stellar mass function. Mitchell et al. ( 2020 ) show a comparison of the ISM-scale mass loadings measured in the EAGLE simulation and those reported by Nelson et al. ( 2019 ) from the TNG50 run of the IllustrisTNG simulation suite. Their results suggest that the gas outflow rates in TNG decline rapidly with increasing radius, such that most of the outflowing ISM material does not make it out of the CGM, while in EAGLE outflows are actually powerful enough to leave the CGM and halo entirelyeven entraining some gas on the way. However, we stress that the comparisons between these simulations were not conducted with the same methodology - Nelson et al. ( 2019 ) calculate Eulerian instantaneous flow rates, while Mitchell et al. ( 2020 ) instead use a Lagrangian particle tracking method. In addition, flow quantities have not necessarily been measured at the same physical scales or using the same definitions of halo properties.

To date, a comprehensive examination of the distinct gas cycling paradigms within these simulations has not been carried out. In this paper, we aim to fill this gap by conducting a detailed comparison of the baryon cycle in the EAGLE, TNG, and SIMBA simulations using the same definitions and methodology. Our primary focus is to gain insight into the physical extent of feedback-dri ven outflo ws across a broad range of halo masses. By undertaking such an analysis, we lay the foundation for developing targeted observational tests that can help to assess specific simulation methodologies, and mo v e towards a more constrained understanding of the role of feedback processes in the baryon cycle.

This paper is arranged as follows: in Section 2 , we outline each of the simulations -EAGLE, TNG, and SIMBA -that we analyse in this study, and the uniform methodology we employ in each to measure gas flow rates. In Section 3 , we outline the similarities and differences between statically measured baryon content and baryon distribution in haloes in these simulations. In Section 4 , we present a detailed analysis of the gas flow rates at different scales surrounding galaxies in these simulations, and how these differences can be used to explain the findings in Section 3 . In Section 5 , we discuss the implications of our findings for the baryon cycle, and discuss future targeted observational tests for different feedback models. Finally, we conclude with Section 6 , where we summarize the main findings of the paper.

In this paper, we make use of three modern cosmological hydrodynamical simulations with publicly available data -namely EAGLE, TNG, and SIMBA. These simulations are all based on a CDM (lambda-cold dark matter) cosmology -in which hierarchically assembled CDM haloes provide the formation site of galaxies; all within the context of an expanding, dark energy dominated universe. As noted abo v e, all of these simulations contain phenomenological subgrid recipes, which are discussed in more detail below, containing parameters that are calibrated via comparison to observational data or quantities derived from observations. These calibrations are anchored using observations of the stellar mass and/or distribution of galaxies, of the models (see Crain & van de Voort 2023 for an o v erview). 1 In Section 2.1 , we introduce the simulations we employ for this study and describe their respective subgrid feedback implementations. Subsequently, in Section 2.2 , we outline the method with which we self-consistently measure gas flows rates in each of the simulations.

In Table 1 , we compare some key features of the simulations used in this paper. We then discuss these features in further detail in Sections 2. 1.1 , 2.1.2 , and 2.1.3 for the EAGLE, TNG, and SIMBA simulations respectively. Each of these simulations share a base CDM cosmology, with the choice of cosmological parameters differing by only of order ≈ 1 per cent . While the simulations take slightly different approaches to solving the equations describing the hydrodynamics of the gas, the y hav e man y characteristics in common: each has a spatial resolution of order 1 -10 kpc , make use of very similar tabulated gas cooling prescriptions, and utilize an ef fecti ve equation of state (EoS) for cool gas within the ISM. The critical difference between these simulations, for the purposes of this study, is in their approach to modelling the unresolved feedback from star formation (stellar feedback) and SMBHs (AGN feedback), which we focus on below. 1 An exception to this is the TNG AGN feedback model, which is calibrated to reproduce z ≈ 0 halo gas fractions at high mass, and SFR densities o v er cosmic time. Even so, there is considerable degeneracy in how the exact interplay between gas flows and the baryon cycle can produce these results.

The EAGLE simulations, as described by Schaye et al. ( 2015 ) and Crain et al. ( 2015 ), utilized a modified version of the Nbody Tree-ParticleMesh ( TREEPM ) smoothed particle hydrodynamics (SPH) solver called GADGET -3, initially presented in Springel ( 2005 ). The set of modifications ( ANARCHY ) are outlined in Schaller et al. ( 2015 ); and include enhancements to the hydrodynamics solver and subgrid physics modules. The initial conditions for this simulation were generated using the method described in Jenkins ( 2013 ), assuming the cosmological parameters from Planck Collaboration XVI ( 2014 ): = 0 . 693, m = 0 . 307, b = 0 . 04825, H 0 = 67 . 77 km s -1 Mpc -1 , σ 8 = 0 . 8288, and n s = 0 . 9611.

The flagship box of the suite, referred to as Ref-L100N1504, has a periodic side length of 67 . 8 h -1 Mpc . It was e x ecuted with 1 , 504 3 DM and baryonic particles, resulting in particle masses of 1 . 8 × 10 6 and 9 . 7 × 10 6 M for baryons and DM, respectively. The Plummer-equi v alent gravitational softening length is set to 2 . 66 ckpc , limited to a maximum proper length of 0 . 7 pkpc . Additionally, we present findings from the EAGLE-Recal simulation (Recal-L25N752 in Schaye et al. 2015 ) in Appendices B and C . EAGLE-Recal has a box side length of 16 . 9 h -1 Mpc and achieves twice the spatial resolution with an initial baryonic particle mass of 2 . 3 × 10 5 M . We present results from the EAGLE-Recal simulation at z ≈ 0 to compare with the behaviour of the fiducial EAGLE model, though the findings we will present in this paper are qualitatively similar for the different resolutions.

Photoheating and radiative cooling are implemented in EAGLE based on the work of Wiersma et al. ( 2009 ), including the influence of eleven elements: H, He, C, N, O, Ne, Mg, Si, S, Ca, and Fe (Schaller et al. 2015 ); with the ultraviolet (UV) and X-ray background described by Haardt & Madau ( 2001 ). Metals are not passively diffused between gas particles in proximity, ho we ver cooling calculations MNRAS 532, 3417-3440 (2024) are conducted using SPH-kernel weighted metallicities. Since the EAGLE simulations do not have the resolution to model cold, interstellar gas, a density-dependent temperature floor is imposed (normalized to T = 8000 K at n H = 10 -1 cm -3 ). To model star formation, a metallicity-dependent density threshold is set, abo v e which star formation is locally permitted (Schaye 2004 ), and gas particles meeting this threshold are converted to star particles stochastically. The SFR is set by a tuned pressure law (Schaye, Carswell & Kim 2007 ), calibrated to the observations of Kennicutt ( 1983 ) at z = 0.

The stellar feedback subgrid model in EAGLE accounts for energy deposition into the ISM from radiation, stellar winds, and superno va e xplosions (types Ia and II). This is implemented based on the prescription outlined in Dalla Vecchia & Schaye ( 2012 ), via a stochastic thermal energy injection to gas particles in the form of a fixed temperature boost, T SF = 10 7 . 5 K. The average energy injection rate from young stars is given by f th × 8 . 73 × 10 15 erg per gram of stellar mass formed, assuming simple stellar population with a Chabrier ( 2003 ) stellar initial mass function (IMF) and that 10 51 erg is liberated per supernov a e vent. The v alue of f th is a function of local gas density and metallicity, ranging between 0 . 3 -3 in the fiducial EAGLE model. The exact parameter choices for the stellar feedback model have been calibrated to observations of the z ≈ 0 galaxy stellar mass function and M -R 50 relation (see Crain et al. 2015 ).

SMBHs of mass 10 5 h -1 M are seeded when a halo exceeds a virial mass of 10 10 h -1 M . Subsequently, SMBHs can grow via Eddington-limited accretion (utilizing a Bondi-Hoyle parametrization), as well as mergers with other SMBHs (Bondi 1952 ;Springel, Di Matteo & Hernquist 2005 ;Schaye et al. 2015 ). Similar to stellar feedback, AGN feedback in EAGLE also involves the injection of thermal energy into surrounding particles in the form of a temperature boost of T AGN = 10 8 . 5 K (in the reference physics run; Schaye et al. 2015 ). The rate of energy injection from AGN feedback is determined using the SMBH accretion rate, and a fix ed conv ersion efficienc y of accreted rest mass to energy . Importantly , we also note that in EAGLE, particles influenced by stellar or AGN feedback are not decoupled from the hydrodynamics when they receive a temperature boost. The choice of temperature boost values are relatively high in order to a v oid rapid cooling and dissipation of the feedback energy. The exact parameter choices in the EA GLE A GN feedback sub-grid model are informed by observations of the M -M BH relation (Crain et al. 2015 ).

The Next Generation Illustris simulations, known as IllustrisTNG (often shortened to TNG), are a collection of cosmological magnetohydrodynamical simulations conducted using the moving-mesh refinement code AREPO (Springel 2010 ;Pakmor, Bauer & Springel 2011 ). Several physics modules in IllustrisTNG were modified relative to its predecessor, Illustris, to achieve better agreement with observations (e.g. Weinberger et al. 2017 ;Springel et al. 2018 ). All of the TNG boxes assume a Planck Collaboration XIII ( 2016 ) cosmology, with = 0 . 692, m = 0 . 31, b = 0 . 0486, H 0 = 67 . 7 km s -1 Mpc -1 , σ 8 = 0 . 8159, and n s = 0 . 97.

The IllustrisTNG suite comprises three primary simulation sets: TNG50, TNG100, and TNG300 (Pillepich et al. 2018 ;Nelson et al. 2019 ). Each simulation set adopts a different box size and consists of three runs with varying mass resolutions. For our work in this paper, we make use of the TNG100 run with side-length 75 . 0 h -1 Mpc to balance mass resolution and volume. The highest resolution TNG100 box, which we adopt for this study, has an initial gas cell mass of 1 . 40 × 10 6 M , and DM particle mass of 7 . 5 × 10 6 M .

Following the original Illustris simulation (Vogelsberger et al. 2014 ), TNG adopts the model of Springel & Hernquist ( 2003 ) to model star formation and the pressurization of the ISM. TNG tracks stellar population evolution and chemical enrichment from supernovae Type Ia, II, as well as asymptotic giant branch stars, with individual accounting for nine elements (H, He, C, N, O, Ne, Mg, Si, and Fe). Metal-enriched gas undergoes radiative cooling in a redshift-dependent, spatially uniform ionizing UV background radiation field, with corrections for self-shielding in the ISM (Katz, Hernquist & Weinberg 1992 ;Faucher-Gigu ère et al. 2009 ). The cooling contribution from metal lines is included from density, temperature, metallicity, and redshift, following Wiersma et al. ( 2009 ) and Smith, Sigurdsson & Abel ( 2008 ). Cooling is further influenced by the radiation field of nearby AGNs, combining the UV background with the AGN radiation field (Vogelsberger et al. 2014 ). In TNG, diffusion of metals between adjacent gas elements is implemented, modelled based on a gradient extrapolation method (Vogelsberger et al. 2014 ;Pillepich et al. 2018 ).

Stellar-driven winds are modelled by isotropically injecting thermal and kinetic energy into wind particles, which are temporarily decoupled from hydrodynamic forces. With 10 per cent of the feedback energy injected thermally, the kinetic component is parametrized by following prescriptions for mass loading and injection velocity that depend on redshift, gas phase metallicity and the velocity dispersion within a weighted kernel o v er the 64 nearest DM particles (denoted as σ DM ). A wind particle is recoupled to the gas cell it currently occupies when it either drops below a specific density (5 per cent of the density threshold for star formation) or reaches a maximum traveltime (2.5 per cent of the current Hubble time). The former criterion is the dominant re-coupling mode, and normally corresponds to a wind particle travelling a distance of a few kpc (Pillepich et al. 2018 ).

The black hole growth and feedback physics in TNG is described in detail in Weinberger et al. ( 2017 ). When haloes exceed a critical mass of 5 × 10 10 M , a black hole is seeded with an initial mass of 8 × 10 5 h -1 M . After seeding, black holes can grow via mergers, or two accretion modes: (i) low accretion state, for which feedback is implemented in the form of kinetic winds; or (ii) the high accretion mode, for which feedback is implemented by depositing thermal energy. The accretion and feedback mode is determined by comparing the accretion rate onto the BH with a critical rate that is black hole mass-dependent. In both of these accretion modes, feedback is injected into surrounding gas cells (the 'feedback region') with no preferential direction, and there is no hydrodynamic decoupling of feedback-affected gas elements.

The calibration process for the TNG model, discussed in Pillepich et al. ( 2018 ), uses various observations such as the stellar mass function, cosmic SFR density as a function of redshift, BH mass versus stellar mass relation at z = 0, hot gas fraction in galaxy clusters at z = 0, and the galaxy mass n . size relation at z = 0.

2. 1.3 The SIMBA simulations SIMBA (Dav é et al. 2019 ) is the next generation of the MUFASA cosmological galaxy formation simulations (Dav é, Thompson & Hopkins 2016 ), and use the meshless finite mass (MFM) mode of the GIZMO hydrodynamics code (Hopkins 2015(Hopkins , 2017 ) ), with gravity solver based on GADGET -3 (Springel 2005 ). The SIMBA MNRAS 532, 3417-3440 (2024) simulations adopt cosmological parameters = 0 . 70, m = 0 . 30, b = 0 . 048, H 0 = 68 km s -1 Mpc -1 , σ 8 = 0 . 82, and n s = 0 . 97. The flagship SIMBA simulation box has a side length of 100 h -1 Mpc , a gas element mass of 1 . 8 × 10 7 M , and DM particle mass of 9 . 7 × 10 7 M . Due to data storage constraints, we make use of a smaller 50 h -1 Mpc SIMBA box (with the same element resolution), which we find to provide sufficient statistical power for the current study (m50n1024 in Dav é et al. 2019 ).

The simulations track the production of eleven elements (H, He, C, N, O, Ne, Mg, Si, S, Ca, and Fe) originating from Type II and Ia supernovae as well as stellar evolution, with some metals lock ed aw ay in the process of dust formation. Radiative cooling and photoionization heating are implemented using the GRACKLE -3.1 library (Smith et al. 2017 ). Star formation in SIMBA follows the Kennicutt-Schmidt Law (Kennicutt 1998 ), scaled by the H 2 fraction calculated for each particle based on local column density and metallicity (Krumholz & Gnedin 2011 ). Star-formationdri ven outflo ws are implemented as decoupled two-phase winds, characterized by updated mass loading factors derived from particle tracking in the Feedback in Realistic Environments (FIRE) zoomin simulations (Angl és-Alc ázar et al. 2017b ). Wind velocities are deriv ed from Murato v et al. ( 2015 ). 30 per cent of wind elements are injected 'hot', by coupling thermal energy to the surrounding gas. Wind elements are recoupled when their velocity relative to surrounding gas elements is less than 50 per cent of the local sound speed. The exact parameter choices in the SIMBA stellar feedback subgrid model are informed by observations of the z ≈ 0 stellar mass function.

The SIMBA simulations employ two distinct modes for SMBH accretion. One mode involves cold gas supported by rotation, following a gravitational torque model based on Hopkins & Quataert ( 2011 ) and Angl és-Alc ázar et al. ( 2017a ). Accretion in this mode is capped at three times the Eddington limit. The second mode applies to hot gas supported by pressure following the standard Bondi-Hoyle-Lyttleton prescription Bondi ( 1952 ), capped at the Eddington limit. A unique feature of the Angl és-Alc ázar et al. ( 2017a ) torquebased accretion model is that there is no need to calibrate the AGN feedback efficiency to reproduce the M -M BH relation, rather, only an assumption about the accretion efficiency is required.

SIMBA also incorporates feedback from AGN through radiative and jet modes. The radiative mode drives winds with velocities from 500 -1500 km s -1 , while the jet mode occurs for SMBHs with Eddington ratios below 0.2 and masses of at least 10 7 . 5 M . The wind velocity depends on black hole mass (where v wind = 1000 km s -1 for M BH = 10 9 M ), and the jet mode can add a velocity boost of up to 7000 km s -1 , meaning outflows can reach speeds of up to 8000 kms -1 . Gas ejected by the jets remains decoupled for a duration that scales with the Hubble time, allowing jets to travel up to ≈ 10 kpc before the energy is deposited into the surrounding gas. Gas in the jets is injected at the halo virial temperature.

Gas flows across a boundary can be measured in either a Lagrangian or Eulerian manner. In hydrodynamical simulations, the Lagrangian method often involves a gas element tracking scheme whereby the mass of gas elements that are found to have crossed a boundary between two outputs, t i and t f , is summed and normalized by t = t ft i (e.g. Mitchell et al. 2020 ;Wright et al. 2020 ). This method is simple to implement in codes where the gas is discretized such that a single gas element represents the same 'parcel' of gas between snapshots; for instance in the case of EAGLE (SPH) and SIMBA (MFM).

The other method for measuring gas flux at a boundary is to take an Eulerian approach, where the instantaneous flow of mass at a boundary can be calculated based on the velocity and mass of boundary elements (e.g. Nelson et al. 2019 ). This method does not require the gas elements to represent the same gas parcel across snapshots as the former method does, it requires only one simulation output as opposed to two, and is more comparable to instantaneous gas flow measurements inferred from observations. As such, we elect to use the Eulerian method for our gas flow calculations, which can be directly applied to the EAGLE SPH outputs, the SIMBA/ GIZMO MFM outputs, and the TNG/ AREPO outputs (without the need to use tracer particles in the latter case, Genel et al. 2013 ).

In this work, we focus on calculating gas flow rates across spherical apertures surrounding central galaxies in each of the simulations. We thus make halo catalogues which are derived using a friends-offriends method combined with (i) SUBFIND in the case of EAGLE and TNG (Springel, Yoshida & White 2001 ;Dolag et al. 2009 ), or (ii) ROCKSTAR in the case of SIMBA, to identify central versus satellite subhaloes in all of the simulations. The only characteristics extracted from these catalogues are halo centre of potential (which, by definition, is centred on the central galaxy), velocity, and halo mass/radius based on a spherical o v erdensity of 200 × ρ c ( M 200c and R 200c , respectively). The rest of the derived baryonic characteristics of haloes are calculated with our own analysis, and we do not expect the small differences in halo finding techniques to affect our results.

At each of the spherical apertures of interest centred at radius r = R, we identify boundary gas elements as those residing in a spherical shell of r = 0 . 2 × R, or between r = 0 . 9 × R and 1 . 1 × R. We see very little quantitative differences between gas flow rates measured for reasonable choices of r (0 . 05 -0 . 5 × R, or fixed r = 1 -10 kpc ). These boundary elements are categorized as being either outflow or inflow depending on the sign of their radial velocity, where the radial velocity of a gas element i relative to a halo centre j can be calculated as v r, ij = v i j • r i j / | r i j | . The gas flow rates at shell r = R around halo j for each of these subsets of boundary gas elements, i ∈ k, can be calculated with equation ( 1 ) as follows:

where Ṁ k ( r = R) is the mass flow rate across the spherical boundary at r = R. For completeness, in Appendix A we demonstrate that our method reproduces previous Lagrangian-based gas flow rates in EAGLE (from Mitchell et al. 2020 ; Fig. A1 ), and previous Eulerianbased gas outflow rates in TNG (from Nelson et al. 2019 ; Fig. A2 ).

In this section, before presenting our analysis of gas flows between the simulations in Section 2.2 , we compare static properties of the baryons within haloes in EAGLE, TNG, and SIMBA at z = 0. Here, we compare our measurements the baryon content of haloes between simulations, and its breakdown between different components. We note that we reserve a detailed discussion of the link between these findings and previous studies for Section 5.1 . Fig. 1 shows the stellar mass relative to the halo mass as a function of halo mass, for central galaxies in each of the three simulations. For comparison, we show the relation derived by Moster et al. ( 2018 ), who use an empirical model that carefully connects observed MNRAS 532, 3417-3440 (2024) galaxy properties to simulated DM haloes. In EAGLE, TNG, and SIMBA, the stellar-halo mass ratios are calculated directly from the simulations, including star particles within a 30 kpc spherical aperture of the halo centre of potential, and normalizing by the halo mass within the R 200c radius, M 200c . We convert the virial masses used as the x -values in Moster et al. ( 2018 ) from the Bryan & Norman ( 1998 ) o v erdensity criterion to the M 200c definition we use here with HMFCALC (Murray, Power & Robotham 2013 ), ho we ver note that we do not adjust the Moster et al. ( 2018 ) f values at a given halo mass. We also note that the observations that constrain the Moster et al. ( 2018 ) empirical model are diverse, and subject to systematic uncertainties and projection effects which we do not attempt to mimic with our simulation-based measurements. Despite these differences in methodology, Fig. 1 illustrates that each simulation reproduces the reduced galaxy formation efficiency at low and high masses required to match observationally derived stellar mass functions (see also Somerville & Dav é 2015 ;Dav é et al. 2020 ). As noted abo v e, and discussed in Somerville & Dav é ( 2015 ), the relatively good agreement of the different simulations in spite of the very different implementation of subgrid physics is largely due to the fact that a high weight is placed on reproducing these observational constraints when calibrating the subgrid parameters.

In Fig. 2 , we compare the total baryon content of haloes at z = 0 as a function of halo mass (top left hand panel), and in the remaining columns show how this baryon content is distributed between stars, cold gas, and hot gas in each of the simulations (as a fraction of total halo mass in the top panels, and as a breakdown of the total baryon content in the bottom panels). In each simulation, we define ' cold gas' as the gas that is either considered star-forming, or below 5 × 10 4 K, and 'hot gas' as the remaining gas within R 200c that does not meet this criteria. We do not illustrate the contribution from black holes in the breakdown panels as this component is typically very small. In the case of haloes with several galaxies within R 200c , we include the contribution from satellite galaxies in calculating the total baryonic mass in each reservoir. Also, unlike Fig. 1 , in this plot the stellar content of the halo is determined using all stellar particles within R 200c , as opposed to the central 30 kpc . This means that in Fig. 2 , the stellar component also includes contributions from intrahalo stars and satellite galaxies within R 200c . The 'cold gas' component is not decomposed further or post-processed to provide H I and H 2 content -for a direct comparison of galaxy neutral hydrogen abundances in each of the simulations, we refer the reader to Dav é et al. ( 2020 ). Furthermore, a phase decomposition of CGM gas as a function of halo mass is also presented for SIMBA in Sorini et al. ( 2022 ).

One of the striking differences between the simulations in terms of halo baryon fraction is their behaviour in the range log 10 ( M 200c / M ) ∈ [ 11 . 5 , 12 . 5]. SIMBA and TNG roughly exhibit the same functional form in this mass range -halo baryon fractions increase with halo mass up to a turno v er mass which we denote as M T 1 . Abo v e M T 1 , halo baryon fractions decrease with halo mass until a second turno v er mass, denoted by M T 2 , is reached. Abo v e M T 2 , halo baryon fractions steadily rise again to the highest halo masses presented here at log 10 ( M 200c / M ) = 14. Physically, these transitions produce three distinct 'regimes' in terms of the M 200cf b relation: (i) sub-M T 1 , where baryon fractions increase with halo mass as a result of the decreasing efficiency of stellar feedback in removing gas from a deeper potential well; (ii) the transition range M T 1 to M T 2 , where this trend reverses as AGN feedback grows in importance, and (iii) high-mass > M T 2 haloes where similar to regime (i), the efficiency of AGN feedback in removing baryonic material decreases in very massive group-cluster mass haloes. We note that in the case of EAGLE, in the M T 1 -M T 2 range, there is only a slight flattening in halo baryon fractions with increasing mass, as opposed to a clear turno v er. These halo baryon fractions, as well as the behaviour of simulation runs without AGN feedback (shown in fig. 8 in Pillepich et al. 2018 for TNG, and fig. 7 in Correa et al. 2018a for EAGLE) indicate that AGN activity is responsible for this turno v er/flattening behaviour, respectively, between M T 1 -M T 2 .

At halo masses below log 10 ( M 200c / M ) ≈ 11 . 5, corresponding to sub-M T 1 masses, the left-hand panel of Fig. 2 demonstrates that all of the simulations presented have total baryon fractions that increase with halo mass. While this behaviour with halo mass is the same, the baryon fractions themselves are markedly different. In this mass regime, the fiducial L100-REF EAGLE model predicts the lowest baryon fractions of all the simulations presented. At a halo mass of log 10 ( M 200c / M ) ≈11, the L100-REF model predicts a median halo baryon fraction of ≈ 15 per cent of the universal baryon fraction ( b / m ≈ 0 . 16). While we do not display the results here, we note that in the L25-RECAL model at this mass, halo baryon fractions are slightly higher at ≈ 30 per cent , though still lower than the corresponding baryon fraction predictions from SIMBA at ≈ 40 per cent of the expected cosmological baryon fraction. TNG predicts the highest halo baryon fractions at this mass, at ≈ 60 per cent of the universal value (see also Davies et al. 2019a ;Davies et al. 2019b ).

In this sub-M T 1 mass regime, we also note significant differences in the way that this baryonic matter is partitioned between the simulations. In the fiducial EAGLE model, the low o v erall baryon fractions mean that in order to match the expected stellar content of haloes, stars make up o v er 25 per cent of the baryonic content, compared to SIMBA and TNG where the stellar component only represents 10 per cent of the baryonic matter. For halo masses below log 10 ( M 200c / M ) ≈ 11 . 5, we note that the EAGLE L100-REF run predicts very low halo cold gas content -al w ays less than half of the predicted stellar content. This is not the case in TNG and SIMBA, where the stellar to cold gas mass ratio is near unity in this Downloaded from https://academic.oup.com/mnras/article/532/3/3417/7713461 by guest on 01 September 2024 mass range . Dav é et al. ( 2020 ) show that the EAGLE fiducial run underpredicts the observed z = 0 H I MF from ALF ALF A (Arecibo Le gac y F ast ALFA) observations by ≈ 0 . 5 dex (Jones et al. 2018 ); but matches the local CO luminosity function from Saintonge et al. ( 2017 ) fairly well (see also Lagos et al. 2015 ). This suggests that the fiducial EAGLE run lacks the gas associated with the more diffuse atomic phase at z = 0, ho we ver this result is heavily dependent on how gas phases are partitioned in post-processing. While not directly shown in this figure, we remark that the L25-RECAL EAGLE model produces a better match to the local H I MF; and that for z ≥ 1, that the EAGLE L100-REF and L25-RECAL models are consistent with regard to predicted H I and H 2 mass functions.

While SIMBA and TNG predict a similar functional form for baryon fractions across halo mass, their normalization and the location of transition masses M T 1 and M T 2 are markedly different (as also found in Ni et al. 2023 ;Delgado et al. 2023 ;Crain & van de Voort 2023 ). In TNG, the first transition mass occurs at log 10 ( M T 1 / M ) ≈ 11 . 8; while in SIMBA, the transition occurs at much lower halo mass, log 10 ( M T 1 / M ) ≈ 11 . 3. Furthermore, in TNG the second transition mass occurs at log 10 ( M T 2 / M ) = 12 . 6; while in SIMBA, halo baryon fractions do not increase until log 10 ( M T 2 / M ) ≈ 12 . 8. In this transition range, halo baryon fractions in TNG are, on average, 20 per cent-50 per cent higher than those measured in SIMBA. Furthermore, at log 10 ( M 200c / M ) = 14, baryon fractions in SIMBA only reach ≈70 per cent of b / m , while haloes in TNG haloes approach this uni versal v alue very closely. Instead of a turno v er , EA GLE predicts a flattening in halo baryon fractions between log 10 ( M T 1 / M ) ≈ 12 . 1 and log 10 ( M T 2 / M ) ≈ 12 . 7. In this transition range, the drop in TNG baryon fractions meets the slowly rising EAGLE baryon fractions at the M T 2 transition point.

Abo v e M T 2 , EAGLE and TNG agree very well in terms of predicted baryon fractions as a function of halo mass.

In the high-mass regime, it is possible to compare predictions for halo gas content with observations. In Kugel et al. ( 2023 ), a number of carefully compiled X-ray and weak lensing (WL) observations were used to calibrate a large-volume cosmological simulation, FLAMINGO. Such observational data sets are nominally cast as (hot) gas fractions within R 500c , as a function of M 500c , representing a slightly shrunken aperture relative to the rest of our measurements. In the bottom left panel of Fig. 2 , we present like-forlike measurements of hot halo gas content in each of the simulations within R 500c , and compare with this set of observations compiled in Kugel et al. ( 2023 ). We refer the reader to this study for a detailed explanation of the curated data set, noting that we only make use of the WL data from Akino et al. ( 2022 ), and that the rele v ant X-ray data sets come from the following studies: Vikhlinin et al. ( 2006 ), Maughan et al. ( 2008 ), Rasmussen & Ponman ( 2009 ), Sun et al. ( 2009 ), Pratt et al. ( 2010), Lin et al. ( 2012 ), Lagan á et al. ( 2013), Sanderson et al. ( 2013 ), Gonzalez et al. ( 2013), Lovisari, Reiprich & Schellenberger ( 2015 ), Pearson et al. ( 2017 ), andLovisari et al. ( 2020 ). As abo v e, we define 'hot' gas as an y gas elements with a temperature T ≥ 5 × 10 4 . We note that this choice of threshold does not significantly influence gas fractions in this mass range, as the majority of gas is at or abo v e halo virial temperature (see top panels). Regardless, we stress that this selection is unlikely to exactly mimic the gas accounted for in the observ ations, ho we ver allo ws us to understand the difference between simulation predictions and their broad agreement with observations. Our measurements between M 500c = 10 13 . 5 M -M 500c = 10 14 . 5 M indicate that median halo gas fractions in TNG and EAGLE are slightly higher than the MNRAS 532, 3417-3440 (2024) compiled observations (more so in EAGLE), and slightly lower than observations in SIMBA. The range of predictions from the simulations is largest between M 500c = 10 13 M -M 500c = 10 14 M , while at higher masses ( M 500c ≈ 10 14 . 5 M ), the gas fractions tend towards a tighter range. For reference, we note that the offset between M 200c and M 500c halo mass measurements is relatively constant at ≈ 0 . 15 dex in this mass range. While the X-ray data sets extend to higher masses, there is a scarcity of cluster-mass haloes in the simulations to used for this study due to their box size.

In Appendix D , we show a breakdown of the baryon content of haloes at z = 2, finding that the qualitative behaviour of the simulations in relation to baryon fraction across halo mass remains similar to the results presented here at z = 0, with a universal shift to higher halo baryon fractions. The transition masses M T 1 and M Tfoot_2 remain similar in EAGLE and TNG at z = 2 compared to z ≈ 0 galaxies; though this is not the case for SIMBA. In SIMBA, the transition masses are higher than at z = 0, with log 10 ( M T 1 / M ) = 12 . 2 and log 10 ( M T 2 / M ) = 12 . 8 at z = 2 (see also fig. 7 in Sorini et al. 2022 ).

To aid in the interpretations of the gas flow measurements in the following section (Section 4 ), in Appendices B and C we provide radial density and temperature profiles of central galaxies in each of the simulations. Together with the figures contained in this section, these radial profiles provide additional information about the spatial distribution and properties of gas within haloes in each of the simulations.

In this section, we aim to identify the physical mechanisms behind the differences in the baryon budget predicted by the EAGLE, TNG, and SIMBA simulations as explored in Section 3 . Our method for measuring gas flow rates is detailed in Section 2.2 . We note that for the purposes of this analysis, we do not enforce any additional selection cuts to classify elements as inflow or outflow other than the sign of their radial velocity.

In Fig. 3 , we display measurements of galaxy outflow rates at the ISM scale in a form commonly presented in previous literature: the ISM scale mass loading factor, η = Ṁ out / SFR , where Ṁ out is the mean gas outflow rate, and Ṁ is the mean SFR. We begin by comparing the mass loading results between simulations at z ≈ 2, which builds a picture for the operation of gas flows at the epoch of maximum SFR density (Madau & Dickinson 2014 ). In the bottom panel of Fig. 3 , we also compare the median specific SFR (sSFR, Ṁ /M ) of galaxies. This assists in understanding the influence of SFR on the mass loading measurements provided in the top panel. In the bottom panel, we also include a hatched grey region which represents an indicative 'quenched' region where the sSFR drops below log 10 ( sSFR / yr -1 ) = -11 + 0 . 5 z (Furlong et al. 2015 ).

We define the stellar mass of a galaxy in the same way as Fig. 1 , including all star particles within 30 kpc of the halo centre of potential. Similarly, the SFRs we use in calculating mass loading are e v aluated with g as within 30 kpc of the g alaxy centre, which is fixed regardless of the scale that we measure outflow rates. We define the 'ISM' scale by measuring gas flows across a spherical boundary at r = 0 . 25 × R 200c , which we consider to be the outer scale of the ISM. We stress that we do not enforce a minimum radial velocity for outflow gas; meaning that while we do select slowmoving outflow gas elements/particles, these outflow values are fair to compare between the simulations and represent the true mass flux at each scale. 2 Furthermore, for the mass flow rate measurements at the scale of the ISM, we do not require the gas to be 'cold' to be included in the calculation. This choice is made to simplify subsequent comparisons with gas flows at larger scales, where in many cases the majority of the gas has been heated to near the virial temperature of the halo (see Appendix C ).

We follow the approach of Mitchell et al. ( 2020 , based on arguments presented in Neistein et al. 2012 ) in calculating mass loading: instead of comparing η ≡ Ṁ out / SFR for each halo individually, we calculate η for each halo mass as the mean outflow rate in each halo mass bin divided by the mean SFR in the bin. Neistein et al. ( 2012 ) demonstrate that this method correctly predicts the average rate of mass exchange when inte grated o v er time, and also mitigates the impact of discreteness noise that would otherwise influence the results. We note that the percentile ranges displayed in the top panel MNRAS 532, 3417-3440 (2024) of Fig. 3 correspond to the range of mass loading values predicted based on the 16th and 84th percentiles of mass outflow rate in a given stellar mass bin, without including variance in SFRs at a given mass.

For reference, in Fig. 3 , we also include indicative η scalings with

200c as dashed and dasheddotted lines respectively (see Murray, Quataert & Thompson 2005 for deri v ation).

At the ISM scale of r = 0 . 25 R 200c in the top panel of Fig. 3 , it is clear that each of the simulations predicts a similar ne gativ e scaling in η = Ṁ out / SFR as a function of M , until a turno v er mass where outflows begin to pick up relative to SFRs (also demonstrated in Nelson et al. 2019 andMitchell et al. 2020 for TNG and EAGLE respectively, though using different scales and methodologies). The turno v er mass differs slightly between simulations -in EAGLE and TNG, the turno v er occurs at M ≈ 10 10 and ≈ 10 10 . 3 M , respecti vely; with η v alues at the turno v er masses reaching a minimum of η ≈ 10 0 . 5 and ≈ 10 0 . 7 , respectively. In SIMBA, the turnover is more gradual, with the minimum average mass loading occurring closer to M ≈ 10 11 M . Folding in the stellar-halo mass relation in each of the simulations, in EAGLE and TNG the turno v er mass occurs at log 10 ( M 200c / M ) ≈ 11 . 8 -12, while in SIMBA, the turno v er occurs at a slightly higher mass of log 10 ( M 200c / M ) ≈ 12 . 4 at z ≈ 2 (as is also the case for halo baryon content).

Below the turno v er mass, the scaling of η with mass for each of the simulations falls between the expected 'momentum conserving' ( η ∝ M . 6 at M ≈ 10 9 M , followed by TNG predicting approximately one-third of this outflow rate at η = 10 1 . 2 , and EAGLE predicting approximately half of the TNG outflow rates at η = 10 0 . 9 , for the same stellar mass.

The higher mass loading factors predicted in TNG and SIMBA at this scale are likely a result of the combination of input mass loading prescriptions, velocities, and the temporarily decoupled nature of these wind gas elements. Even though the wind gas elements in both simulations are likely to have recoupled before reaching the scale where we measure outflow rates (0 . 25 × R 200c ), the ability of these gas elements to escape the densest regions of the ISM without losing as much momentum means that the mass flux at this scale can remain quite high. Nelson et al. ( 2019 ) show in relation to the TNG50 simulation that while the velocity and η of SF-driven feedback are prescribed at injection, the way in which these outflows propagate to larger scales is non-trivial. The mass loading factors presented at the ISM scale agree with the measurements at 10 kpc presented in Nelson et al. ( 2019 ), and the SIMBA mass loading at this scale agrees with the input mass loading prescription from Angl és-Alc ázar et al. ( 2017b ) taken from the FIRE simulations. In comparison, in EAGLE, thermal SF-driven outflows do not have prescribed velocities or mass loading. The feedback energy is deposited purely thermally, and the large temperature boosts ( T SF = 10 7 . 5 K) and associated pressure produce strong outflows.

Abo v e the mass turno v er where AGN feedback becomes the dominant ejection mechanism, there are further clear differences between mass loading in the simulations. TNG demonstrates the steepest upturn in ISM-scale mass loading with halo mass, reaching η = 10 1 . 7 at log 10 ( M 200c / M ) = 13. At this same mass, EAGLE predicts a mass loading of approximately half the TNG value at 10 1 . 4 , and SIMBA predicts mass loading at one tenth of the EAGLE value at 10 0 . 6 (corresponding to the increased transition mass at which AGN feedback becomes efficient as per Fig. 2 ). The lower η values at high mass in SIMBA are likely linked to the onset of jet-mode feedback, which injects more energy than the QSO-mode feedback with the same momentum flux, corresponding to a much lower mass loading factor relative to the BH accretion rate.

Comparing the median sSFRs in the bottom panel of Fig. 3 , we note that there is broad agreement between the simulations in terms of the shape of the star formation main sequence. There are, ho we ver, small quantitati ve dif ferences between the simulations -for instance, the median sSFR in SIMBA at M ≈ 10 9 M sits 0 . 1 -0 . 2 dex below EAGLE and TNG, and sits 0 . 3 -0 . 4 dex above EAGLE and TNG for M ≈ 10 10 . 2 -10 10 . 7 M . Together, these slight differences in star formation activity complicate comparisons between outflow behaviour in the simulations using η. Furthermore, instantaneous measurements of mass loading do not account for the delay between star formation and the associated feedback. As such, for the remainder of this paper, we quote flow rates without normalizing by SFR. Additionally, citing the (slight) differences in the stellar-halo mass relation between simulations as presented in Fig. 1 , we choose to compare the baryon cycle between simulations as a function of halo mass rather than stellar mass for the remainder of this work. Given the very similar underlying cosmological parameters, the form of the halo mass functions in the simulations are very closely matched; meaning that using halo mass as the dependent variable allows us to fairly compare the simulations at a set depth of gravitational potential well.

In Figs 4 and 5 , we show the raw outflow rates for central galaxies in each simulation as a function of halo mass, normalized by b / m × M 200c at z = 2 and 0, respectively ('outflow efficiency').

In the left-hand panels, we show these outflow rates at the ISM scale (0 . 25 × R 200c ). In the remaining panels for each simulation, we show outflow rates at the ISM scale, the halo scale (defined at 1 . 00 × R 200c ), and at the IGM scale (defined at 2 . 50 × R 200c ).

The values of outflow efficiency indicate the inverse of the timescale that would be required for a gas mass of M 200c × b / m to be remo v ed at this scale. In the left-hand panels of Figs 4 and 5 , we compare median outflow rates in each of the simulations at the ISM scale (0 . 25 × R 200c ). In each of the remaining panels, for each of the simulations we compare outflow rates at three scales: and 2 . 50 × R 200c . First, focusing on the left-hand panels which compare the simulations at the ISM scale, we note that at both redshifts, below log 10 ( M 200c / M ) ≈ 12, total outflow rates in TNG are actually higher than in SIMBA, in contrast to the mass loading factors which are higher in SIMBA. This can be explained by the slightly lower normalization of the M -sSFR relation in SIMBA compared to TNG (see Fig. 3 for sSFR measurements at z ≈ 2 and Dav é et al. 2019 for measurements at z ≈ 0) -meaning that galaxies in SIMBA eject more gas per unit SFR compared to TNG galaxies, but have lower absolute total outflow rates. This highlights the fact that if a galaxy is gas poor, the amount of gas available for removal is proportionally reduced.

It is when we analyse the behaviour of outflow rates as a function of scale as per the right-hand panels in Figs 4 and 5 , that dramatic differences between gas flows in the simulations begin to emerge. We first focus on the z ≈ 2 case, as presented in MNRAS 532, 3417-3440 (2024) the assumed wind velocities at injection, which in SIMBA are taken from Muratov et al. ( 2015 ). Abo v e log 10 ( M 200c / M ) ≈ 12 . 5 at z ≈ 2, where AGN feedback begins to dominate the driving of outflows, the behaviour between the simulations is also different. In EAGLE, the outflow rate at the ISM and halo scales remain similar in this mass range, but drop dramatically at the 2 . 5 × R 200c scale, with Ṁ out at 2 . 5 × R 200c approximately one-tenth of the value at R 200c in haloes with mass log 10 ( M 200c / M ) ≈ 13. This indicates that AGN-dri ven outflo ws in EAGLE tend to stall near or just outside the virial radius at z ≈ 2, and seldom reach larger scales. TNG shows qualitatively similar beha viour, b ut a less dramatic drop off in outflow rate with increasing scale. SIMBA actually shows a monotonically increasing outflow rate with increasing scale at log 10 ( M 200c / M ) ≈ 11 . 5 -13, and a constant outflow rate from R 200c to 2 . 5 × R 200c at higher halo masses.

This indicates that the outflows driven by AGN in TNG and SIMBA are more powerful in terms of spatial reach, likely related to the fact that in both of these models, AGN feedback is implemented with a significant kinetic component rather than just thermally as in the case of EAGLE.

MNRAS 532, 3417-3440 (2024) Moving to the z ≈ 0 universe in Fig. 5 , we note that the qualitative differences between simulations remain similar to the z ≈ 2 picture at the ISM scale, ho we ver the disparities between how the outflows propagate across scales become even more pronounced. In the lefthand panel, we can see that the normalization of gas flow rates relative to halo mass are much lower than at z ≈ 2, with longer dynamical time-scales and lower gas densities at this epoch. One exception to this is the strength of AGN-dri ven outflo ws in SIMBAwhich are much stronger at z ≈ 0, and also become efficient at lower halo masses. This is likely a result of jet-mode feedback switching on more frequently at low redshift due to the requirement of low Eddington ratios (typical Eddington ratios decline at low redshift; see Thomas et al. 2019 , Fig. 5 ).

One very interesting feature at z ≈ 0 in the EAGLE simulations is the apparent entrainment of outflowing gas with increasing scale. We find that at log 10 ( M 200c / M ) ≈ 11, the mass outflow rate at 2 . 5 × R 200c is o v er 10 times higher than that at 0 . 25 × R 200c in EAGLE at z ≈ 0. This means that the stellar-feedback-driven outflows in EAGLE are not only able to reach scales beyond the halo, but also entrain mass along the way. Mitchell et al. ( 2020 ) have previously investigated this effect, finding that the outflowing gas is significantly o v erpressurized relativ e to the ambient CGM at z = 0. This causes the tenuous material in the CGM to be swept up with the bulk flow.

As demonstrated in Fig. 4 , in the regime dominated by stellar feedback, M halo 10 12 M , we find that in EAGLE there is minimal mass entrainment in outflows at z ≈ 2 compared to z ≈ 0 (with the exception of haloes below 10 11 M , where outflows at 2 . 5 × R 200c are enhanced relative to the ISM and halo scale by a factor of a few). This means that while the outflows do not significantly stall for halo masses of log 10 ( M 200c / M ) ≈ 11 -11 . 5 at z ≈ 2, we do not find any strong evidence for an increase in mass outflow rates at larger scales as is predicted at later times in EAGLE. This result can be interpreted with the aid of the density and temperature profiles of haloes that we include in Appendices B and C . Comparing the radial density profiles of gas in log 10 ( M 200c / M ) ≈ 11 at z ≈ 2 and ≈ 0 in Figs B1 and B2 , we find that in all simulations, gas densities in the CGM universally decrease from z ≈ 2 to ≈ 0 . In TNG and SIMBA, the difference between gas densities at z ≈ 2 and ≈ 0 is ≈ 1 . 3 -1 . 5 dex at r = 0 . 50 × R 200c . This is not an unexpected result -this decrease in density roughly corresponds to the cubed change in scale factor between these redshifts. In the case of EAGLE, ho we ver, the decrease in CGM gas density is much more dramatic, closer to ≈ 2 -2 . 5 dex from z ≈ 2 to ≈ 0. This means that at z ≈ 0, outflows are o v erpressurized in a relativ e sense to the CGM (as per the findings in Mitchell et al. 2020 ). Comparing this picture to z ≈ 2, where CGM densities between all simulations are similar, the lack of mass entrainment at z ≈ 2 in EAGLE is likely the result of a much denser CGM at this redshift compared to z ≈ 0, meaning that outflows are not o v erpressurized relativ e to the surrounding gas, and not able to sweep up the surrounding material with the bulk flow.

An additional factor which is helpful in interpreting different levels of gas entrainment in outflows is the temperature of the gas in the CGM, as presented in Appendix C . Inspecting the temperature profiles of haloes at z ≈ 2 (Fig. C1 ), in the log 10 ( M 200c / M ) ≈ 11 and ≈ 12 bins, median CGM gas temperatures in EAGLE are substantially lower than T vir (estimated from van de Voort 2017 ); particularly where r/R 200c 0 . 75 (the 'inner' CGM). This aligns with the commonly accepted picture that 'cold-mode' accretioncool, filamentary inflow from the cosmic web -is dominant at higher redshift and lower halo masses, and that these inflows can constitute a significant portion of the CGM without being shock-heated (e.g. Kere š et al. 2005 ). Comparatively, at z ≈ 0 in EAGLE, gas temperatures in the CGM are significantly higher -within 0 . 2 -0 . 3 dex of T vir at log 10 ( M 200c / M ) ≈ 12. The low co v ering fractions typical of gas inflows at high redshift (e.g. Wright et al. 2021 ) mean that feedback-driven outflows are less likely to encounter (and thus entrain) such gas in the CGM.

As discussed abo v e, we find that in TNG, outflow rates monotonically decrease with increasing scale below log 10 ( M 200c / M ) ≈ 12 at both z ≈ 0 and ≈ 2, indicating stalling of outflows as they propagate to larger distances. Interestingly, at z ≈ 2, Fig. C1 shows that for TNG haloes in bins log 10 ( M 200c / M ) ≈ 11 and ≈ 12, gas in the inner CGM is significantly hotter than in EAGLE. The inner CGM likely corresponds to the scale at which much of the stellar feedback energy is thermalized. Unlike EAGLE, we find that gas densities in the CGM remain relatively high in this mass range at z ≈ 0 in TNG, which in addition to the relatively small thermal component of stellar feedback in TNG (10 per cent of feedback energy), contributes to the outflowing material being underpressurized relative to the surrounding medium, and thus unable to reach larger scales or entrain mass.

Considering larger haloes at z ≈ 0, an interesting feature in SIMBA is the efficient entrainment of AGN-driven gas outflows abo v e log 10 ( M 200c / M ) ≈ 12. At log 10 ( M 200c / M ) ≈ 12 . 5, the gas outflow rate at 2 . 5 × R 200c is roughly 5-10 times the ouflow rate at r = 0 . 25 × R 200c -an effect not seen in EAGLE, nor TNG. At z ≈ 0, the ambient density of the CGM at this mass scale is quite low (as demonstrated in Fig. 2 and Fig. B2 ). This, together with the direct heating of gas elements surrounding AGN, means that the AGN-dri ven outflo ws are o v erpressurized relativ e to the surrounding medium, and thus able entrain tenuous gas in both the CGM, and the region surrounding the halo.

In this section, we build on the results of Section 4.1 to study the gas inflow rates onto haloes as a function on halo mass, scale, and redshift in the EAGLE, TNG, and SIMBA simulations. This helps us understand feedback not only from an 'ejective' perspective, but also a 'pre ventati ve' perspecti ve; where feedback actively prevents the introduction of baryons to haloes, in addition to directly removing gas via outflows. This is very important for semi-analytic and semiempirical models of galaxy formation, which typically assume that the amount of gas accretion traces the growth of the DM halo (Moster et al. 2018 ;Behroozi et al. 2019 ;Pandya et al. 2020 ;Wright et al. 2020 ).

In Figs 6 and 7 , we show the relationship between inflow rate and halo mass in each of the simulations at z ≈ 2 and ≈ 0. In the lefthand panel, we plot inflow rates (normalized by M 200c × b / m ; ' inflo w ef ficiency') at the R 200c halo boundary as a function of M 200c for EAGLE, TNG, and SIMBA, compared with the 'expected' gas accretion rates from DM accretion rates, computed by multiplying the DM accretion rate times the cosmological baryon fraction (grey lines, Wright et al. 2020 ). The 'inflow efficiency' measurements plotted here correspond to the inverse of the time-scale that would be required for a gas mass of M 200c × b / m to be accreted at this scale. In the second-fourth columns, we compare inflow rates within the same simulations between three scales: 0 . 25 × R 200c (within the ISM), R 200c , and at 2 . 5 × R 200c . This comparison allows us to see ho w inflo w rates change from cosmological (IGM) scales do wn to the ISM scale in each of the simulations.

First focusing on the left-hand panel for the z ≈ 0 case in Fig. 7 , we remark that the respective gas accretion rates between the simulations MNRAS 532, 3417-3440 (2024) 2020) indicating the gas accretion rate expected if gas inflows perfectly traced the infall of DM in grey. In all panels, error bars correspond to the 16th-84th percentile range in inflow rates at a given mass, and hatched regions correspond to the bootstrap-generated 90 per cent confidence interval on the medians at a given mass. All simulations predict that AGN feedback has a strong pre ventati ve impact on gas accretion at the ISM scale, ho we ver there is considerable variability between the prev entativ e impact of stellar feedback. Panels 2-4: inflo w ef ficiencies at dif ferent scales in EAGLE, TNG, and SIMBA, respectively. In each of these panels, ISM-scale inflows are plotted with lighter colours, halo-scale inflows are plotted in medium colours, and IGM-scale inflows are plotted in darker colours. We also include measurements from Wright et al. ( 2020 ) indicating the gas accretion rate expected if gas inflows perfectly traced the infall of DM in grey. In all panels, error bars correspond to the 16th-84th percentile range in inflow rates at a given mass, and hatched regions correspond to the bootstrap-generated 90 per cent confidence interval on the medians at a given mass. The prevention of gas inflow at the halo scale due to stellar feedback is strongest in EAGLE, followed by SIMBA, then TNG. The prevention of gas inflow at the halo scale due to AGN feedback is strongest in SIMBA, followed by TNG, then EAGLE.

at R 200c follow a similar functional form to the halo-wide baryon fractions presented in Fig. 2 . In the stellar feedback-dominated regime, log 10 ( M 200c / M ) 11 . 5, the lowest specific accretion rates are seen in EAGLE, where gas inflow rates are just ≈ 10 per cent of the expectation based on the halo growth rate, due to the farreaching influence of stellar feedback (as directly demonstrated in Section 4.1 ). Accretion rates onto haloes of the same mass are much higher in TNG, with gas accretion rates ≈ 50 per cent of the expectation based on the DM halo growth rate. Halo-scale gas accretion in TNG remains fairly high due to the reduced efficiency of outflows beyond R 200c in this mass range, as was demonstrated in Figs 4 and 5 . The predictions from SIMBA regarding gas inflow rates in the stellar-feedback dominated regime fall between those of EAGLE and TNG -where stellar feedback has a stronger pre ventati ve ef fect than in TNGat the halo scale, but not as dramatic as in the case of EAGLE.

In the AGN-dominated regime, log 10 ( M 200c / M ) 12 . 5, the simulations predict very dif ferent le vels of pre ventati ve feedback. This is most dramatic in SIMBA at log 10 ( M 200c / M ) ≈ 12 -12 . 5, where the amount of gas inflow to haloes is less than 10 per cent of the expectation based on the total halo growth rate times the cosmological baryon fraction. The halo-scale gas inflow efficiency in EAGLE continues to rise with halo mass, indicating that AGN feedback plays little role in preventing gas inflow at this scale (see Wright et al. 2020 , who report a decrease in halo gas accretion of only ≈ 20 per cent compared to the EAGLE run with no AGN feedback).

In TNG, we note a small drop in inflow efficiency due to AGN, where at log 10 ( M 200c / M ) ≈ 12 . 5, the inflow efficiency drops to ≈ 25 per cent of that expected from the total halo growth rate -less dramatic than seen in SIMBA, but about twice as much suppression as predicted by EAGLE. In all simulations, approaching halo masses of log 10 ( M 200c / M ) ≈ 13 . 5, the ability of AGN to suppress halo-scale gas accretion decreases; related to the reduced influence of outflows at larger scales which are unable to o v ercome the deeper potential well (see the decreases in outflow rates at 2 . 5 × R 200c at this mass scale in Figs 4 and 5 , with the exception of SIMBA at z ≈ 0).

In the three right-hand panels of each figure, we explore how inflow rates in EAGLE, TNG, and SIMBA vary across scales. In EAGLE, at z ≈ 0, we clearly see that the rate of gas inflow decreases with radius across all halo masses. Below log 10 ( M 200c / M ) ≈ 11 . 5, the biggest drop in inflow efficiency occurs between r = 2 . 5 × R 200c and 1 × R 200c , indicating that this is where a large portion of the cosmological inflow begins to encounter the outflowing and/or recycling material, in agreement with our findings discussed in Section 4.1 regarding the scale at which outflows stall in EAGLE. At log 10 ( M 200c / M ) 12, we note that while cosmological inflow is not greatly decreased at the halo-scale in EAGLE, the rate of accretion is significantly suppressed at r = 0 . 25 R 200c (as also reported in Correa et al. 2018a , b ). This agrees with our mass loading results, which suggest that the strength of AGN-dri ven outflo ws at the scale of the ISM are similar between EAGLE and TNG, but that they stall at smaller radius in EAGLE.

Inspecting inflow rates across different scales in TNG, we find that for haloes with log 10 ( M 200c / M ) 11 . 5, inflow rates at the ISM scale (0 . 25 × R 200c ) are similar to inflow rates at the halo boundary at z ≈ 2; but at z ≈ 0, they are actually twice as high at the scale of the ISM as those recorded at the scale of the halo. With the relatively dense CGM in TNG in this mass range (see Fig. B2 , partly contributed to by the build-up of ejected, metal enriched gas in this region which does not leave the halo), we argue that the ISM-scale inflow rates are higher than halo-scale inflow rates due to a combination of efficient recycling of cold gas, and efficient metal cooling of warm CGM gas at this scale. Abo v e log 10 ( M 200c / M ) ≈ 12 . 5, we find the opposite effect -ISM-scale inflow rates are suppressed relative to halo-scale inflow rates (as is true for all simulations in this mass regime).

Lastly, focusing on inflow rates in SIMBA, we find that below log 10 ( M 200c / M ) 11 . 5, accretion rates are suppressed equally at r = 0 . 25 R 200c , R 200c , and 2 . 5 R 200c . Based on Figs 4 and 5 , this aligns with a scenario where the spatial range of stellar feedback sits between that of EAGLE (far-reaching, beyond the halo) and TNG (within the inner CGM). For log 10 ( M 200c / M ) 12, SIMBA predicts a steadily increasing inflow prevention effect moving inwards from r = 2 . 5 × R 200c to 0 . 25 × R 200c . Ev en though the prev entativ e effect is greatest at 0 . 25 × R 200c (with gas inflow rates at just 1 per cent-5 per cent of those expected based on the halo growth rate at 10 12 . 5 M ), SIMBA predicts that even gas at 2 . 5 × R 200c experiences significant pre ventati ve feedback due to AGN. While we do not present the results here, we note that in the SIMBA run with no jetmode feedback, the pre ventati ve ef fect is greatly reduced -indicating that it is the low-accretion rate, high kinetic energy jet-mode AGN feedback in SIMBA that is responsible for the strong suppression of gas inflow.

To summarize -in EAGLE, TNG, and SIMBA, the level of inflow prevention at a given scale can be clearly linked with the presence, or lack thereof, of feedback-driven outflows at the same scale. For log 10 ( M 200c / M ) 11 . 5, EAGLE and SIMBA predict that stellar feedback produces a strong pre ventati ve ef fect on halo-scale inflo w rates, while in TNG, halo-scale gas inflow is suppressed by a smaller amount. For halo masses log 10 ( M 200c / M ) 12, SIMBA predicts very strong AGN-driven suppression in halo-scale gas inflow rates, relative to a more mild prev entativ e effect in TNG, and very minimal pre ventati ve feedback at the halo scale in EAGLE.

The EAGLE, TNG, and SIMBA simulations have all been largely successful in producing a reasonable match to the observed galaxy stellar mass function at z ≈ 0, albeit with careful tuning required to calibrate their respective subgrid feedback models that cannot be implemented ab initio . We have demonstrated that different implementations of feedback processes -and thus, the manner in which the baryon cycle operates o v er cosmic time -can be degenerate in producing galaxies at z ≈ 0 with very similar stellar masses and SFRs. In Section 5.1 , we use our measurements of gas flows rates across the simulations in Section 4 to make the causal link between the baryon cycle and static halo properties. In Section 5.2 , we compare our results to previous literature investigating the baryon cycle in simulated galaxies. In Section 5.4 , we outline a number of promising observational avenues by which the degeneracy between different feedback models could be broken. Lastly, in Section 5.5 , we discuss how the finite resolution of the simulations used in this study may influence our results.

In Fig. 8 , we summarize the main findings of this study in visual form -indicating the scale reached by feedback-dri ven outflo ws in each of the simulations and the consequent gas content of the CGM in each of the simulations. We split this diagram into the low-mass ( M 200c 10 12 M ) and high-mass ( M 200c 10 12 M ) regimes, where stellar and AGN feedback respectively dominate. Our findings in Section 4 aid in the interpretation of our findings in Section 3 -in particular Fig. 2 -where we compare the total baryon content of haloes in each simulation, and its breakdown into stellar and gaseous components. First, we discuss the low-mass regime, log 10 ( M 200c / M ) 12, where stellar feedback dominates the behaviour of the baryon cycle. In EAGLE, the efficient removal of gas to large scales -where outflows reach several times the virial radius at both z ≈ 2 and ≈ 0 -can be clearly linked with the reduced halo gas content below log 10 ( M 200c / M ) ≈ 11 . 5 in Fig. 2 . This is maintained by the reduced halo gas inflow rates as a consequence of direct prevention of cosmological gas accretion, as shown in Section 4.2 .

In contrast, in TNG, stellar feedback-dri ven outflo ws typically stall before reaching R 200c . This deposits a significant amount of enriched gas in the CGM that can subsequently be efficiently recycled into the ISM, resulting in increased ISM-scale inflow rates relative to even halo-scale inflow in TNG at z ≈ 0. This efficient intrahalo recycling scenario means that in low-mass haloes, the baryon content of the CGM in TNG is much higher than that recorded in EAGLE at z ≈ 0. SIMBA predicts a picture in between that of EAGLE and TNG, where some gas is fully remo v ed from the halo via stellar feedback (particularly below log 10 ( M 200c / M ) ≈ 11), and some portion of these outflows stall within the CGM. This leads to halo-wide baryon fractions and CGM gas content between that of EAGLE and TNG at z ≈ 0.

Next, we discuss the high-mass regime -log 10 ( M 200c / M ) 12 -where AGN feedback becomes dominant in these simulations, and plays a driving role in the baryon cycle. In EAGLE, we find that the purely thermal implementation of AGN feedback can remo v e gas from the ISM, but is seldom powerful enough to remo v e gas to large scales and significantly influence cosmological inflow. This is likely a consequence of rapid thermal dissipation as the gas travels through the CGM, meaning that the outflows stall before travelling far beyond the halo. Even so, this thermal energy (in combination MNRAS 532, 3417-3440 ( 2024) with the transition to a hot, virialized halo atmosphere -e.g. Correa et al. 2018a , b ) is sufficient to prevent cosmological inflow from entering the ISM (see Figs 6 and 7 ), and ultimately curtail SFRs in high-mass galaxies. Since little gas is expelled from the halo and the pre ventati ve impact of AGN works within the CGM, this leads to the result shown in Fig. 2 where galaxy baryon fractions monotonically increase with halo mass.

In Section 4 , we show in Figs 4 and 5 that AGN-driven outflows in TNG and SIMBA reach much larger scales before stalling, presumably because of the kinetic implementation of feedback in the 'jet mode' in both models. At z ≈ 0, we show that for haloes with log 10 ( M 200c / M ) ≈ 12 . 5, gas outflow rates at 2 . 5 × R 200c are just as high as they are at R 200c in TNG, and even higher in the case of SIMBA. Unlike in EAGLE (where the pre ventati ve ef fect of AGN is confined to within the CGM), these outflows also act to prevent cosmological gas accretion at the halo scale. Together, particularly in the case of SIMBA, this causes the downturn in halo baryon fractions at intermediate halo mass which we demonstrate in Fig. 2 .

In this section, we outline how our results compare with previous theoretical studies relating to the baryon cycle of galaxies. In particular, we focus on comparing how outflow and inflow rates change with scale relative to the simulations we study here.

As outlined in Section 2.2 and Appendix A , we have demonstrated that our methodology for measuring gas flow rates produces results that are consistent with previous studies utilizing the EAGLE simulations (Mitchell et al. 2020 ) and the TNG suite (Nelson et al. 2019 ). Fig. 14 in Mitchell et al. ( 2020 ) shows the difference between z ≈ 2 outflow mass loading reported by Nelson et al. ( 2019 ) in the TNG50 simulations, compared with those measured in EAGLE. For galaxies with stellar mass below ≈ 10 10 M , where stellar feedback is ef ficient, gas outflo w rates in TNG50 tend to decrease from smaller (20 kpc ) to larger (50 kpc ) radii. Comparatively, in EAGLE the mass loading is identical at both of these scales. Our work, which uniquely applies a lik e-for-lik e methodology to measure gas flow rates between the simulations, provides a comprehensive explanation for this result -showing that SF-driven outflows produce a larger-scale impact on the baryon cycle in EAGLE compared to TNG.

We can also compare our results with those from zoom-in simulations, where impro v ements in resolution enable more detailed subgrid modelling for gas cooling, star formation and feedback within the ISM. 3 3 . This is, perhaps, unsurprising, given that SIMBA uses the FIRE power-law fits as the prescriptive mass loading factor and wind velocity. The mass loading values at fixed stellar mass in TNG and EAGLE are, thus, lower than the predictions from FIRE (for the mass range considered here, M 10 9 M ). Focusing on an individual system (m12i), Muratov et al. ( 2015 ) show that in a halo with M 200c ≈ 10 12 M at z ≈ 2, mass loading values at the ISM scale (0 . 25 × R 200c ) are stochastic (like the outflows themselves), ho we ver typically take a value of η ≈ 10 1 in episodes of star formation. In this same system, Muratov et al. ( 2015 ) show that outflow rates at R 200c reach only ≈ 10 per cent -20 per cent of the outflow rate at the ISM in episodes of star formation. This aligns with the predictions we outline from TNG and SIMBA, where MNRAS 532, 3417-3440 (2024) SF-dri ven outflo ws typically recycle within the CGM at this mass scale. Comparatively, in EAGLE, mass flow rates at both radii are very similar at z ≈ 2 (on average). Pandya et al. ( 2021 ) investigate the mass, momentum, and energy loading in the subsequent FIRE-2 simulations (Hopkins et al. 2018 ). When using the same radial shell to define the ISM and no minimum outgoing velocity requirement (as per the method in Muratov et al. 2015 ), they find a very similar mass loading scaling with stellar mass. At z ≈ 0 in a halo of mass M 200c ≈ 10 12 M (m12f), they show that halo-scale gas outflow rates rarely exceed those from the ISM escape -again, in alignment with the picture painted by TNG and SIMBA at this mass. Interestingly, for dwarf haloes with M 200c ≈ 10 10 M , Pandya et al. ( 2021 ) find that mass outflow rates at the halo scale nominally exceed that measured at the ISM scale, and at M 200c ≈ 10 11 M , that outflow rates at the two scales can be roughly equal (if the time delay as the outflow traverses the ISM is taken into account). The former case is similar to the picture painted by EAGLE at M 200c ≈ 10 11 M , where mass is entrained in outflows through the CGM. This indicates that in FIRE-2, the scale reached by outflows is strongly influenced by halo mass. Christensen et al. ( 2016 ) use a particle-tracking technique to study gas flows in and around galaxies in a suite of GASOLINE -based simulations (Wadsley, Stadel & Quinn 2004 ). These simulations use a 'blast-wave' technique (where cooling is temporarily disabled) to distribute energy from stellar feedback, as per Stinson et al. ( 2006 ). As demonstrated in Mitchell et al. ( 2020 ), Christensen et al. ( 2016 ) find lower ISM-scale mass loading factors than those measured in EAGLE, and a steeper decline in η with mass. Since the ISM-scale mass loading values are lowest in EAGLE of the three simulations we discuss here (see Fig. 3 ), this indicates that the mass loading values in Christensen et al. ( 2016 ) are also lower than measured in SIMBA and TNG. While direct measurements of halo-scale outflow rates are not presented, they find that a substantial fraction of the gas that leaves the ISM will also leave the halo virial radius. Tollet et al. ( 2019 ) present an analysis of gas flows in the NIHAO simulations using a Lagrangian particle tracking method (Wang et al. 2015 ). NIHAO also uses the blast-wave technique to distribute energy from stellar feedback, as per Stinson et al. ( 2006 ). As a function of stellar mass, mass loading values in NIHAO are similar to TNG and SIMBA at M ≈ 10 9 M ( η ≈ 10 1 . 5 ), ho we ver decrease sharply with stellar mass to η ≈ 10 -0 . 5 at M ≈ 10 10 . 5 M . At the same stellar mass, predictions from EAGLE, TNG, and SIMBA are one full decade higher at η ≈ 10 0 . 5 . For haloes with M 200c 10 12 M , they show that the majority of gas does not exit R 200c after being ejected from the ISM.

The total baryon fraction of haloes within R 200c we measure align with the findings presented in Ayromlou, Nelson & Pillepich ( 2023 ), who also analyse the EAGLE, TNG, and SIMBA simulations to sho w ho w feedback redistributes the baryons in and around haloes. This is measured in terms of the 'closure radius' of haloes -the minimum radius at which the enclosed baryon content reaches the uni versal v alue. They find that for low-mass haloes where stellar feedback dominates ( M 200c 10 11 . 5 M ), the closure radius is the largest for EAGLE galaxies, followed by SIMBA galaxies, and then TNG galaxies. In the case of higher mass haloes where AGN feedback dominates ( M 200c 10 12 M ), the closure radius is the largest for SIMBA galaxies, followed by TNG galaxies, and then EAGLE galaxies.

Our findings also provide insight into the measurements presented in Dav é et al. ( 2020 ), who investigate the global cold gas (H I and H 2 ; decomposed in post-processing) content of galaxies in the same simulations and compare with available observations. At z ≈ 0, Dav é et al. ( 2020 ) show that the fiducial resolution EAGLE run underpredicts the H I content of galaxies compared to local ALF ALF A observations of the H I -MF (Jones et al. 2018 ), as well as the H I fraction of galaxies as a function of M . Given that the H I identified by this decomposition algorithm tends to trace the gas that is slightly more diffuse and extended than H 2 gas, the low H I content of EAGLE galaxies could be explained by the stellar feedback scenario we outline abo v e, which e v acuates much of the CGM -and, likely, some of the gas identified as H I which resides the edge of the ISM. We do note, ho we ver, that predictions from the recalibrated, higher resolution EAGLE run and H I observations and are in closer agreement (see also Figs B1 and B2 ).

Our results highlight the fact that differences in baryon cycling can indeed leave signatures even at the scale of the ISM, particularly when considering multiphase gas. With the resolution and physics impro v ements in the next generation of cosmological simulations, it will become possible to better resolve the phase structure of the ISM in a self-consistent manner -and thus, to produce statistical predictions for the spatially resolved properties of the ISM, as well as the influence of gas flows on these distributions. Such advancements will also provide the statistical context necessary to generalize the results of high-resolution zoom simulations (some of which we summarize abo v e), which pro vide more detailed insight into the link between the ISM and CGM.

Semi-analytic models (SAMs) have been used extensively as a complement to numerical hydrodynamic simulations to make predictions for how galaxies form and evolve and to interpret observations (for a re vie w, see Somerville & Dav é 2015 ). Instead of explicitly solving the equations of gravity, hydrodynamics, etc. numerically, SAMs instead track bulk flows of mass and metals into and out of reservoirs that typically include gas (and gas phase metals) in the IGM, CGM, and ISM, as well as stars within galaxies. The models are set within cosmological 'merger trees' which represent the formation of structure and growth of haloes via mergers and cosmological accretion with the CDM paradigm. Phenomenological, parametrized scaling relations describe processes such as the rates of cooling of hot halo gas and its flow into the ISM, conversion of cold ISM gas into stars, ejection of gas by stellar feedback, and accretion onto and feedback from SMBHs. These free parameters are typically tuned to match a set of galaxy observables or quasi-observables, which have traditionally focused on stellar and ISM properties such as stellar mass functions (or stellar mass versus halo mass relations), cold ISM gas fraction, stellar mass versus stellar phase metallicity, SFR distrib utions, etc. Different SAM groups ha v e adopted v ery different mass loading factors for stellar driven winds (see e.g. the comparison presented in Fig. 9 of Mitchell et al. 2020 ), in spite of predicting very similar stellar mass functions and stellar-to-halo mass ratios. The results presented here provide important insights for future SAMs.

Some SAMs do not have a direct channel for stellar feedback to eject CGM material from the halo (e.g. Santa Cruz, GALFORM SAMs; Somerville et al. 2008 ;Bower, Benson & Crain 2012 ;Lacey et al. 2016 ). Gas and metals can be ejected from the ISM and are either deposited in the CGM or into an 'ejected' reservoir. Some SAMs compute the total energy budget produced by stellar populations and parametrize the fraction that drives an ISM-scale outflow, and assume that any remaining energy can be used to eject gas from the CGM (e.g. L-galaxies, SHARK, and related models; Henriques et al. 2015 ;Lagos et al. 2018 ). AGN feedback in some models can also help eject MNRAS 532, 3417-3440 (2024) material out of the ISM to the CGM, and potentially from the CGM to the IGM (e.g. Lagos et al. 2024 ).

Our study shows that in some hydrodynamical simulations (e.g. EAGLE), the ejected mass on CGM scales ( R 200c ) can be higher than that on ISM scales (0.25 R 200c ), indicating that mass is entrained in the winds. This suggests that SAMs could include this additional channel for removal of material from the CGM by stellar-driven winds, and also provides some guidance for priors on the coupling factors that are currently treated as free parameters.

In most SAMs, gas and metals in the ejected reservoir can 'reaccrete' into the CGM on a specified time-scale. This time-scale, and how it scales with halo mass and time, are uncertain and different choices are adopted in different SAMs. Though not studied explicitly here, our results demonstrate that different hydrodynamic simulations are likely to predict very different time-scales for this 'reaccreted' gas, since they show very different rates of gas accretion at halo scales ( R 200c ).

SAMs, as well as empirical models (Behroozi et al. 2019 ;Moster et al. 2018 ) typically assume that the mass inflow rate into the CGM is equal to f b Ṁ tot , where Ṁ tot is the total rate of mass accretion into the halo. 4 As discussed e xtensiv ely in Section 5.1 , 'prev entativ e' feedback (where inflows are considerably smaller than f b Ṁ tot ) can be quite important in simulations, but its strength at different scales varies depending on the subgrid physics. Only a few SAMs have included pre ventati ve feedback on halo scales due to stellar feedback processes (Lu, Mo & Wechsler 2015 ;Hirschmann et al. 2014 ;Pandya et al. 2023 ).

Most SAMs model AGN feedback by allowing thermal energy from a notional radio jet to offset some of the CGM cooling, reducing the rate of cooling from the CGM to the ISM. The AGN is not allowed to eject gas from the CGM to the IGM, or to change the density profile or temperature of the CGM. The Santa Cruz and SHARK SAMs (Somerville et al. 2008 ;Lagos et al. 2024 ) include mass ejection from the ISM due to radiatively efficient AGN winds. This material is assumed to leave the halo and not to return, due to the high assumed launch velocities in the case of the Santa Cruz SAM, while in SHARK the material may be trapped in the CGM or escape to the IGM depending on the wind energy. In both models, ho we ver, this was found to have little impact on the results. It is clear from the results presented here and elsewhere (e.g. Ayromlou et al. 2023 ) that in order to more faithfully depict the baryon cycle in numerical hydrodynamic simulations, SAMs must include recipes for how AGN winds and jets remo v e gas and metals from the CGM (e.g. Lagos et al. 2024 ), as well as for the associated pre ventati ve feedback.

A promising approach for modelling ejection from the CGM and pre ventati ve feedback associated with either stars or an AGN has been proposed by Carr et al. ( 2023 ) and Pandya et al. ( 2023 ): in addition to tracking inflows and outflows of mass and metals, these SAMs also include tracking of energy sources (e.g. from stellar and AGN feedback) and sinks (e.g. from cooling and turbulent dissipation). Excess energy deposited into the CGM can 'lift' it out of the halo. An extension of these ideas is presented by Voit et al. ( 2024aVoit et al. ( , 2024b ) ) in which the CGM is in generalized energy balance within the gravitational potential of the halo, and can expand or contract in response to energy deposition or loss.

As discussed in Introduction, direct observation and quantification of the gas flows surrounding galaxies is very challenging -particularly when seeking to compare the baryon cycle at a statistical level. This means that currently, it is very difficult to place constraints on which feedback scenarios are fa v oured by observations. Here, we briefly outline some avenues against which these models could be tested with current and/or future observations.

Comparing the average gas density profiles between simulations at z ≈ 2 and ≈ 0 in Figs B1 and B2 , it is clear that the differences in total halo gas content arising from disparate baryon cycling paradigms are clearest at z ≈ 0, after the variability between models has had adequate time to consolidate in terms of static halo properties. Presently, a combination of X-ray observations and WL studies have been able to place constraints on the gas fraction of high-mass galaxy groups and clusters ( M 500c 10 13 M , e.g. Hoekstra et al. 2015 ;Pearson et al. 2017 ;Mulroy et al. 2019 ;Lovisari et al. 2020 ;Akino et al. 2022 ); which can be used to calibrate the strength of AGN feedback in the next generation of cosmological hydrodynamical models (e.g Schaye et al. 2023 ;Kugel et al. 2023 ). As demonstrated in Fig. 2 , such observations indicate that the baryon fractions of haloes of mass M 500c ≈ 10 14 M within R 500c are slightly too high in EAGLE and TNG (very close to the universal value), and slightly too low in SIMBA.

As demonstrated in Fig. 2 , the most significant AGN-driven difference in halo baryon content is induced below the mass scale of clusters, between M 200c ≈ 10 12 -10 13 M . Measurements of the baryon content of haloes in this mass range would provide significant constraining power on the strength of AGN feedback, as well as the mass scale at which SMBH accretion and feedback becomes efficient. Preliminary studies using the eROSITA telescope (Predehl et al. 2021 ) and stacked soft X-ray observations indicate that for galaxies with stellar mass M ≈ 10 11 M , forward-modelled X-ray luminosity profiles (up to ≈ R 200c ) from EAGLE and TNG are consistent with observational constraints (Oppenheimer et al. 2020 ;Chadayammuri et al. 2022 ). Further constraints on the baryon fraction of haloes in this mass range will become possible with observations utilizing the Sun yaev-Zeldo vich (SZ) effect -for instance, Bregman et al. ( 2022 ) demonstrate the feasibility of placing such constraints in a small sample of local ≈ L * haloes. Yang et al. ( 2022 ) demonstrate that the disparate baryon content of haloes between TNG and SIMBA abo v e log 10 ( M 200c / M ) ≈ 12 . 5 produces an observationally detectable difference in forward-modelled SZ measurements, provided that an instrument with adequate angular resolution is utilized (such as the upcoming Simons Observatory, Ade et al. 2019 ). In this mass range, studies investigating the H I content of groups using stacking techniques may also pro v e useful to constraining baryon content (e.g. Dev et al. 2023 ), particularly with radio surv e ys such as DINGO (Deep Investigation of Neutral Gas Origins) with the ASKAP (Australian Square Kilometre Array Pathfinder) telescope (Hotan et al. 2021 ) and MIGHTEE-HI (MeerKAT International GigaHertz Tiered Extragalactic Exploration) (Maddox et al. 2021 ).

Below mass scales of M 200c ≈ 10 12 M , it is very difficult to observationally constrain the baryon content of galaxies and their CGM. In this mass range, the CGM is expected to be multiphase and anisotropic (e. g. McCourt et al. 2018 ;Fielding et al. 2020 ), constituted of a combination of cool filamentary inflow structures together with enriched recycling gas. Observations have shown that sightline measurements of CGM metallicity can vary by up to 2 dex within the same halo (e.g. Lehner et al. 2013 ;Prochaska et al. 2017 ;Zahedy et al. 2019 ). The resolution of the cosmological MNRAS 532, 3417-3440 ( 2024) hydrodynamical simulations we study here is inadequate to resolve detailed phase structure, and more detailed subgrid models are also likely required to accurately model the CGM (e.g. Crain & van de Voort 2023 ). We can, ho we ver, make general comments regarding the expected impact of different feedback regimes on the properties of the CGM.

In the case of EAGLE, where the CGM is predicted to be nearly entirely devoid of gas at low halo masses, Wright et al. ( 2021 ) show that the integrated metallicity of the CGM is extremely sensitive to the rate of pristine halo-scale gas accretion. This is a result of newly introduced gas being able to efficiently dilute the metallicity of the already small gaseous reservoir, which is only possible due to the pre vious e v acuation of the CGM via far-reaching, strong stellar feedback. Thus, in a model like EAGLE, we might expect to find significant amounts of very metal poor to pristine gas in the CGM, which would not be the case in a model like TNG where the metalenriched gas from SF-driven outflows is directly deposited into the CGM. Despite this, we refrain from comparing the models in terms of the distribution of metals in the CGM in this work, as it is nontrivial to disentangle differences in the modelling of metal diffusion (included in TNG, but not in EAGLE or SIMBA) from changes induced by differences in gas flows.

In any case, it is clear that the scale and strength of stellar feedback will leave signatures in the distribution of metals in the CGM. P éroux et al. ( 2020 ) demonstrate that TNG and EAGLE predict similar net gas flows rates as a function of galactocentric azimuthal angle for galaxies with M ≈ 10 9 . 5 and ≈ 10 10 . 5 M . Investigating the scale of metal-enrichment at different impact parameters, particularly studying how angular and radial variations would manifest in discrete sightline observations, could constitute a promising test of model accuracy. This is potentially possible with carefully constructed comparisons with current surv e ys that probe different phases of the CGM, such as COS-Halos (Cosmic Origins Spectrograph -Halos; Tumlinson et al. 2013 ;Werk et al. 2014 ;Prochaska et al. 2017 ) and MEGAFLOW (MusE GAs FLOw and Wind;Schroetter et al. 2016 ;Zabl et al. 2019 ).

Alone this line, DeFelippis et al. ( 2021 ) compare results from the MEGAFLOW surv e y with predictions from the TNG50 simulation (several fold higher resolution than the TNG100 simulation we investigate here, but the same feedback physics); indicating that TNG50 can reasonably reproduce observations of the kinematic diversity of strong Mg II absorbers (tracing the cool CGM) in the halo mass range 10 11 . 5 -10 12 M . Appleby et al. ( 2021 ) compare results from SIMBA and the COS-Halos and COS-Dwarfs surv e y, finding general agreement with observations regarding the predicted abundance of H I and metal absorbers in the CGM of star-forming SIMBA galaxies, but find slight tension in the quenched population. Appleby et al. ( 2023 ) also demonstrate the feasibility of using machine learning to infer CGM conditions from absorption line properties, finding that their model trained on SIMBA can accurately predict absorber o v erdensity, temperature, and metallicities.

The Habitable Worlds Observatory, a large-space-based UVoptical-IR mission concept endorsed by the US Astro2020 Decadal report, will provide greatly enhanced spectroscopic sensitivity in the UV relative to the Hubble Space Telescope Cosmic Origins Spectrograph, enabling significantly larger samples of z 1 galaxies with absorption line studies of their CGM (NAS 2021 ). Furthermore, new ways of probing the CGM continue to be proposed and tested -for instance using fast radio bursts (Macquart et al. 2020 ), and the characterization of of Ly α haloes in both emission (e.g. Lokhorst et al. 2019 ;Augustin et al. 2019 ) and absorption (i.e. Ly α forest, see Tillman et al. 2023a , b ).

Observations of gas flow rates -as opposed to static measurements of halo baryon content -will also pro v e useful in constraining feedback models in simulations. Different techniques can provide measurements on the multiphase nature of galactic outflows in individual systems; for instance CO tracing cool molecular outflows (e.g. Leroy et al. 2015 using Atacama Large Millimetre Array observations), or with metal absorption lines and detailed ionization modelling (e.g. see discussions in Chisholm et al. 2016 ;Veilleux et al. 2020 ). Combining careful forward modelling from simulations with observational approaches across wavelength may be able to place constraints on total gas outflow rates as the sample size of outflow measurements continues to grow.

It is important to note that the resolution of each of these simulations we discuss is inadequate to resolve the multiphase structure of the CGM -in particular, any interfaces between cold and hot CGM gas (see e.g. McCourt et al. 2018 ;Fielding et al. 2020 ). Creasey et al. ( 2011 ) show that a mass resolution of < 10 6 M is desirable to a v oid numerical o v ercooling of accretion shocks onto haloes with SPH -a criterion which is not met in any of the simulations, particularly SIMBA. At even finer resolutions, Rey et al. ( 2024 ) show in AMR simulations of a galaxy with mass M = 10 8 M that increasing spatial resolution from > 200 to ≈ 20 pc can lead to twofold increases in mass loading factors at r = 5 kpc . They show this is due to the better resolved 'hot' outflowing phase ( T 10 5 K) being heated to higher temperatures, and staying hot for longer periods. Unfortunately, without any subgrid refinement schemes, the simulations we analyse here simply cannot simultaneously resolve scales from Mpc down to the ≈ 0 . 1 -10 pc required to model cold gas in the ISM or CGM.

These considerations beg the question as to whether improved resolution could influence the results presented in Sections 3 and 4 . As shown in Table 1 , the mass resolution of the EAGLE and TNG simulations are relatively similar, but there is a notable difference in SIMBA, where the gas element mass is roughly 10 times higher than the two aforementioned simulations. While we do not show the results here, we have conducted a comparison between each simulation in this paper with a higher resolution counterpart (EAGLE-L25N752-Ref in the case of EAGLE, TNG50 in the case of TNG100, and a higher resolution 25 Mpc SIMBA box in the case of SIMBA).

Upon investigation, we noted that there is very little difference between the propagation of feedback-driven outflows in TNG100 and TNG50. In the case of stellar feedback, this likely due to the fact that feedback affected gas is 'decoupled' from the hydrodynamics at injection (directly a v oiding any immediate overcooling), and the mass loading is set based on local DM dispersion as opposed to baryonic properties. At higher halo masses, TNG50 haloes are very slightly more baryon-rich than those in TNG100, ho we ver the scale reached by AGN-driven outflows remain identical -indicating that resolution does not qualitatively influence the conclusions we make in this paper.

In the case of EAGLE and SIMBA, we found slight differences in the propagation of outflows in low-mass haloes 10 11 . 5 M with enhanced resolution. ISM-scale outflows in both higher resolution EAGLE and SIMBA runs are enhanced relative to the standard resolution runs, where numerical o v ercooling could inhibit the escape of outflows. These slightly enhanced ISM-scale outflow rates relative to the standard resolution runs leads to the deposition of more gas in MNRAS 532, 3417-3440 (2024) the CGM, and an associated increase in ISM-scale inflow rates that offset any differences that could otherwise arise in terms of SFR.

Despite these small quantitative differences, we stress that the central findings of our study regarding the differences between subgrid prescriptions and how this manifests in the larger scale baryon cycle do remain qualitatively sound when comparing with higher resolution counterparts of each simulation.

In this paper, we have conducted a lik e-for-lik e comparison of the baryon cycle in three modern cosmological hydrodynamical simulations: EAGLE, TNG, and SIMBA. We use a common methodology to analyse the gas flows in and around central galaxies, from the scale of the ISM to several halo virial radii. While galaxies in these simulations share very similar stellar mass content and SFRs at z ≈ 0, our work has highlighted that this agreement is achieved for very different physical reasons.

Concerning haloes with mass M 200c 10 11 . 5 M , where stellar feedback is dominant, we summarize our results as follows:

(i) In EAGLE, SF-driven outflows are able to reach several times R 200c before stalling. At z ≈ 0, mass outflow rates increase with scale with the gas surrounding galaxies being very tenuous and underpressurized relative to the outflows. At z ≈ 2 we do not see mass entrainment in outflows with increasing scale, where the CGM is still relatively dense and more cool/filamentary in nature.

(ii) The extended scale of SF-driven outflows in EAGLE leads to a significant pre ventati ve ef fect on gas inflo w at the halo scale. Overall, this leads to haloes with a low density CGM and very low baryon fractions (just 10 per cent of the cosmological value) in this mass range at z ≈ 0.

(iii) In TNG, SF-driven outflows are very strong at the scale of the ISM, ho we ver, tend to stall within the CGM before reaching R 200c . This is also true in SIMBA at z ≈ 2, ho we ver at z ≈ 0, SF-dri ven outflows in SIMBA can reach a similar scale as in EAGLE (albeit without mass entrainment).

(iv) In TNG, halo baryon fractions remain abo v e 50 per cent of the cosmological value in this mass regime at z ≈ 0. SF-driven outflows are recycled within the scale of the halo, with inflow rates at the ISM scale higher than at the halo scale due to very efficient cooling within the dense, enriched CGM.

(v) In SIMBA, halo baryon fractions in this mass regime are between those measured in EAGLE and TNG, reflecting moderate pre ventati ve feedback on halo-scale gas inflows.

Furthermore, in relation to haloes with mass M 200c 10 12 M , where AGN feedback is dominant:

(i) In EA GLE, A GN-dri ven outflo ws do not reach far beyond R 200c at either z ≈ 0 nor ≈ 2. This leads to minimal pre ventati ve impact on gas inflow at the halo scale, but significantly prevents gas inflow at the scale of the ISM.

(ii) In TNG, AGN-driven outflows are more powerful than in EA-GLE, with at least some gas in these outflows reaching 2 -3 × R 200c . This leads to a moderate pre ventati ve impact on cosmological inflow at the halo scale, with further inflo w pre vention within the heated CGM. As a result, there is a significant turno v er in halo baryon fractions at z ≈ 0, reaching values as low as 50 per cent of the cosmological value at M 200c ≈ 10 12 . 5 M .

(iii) In SIMBA, AGN-driven outflows are significantly stronger than in both EAGLE and TNG, particularly at z ≈ 0. These outflows continually entrain mass from the ISM scale to several times R 200c .

This causes a very significant drop in halo-scale cosmological gas inflow, and a corresponding drop in halo baryon fractions to as low as 20 per cent of the universal value between M 200c = 10 12 -10 13 M .

This work lays the foundation for developing targeted observational tests that can fa v our or disfa v our certain feedback scenarios, and mo v e towards a more constrained understanding of the role of feedback processes in the baryon cycle. Such an understanding will be imperative to inform the next generation of cosmological simulations, and to maximize their scientific utility . Additionally , a precise understanding of the baryonic effects on the large-scale matter distribution (see Delgado et al. 2023 ;Gebhardt et al. 2024 ) will be imperative in order to extract meaningful constraints from future WL studies with regard to cosmological parameters, cosmic shear, and the growth of cosmic structure.

In this section, we compare our measurements of gas flow rates with those measured from the same simulations (where applicable) in previous literature. As outlined in Section 2.2 , we elect to use an Eulerian or 'instantaneous' approach to calculating gas flow rates, rather than a Lagrangian particle tracking method. 2020) take a Lagrangian particle-tracking approach to measure gas flow rates, we reco v er good agreement with their results (to within ≈ 0 . 2 de x across the halo mass range). Our instantaneously measured mass loading factors are systematically (very) slightly higher than the Mitchell et al. ( 2020 ) results, likely due to the latter not including the subset of gas that is very quickly recycled within the t used for the Lagrangian calculation. This offset is highest at M 200c ≈ 10 13 M , indicating that some AGN-driven outflows in EAGLE may be recycled quite quickly at the ISM scale.

Fig. A2 shows a comparison of z ≈ 2 10 kpc -scale mass loading factors as a function of M in TNG for two different choices of outflow velocity cut: 50 and 150 km s -1 . Utilizing a very similar method, our measurements agree very well, to within ≈ 0 . 1 dex for both velocity thresholds, across the range of stellar masses shown. As discussed in Section 2.2 , in the main body of this paper we do not enforce a minimum velocity threshold for the outflows. With experimentation, we find that imposing a velocity threshold does not qualitatively influence our results. The data presented in orange correspond to outflow measurements with a minimum outflow velocity of 150 km s -1 , while the data presented in blue correspond to outflow measurements with a minimum outflow velocity of 50 km s -1 . In our calculations, hatched shaded regions correspond to the bootstrap-generated 90 per cent confidence interval on the medians at a given mass. For both velocity thresholds, our measurements agree very well, to within ≈ 0 . 1 dex .

. Radial density profiles of galaxies at z = 2 in EAGLE (blue), EAGLE-RECAL (teal), TNG (purple), and SIMBA (pink). Each column represents a different bin in halo mass, increasing top to bottom. Error bars correspond to the 16th-84th percentile range in inflow rates at a given mass, and hatched regions correspond to the bootstrap-generated 90 per cent confidence interval on the medians at a given mass. Gas density profiles are relatively similar between the simulations at this redshift, corresponding to the universally higher halo baryon fractions presented in Fig. D1 .

In Figs B1 and B2 , we show the density profiles of galaxies in each of the simulations at z = 2 and 0, respectively. In general, density profiles at z = 2 agree fairly well between the simulations. This corresponds to the universally higher baryon content of galaxies and their CGM at this redshift, as demonstrated in Fig. D1 . Considering the lowest mass bin at z = 0, we note that the density of the CGM in EAGLE is quite low -≈ 1 dex lower than the predictions from TNG and SIMBA at the edge of the ISM. This low density allows for the entrainment of gas with outflows in EAGLE, with the outflows being o v erpressurized relativ e to the surrounding medium. We observe a similar effect in SIMBA in the two higher mass bins, with CGM densities in the intermediate halo mass bin being ≈ 1 dex lower than the predictions from EAGLE and TNG. This signifies the ef ficient e v acuation of the CGM via AGN feedback, which also allows subsequent outflows to be o v erpressurized relativ e to the CGM and entrain gas in the outflows towards larger scales. MNRAS 532, 3417-3440 (2024) A P P E N D I X C : G A L A X Y T E M P E R AT U R E P RO F I L E S In Figs C1 and C2 , we show the temperature profiles of galaxies in each of the simulations at z = 2 and 0, respectively. Each row represents a different bin in halo mass -top: log 10 ( M 200c / M ) ∈ [10 . 75 , 11 . 25]; middle: log 10 ( M 200c / M ) ∈ [11 . 75 , 12 . 25]; and bottom: log 10 ( M 200c / M ) ∈ [12 . 75 , 13 . 25]). Radial x-values are quoted as fractions of R 200c , with a broken scale abo v e 1 × R 200c to include radial bins up to 3 × R 200c . For reference, in each panel we include a horizontal line indicating the expected virial temperature of haloes in this mass range, with T vir ≈ 1 . 1 × 10 6 K ( μ/ 0 . 59 ) M halo / 10 12 M 2 / 3 (1 + z), from van de Voort ( 2017 ).

As discussed in Section 4.1 , considering haloes with mass M 200c 10 13 M at z ≈ 2, predictions for temperatures in the inner CGM are quite different between the simulations. In particular, in EAGLE (and to some extent, SIMBA) there remains a significant presence of cold gas between 0 . 3 -0 . 6 × R 200c . Comparatively, in TNG, temperatures in this region typically exceed 10 5 K. At lower redshift, the predictions between the simulations are relatively similar in the two lower mass bins (signifying a universal shift to a hotter CGM, and the dominance of hot-mode accretion).

In this section, we add to the z ≈ 0 results shown in Section 3 to compare the total baryon content of haloes at as a function of halo mass, and how this baryon content is distributed between stars, cold gas, and hot gas in each of the simulations. As in the main body of this paper, in each simulation, we define 'cold gas' as the gas that is either considered star-forming, or below 5 × 10 4 K, and 'hot gas' as the remaining gas within R 200c that does not meet this criteria.

At this redshift, the baryon content of haloes is universally higher in the simulations -for M 200c 10 11 . 5 M , all simulations predict a baryon fraction abo v e 50 per cent of the cosmic mean. Similar to z ≈ 0, the baryon content of EAGLE galaxies at low halo mass -M 200c 10 12 M -is lower than the predictions from SIMBA and TNG, indicating that the ejective and prev entativ e influence of stellar feedback is already at play by z ≈ 2. The transition masses M T 1 and M T 2 are similar to those reported at z ≈ 0, with the exception of SIMBA, where the M T 1 mass at z ≈ 2 increases to ≈ 10 12 M (as opposed to ≈ 10 11 . 3 M at z ≈ 0, with M T 2 remaining at ≈ 10 12 . 8 M ). For SIMBA, this indicates that AGN feedback becomes more efficient at lower masses towards z ≈ 0.