In recent years cosmological hydrodynamical simulations have proven their utility as key interpretative tools in the study of galaxy formation and evolution. In this work we present a comparative analysis of the baryon cycle in three publicly available leading cosmological simulation suites EAGLE IllustrisTNG and SIMBA. While these simulations broadly agree in terms of their predictions for the stellar mass content and star formation rates of galaxies at zapprox 0 they achieve this result for markedly different reasons. In EAGLE and SIMBA we demonstrate that at low halo masses Mrm 200clesssim 1011.5 mathrmModot  stellar feedback SFdriven outflows can reach far beyond the scale of the halo extending up to 23times Rrm 200c. In contrast in TNG SFdriven outflows while stronger at the scale of the interstellar medium recycle within the circumgalactic medium within Rrm 200c. We find that active galactic nucleus AGNdriven outflows in SIMBA are notably potent reaching several times Rrm 200c even at halo masses up to Mrm 200capprox 1013.5 mathrmModot . In both TNG and EAGLE AGN feedback can eject gas beyond Rrm 200c at this mass scale but seldom beyond 23times Rrm 200c. We find that the scale of feedbackdriven outflows can be directly linked with the prevention of cosmological inflow as well as the total baryon fraction of haloes within Rrm 200c. This work lays the foundation to develop targeted observational tests that can discriminate between feedback scenarios and inform subgrid feedback models in the next generation of simulations.
['Astrophysics - Astrophysics of Galaxies']
The baryon cycle in modern cosmological hydrodynamical simulations
{'The baryon cycle in modern cosmological hydrodynamical simulations': 'Ruby J. Wright 1 , 2 , 3 , 4 ★ , Rachel S. Somerville 2 , Claudia del P. Lagos 3 , 4 , Matthieu Schaller 5 , 6 , Romeel Davé 7 , 9 , Daniel Anglés-Alcázar 8 , 2 , Shy Genel 2 , 10  \n- 1 Department of Physics, University of Helsinki, Gustaf Hällströmin katu 2, FI-00014 Helsinki, Finland\n- 2 Center for Computational Astrophysics, Flatiron Institute, 162 5 th Avenue, New York, NY 10010, USA\n- 3 International Centre for Radio Astronomy Research, University of Western Australia, 7 Fairway, Crawley, WA 6009, Australia\n- 4 ARC Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D)\n- 5 Lorentz Institute for Theoretical Physics, Leiden University, PO Box 9506, NL-2300 RA Leiden, The Netherlands\n- 6 Leiden Observatory, Leiden University, PO Box 9513, NL-2300 RA Leiden, The Netherlands\n- 7 Institute for Astronomy, University of Edinburgh, Royal Observatory, Edinburgh EH9 3HJ, UK\n- 8 Department of Physics, University of Connecticut, 196 Auditorium Road, U-3046, Storrs, CT, 06269, USA\n- 9 Department of Physics and Astronomy, University of the Western Cape, Bellville, Cape Town 7535, South Africa\n- 10 Columbia Astrophysics Laboratory, Columbia University, 550 West 120th Street, New York, NY 10027, USA  \nAccepted XXX. Received YYY; in original form ZZZ', 'ABSTRACT': 'In recent years, cosmological hydrodynamical simulations have proven their utility as key interpretative tools in the study of galaxy formation and evolution. In this work, we present a comparative analysis of the baryon cycle in three publicly available, leading cosmological simulation suites: EAGLE, IllustrisTNG, and SIMBA. While these simulations broadly agree in terms of their predictions for the stellar mass content and star formation rates of galaxies at 𝑧 ≈ 0, they achieve this result for markedly different reasons. In EAGLE and SIMBA, we demonstrate that at low halo masses ( 𝑀 200c ≲ 10 11 . 5 𝑀 ⊙ ), stellar feedback (SF)driven outflows can reach far beyond the scale of the halo, extending up to 2 -3 × 𝑅 200c . In contrast, in TNG, SF-driven outflows, while stronger at the scale of the ISM, recycle within the CGM (within 𝑅 200c ). We find that AGN-driven outflows in SIMBA are notably potent, reaching several times 𝑅 200c even at halo masses up to 𝑀 200c ≈ 10 13 . 5 𝑀 ⊙ . In both TNG and EAGLE, AGN feedback can eject gas beyond 𝑅 200c at this mass scale, but seldom beyond 2 -3 × 𝑅 200c . We find that the scale of feedback-driven outflows can be directly linked with the prevention of cosmological inflow, as well as the total baryon fraction of haloes within 𝑅 200c . This work lays the foundation to develop targeted observational tests that can discriminate between feedback scenarios, and inform sub-grid feedback models in the next generation of simulations.  \nKey words: galaxies: formation - galaxies: evolution - galaxies: haloes - methods: numerical', '1 INTRODUCTION': 'The baryon cycle represents a complex and interconnected set of processes, including gas accretion, star formation, chemical enrichment, and large scale outflows, that collectively drive the formation and evolution of galaxies (e.g., Oppenheimer & Davé 2006; Ford et al. 2014; Somerville & Davé 2015; Tumlinson et al. 2017; AnglésAlcázar et al. 2017b; Oppenheimer et al. 2018). The interplay of these processes establishes a causal link between different scales: the interstellar medium (ISM), circum-galactic medium (CGM), and intergalactic medium (IGM); as well as a connection between the different baryonic components involved at each of these scales: stars, gas, and black holes.  \nThe cosmological inflow of gas from the IGM, through the CGM, and eventually into the ISM provides the requisite fuel for the process of star formation (Kereš et al. 2005; Sánchez Almeida et al. 2014). Star formation occurs within the dense regions of the ISM, where  \nefficient 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 (SMBH), which are known to inhabit the dense central regions of many galaxies (see reviews by Ferrarese & Ford 2005; Kormendy & Ho 2013; Heckman & Best 2014).  \nBoth star formation and accretion onto SMBH drive large scale outflows, which play a crucial role in regulating the baryon cycle (Veilleux et al. 2005, 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 supermassive black holes at galactic centers, 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).  \nOn the scale of the ISM, various observational techniques provide  \nrelatively well-constrained measurements of the stellar mass, galaxy cold gas (H/i.pc and H 2 ) content, and star formation rates. 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 review see Rubin 2017); and in a number of studies where it is possible to use a fortuitouslyaligned background quasar to probe circumgalactic gas flows (e.g. Lehner et al. 2013; Bouché et al. 2016; Martin et al. 2019). Even with these discoveries, the rate at which the gas accretes and its associated properties have yet to be well-quantified on a statistical basis; and are associated with considerable uncertainties.  \nThere is a significant body of literature presenting supporting evidence for the presence of feedback-driven gas outflows from galaxies (e.g., Heckman et al. 2000; Martin 2005; Veilleux et al. 2005; Feruglio et al. 2010; Newman et al. 2012; Rubin et al. 2014; Heckman et al. 2015; Schroetter et al. 2016; Chisholm et al. 2016; Rupke et al. 2019; for a recent review see Veilleux et al. 2020). However, 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 gas rich galaxies exhibiting particularly strong outflows - either those at high redshift, or local galaxies experiencing a starburst.  \nModern 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. 2010; Dubois et al. 2014; Vogelsberger et al. 2014; Christensen et al. 2016; Schaye et al. 2015; Pillepich et al. 2018; Hopkins et al. 2018; Davé et al. 2019). However, 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 "sub-grid" techniques, which effectively account for their influence on galaxy evolution in a physically-motivated, but necessarily coarse-grained manner (for reviews, see Somerville & Davé 2015; Vogelsberger et al. 2020). A wide range of sub-grid 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 §3).  \nDespite these promising results, recent studies have suggested that the behavior 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 over cosmic time (Naab & Ostriker 2017; Mitchell et al. 2018, 2020; Pandya et al. 2020; Davé et al. 2020; VillaescusaNavarro et al. 2021; Ni et al. 2023; Tillman et al. 2023a). For example, 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 overview of the sub-grid prescriptions in each simulation, see §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) over cosmic time, within a range of ≲ 0 . 3 dex. This phase is strongly correlated with star formation, which also agrees fairly well between  \nmodels as also established by the aforementioned stellar mass function. Predictions for the slightly more tenuous atomic phase, however, show greater divergence.  \nSimilarly, the emergent gas mass outflow rates and mass loading factors ( 𝜂 ≡ / 𝑀 out ( 𝑟 )/ / 𝑀 ★ ) at ISM scales that arise from different sub-grid implementations of stellar-driven winds appear to differ significantly even across simulations that are all calibrated to reproduce the observed 𝑧 = 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 entirely - even 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.  \nTo 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-driven outflows 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 move towards a more constrained understanding of the role of feedback processes in the baryon cycle.  \nThis paper is arranged as follows: in §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 §3, we outline the similarities and differences between statically measured baryon content and baryon distribution in haloes in these simulations. In §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 §3. In §5, we discuss the implications of our findings forthe baryon cycle, and discuss future targeted observational tests for different feedback models. Finally, we conclude with §6, where we summarise the main findings of the paper.', '2 METHODS': "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 (cold dark matter) cosmology - in which hierarchically assembled cold dark matter haloes provide the formation site of galaxies; all within the context of an expanding, dark energy dominated universe. As noted above, all of these simulations contain phenomenological sub-grid 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, otherwise leaving the behaviour of the gaseous phase as  \na 'prediction' of the models (see Crain & van de Voort 2023 for an overview). 1 .  \nIn §2.1 we introduce the simulations we employ for this study and describe their respective sub-grid feedback implementations. Subsequently, in §2.2 we outline the method with which we selfconsistently measure gas flows rates in each of the simulations.", '2.1 Cosmological 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 §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 Λ CDMcosmology, with the choice of cosmological parameters differing by only of order ≈ 1%. While the simulations take slightly different approaches to solving the equations describing the hydrodynamics of the gas, they have many characteristics in common: each has a spatial resolution of order 1 -10 kpc, make use of very similar tabulated gas cooling prescriptions, and utilise an effective 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 modeling the unresolved feedback from star formation (stellar feedback) and SMBHs (AGN feedback), which we focus on below.', '2.1.1 The EAGLE simulations': 'The EAGLE simulations, as described by Schaye et al. (2015) and Crain et al. (2015), utilised a modified version of the N-body TreeParticleMesh (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 sub-grid physics modules. The initial conditions for this simulation were generated using the method described in Jenkins 2013, assuming the cosmological parameters from Planck Collaboration et al. (2014): Ω Λ = 0 . 693, Ω m = 0 . 307, Ω 𝑏 = 0 . 04825, 𝐻 0 = 67 . 77kms -1 Mpc -1 , 𝜎 8 = 0 . 8288, and 𝑛 s = 0 . 9611.  \nThe flagship box of the suite, referred to as Ref-L100N1504, has a periodic side length of 67 . 8 ℎ -1 Mpc. It was executed with 1 , 504 3 dark matter and baryonic particles, resulting in particle masses of 1 . 8 × 10 6 M ⊙ and 9 . 7 × 10 6 M ⊙ for baryons and dark matter respectively. The Plummer-equivalent gravitational softening length is set to 2 . 66ckpc, limited to a maximum proper length of 0 . 7pkpc. 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 ℎ -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 𝑧 ≈ 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.  \nPhotoheating 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 UV and X-ray background described  \n1 An exception to this is the TNG AGN feedback model, which is calibrated to reproduce 𝑧 ≈ 0 halo gas fractions at high mass, and SFR densities over cosmic time. Even so, there is considerable degeneracy in how the exact interplay between gas flows and the baryon cycle can produce these results.  \nby Haardt & Madau (2001). Metals are not passively diffused between gas particles in proximity, however cooling calculations are conducted using SPH-kernel weighted metallicities. Since the EAGLEsimulations do not have the resolution to model cold, interstellar gas, a density-dependent temperature floor is imposed (normalised to 𝑇 = 8 000 K at 𝑛 H = 10 -1 cm -3 ). To model star formation, a metallicity-dependent density threshold is set, above which star formation is locally permitted (Schaye 2004), and gas particles meeting this threshold are converted to star particles stochastically. The star formation rate is set by a tuned pressure law (Schaye et al. 2007), calibrated to the observations of Kennicutt (1983) at 𝑧 = 0.  \nThe stellar feedback sub-grid model in EAGLE accounts for energy deposition into the ISM from radiation, stellar winds, and supernova explosions (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, Δ 𝑇 SF = 10 7 . 5 K. The average energy injection rate from young stars is given by 𝑓 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 supernova event. The value of 𝑓 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 𝑧 ≈ 0 galaxy stellar mass function and 𝑀 ★ -𝑅 50 relation (see Crain et al. 2015).  \nSMBHs of mass 10 5 ℎ -1 M ⊙ are seeded when a halo exceeds a virial mass of 10 10 ℎ -1 M ⊙ . Subsequently, SMBHs can grow via Eddington-limited accretion (utilising a Bondi-Hoyle parameterisation), as well as mergers with other SMBHs (Bondi 1952; Springel et al. 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 Δ 𝑇 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 fixed conversion efficiency 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 avoid rapid cooling and dissipation of the feedback energy. The exact parameter choices in the EAGLE AGN feedback sub-grid model are informed by observations of the 𝑀 ★ -𝑀 BH relation (Crain et al. 2015).', '2.1.2 The IllustrisTNG simulations': '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 et al. 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 et al. (2016) cosmology, with Ω Λ = 0 . 692, Ω m = 0 . 31, Ω 𝑏 = 0 . 0486, 𝐻 0 = 67 . 7kms -1 Mpc -1 , 𝜎 8 = 0 . 8159, and 𝑛 s = 0 . 97.  \nThe 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 ℎ -1 Mpc to  \nTable 1. Description of the simulations used for this study: EAGLE, TNG, and SIMBA. Key parameters include the simulation box size, resolution, hydrodynamics scheme, gas cooling & star formation prescription, and the implementation of stellar and AGN feedback.  \nbalance 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 dark matter particle mass of 7 . 5 × 10 6 M ⊙ .  \nFollowing 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 AGB 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 et al. 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 et al. (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).  \nStellar-driven winds are modeled by isotropically injecting thermal and kinetic energy into wind particles, which are temporarily decoupled from hydrodynamic forces. With 10% of the feedback energy injected thermally, the kinetic component is parameterised by following prescriptions for mass loading and injection velocity that depend on redshift, gas phase metallicity and the velocity dispersion  \nwithin a weighted kernel over the 64 nearest dark matter 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% of the density threshold for star formation) or reaches a maximum travel time (2.5% 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).  \nThe 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 ℎ -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.  \nThe calibration process for the TNG model, discussed in Pillepich et al. (2018), uses various observations such as the stellar mass function, cosmic star formation rate density as a function of redshift, BH mass versus stellar mass relation at 𝑧 = 0, hot gas fraction in galaxy clusters at 𝑧 = 0, and the galaxy mass-size relation at 𝑧 = 0.  \nSIMBA (Davé et al. 2019) is the next generation of the M/u.pc/f.pc/a.pc/s.pc/a.pc cosmological galaxy formation simulations (Davé et al. 2016), and use the meshless finite mass (MFM) mode of the GIZMO hydrodynamics code (Hopkins 2015, 2017), with gravity solver based on GADGET-3 (Springel 2005). The SIMBA simulations adopt cosmological parameters Ω Λ = 0 . 70, Ω m = 0 . 30, Ω 𝑏 = 0 . 048, 𝐻 0 = 68kms -1 Mpc -1 , 𝜎 8 = 0 . 82, and 𝑛 s = 0 . 97.  \nThe flagship SIMBA simulation box has a side length of 100 ℎ -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 ℎ -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).  \nThe 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 locked away in the process of dust formation. Radiative cooling and photoionisation 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-formation driven outflows are implemented as decoupled two-phase winds, characterized by updated mass-loading factors derived from particle tracking in the Feedback in Realistic Environments (FIRE) zoom-in simulations (Anglés-Alcázar et al. 2017b). Wind velocities are derived from Muratov et al. (2015). 30% 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% of the local sound speed. The exact parameter choices in the SIMBA stellar feedback sub-grid model are informed by observations of the 𝑧 ≈ 0 stellar mass function.  \nThe 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-HoyleLyttleton prescription Bondi (1952), capped at the Eddington limit. A unique feature of the Anglés-Alcázar et al. (2017a) torque-based accretion model is that there is no need to calibrate the AGN feedback efficiency to reproduce the 𝑀 ★ -𝑀 BH relation, rather, only an assumption about the accretion efficiency is required.  \nSIMBA also incorporates feedback from AGN through radiative and jet modes. The radiative mode drives winds with velocities from 500 -1500 kms -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 𝑣 wind = 1000 km / s for 𝑀 BH = 10 9 M ⊙ ), and the jet-mode can add a velocity boost of up to 7000 km / s, 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 kiloparsecs before the energy is deposited into the surrounding gas. Gas in the jets is injected at the halo virial temperature.', '2.2 Measuring gas flows': "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  \nmass of gas elements that are found to have crossed a boundary between two outputs, 𝑡 𝑖 and 𝑡 f , is summed and normalised by Δ 𝑡 = 𝑡 f -𝑡 i (e.g. Mitchell et al. 2020; Wright et al. 2020). This method is simple to implement in codes where the gas is discretised 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).  \nThe 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).  \nIn 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 friendsof-friends (FOF) method combined with (i) SUBFIND in the case of EAGLE and TNG (Springel et al. 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 center of potential (which, by definition, is centered on the central galaxy), velocity, and halo mass/radius based on a spherical overdensity of 200 × 𝜌 𝑐 ( 𝑀 200c and 𝑅 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.  \nAt each of the spherical apertures of interest centered at radius 𝑟 = 𝑅 , we identify boundary gas elements as those residing in a spherical shell of Δ 𝑟 = 0 . 2 × 𝑅 , or between 𝑟 = 0 . 9 × 𝑅 and 𝑟 = 1 . 1 × 𝑅 . We see very little quantitative differences between gas flow rates measured for reasonable choices of Δ 𝑟 (0 . 05 -0 . 5 × 𝑅 , or fixed Δ 𝑟 = 1 -10 kpc). These boundary elements are categorised as being either outflow or inflow depending on the sign of their radial velocity, where the radial velocity of a gas element 𝑖 relative to a halo center 𝑗 can be calculated as 𝑣 r , 𝑖 𝑗 = fi 𝑣 𝑖 𝑗 · fi 𝑟 𝑖 𝑗 /| fi 𝑟 𝑖 𝑗 | . The gas flow rates at shell 𝑟 = 𝑅 around halo 𝑗 for each of these subsets of boundary gas elements, 𝑖 ∈ 𝑘 , can be calculated with Equation 1 as follows:  \n/ 𝑀 𝑘 ( 𝑟 = 𝑅 ) = 1 Δ 𝑟 × ∑︁ 𝑖 ∈ 𝑘 GLYPH<18> 𝑚 𝑖 fi 𝑣 𝑖 𝑗 · fi 𝑟 𝑖 𝑗 | fi 𝑟 𝑖 𝑗 | GLYPH<19> , (1)  \nwhere / 𝑀 𝑘 ( 𝑟 = 𝑅 ) is the mass flow rate across the spherical boundary at 𝑟 = 𝑅 . 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).", '3 THE BARYON BUDGET BETWEEN SIMULATIONS': "In this section, before presenting our analysis of gas flows between the simulations in §2.2, we compare static properties of the baryons within haloes in EAGLE, TNG, and SIMBA at 𝑧 = 0. Here, we compare our measurements the baryon content of haloes between simulations, and its breakdown between different components. We  \nFigure 1. Median stellar to halo mass fractions as a function of 𝑀 200c mass in EAGLE (blue), TNG (purple), and SIMBA (pink) at 𝑧 = 0. The inferred relation presented in Moster et al. (2018) ± 1 𝜎 from their fit is included for reference as a grey shaded region. In each of the simulations, the hatched shaded regions represent the 90% confidence interval on the median from 100 bootstrap re-samplings in each bin. There is general agreement between the simulations and the Moster et al. (2018) relation, despite the simulations each utilising very different sub-grid physics.  \n<!-- image -->  \nnote that we reserve a detailed discussion of the link between these findings and previous studies for §5.1.  \nFig. 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 galaxy properties to simulated dark matter 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 normalising by the halo mass within the 𝑅 200c radius, 𝑀 200c . We convert the virial masses used as the x-values in Moster et al. (2018) from the Bryan & Norman (1998) overdensity criterion to the 𝑀 200c definition we use here with /h.pc/m.pc/f.pc/c.pc/a.pc/l.pc/c.pc (Murray et al. 2013), however note that we do not adjust the Moster et al. (2018) 𝑓 ★ 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 above, and discussed in (Somerville & Davé 2015), the relatively good agreement of the different simulations in spite of the very different implementation of sub-grid physics is largely due to the fact that a high weight is placed on reproducing these observational constraints when calibrating the sub-grid parameters.  \nIn Fig. 2, we compare the total baryon content of haloes at 𝑧 = 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  \n' 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 𝑅 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 𝑅 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 𝑅 200c , as opposed to the central 30kpc. This means that in Fig. 2, the stellar component also includes contributions from intra-halo stars and satellite galaxies within 𝑅 200c . The 'cold gas' componentisnotdecomposedfurtheror post-processed to provide H/i.pc 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).  \nOne of the striking differences between the simulations in terms of halo baryon fraction is their behaviour in the range log 10 ( 𝑀 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 turnover mass which we denote as M T 1 . Above M T 1 , halo baryon fractions decrease with halo mass until a second turnover mass, denoted by M T 2 , is reached. Above M T 2 , halo baryon fractions steadily rise again to the highest halo masses presented here at log 10 ( 𝑀 200c / M ⊙ ) = 14. Physically, these transitions produce three distinct 'regimes' in terms of the 𝑀 200c -𝑓 b relation: (i) subM 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 turnover. 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 turnover/flattening behaviour respectively between M T 1 -M T 2 .  \nAt halo masses below log 10 ( 𝑀 200c / M ⊙ ) ≈ 11 . 5, corresponding to subM 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 ( 𝑀 200c / M ⊙ ) ≈ 11, the L100-REF model predicts a median halo baryon fraction of ≈ 15% 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%, though still lower than the corresponding baryon fraction predictions from SIMBA at ≈ 40% of the expected cosmological baryon fraction. TNG predicts the highest halo baryon fractions at this mass, at ≈ 60% of the universal value (see also Davies et al. 2019a; Davies et al. 2019b).  \nIn this subM 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 overall baryon fractions mean that in order to match the expected stellar content of haloes, stars make up over 25% of the baryonic content, compared to SIMBA  \nFigure 2. The baryon content of haloes at 𝑧 ≈ 0 in EAGLE (blue), TNG (purple) and SIMBA (pink). Top panel, 1 st column : The median total baryon fraction within 𝑅 200c (stars, gas, black holes) of haloes as a function of 𝑀 200c mass. Bottom panel, 1 st column : The gas fraction of haloes within 𝑅 500c as a function of 𝑀 500c , shown in comparison to weak lensing observations from Akino et al. (2022) (diamond point); and X-ray data compiled in Kugel et al. (2023) (square points). Top panels, 2 nd -4 th columns : Total baryon content for each simulation within 𝑅 200c broken down into stars (hatched bars), 'cold' gas (coloured shaded bars), and 'hot gas' (grey bars) in each simulation. In the top panels, this baryon content is normalised by 𝑓 b × 𝑀 200c . Bottom panels, 2 nd -4 th columns : These panels shows the same breakdown, but re-normalised by the actual baryonic mass in each halo, 𝑀 bar ( 𝑟 < 𝑅 200c ) . In each of the simulations in the left-hand panel, error-bars represent the 16 th -84 th percentiles in halo baryon fractions for a given mass bin; the hatched regions represent the 90% confidence interval on the median from 100 bootstrap re-samplings in each bin. While the simulations produce similar results for the stellar content of haloes at 𝑧 ≈ 0, the same cannot be said for the gaseous phase.  \n<!-- image -->  \nand TNG where the stellar component only represents ≲ 10% of the baryonic matter. For halo masses below log 10 ( 𝑀 200c / M ⊙ ) ≈ 11 . 5, we note that the EAGLE L100-REF run predicts very low halo cold gas content - always 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 mass range. Davé et al. (2020) show that the EAGLE fiducial run underpredicts the observed 𝑧 = 0 H/i.pcMF from ALFALFA 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 𝑧 = 0, however 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.pcMF; and that for 𝑧 ≥ 1, that the EAGLE L100-REF and L25-RECAL models are consistent with regard to predicted H/i.pc and H 2 mass functions.  \nWhile SIMBA and TNG predict a similar functional form for baryon fractions across halo mass, their normalisation 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%-50% higher than those measured in SIMBA. Furthermore,  \nat log 10 ( 𝑀 200c / M ⊙ ) = 14, baryon fractions in SIMBA only reach ≈ 70% of Ω b / Ω m , while haloes in TNG haloes approach this universal value very closely. Instead of a turnover, EAGLE 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. Above M T 2 , EAGLE and TNG agree very well in terms of predicted baryon fractions as a function of halo mass.  \nIn 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 datasets are nominally cast as (hot) gas fractions within 𝑅 500c , as a function of 𝑀 500c , representing a slightly shrunken aperture relative to the rest of our measurements. In the bottom left panel of Figure 2, we present like-for-like measurements of hot halo gas content in each of the simulations within 𝑅 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 dataset, noting that we only make use of the WL data from Akino et al. (2022), and that the relevant X-ray datasets 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 et al. (2015); Pearson et al. (2017); Lovisari et al. (2020). As above, we define 'hot' gas as any gas elements with a temperature 𝑇 ≥ 5 × 10 4 . We note that this choice  \nof threshold does not significantly influence gas fractions in this mass range, as the majority of gas is at or above halo virial temperature (see top panels). Regardless, we stress that this selection is unlikely to exactly mimic the gas accounted for in the observations, however allows us to understand the difference between simulation predictions and their broad agreement with observations. Our measurements between 𝑀 500c = 10 13 . 5 M ⊙ -𝑀 500c = 10 14 . 5 M ⊙ indicate that median halo gas fractions in TNG and EAGLE are slightly higher than the compiled observations (more so in EAGLE), and slightly lower than observations in SIMBA. The range of predictions from the simulations is largest between 𝑀 500c = 10 13 M ⊙ -𝑀 500c = 10 14 M ⊙ , while at higher masses ( 𝑀 500c ≈ 10 14 . 5 M ⊙ ), the gas fractions tend towards a tighter range. For reference, we note that the offset between 𝑀 200c and 𝑀 500c halo mass measurements is relatively constant at ≈ 0 . 15 dex in this mass range. While the X-ray datasets 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.  \nIn Appendix D, we show a breakdown of the baryon content of haloes at 𝑧 = 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 𝑧 = 0, with a universal shift to higher halo baryon fractions. The transition masses M T 1 and M T 2 remain similar in EAGLE and TNG at 𝑧 = 2 compared to 𝑧 ≈ 0 galaxies; though this is not the case for SIMBA. In SIMBA, the transition masses are higher than at 𝑧 = 0, with log 10 (M T 1 / M ⊙ ) = 12 . 2 and log 10 (M T 2 / M ⊙ ) = 12 . 8 at 𝑧 = 2 (see also Fig. 7 in Sorini et al. 2022).  \nToaidintheinterpretations of the gas flow measurements in the following 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.", '4 GAS FLOWS THROUGH THE ISM, CGM, AND IGM': '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 §3. Our method for measuring gas flow rates is detailed in §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.', '4.1 Ejective feedback: gas outflow rates as a function of scale': "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 star formation rate. We begin by comparing the mass loading results between simulations at 𝑧 ≈ 2, which builds a picture for the operation of gas flows at the epoch of maximum star formation rate density (Madau & Dickinson 2014). In the bottom panel of Fig. 3, we also compare the median specific star formation rate ( / 𝑀 ★ / 𝑀 ★ ) of galaxies. This assists in understanding the influence of star formation rate 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 specific star formation rate drops below log 10 ( sSFR / yr -1 ) = -11 + 0 . 5 𝑧 (Furlong et al. 2015).  \nWe define the stellar mass of a galaxy in the same way as Fig.  \n1, including all star particles within 30 kpc of the halo centre of potential. Similarly, the star formation rates we use in calculating mass loading are evaluated with gas within 30kpc of the galaxy 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 𝑟 = 0 . 25 × 𝑅 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).  \nWe 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 star formation rate in the bin. Neistein et al. (2012) demonstrate that this method correctly predicts the average rate of mass exchange when integrated over 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 of Fig. 3 correspond to the range of mass loading values predicted based on the 16 th and 84 th percentiles of mass outflow rate in a given stellar mass bin, without including variance in star formation rates at a given mass.  \nFor reference, in Fig. 3, we also include indicative 𝜂 scalings with 𝑉 c ≡ √︁ 𝐺𝑀 200c / 𝑅 200c (halo circular velocity): a 'momentum conserving' scenario with 𝜂 ∝ 𝑉 -1 c ∝ 𝑀 -1 / 3 200c ; and an 'energy conserving' scenario with 𝜂 ∝ 𝑉 -2 c ∝ 𝑀 -2 / 3 200c as dashed and dash-dotted lines respectively (see Murray et al. 2005 for derivation).  \nAt the ISM-scale of 𝑟 = 0 . 25 𝑅 200c in the top panel of Fig. 3, it is clear that each of the simulations predicts a similar negative scaling in 𝜂 = / 𝑀 out / SFR as a function of 𝑀 ★ , until a turnover mass where outflows begin to pick up relative to star formation rates (also demonstrated in Nelson et al. 2019 and Mitchell et al. 2020 for TNG and EAGLE respectively, though using different scales and methodologies). The turnover mass differs slightly between simulations in EAGLE and TNG, the turnover occurs at 𝑀 ★ ≈ 10 10 M ⊙ and 𝑀 ★ ≈ 10 10 . 3 M ⊙ respectively; with 𝜂 values at the turnover 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 𝑀 ★ ≈ 10 11 M ⊙ . Folding in the stellar-halo mass relation in each of the simulations, in EAGLE and TNG the turnover mass occurs at log 10 ( 𝑀 200c / M ⊙ ) ≈ 11 . 8 -12, while in SIMBA, the turnover occurs at a slightly higher mass of log 10 ( 𝑀 200c / M ⊙ ) ≈ 12 . 4 at 𝑧 ≈ 2 (as is also the case for halo baryon content).  \nBelow the turnover mass, the scaling of 𝜂 with mass for each of the simulations falls between the expected 'momentum conserving' ( 𝜂 ∝ 𝑀 -1 / 3 200c ) and 'energy conserving' ( 𝜂 ∝ 𝑀 -2 / 3 200c ) scaling relations. The relation is steepest in SIMBA, for which 𝜂 approximately follows 𝑀 -2 / 3 200c , and is closer to 𝑀 -1 / 2 200c in EAGLE and TNG. The mass  \nFigure 3. Mass loading ( 𝜂 ) and specific star formation rates (sSFR) at 𝑧 ≈ 2 as a function of stellar mass in EAGLE (blue), TNG (purple) and SIMBA (pink). Top panel : averaged ISM-scale ( 𝑟 = 0 . 25 𝑅 200c ) mass loading (see text for definition) as a function of stellar mass. Bottom panel : mean specific star formation rate as a function of stellar mass in each of the simulations. In the top panel, the error-bars displayed correspond to the range of measured mass loading values based on the 16 th and 84 th percentiles of mass outflow rate in a given stellar mass bin. In the bottom panel, the error-bars correspond to the 16 th -84 th percentile range of sSFR at a given mass. In both panels, coloured hatched regions correspond to the bootstrap-generated 90% confidence interval on the respective averages at a given mass. 𝜂 scales negatively with 𝑀 ★ for all simulations below 𝑀 ★ ≈ 10 10 M ⊙ , however the 𝜂 values in this mass range can vary between simulations by ≳ 0 . 5 dex on average.  \n<!-- image -->  \noutflow rate from the ISM per unit star formation rate in this mass range is highest in SIMBA, which predicts 𝜂 = 10 1 . 6 at 𝑀 ★ ≈ 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.  \nThe 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 × 𝑅 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  \nwith 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 ( Δ 𝑇 SF = 10 7 . 5 K) and associated pressure produce strong outflows.  \nAbove the mass turnover 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 ( 𝑀 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.  \nComparing the median specific star formation rates 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, however, small quantitative differences between the simulations - for instance, the median specific star formation rate in SIMBA at 𝑀 ★ ≈ 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 𝑀 ★ ≈ 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 star formation rate. 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.  \nIn Fig. 4 and 5, we show the raw outflow rates for central galaxies in each simulation as a function of halo mass, normalised by Ω b / Ω m × 𝑀 200c at 𝑧 = 2 and 𝑧 = 0 respectively ('outflow efficiency'). In the left hand panels, we show these outflow rates at the ISM scale (0 . 25 × 𝑅 200c ). In the remaining panels for each simulation, we show outflow rates at the ISM scale, the halo scale (defined at 1 . 00 × 𝑅 200c ), and at the IGM scale (defined at 2 . 50 × 𝑅 200c ).  \nThe values of outflow efficiency indicate the inverse of the timescale that would be required for a gas mass of 𝑀 200c × Ω b / Ω m to be removed at this scale. In the left-hand panels of Fig. 4 and 5, we compare median outflow rates in each of the simulations at the ISM scale (0 . 25 × 𝑅 200c ). In each of the remaining panels, for each of the simulations we compare outflow rates at three scales: 𝑟 = 0 . 25 × 𝑅 200c , 𝑟 = 1 . 00 × 𝑅 200c , and 𝑟 = 2 . 50 × 𝑅 200c . Firstly, focusing on the left-hand panels which compare the simulations at the ISM scale, we note that at both redshifts, below log 10 ( 𝑀 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 normalisation of the 𝑀 ★ -sSFR relation in SIMBA compared to TNG (see Fig. 3 for sSFR measurements at 𝑧 ≈ 2 and Davé et al. 2019 for measurements  \nFigure 4. Mass outflow rates as a function of 𝑀 200c at 𝑧 ≈ 2 in EAGLE (blue), TNG (purple) and SIMBA (pink). Mass outflow rates are normalised by 𝑓 b × 𝑀 200c in each case to indicate gas outflow 'efficiency' across different mass scales. Left panel : Outflow efficiency at the ISM scale, compared between the simulations. Panels 2-4 : Outflow efficiencies at different scales in EAGLE, TNG, and SIMBA respectively. In each of these panels, ISM-scale outflows are plotted with lighter colours, halo-scale outflowsare plotted in medium colours, and larger-scale outflows are plotted in darker colours. In all panels, error-bars correspond to the 16 th -84 th percentile range in outflow rates at a given mass, and hatched regions correspond to the bootstrap-generated 90% confidence interval on the medians at a given mass. SF-driven outflows reach the greatest scale in EAGLE, while AGN-driven outflows reach the greatest scale in SIMBA.  \n<!-- image -->  \nFigure 5. Mass outflow rates as a function of 𝑀 200c at 𝑧 ≈ 0 in EAGLE (blue), TNG (purple) and SIMBA (pink). Mass outflow rates are normalised by 𝑓 b × 𝑀 200c in each case to indicate gas outflow 'efficiency' across different mass scales. Left panel : Outflow efficiency at the ISM scale, compared between the simulations. Panels 2-4 : Outflow efficiencies at different scales in EAGLE, TNG, and SIMBA respectively. In each of these panels, ISM-scale outflows are plotted with lighter colours, halo-scale outflowsare plotted in medium colours, and larger-scale outflows are plotted in darker colours. In all panels, error-bars correspond to the 16 th -84 th percentile range in outflow rates at a given mass, and hatched regions correspond to the bootstrap-generated 90% confidence interval on the medians at a given mass. EAGLE predicts efficient SF-driven ejection of gas from haloes, while TNG and SIMBA predict an intra-halo 'recycling' scenario for SF-affected gas. AGN-driven outflows reach the greatest scale in SIMBA.  \n<!-- image -->  \nat 𝑧 ≈ 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.  \nIt is when we analyse the behaviour of outflow rates as a function of scale as per the right-hand panels in Fig. 4 and 5, that dramatic differences between gas flows in the simulations begin to emerge. We first focus on the 𝑧 ≈ 2 case, as presented in Fig. 4. Below log 10 ( 𝑀 200c / M ⊙ ) ≈ 12 in EAGLE (where stellar feedback is dominant in driving outflows), mass outflow rates at the ISM (0 . 25 × 𝑅 200c ), halo (1 × 𝑅 200c ), and IGM (2 . 5 × 𝑅 200c ) scales remain very similar. This indicates that the galactic winds in EAGLE can travel to at least 2 . 5 × 𝑅 200c at 𝑧 ≈ 2 before stalling and falling back.  \nUnlike the picture in EAGLE, below log 10 ( 𝑀 200c / M ⊙ ) ≈ 12 outflow rates decrease with increasing scale in TNG. At log 10 ( 𝑀 200c / M ⊙ ) ≈ 11 . 5 -12 in TNG, outflow rates at 𝑟 = 2 . 5 × 𝑅 200c are just 10% of those measured at 𝑟 = 0 . 25 × 𝑅 200c .  \nWe do note, however, that this stalling is less obvious in lower-mass systems with log 10 ( 𝑀 200c / M ⊙ ) ≈ 11 - at this mass, average outflow rates at 𝑅 200c are ≈ 50% of that measured at the ISM-scale. In general, the stalling of SF-driven outflows in TNG is related to the velocities imparted to wind particles at injection, which scale with 𝜎 DM , which is almost always below the escape velocity of the halo (Nelson et al. 2019). This means that SF-driven outflows in TNG are ultimately confined to recycle within the halo, except for perhaps in the context of very low-mass haloes ( ≲ 10 10 . 5 M ⊙ ). In SIMBA, the picture for SF-driven outflows is somewhere between EAGLE and TNG. At low halo masses ( 𝑀 200c ≈ 10 10 . 5 M ⊙ ) outflows are able to reach a similar scale as predicted in EAGLE (2 -3 × 𝑅 200c ), while towards larger masses ( 𝑀 200c ≈ 10 12 M ⊙ ) outflows largely stall before reaching 𝑅 200c . This result is again a direct consequence of the assumed wind velocities at injection, which in SIMBA are taken from Muratov et al. (2015).  \nAbove log 10 ( 𝑀 200c / M ⊙ ) ≈ 12 . 5 at 𝑧 ≈ 2, where AGN feedback  \nbegins 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 × 𝑅 200c scale, with / 𝑀 out at 2 . 5 × 𝑅 200c approximately one tenth of the value at 𝑅 200c in halos with mass log 10 ( 𝑀 200c / M ⊙ ) ≈ 13. This indicates that AGN-driven outflows in EAGLE tend to stall near or just outside the virial radius at 𝑧 ≈ 2, and seldom reach larger scales. TNG shows qualitatively similar behaviour, but 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 ( 𝑀 200c / M ⊙ ) ≈ 11 . 5 -13, and a constant outflow rate from 𝑅 200c to 2 . 5 × 𝑅 200c at higher halo masses.  \nThis 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.  \nMoving to the 𝑧 ≈ 0 Universe in Fig. 5, we note that the qualitative differences between simulations remain similar to the 𝑧 ≈ 2 picture at the ISM scale, however the disparities between how the outflows propagate across scales become even more pronounced. In the lefthand panel, we can see that the normalisation of gas flow rates relative to halo mass are much lower than at 𝑧 ≈ 2, with longer dynamical time-scales and lower gas densities at this epoch. One exception to this is the strength of AGN-driven outflows in SIMBA - which are much stronger at 𝑧 ≈ 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).  \nOne very interesting feature at 𝑧 ≈ 0 in the EAGLE simulations is the apparent entrainment of outflowing gas with increasing scale. We find that at log 10 ( 𝑀 200c / M ⊙ ) ≈ 11, the mass outflow rate at 2 . 5 × 𝑅 200c is over 10 times higher than that at 0 . 25 × 𝑅 200c in EAGLEat 𝑧 ≈ 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 out-flowing gas is significantly over-pressurised relative to the ambient CGM at 𝑧 = 0. This causes the tenuous material in the CGM to be swept up with the bulk flow.  \nAs demonstrated in Fig. 4, in the regime dominated by stellar feedback, 𝑀 halo ≲ 10 12 M ⊙ , we find that in EAGLE there is minimal mass entrainment in outflows at 𝑧 ≈ 2 compared to 𝑧 ≈ 0 (with the exception of haloes below 10 11 M ⊙ , where outflows at 2 . 5 × 𝑅 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 ( 𝑀 200c / M ⊙ ) ≈ 11 -11 . 5 at 𝑧 ≈ 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.  \nThis 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 ( 𝑀 200c / M ⊙ ) ≈ 11 at 𝑧 ≈ 2 and 𝑧 ≈ 0 in Fig. B1 and B2, we find that in all simulations, gas densities in the CGM universally decrease from 𝑧 ≈ 2 to 𝑧 ≈ 0. In TNG and SIMBA, the difference between gas densities at 𝑧 ≈ 2 and 𝑧 ≈ 0 is ≈ 1 . 3 -1 . 5 dex at 𝑟 = 0 . 50 × 𝑅 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, however, the decrease in CGM gas density is much more dramatic, closer to ≈ 2 -2 . 5 dex from 𝑧 ≈ 2 to 𝑧 ≈ 0. This means that at 𝑧 ≈ 0, outflows are over-pressurised in a relative sense to the CGM (as per the findings in Mitchell et al. 2020). Comparing this picture  \nto 𝑧 ≈ 2, where CGM densities between all simulations are similar, the lack of mass entrainment at 𝑧 ≈ 2 in EAGLE is likely the result of a much denser CGM at this redshift compared to 𝑧 ≈ 0, meaning that outflows are not over-pressurised relative to the surrounding gas, and not able to sweep up the surrounding material with the bulk flow.  \nAn 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 𝑧 ≈ 2 (Fig. C1), in the log 10 ( 𝑀 200c / M ⊙ ) ≈ 11 and log 10 ( 𝑀 200c / M ⊙ ) ≈ 12 bins, median CGM gas temperatures in EAGLE are substantially lower than 𝑇 vir (estimated from van de Voort 2017); particularly where 𝑟 / 𝑅 200c ≲ 0 . 75 (the 'inner' CGM). This aligns with the commonly accepted picture that 'cold-mode' accretion - cool, 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 𝑧 ≈ 0 in EAGLE, gas temperatures in the CGM are significantly higher - within 0 . 2 -0 . 3 dex of 𝑇 vir at log 10 ( 𝑀 200c / M ⊙ ) ≈ 12. The low covering 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.  \nAs discussed above, we find that in TNG, outflow rates monotonically decrease with increasing scale below log 10 ( 𝑀 200c / M ⊙ ) ≈ 12 at both 𝑧 ≈ 0 and 𝑧 ≈ 2, indicating stalling of outflows as they propagate to larger distances. Interestingly, at 𝑧 ≈ 2, Fig. C1 shows that for TNG haloes in bins log 10 ( 𝑀 200c / M ⊙ ) ≈ 11 and log 10 ( 𝑀 200c / M ⊙ ) ≈ 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 thermalised. Unlike EAGLE, we find that gas densities in the CGM remain relatively high in this mass range at 𝑧 ≈ 0 in TNG, which in addition to the relatively small thermal component of stellar feedback in TNG (10% of feedback energy), contributes to the outflowing material being under-pressurised relative to the surrounding medium, and thus unable to reach larger scales or entrain mass.  \nConsidering larger haloes at 𝑧 ≈ 0, an interesting feature in SIMBA is the efficient entrainment of AGN-driven gas outflows above log 10 ( 𝑀 200c / M ⊙ ) ≈ 12. At log 10 ( 𝑀 200c / M ⊙ ) ≈ 12 . 5, the gas outflow rate at 2 . 5 × 𝑅 200c is roughly 5-10 times the ouflow rate at 𝑟 = 0 . 25 × 𝑅 200c - an effect not seen in EAGLE, nor TNG. At 𝑧 ≈ 0, the ambient density of the CGM at this mass scale is quite low (as demonstrated in Fig. 2 and B2). This, together with the direct heating of gas elements surrounding AGN, means that the AGN-driven outflows are over-pressurised relative to the surrounding medium, and thus able entrain tenuous gas in both the CGM, and the region surrounding the halo.", '4.2 Preventative feedback: gas inflow rates as a function of scale': "In this section, we build on the results of §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 'preventative' perspective; 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).  \nIn Fig. 6 and 7, we show the relationship between inflow rate and  \nFigure 6. Mass inflow rates as a function of 𝑀 200c at 𝑧 ≈ 2 in EAGLE (blue), TNG (purple) and SIMBA (pink). Inflow rates are normalised by 𝑓 b × 𝑀 200c in each case to represent gas inflow 'efficiency' across different mass scales. Left panel : Inflow efficiency at the halo scale, compared between the simulations. Panels 2-4 : Inflow efficiencies at different 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 16 th -84 th percentile range in inflow rates at a given mass, and hatched regions correspond to the bootstrap-generated 90% confidence interval on the medians at a given mass. All simulations predict that AGN feedback has a strong preventative impact on gas accretion at the ISM scale, however there is considerable variability between the preventative impact of stellar feedback.  \n<!-- image -->  \nFigure 7. Mass inflow rates as a function of 𝑀 200c at 𝑧 ≈ 0 in EAGLE (blue), TNG (purple) and SIMBA (pink). Inflow rates are normalised by 𝑓 b × 𝑀 200c in each case to represent gas inflow 'efficiency' across different mass scales. Left panel : Inflow efficiency at the halo scale, compared between the simulations. Panels 2-4 : Inflow efficiencies at different 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 16 th -84 th percentile range in inflow rates at a given mass, and hatched regions correspond to the bootstrap-generated 90% 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.  \n<!-- image -->  \nhalo mass in each of the simulations at 𝑧 ≈ 2 and 𝑧 ≈ 0. In the left hand panel, we plot inflow rates (normalised by 𝑀 200c × Ω b / Ω m ; ' inflow efficiency') at the 𝑅 200c halo boundary as a function of 𝑀 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 𝑀 200c × Ω b / Ω m to be accreted at this scale. In the 2 nd -4 th columns, we compare inflow rates within the same simulations between three scales: 0 . 25 × 𝑅 200c (within the ISM), 𝑅 200c , and at 2 . 5 × 𝑅 200c . This comparison allows us to see how inflow rates change from cosmological (IGM) scales down to the ISM scale in each of the simulations.  \nFirst focusing on the left-hand panel for the 𝑧 ≈ 0 case in Fig. 7, we remark that the respective gas accretion rates between the simulations at 𝑅 200c follow a similar functional form to the halo-wide baryon fractions presented in Fig. 2. In the stellar feedback-dominated regime, log 10 ( 𝑀 200c / M ⊙ ) ≲ 11 . 5, the lowest specific accretion rates are seen in EAGLE, where gas inflow rates are just ≈ 10% of the expectation based on the halo growth rate, due to the far-reaching influence of stellar feedback (as directly demonstrated in §4.1). Accretion rates onto haloes of the same mass are much higher in TNG, with gas accretion rates ≈ 50% 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 𝑅 200c in this mass range, as was demonstrated in Fig. 4 and 5. The predictions from SIMBA regarding gas inflow rates in the stellar-feedback dominated  \nregime fall between those of EAGLE and TNG- where stellar feedback has a stronger preventative effect than in TNGat the halo-scale, but not as dramatic as in the case of EAGLE.  \nIn the AGN-dominated regime, log 10 ( 𝑀 200c / M ⊙ ) ≳ 12 . 5, the simulations predict very different levels of preventative feedback. This is most dramatic in SIMBA at log 10 ( 𝑀 200c / M ⊙ ) ≈ 12 -12 . 5, where the amount of gas inflow to haloes is less than 10% 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% 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 ( 𝑀 200c / M ⊙ ) ≈ 12 . 5, the inflow efficiency drops to ≈ 25% 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 ( 𝑀 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 overcome the deeper potential well (see the decreases in outflow rates at 2 . 5 × 𝑅 200c at this mass scale in Fig. 4 and 5, with the exception of SIMBA at 𝑧 ≈ 0).  \nIn the three right-hand panels of each Figure, we explore how inflow rates in EAGLE, TNG, and SIMBA vary across scales. In EAGLE,at 𝑧 ≈ 0, we clearly see that the rate of gas inflow decreases with radius across all halo masses. Below log 10 ( 𝑀 200c / M ⊙ ) ≈ 11 . 5, the biggest drop in inflow efficiency occurs between 𝑟 = 2 . 5 × 𝑅 200c and 𝑟 = 1 × 𝑅 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 §4.1 regarding the scale at which outflows stall in EAGLE. At log 10 ( 𝑀 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 𝑟 = 0 . 25 𝑅 200c (as also reported in Correa et al. 2018a,b). This agrees with our mass-loading results, which suggest that the strength of AGN-driven outflows at the scale of the ISM are similar between EAGLE and TNG, but that they stall at smaller radius in EAGLE.  \nInspecting inflow rates across different scales in TNG, we find that for haloes with log 10 ( 𝑀 200c / M ⊙ ) ≲ 11 . 5, inflow rates at the ISM scale (0 . 25 × 𝑅 200c ) are similar to inflow rates at the halo boundary at 𝑧 ≈ 2; but at 𝑧 ≈ 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 CGMinTNGinthismassrange(seeFig.B2,partlycontributed 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. Above log 10 ( 𝑀 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).  \nLastly, focusing on inflow rates in SIMBA, we find that below log 10 ( 𝑀 200c / M ⊙ ) ≲ 11 . 5, accretion rates are suppressed equally at 𝑟 = 0 . 25 𝑅 200c , 𝑅 200c , and 2 . 5 𝑅 200c . Based on Fig. 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 ( 𝑀 200c / M ⊙ ) ≳ 12, SIMBA predicts a steadily increasing inflow prevention effect moving inwards from 𝑟 = 2 . 5 × 𝑅 200c to 𝑟 = 0 . 25 × 𝑅 200c . Even though the preventative effect is greatest at 0 . 25 × 𝑅 200c (with gas inflow rates at just 1-5% of those expected based on the halo growth rate at 10 12 . 5 M ⊙ ),  \nSIMBA predicts that even gas at 2 . 5 × 𝑅 200c experiences significant preventative feedback due to AGN. While we do not present the results here, we note that in the SIMBA run with no jet-mode feedback, the preventative effect 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.  \nTo summarise - 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 ( 𝑀 200c / M ⊙ ) ≲ 11 . 5, EAGLE and SIMBA predict that stellar feedback produces a strong preventative effect on halo-scale inflow rates, while in TNG, halo-scale gas inflow is suppressed by a smaller amount. For halo masses log 10 ( 𝑀 200c / M ⊙ ) ≳ 12, SIMBA predicts very strong AGN-driven suppression in halo-scale gas inflow rates, relative to a more mild preventative effect in TNG, and very minimal preventative feedback at the halo scale in EAGLE.", '5 DISCUSSION': 'The EAGLE, TNG, and SIMBA simulations have all been largely successful in producing a reasonable match to the observed galaxy stellar mass function at 𝑧 ≈ 0, albeit with careful tuning required to calibrate their respective sub-grid 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 over cosmic time - can be degenerate in producing galaxies at 𝑧 ≈ 0 with very similar stellar masses and star formation rates. In §5.1, we use our measurements of gas flows rates across the simulations in §4 to make the causal link between the baryon cycle and static halo properties. In §5.2, we compare our results to previous literature investigating the baryon cycle in simulated galaxies. In §5.4 we outline a number of promising observational avenues by which the degeneracy between different feedback models could be broken. Lastly, in §5.5, we discuss how the finite resolution of the simulations used in this study may influence our results.', '5.1 The influence of gas flows on halo baryon content': 'In Figure 8, we summarise the main findings of this study in visual form - indicating the scale reached by feedback-driven outflows 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 ( 𝑀 200c ≲ 10 12 M ⊙ ) and high mass ( 𝑀 200c ≳ 10 12 M ⊙ ) regimes, where stellar and AGN feedback respectively dominate. Our findings in §4 aid in the interpretation of our findings in §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. Firstly, we discuss the low-mass regime, log 10 ( 𝑀 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 𝑧 ≈ 2 and 𝑧 ≈ 0 - can be clearly linked with the reduced halo gas content below log 10 ( 𝑀 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 §4.2.  \nIn contrast, in TNG, stellar feedback driven outflows typically stall before reaching 𝑅 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 𝑧 ≈ 0. This efficient intra-halo recycling scenario means that in low-mass haloes, the baryon content of the  \nFigure 8. Schematic diagram summarising the different gas flow paradigms found in EAGLE, TNG and SIMBA, split into low mass haloes where stellar feedback dominates (left), and higher mass haloes where AGN feedback dominates (right).  \n<!-- image -->  \nCGMin TNG is much higher than that recorded in EAGLE at 𝑧 ≈ 0. SIMBA predicts a picture in between that of EAGLE and TNG, where some gas is fully removed from the halo via stellar feedback (particularly below log 10 ( 𝑀 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 𝑧 ≈ 0.  \nNext, we discuss the high-mass regime - log 10 ( 𝑀 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 remove gas from the ISM, but is seldom powerful enough to remove 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 with the transition to a hot, virialised halo atmosphere - e.g. Correa et al. 2018a,b) is sufficient to prevent cosmological inflow from entering the ISM (see Fig. 6 and 7), and ultimately curtail star formation rates in high mass galaxies. Since little gas is expelled from the halo and the preventative 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.  \nIn §4, we show in Fig. 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 𝑧 ≈ 0 we show that for haloes with log 10 ( 𝑀 200c / M ⊙ ) ≈ 12 . 5, gas outflow rates at 2 . 5 × 𝑅 200c are just as high as they are at 𝑅 200c in TNG, and even higher in the case of SIMBA. Unlike in EAGLE (where the preventative effect 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.', '5.2 Comparison to previous literature': "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.  \nAs outlined in §2.2 and Appendix A, we have demonstrated that our methodology for measuring gas flow rates produces results that are consistent with previous studies utilising 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 𝑧 ≈ 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 efficient, gas outflow 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 like-for-like 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.  \nWe can also compare our results with those from zoom-in simulations, where improvements in resolution enable more detailed subgrid modelling for gas cooling, star formation and feedback within the ISM. 3 Muratov et al. (2015) and Anglés-Alcázar et al. (2017b) investigate mass loading and baryon cycle predictions from the FIRE simulations, introduced in Hopkins et al. (2014). FIRE implements stellar feedback by injecting energy, momentum, mass, and metals from stellar radiation pressure, H/i.pc/i.pc photo-ionization and photoelectric heating, Type I and Type II supernovae, and stellar winds - while not hydrodynamically decoupling or disabling cooling for affected gas elements. The 𝑀 ★ -𝜂 relations presented in Muratov  \net al. (2015) and Anglés-Alcázar et al. (2017b) (with an Eulerian, shell-based technique and Lagrangian, particle tracking technique respectively) are very similar to those we present for SIMBA in Fig. 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, 𝑀 ★ ≳ 10 9 M ⊙ ).  \nFocusing on an individual system (m12i), Muratov et al. (2015) show that in a halo with 𝑀 200c ≈ 10 12 M ⊙ at 𝑧 ≈ 2, mass loading values at the ISM scale (0 . 25 × 𝑅 200c ) are stochastic (like the outflows themselves), however 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 𝑅 200c reach only ≈ 10 -20% of the outflow rate at the ISM in episodes of star formation. This aligns with the predictions we outline from TNG and SIMBA, where SF-driven outflows typically recycle within the CGM at this mass scale. Comparatively, in EAGLE, mass flow rates at both radii are very similar at 𝑧 ≈ 2 (on average).  \nPandya 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 𝑧 ≈ 0 in a halo of mass 𝑀 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 𝑀 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 𝑀 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 𝑀 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.  \nChristensen et al. (2016) use a particle-tracking technique to study gas flows in and around galaxies in a suite of /g.pc/a.pc/s.pc/o.pc/l.pc/i.pc/n.pc/e.pc-based simulations (Wadsley et al. 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 3 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.  \nTollet et al. (2019) present an analysis of gas flows in the NIHAO simulations using a Lagrangian particle tracking method (Wang et al. 2015). NIHAOalso 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 𝑀 ★ ≈ 10 9 M ⊙ ( 𝜂 ≈ 10 1 . 5 ), however decrease sharply with stellar mass to 𝜂 ≈ 10 -0 . 5 at 𝑀 ★ ≈ 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 𝑀 200c ≲ 10 12 M ⊙ , they show that the majority of gas does not exit 𝑅 200c after being ejected from the ISM.  \nThe total baryon fraction of haloes within 𝑅 200c we measure align  \nwith the findings presented in Ayromlou et al. (2023), who also analyse the EAGLE, TNG, and SIMBA simulations to show how 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 universal value. They find that for low mass haloes where stellar feedback dominates ( 𝑀 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 ( 𝑀 200c ≳ 10 12 M ⊙ ), the closure radius is the largest for SIMBAgalaxies, followed by TNG galaxies, and then EAGLE galaxies.  \nOur findings also provide insight into the measurements presented in Davé et al. (2020), who investigate the global cold gas (H/i.pc and H 2 ; decomposed in post-processing) content of galaxies in the same simulations and compare with available observations. At 𝑧 ≈ 0, Davé et al. (2020) show that the fiducial resolution EAGLE run underpredicts the H/i.pc content of galaxies compared to local ALFALFA observations of the H/i.pc-MF (Jones et al. 2018), as well as the H/i.pc fraction of galaxies as a function of 𝑀 ★ . Given that the H/i.pc 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.pc content of EAGLE galaxies could be explained by the stellar feedback scenario we outline above, which evacuates much of the CGM - and, likely, some of the gas identified as H/i.pc which resides the edge of the ISM. We do note, however, that predictions from the recalibrated, higher resolution EAGLE run and H/i.pc observations and are in closer agreement (see also Fig. B1 and B2).  \nOur results highlight the fact that differences in baryon cycling can indeed leave signatures even at the scale of the ISM, particularly when considering multi-phase gas. With the resolution and physics improvements 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 generalise the results of high-resolution zoom simulations (some of which we summarise above), which provide more detailed insight into the link between the ISM and CGM.", '5.3 Implications for semi-analytic models': "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 review 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 halos via mergers and cosmological accretion with the Λ CDM paradigm. Phenomenological, parameterized 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 supermassive black holes. 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 vs. halo mass relations), cold ISM gas fraction, stellar mass vs. stellar phase metallicity, SFR  \ndistributions, etc. Different SAM groups have adopted very 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.  \nSome 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 et al. 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 parameterize 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. Lgalaxies, SHARK, and related models; Henriques et al. 2015; Lagos et al. 2018). AGN feedback in some models can also help eject material out of the ISM to the CGM, and potentially from the CGM to the IGM (e.g. Lagos et al. 2024).  \nOur study shows that in some hydrodynamical simulations (e.g. EAGLE), the ejected mass on CGM scales ( 𝑅 200c ) can be higher than that on ISM scales (0.25 𝑅 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.  \nIn most SAMs, gas and metals in the ejected reservoir can 'reaccrete' into the CGM on a specified timescale. This timescale, 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 timescales for this 're-accreted' gas, since they show very different rates of gas accretion at halo scales ( 𝑅 200c ).  \nSAMs, 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 𝑓 𝑏 / 𝑀 tot , where / 𝑀 tot is the total rate of mass accretion into the halo 4 . As discussed extensively in Section 5.1, 'preventative' feedback (where inflows are considerably smaller than 𝑓 𝑏 / 𝑀 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 preventative feedback on halo scales due to stellar feedback processes (Lu et al. 2015; Hirschmann et al. 2014; Pandya et al. 2023).  \nMost 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, however, 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  \nfor how AGN winds and jets remove gas and metals from the CGM (e.g. Lagos et al. 2024), as well as for the associated preventative feedback.  \nA promising approach for modeling ejection from the CGM and preventative 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. (2024a,b, in prep.), 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.", '5.4 Identifying observational feedback tests': 'As discussed in §1, 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 favoured by observations. Here, we briefly outline some avenues against which these models could be tested with current and/or future observations.  \nComparing the average gas density profiles between simulations at 𝑧 ≈ 2 and 𝑧 ≈ 0 in Fig. B1 and B2, it is clear that the differences in total halo gas content arising from disparate baryon cycling paradigms are clearest at 𝑧 ≈ 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 weak lensing studies have been able to place constraints on the gas fraction of high-mass galaxy groups and clusters ( 𝑀 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 Figure 2, such observations indicate that the baryon fractions of haloes of mass 𝑀 500c ≈ 10 14 M ⊙ within 𝑅 500c are slightly too high in EAGLE and TNG (very close to the universal value), and slightly too low in SIMBA.  \nAs demonstrated in Fig. 2, the most significant AGN-driven difference in halo baryon content is induced below the mass scale of clusters, between 𝑀 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 𝑀 ★ ≈ 10 11 M ⊙ , forward-modelled X-ray luminosity profiles (up to ≈ 𝑅 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 utilising the Sunyaev Zeldovich (SZ) effect - for instance, Bregman et al. (2022) demonstrate the feasibility of placing such constraints in a small sample of local ≈ 𝐿 ∗ haloes. Yang et al. (2022) demonstrate that the disparate baryon content of haloes between TNG and SIMBA above log 10 ( 𝑀 200c / M ⊙ ) ≈ 12 . 5 produces an observationally detectable difference in forward-modelled SZ measurements, provided that an instrument with adequate angular resolution is utilised (such as the upcoming Simons Observatory, Ade et al. 2019). In this mass range, studies investigating the H/i.pc content of groups using stacking tech-  \nstraining baryon content (e.g. Dev et al. 2023), particularly with radio surveys such as DINGO with the ASKAP telescope (Hotan et al. 2021) and MIGHTEE-HI with MeerKAT (Maddox et al. 2021).  \nBelow mass scales of 𝑀 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 multi-phase 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 sight-line 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 hydrodynamical simulations we study here is inadequate to resolve detailed phase structure, and more detailed sub-grid models are also likely required to accurately model the CGM (e.g. Crain & van de Voort 2023). We can, however, make general comments regarding the expected impact of different feedback regimes on the properties of the CGM.  \nIn 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 previous evacuation 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.  \nIn 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 galacto-centric azimuthal angle for galaxies with 𝑀 ★ ≈ 10 9 . 5 M ⊙ 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 sight-line observations, could constitute a promising test of model accuracy. This is potentially possible with carefully constructed comparisons with current surveys that probe different phases of the CGM, such as COS-Halos (Tumlinson et al. 2013; Werk et al. 2014; Prochaska et al. 2017) and MEGAFLOW (Schroetter et al. 2016; Zabl et al. 2019).  \nAlone this line, DeFelippis et al. (2021) compare results from the MEGAFLOW survey 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 MgII absorbers (tracing the cool CGM) in the halo mass range 10 11 . 5 M ⊙ -10 12 M ⊙ . Appleby et al. (2021) compare results from SIMBAandthe COS-Halos and COS-Dwarfs survey, finding general agreement with observations regarding the predicted abundance of H/i.pc 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  \nthat their model trained on SIMBA can accurately predict absorber overdensity, temperature, and metallicities.  \nThe 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 HST Cosmic Origins Spectrograph, enabling significantly larger samples of 𝑧 ≲ 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 (FRBs; Macquart et al. 2020), and the characterisation 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).  \nObservations of gas flow rates - as opposed to static measurements of halo baryon content - will also prove useful in constraining feedback models in simulations. Different techniques can provide measurements on the multi-phase nature of galactic outflows in individual systems; for instance CO tracing cool molecular outflows (e.g. Leroy et al. 2015 using ALMA observations), or with metal absorption lines and detailed ionisation 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.', '5.5 The impact of resolution': "It is important to note that the resolution of each of these simulations we discuss is inadequate to resolve the multi-phase 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 avoid numerical overcooling 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 𝑀 ★ =10 8 M ⊙ that increasing spatial resolution from > 200 pc to ≈ 20 pc can lead to two-fold increases in mass loading factors at 𝑟 = 5 kpc. They show this is due to the better resolved 'hot' outflowing phase ( 𝑇 ≳ 10 5 K) being heated to higher temperatures, and staying hot for longer periods. Unfortunately, without any sub-grid 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.  \nThese considerations beg the question as to whether improved resolution could influence the results presented in §3 & §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 25Mpc SIMBA box in the case of SIMBA).  \nUpon 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 avoiding any immediate over-cooling), and the mass loading is set based on local DM dispersion as opposed to baryonic properties. At higher halo masses, TNG50 haloes are very  \nslightly more baryon-rich than those in TNG100, however the scale reached by AGN-driven outflows remain identical - indicating that resolution does not qualitatively influence the conclusions we make in this paper.  \nIn 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 over-cooling 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 the CGM, and an associated increase in ISM-scale inflow rates that offset any differences that could otherwise arise in terms of star formation rate.  \nDespite 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 higherresolution counterparts of each simulation.", '6 SUMMARY': 'In this paper, we have conducted a like-for-like 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 star-formation rates at 𝑧 ≈ 0, our work has highlighted that this agreement is achieved for very different physical reasons.  \nConcerning haloes with mass 𝑀 200c ≲ 10 11 . 5 M ⊙ , where stellar feedback is dominant, we summarise our results as follows:  \n- · In EAGLE, SF-driven outflows are able to reach several times 𝑅 200c before stalling. At 𝑧 ≈ 0, mass outflow rates increase with scale with the gas surrounding galaxies being very tenuous and underpressurised relative to the outflows. At 𝑧 ≈ 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.\n- · The extended scale of SF-driven outflows in EAGLE leads to a significant preventative effect on gas inflow at the halo-scale.  \nOverall, this leads to haloes with a low density CGM and very low baryon fractions (just 10% of the cosmological value) in this mass range at 𝑧 ≈ 0.  \n- · In TNG, SF-driven outflows are very strong at the scale of the ISM, however, tend to stall within the CGM before reaching 𝑅 200c . This is also true in SIMBA at 𝑧 ≈ 2, however at 𝑧 ≈ 0, SF-driven outflows in SIMBA can reach a similar scale as in EAGLE (albeit without mass entrainment).\n- · In TNG, halo baryon fractions remain above 50% of the cosmological value in this mass regime at 𝑧 ≈ 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.\n- · In SIMBA, halo baryon fractions in this mass regime are between those measured in EAGLE and TNG, reflecting moderate preventative feedback on halo-scale gas inflows.  \nFurthermore, in relation to haloes with mass 𝑀 200c ≳ 10 12 M ⊙ , where AGN feedback is dominant:  \n- · In EAGLE, AGN-driven outflows do not reach far beyond 𝑅 200c  \nat either 𝑧 ≈ 0 nor 𝑧 ≈ 2. This leads to minimal preventative impact on gas inflow at the halo scale, but significantly prevents gas inflow at the scale of the ISM.  \n- · In TNG, AGN-driven outflows are more powerful than in EAGLE, with at least some gas in these outflows reaching 2 -3 × 𝑅 200c . This leads to a moderate preventative impact on cosmological inflow at the halo-scale, with further inflow prevention within the heated CGM. As a result, there is a significant turnover in halo baryon fractions at 𝑧 ≈ 0, reaching values as low as 50% of the cosmological value at 𝑀 200c ≈ 10 12 . 5 M ⊙ .\n- · In SIMBA, AGN-driven outflows are significantly stronger than in both EAGLE and TNG, particularly at 𝑧 ≈ 0. These outflows continually entrain mass from the ISM scale to several times 𝑅 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% of the universal value between 𝑀 200c = 10 12 M ⊙ -10 13 M ⊙ .  \nThis work lays the foundation for developing targeted observational tests that can favour or disfavour certain feedback scenarios, and move 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 maximise 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. 2023) will be imperative in order to extract meaningful constraints from future weak lensing studies with regard to cosmological parameters, cosmic shear, and the growth of cosmic structure.', 'ACKNOWLEDGEMENTS': 'RJW acknowledges support from the European Research Council via ERC Consolidator Grant KETJU (no. 818930), the Fulbright Australia Commission for a Postdoctoral Visiting Scholar Grant, and the Australian Research Council Centre of Excellence for All Sky Astrophysics in 3 Dimensions (ASTRO 3D).  \nDAA acknowledges support by NSF grants AST-2009687 and AST-2108944, CXO grant TM2-23006X, JWST grant GO01712.009-A, Simons Foundation Award CCA-1018464, and Cottrell Scholar Award CS-CSA-2023-028 by the Research Corporation for Science Advancement.  \nFor the purpose of open access, the author has applied a Creative Commons Attribution (CC BY) licence to any Author Accepted Manuscript version arising from this submission.  \nThe authors used the following software tools for the data analysis and visualisation in the paper:  \n- · python3 (van Rossum & Drake Jr 1995)\n- · numpy (Harris et al. 2020)\n- · scipy (Virtanen et al. 2020)\n- · matplotlib (Hunter 2007)', 'DATA AVAILABILITY': 'The EAGLE ( http://icc.dur.ac.uk/Eagle/ ), TNG ( https: //www.tng-project.org/ ), and SIMBA ( http://simba.roe. ac.uk/ ) simulation outputs used for our analysis are all publicly available. The gas flow catalogs generated for this work will be shared by the corresponding author upon reasonable request.', 'APPENDIX A: COMPARISON OF GAS FLOW MEASUREMENTS WITH PREVIOUS LITERATURE': "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 §2.2, we elect to use an Eulerian or 'instantaneous' approach to calculating gas flow rates, rather than a Lagrangian particle tracking method.  \nFigure A1 shows a comparison of 𝑧 ≈ 2 and 𝑧 ≈ 0 halo-scale mass loading factors at the ISM scale (here, 0 . 20 × 𝑅 200c ) in EAGLE as a function of 𝑀 200c mass. While Mitchell et al. (2020) take a Lagrangian particle-tracking approach to measure gas flow rates, we recover good agreement with their results (to within ≈ 0 . 2 dex 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 Δ 𝑡 used for the Lagrangian calculation. This offset is highest at 𝑀 200c ≈ 10 13 M ⊙ , indicating  \nFigure A1. Comparison of 𝑧 ≈ 2 (blue) and 𝑧 ≈ 0 (black) halo-scale mass loading factors ( 𝜂 ≡ ⟨ / 𝑀 out ⟩/⟨ / 𝑀 ★ ⟩ ) as a function of 𝑀 200c mass measured in this work (solid lines, no shading), and the measurements presented in Mitchell et al. (2020) (square markers connected with dot-filled solid lines). In our calculations, hatched shaded regions correspond to the bootstrapgenerated 90% confidence interval on the medians at a given mass. At both redshifts, our measurements agree well, to within ≈ 0 . 2 dex.  \n<!-- image -->  \nthat some AGN-driven outflows in EAGLE may be recycled quite quickly at the ISM scale.  \nFigure A2 shows a comparison of 𝑧 ≈ 2 10kpc-scale mass loading factors as a function of 𝑀 ★ in TNG for two different choices of outflowvelocity cut: 50 km s -1 , and 150 km s -1 . Utilising 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 §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.", 'APPENDIX B: GALAXY DENSITY PROFILES': 'In Fig. B1 and B2, we show the density profiles of galaxies in each of the simulations at 𝑧 = 2 and 𝑧 = 0 respectively. Each row represents a different bin in halo mass - top: log 10 ( 𝑀 200c / M ⊙ ) ∈ [ 10 . 75 , 11 . 25 ] ; middle: log 10 ( 𝑀 200c / M ⊙ ) ∈ [ 11 . 75 , 12 . 25 ] ; bottom: log 10 ( 𝑀 200c / M ⊙ ) ∈ [ 12 . 75 , 13 . 25 ] ). Radial 𝑥 -values are quoted as fractions of 𝑅 200c , with a broken scale above 1 × 𝑅 200c to include radial bins up to 3 × 𝑅 200c . For reference, in each panel we include a dashed line indicating the expected outer slope of an isothermal 𝛽 profile with 𝛽 = 2 / 3; i.e. 𝜌 ( 𝑟 ) ∝ 𝑟 -2 .  \nIn general, density profiles at 𝑧 = 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 𝑧 = 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 overpressurised relative to the surrounding  \nFigure A2. Comparison of 𝑧 ≈ 2 mass loading factors ( 𝜂 ≡ ⟨ / 𝑀 out ⟩/⟨ / 𝑀 ★ ⟩ ) at 10 kpc measured in this work (solid lines, no shading) and those presented in Nelson et al. (2019) (square markers connected with dot-filled solid lines). 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% confidence interval on the medians at a given mass. For both velocity thresholds, our measurements agree very well, to within ≈ 0 . 1 dex.  \n<!-- image -->  \nmedium. 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 efficient evacuation of the CGM via AGN feedback, which also allows subsequent outflows to be overpressurised relative to the CGM and entrain gas in the outflows towards larger scales.', 'APPENDIX C: GALAXY TEMPERATURE PROFILES': 'In Fig. C1 and C2, we show the temperature profiles of galaxies in each of the simulations at 𝑧 = 2 and 𝑧 = 0 respectively. Each row represents a different bin in halo mass - top: log 10 ( 𝑀 200c / M ⊙ ) ∈ [ 10 . 75 , 11 . 25 ] ; middle: log 10 ( 𝑀 200c / M ⊙ ) ∈ [ 11 . 75 , 12 . 25 ] ; bottom: log 10 ( 𝑀 200c / M ⊙ ) ∈ [ 12 . 75 , 13 . 25 ] ). Radial 𝑥 -values are quoted as fractions of 𝑅 200c , with a broken scale above 1 × 𝑅 200c to include radial bins up to 3 × 𝑅 200c . For reference, in each panel we include a horizontal line indicating the expected virial temperature of haloes in this mass range, with 𝑇 vir ≈ 1 . 1 × 10 6 K ( 𝜇 / 0 . 59 ) GLYPH<16> 𝑀 halo / 10 12 M ⊙ GLYPH<17> 2 / 3 ( 1 + 𝑧 ) , from van de Voort (2017).  \nAs discussed in §4.1, considering haloes with mass 𝑀 200c ≲ 10 13 M ⊙ at 𝑧 ≈ 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 × 𝑅 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  \nFigure B1. Radial density profiles of galaxies at 𝑧 = 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 16 th -84 th percentile range in inflow rates at a given mass, and hatched regions correspond to the bootstrap-generated 90% 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.  \n<!-- image -->  \nlower mass bins (signifying a universal shift to a hotter CGM, and the dominance of hot-mode accretion).', 'APPENDIX D: HALO BARYON FRACTIONS AT REDSHIFT 2': "In this section, we add to the 𝑧 ≈ 0 results shown in §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  \nFigure B2. Radial density profiles of galaxies at 𝑧 = 0 in EAGLE (blue), EAGLE-RECAL (teal), TNG (purple) and SIMBA (pink). Each row represents a different bin in halo mass, increasing top to bottom. Error-bars correspond to the 16 th -84 th percentile range in inflow rates at a given mass, and hatched regions correspond to the bootstrap-generated 90% confidence interval on the medians at a given mass. Density profiles diverge between the simulations at this redshift, with particularly low densities in EAGLE in the lowest halo mass bin, and SIMBA for the two higher halo mass bins.  \n<!-- image -->  \nin 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 𝑅 200c that does not meet this criteria.  \nAt this redshift, the baryon content of haloes is universally higher in the simulations - for 𝑀 200c ≳ 10 11 . 5 M ⊙ , all simulations predict a baryon fraction above 50% of the cosmic mean. Similar to 𝑧 ≈ 0, the baryon content of EAGLEgalaxies at low halo mass 𝑀 200c ≲ 10 12 M ⊙ - is lower than the predictions from SIMBA and TNG, indicating that the ejective and preventative influence of stellar  \nFigure C1. Radial temperature profiles of galaxies at 𝑧 = in EAGLE (blue), EAGLE-RECAL (teal), TNG (purple) and SIMBA (pink). Each row represents a different bin in halo mass, increasing top to bottom. Error-bars correspond to the 16 th -84 th percentile range in inflow rates at a given mass, and hatched regions correspond to the bootstrap-generated 90% confidence interval on the medians at a given mass. Temperature predictions diverge between the simulations at this redshift in the two lower mass bins, particularly in the inner CGM.  \n<!-- image -->  \nfeedback is already at play by 𝑧 ≈ 2. The transition masses M T 1 and M T 2 are similar to those reported at 𝑧 ≈ 0, with the exception of SIMBA,wherethe M T 1 mass at 𝑧 ≈ 2 increases to ≈ 10 12 M ⊙ (as opposed to ≈ 10 11 . 3 M ⊙ at 𝑧 ≈ 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 𝑧 ≈ 0.  \nThis paper has been typeset from a T E X/L A T E X file prepared by the author.  \nFigure C2. Radial temperature profiles of galaxies at 𝑧 = 0 in EAGLE (blue), EAGLE-RECAL (teal), TNG (purple) and SIMBA (pink). Each row represents a different bin in halo mass, increasing top to bottom. Error-bars correspond to the 16 th -84 th percentile range in inflow rates at a given mass, and hatched regions correspond to the bootstrap-generated 90% confidence interval on the medians at a given mass. Temperature profiles are in slightly better agreement between the simulations compared to 𝑧 ≈ 2, with generally hotter CGM.  \n<!-- image -->  \nFigure D1. The baryon content of haloes at 𝑧 ≈ 2 in EAGLE (blue), TNG (purple) and SIMBA (pink). Top panel, 1 st column : The median total baryon fraction within 𝑅 200c (stars, gas, black holes) of haloes as a function of 𝑀 200c mass. Top panels, 2 nd -4 th columns : Total baryon content for each simulation within 𝑅 200c broken down into stars (hatched bars), 'cold' gas (coloured shaded bars), and 'hot gas' (grey bars) in each simulation. In the top panels, this baryon content is normalised by 𝑓 b × 𝑀 200c . Bottom panels : The bottom row of panels shows the same breakdown, but re-normalised by the actual baryonic mass in each halo, 𝑀 bar ( 𝑟 < 𝑅 200c ) . In each of the simulations in the left-hand panel, error-bars represent the 16 th -84th percentiles in halo baryon fractions for a given mass bin; the hatched regions represent the 90% confidence interval on the median from 100 bootstrap re-samplings in each bin. The difference between the simulations in terms of 𝑧 ≈ 2 baryon fractions are qualitatively similar compared to that shown at 𝑧 ≈ 0 in Fig. 2, but the baryon content of haloes is universally higher.  \n<!-- image -->"}