Galactic outflow rates in the EAGLE simulations

We present measurements of galactic outﬂow rates from the EAGLE suite of cosmological simulations. We ﬁnd that gas is removed from the interstellar medium (ISM) of central galaxies with a dimensionless mass loading factor that scales approximately with circular velocity as V − 3 / 2 c in the low-mass regime where stellar feedback dominates. Feedback from active galactic nuclei causes an upturn in the mass loading for halo masses > 10 12 M (cid:2) . We ﬁnd that more gas outﬂows through the halo virial radius than is removed from the ISM of galaxies, particularly at low redshifts, implying substantial mass loading within the circumgalactic medium. Outﬂow velocities span a wide range at a given halo mass/redshift, and on average increase positively with redshift and halo mass up to M 200 ∼ 10 12 M (cid:2) . Outﬂows exhibit a bimodal ﬂow pattern on circumgalactic scales, aligned with the galactic minor axis. We present a number of like-for-like comparisons to outﬂow rates from other recent cosmological hydrodynamical simulations, and show that comparing the propagation of galactic winds as a function of radius reveals substantial discrepancies between different models. Relative to some other simulations, EAGLE favours a scenario for stellar feedback where agreement with the galaxy stellar mass function is achieved by removing smaller amounts of gas from the ISM, but with galactic winds that then propagate and entrain ambient gas out to larger radii.


I N T RO D U C T I O N
In the modern cosmological paradigm, galaxies grow within dark matter haloes, which represent collapsed density fluctuations that in turn grow via gravitational instability from a near-homogeneous initial density field. In this picture, galaxies do not form in monolithic formation events, and instead grow gradually via sustained periods of gaseous inflow from the larger scale environment, tracing the hierarchical build-up of dark matter haloes (e.g. Blumenthal et al. 1984). Star formation proceeds with sufficient efficiency to deplete gas from the interstellar medium (ISM) over a time-scale that is comparable or shorter than a Hubble time, and as such galaxy evolution is to zeroth order set by the fluxes of gas into and out of the ISM (e.g. Schaye et al. 2010;Davé, Finlator & Oppenheimer 2012).
Observationally, direct measurements of inflowing gas fluxes have remained elusive, with only a handful of reported detections (e.g. Rubin et al. 2012;Fox et al. 2014;Roberts-Borsani & Saintonge 2019). Detections and evidence for outflowing gas are comparatively plentiful (e.g. Heckman et al. 2000;Strickland & Heckman 2009;Feruglio et al. 2010;Steidel et al. 2010;Rubin et al. 2014;Schroetter et al. 2016), although determinations of the associated mass flux are likely beset by a number of systematic uncertainties (e.g. Chisholm et al. 2016), and a given outflow tracer probes gas over only a subset of the relevant spatial scales and gas phases. E-mail: mitchell@strw.leidenuniv.nl The need for substantial outflowing fluxes has long been recognized, for example in order to explain the form of the observed galaxy luminosity function (e.g. White & Frenk 1991;Benson et al. 2003), the correlation between galaxy mass and metallicity (e.g. Larson 1974), and the presence of metals in the diffuse intergalactic medium (e.g. Aguirre et al. 2001). Feedback in the form of mass, momentum, and energy input from massive stars and supermassive black holes (SMBHs) is thought to be responsible for driving outflows from galaxies (e.g. Larson 1974;Silk & Rees 1998). These feedback mechanisms are a core element of modern phenomenological models and simulations that reproduce the observed properties of the overall galaxy population (e.g. Somerville et al. 2008;Vogelsberger et al. 2014;Schaye et al. 2015).
Determining the efficiency with which galactic winds are driven as a function of the rates at which mass, momentum, and energy are injected into the ISM represents one of the major outstanding challenges of modern astrophysics, both from the observational and theoretical perspectives. Relevant radiative losses occur in principle over an enormous dynamic range in scale, and depend on the properties of the ambient medium over this range. Numerical simulations are routinely used to explore this problem, again over scales ranging from the small-scale ISM (e.g. Chevalier 1974;Walch & Naab 2015), to the entire galaxy population , and scales in between (e.g. Hopkins, Quataert & Murray 2012;Creasey, Theuns & Bower 2013;Kim & Ostriker 2018).
On the large-scale end of this distribution of numerical studies, the EAGLE simulation project simulates the formation and evolution of galaxies within the full Cold Dark Matter ( CDM) context, integrating periodic cubic boxes (up to 100 3 Mpc 3 in volume) down to z = 0 (Crain et al. 2015;Schaye et al. 2015). At the reference resolution of the project, these simulations employ a fiducial baryonic particle mass of 1.81 × 10 6 M , and reach a maximum spatial resolution of about 1 kpc at z = 0, and so do not resolve the physics of the ISM. As with other simulations of this type (e.g. Schaye et al. 2010;Dubois et al. 2014;Vogelsberger et al. 2014;Davé et al. 2017), this means that the EAGLE simulations cannot make accurate predictions for the radiative losses that occur on ISM scales, and a strategy must be adopted to avoid the spurious losses that would occur should the energy injected by feedback be smoothly distributed.
In the case of EAGLE, spurious losses are mitigated by heating relatively few ISM particles to a high temperature (10 7.5 K for stellar feedback; , with the unresolved radiative losses then set by hand with model parameters that are calibrated by comparing to various observational constraints. As discussed by Crain et al. (2015), it is possible to produce an acceptable fit to the galaxy stellar mass function inferred from observations by assuming that 100 per cent of the energy available from Type II supernovae is able to heat gas to high temperatures [in addition to the energy injection provided by active galactic nuclei (AGNs)]. To also reproduce the observed distributions of galaxy sizes as a function of mass, it was found that the energy injected per unit stellar mass had to vary by factors of a few, scaling negatively with gas metallicity and positively with density. EAGLE is therefore differentiated from a number of similar projects (e.g. Vogelsberger et al. 2014;Davé et al. 2017) that instead mitigate spurious losses by temporarily decoupling the particles that are kicked by feedback from the hydrodynamical scheme, while also disabling radiative cooling for these particles. In such alternative schemes, particles are explicitly kicked with a velocity that scales linearly with the circular velocity of the system, and the rate of mass of particles kicked per unit rate of mass of stars formed (defining the dimensionless mass loading factor) is assumed to scale negatively with circular velocity. As no such explicit scaling with galaxy properties is utilized in EAGLE, 1 the mass loading and velocities of galactic winds are instead emergent phenomena, presumably determined (for example) by the escape velocity of system, and the column density of gas that winds must push through to break out of the ISM. This does not guarantee that outflow scalings in the EAGLE simulations are necessarily more realistic than other schemes used at the same resolution; the EAGLE feedback schemes are still approximate in nature, and do not explicitly simulate much of the relevant astrophysics that is studied analytically or numerically at much higher numerical resolution. Furthermore, phenomenological feedback schemes that use decoupled winds only do so temporarily, recoupling wind particles once they leave the ISM, meaning that emergent behaviours such as anisotropic flow patterns develop despite not being imposed at injection , and the interaction between outflows and ambient circumgalactic gas is fully simulated.
We set out in this study to measure the outflow rates of galactic winds from central galaxies in the EAGLE simulations. At a basic level, this allows us to better understand how and why different aspects of galaxy evolution proceed in a given manner within the simulation, adding valuable information that can be used to interpret the myriad of other results already published based on analyses of EAGLE. This work also serves as an introduction to a more complete upcoming study of the network of inflows, outflows, and recycling of gas flows from EAGLE, and we take care to explain our methodology within this context. For a more observations-focused analysis of outflows in the EAGLE simulations, we refer readers to Tescari et al. (2018), who analyse the simulations within the context of recent integral field unit observations. In addition, a preliminary version of our inflow rate measurements is used in Collacchioni et al. (2019) to study the connection between inflows and radial metallicity gradients.
On a broader level, we use our measurements of outflow rates to provide a viable quantitative scenario for how galaxy evolution might proceed across most of the relevant redshift range and galaxy mass scales. We make the effort to show like-for-like comparisons with other simulation projects (both large-volume simulations and zoom-in simulations) to check whether there is yet any consensus emerging from cosmological simulations (the short answer is that there is little quantitative agreement at present, but there is rough qualitative agreement). All of the simulations we compare to achieve (to a greater or lesser extent) at least somewhat reasonable agreement with the observed stellar properties of galaxies, and so the range of outflow rates shown in the comparisons might guide observers as well as smaller scale simulators as to what is likely required from galactic winds in order to explain the observed galaxy stellar mass function.
The layout of this paper as follows: We introduce our methodology for measuring outflow rates in Section 2, we present measurements of outflow rates and velocities from EAGLE in Section 3. We finish by placing our work into the wider context of theoretical models, simulations, and observations in Section 5, and we summarize our results in Section 6.

Rationale
Our objective is to measure the amount of gas that is ejected from galaxies and their associated dark matter haloes in the EAGLE simulations. This is essential in order to understand the emergent relationship between stellar mass, gas mass [in the ISM and also the circumgalactic medium (CGM) out to the virial radius], and total halo mass. Outflow rates can be measured from simulations using either Eulerian or Lagrangian methods. The former involves measuring the instantaneous flux of outflowing gas through a surface (or within a shell) at a given distance from the centre of the galaxy or halo (e.g. Mitchell et al. 2018b;Nelson et al. 2019). The latter method involves measuring the flux of mass that crosses a surface over a discrete time interval (e.g. Neistein et al. 2012;Christensen et al. 2016;Anglés-Alcázar et al. 2017).
We opt to use a Lagrangian method to measure outflow rates. Our primary motivation for this choice is that the method enables accurate measurements of the correct time-integrated outflow rate of a given galaxy. This is particularly pertinent for the EAGLE simulations, where the high heating temperature used in the subgrid model leads to highly time-variable instantaneous outflow rates. The primary drawback of the Lagrangian method is that correct time-integrated fluxes are only obtained if fluid elements cross the surface only once over the finite time interval adopted (fluid elements that cross multiple times cause an underestimate of the true time-integrated flux). In practice, this means that a substantial number of simulation outputs (roughly 200 in our case) are required to achieve converged outflow rates of gas being ejected from the ISM (see Appendix A1), as the time-scale between gas entering and exiting the ISM can be short compared to the halo dynamical time. Note that, when we show average radial velocities, or energy and momentum fluxes, we will switch to Eulerian measurements based on discrete shells; this is because (unlike mass) these quantities are not necessarily conserved after leaving the ISM, and so are more clearly defined at a fixed radius.
Another aspect of measuring gas fluxes from simulations is the choice of surface or shell, and the choice of which subset of the fluid elements flowing through the surface should be selected for the measurement. On the one hand, simple choices for both yield measurements that are easy to reproduce and compare with other simulations, and the same also applies for comparison with observational studies to some extent. On the other hand, adopting an arbitrary choice of surface runs the risk of not capturing the desired quantity, which we take to be the flux of gas being removed from the ISM. In simulations like EAGLE that model the galaxy population across a wide range in mass and redshift, the star-forming gaseous content of a galaxy can vary hugely in structure and spatial scale (both in an absolute sense and relative to the halo), as is ably demonstrated by the two examples shown in appendix C of Mitchell et al. (2018a). Furthermore, non-negligible amounts of the outflowing flux on scales close to the ISM can be associated with gas that is moving past pericentre on orbits that are driven primarily by gravity (rather than by feedback).
For these reasons, we have adopted (and laboriously checked) criteria that select gas that was within the ISM (at the previous simulation output) and has now (at the current simulation output) exited the ISM, and is in the process of moving out over a significant distance into the CGM. A direct comparison of simple Eulerian measurements with our full Lagrangian criteria is shown in Appendix A4, for readers who may be interested to see the impact of our selection criteria on our conclusions. Our methodology is similar to that of Christensen et al. (2016), who measure gas particles that leave an ISM defined in a similar way using phase cuts, and that outflow with kinetic energy exceeding that of the gravitational potential, as well as that of Anglés-Alcázar et al. (2017), who perform similar measurements but instead define the ISM with a Friends-of-Friends (FoF) algorithm, along with a cut in gas density.

Simulations and subgrid physics
The EAGLE project is a suite of hydrodynamical simulations that simulate the formation and evolution of galaxies within the context of the CDM cosmological model , and that have been publically released (McAlpine et al. 2016). The suite was created using a modified version of the GADGET-3 code (last presented in Springel et al. 2005), and features a number of cosmological periodic boxes containing both gas and dark matter, integrated down to z = 0. Cosmological parameters are set following Planck Collaboration XVI (2014), with m = 0.307, = 0.693, b = 0.048 25, h = 0.6777, and σ 8 = 0.8288. The suite employs a state-ofthe-art implementation of smoothed particle hydrodynamics (SPH; see Schaye et al. 2015;Schaller et al. 2015a), and a range of subgrid models that account for important physical processes that are not resolved by the simulation (radiative cooling, star formation, stellar mass-loss and metal enrichment, SMBH growth, and energy injection from stellar and AGN feedback).
Unless otherwise stated, all results presented here are produced using the reference 100 3 cMpc 3 simulation, which includes 1504 3 particles for both gas and dark matter, with particles masses of 1.81 × 10 6 M and 9.70 × 10 7 M for gas and dark matter, respec-tively. This simulation (referred to as L0100N1504 in Schaye et al. 2015) uses the subgrid models and parameters of the EAGLE reference model described by Schaye et al. (2015) (and also discussed in detail by Crain et al. 2015). Hereafter, we refer to this simulation as the 100 Mpc reference run. In some parts we also utilize smaller 25 3 and 50 3 cMpc 3 versions of the reference simulation (with the same physics and resolution), as well as a 50 3 cMpc 3 simulation that was simulated without AGN feedback.
An overview of the salient aspects of the EAGLE reference model within the context of this study is as follows. First, stars are allowed to form above the metallicity-dependent threshold for which the gas is expected to become cold and molecular (Schaye 2004), where Z is the gas metallicity. 2 Gas particles are artificially pressurized up to a minimum pressure floor set proportional to gas density as P ∝ ρ 4/3 g , normalized to a temperature of T = 8 × 10 3 K at a hydrogen density of n H = 0.1 cm −3 . This acts to ensure that the thermal Jeans mass is always at least marginally resolved, but prevents the formation of a cold ISM phase. In addition to equation (1), gas particles are eligible to form stars only if they are within 0.5 dex in temperature from the temperature floor.
Star formation is implemented stochastically as described in , with individual gas particles being converted into collisionless star particles by sampling from a probability distribution such that the star formation rate is given by where m gas is the gas particle mass, P is the local gas pressure, γ = 5/3 is the ratio of specific heats, G is the gravitational constant, and f g is the gas mass fraction (set to unity). A and n are taken from the observed Kennicutt-Schmidt star formation law,˙ = A( g /1 M pc −2 ) n , and are set to A = 1.515 × 10 −4 M yr −1 kpc −2 and n = 1.4 (Kennicutt 1998), with n changed to n = 2 for hydrogen densities greater than n H = 10 3 cm −3 . Stellar feedback is represented by stochastic thermal energy injection, following the methodology introduced by . In this scheme, gas particles are heated by neighbouring star particles by a fixed temperature jump, T = 10 7.5 K, with a probability set such that the average thermal energy injected is f th × 8.73 × 10 15 erg g −1 of stellar mass formed, where f th is a model parameter. For f th = 1, the injected energy per unit stellar mass corresponds to that of a simple stellar population with a Chabrier initial mass function, assuming that 6-100 M stars explode as supernovae and that each supernova injects 10 51 erg of energy. Neighbouring gas particles are heated by stellar feedback 30 Myr after the formation of a star particle.
In order to empirically recover an adequate match to both the galaxy stellar mass function and the galaxy size versus stellar mass distribution inferred from observations (Crain et al. 2015), f th is varied as a function of local gas metallicity, Z, and the gas density, n H, birth , inherited by the star particle from the gas from which it formed, with the parametrization given by where f th, min and f th, max are model parameters that are the asymptotic values of a sigmoid function in metallicity, with a transition scale at a characteristic metallicity, 0.1Z (above which radiative losses are expected to increase due to metal cooling; Wiersma, Schaye & Smith 2009), and with a width controlled by n Z . An additional dependence on local gas density is controlled by model parameters, n H, 0 , and n n . The two asymptotes, f th, min and f th, max , are set to 0.3 and 3, respectively, such that between 0.3 and 3 times the canonical supernova energy is injected. n Z and n n are both set to 2/ln (10), and n H, 0 is set to 0.67 cm −3 . SMBH growth is modelled first by seeding SMBH particles at the position of the highest density gas particle within dark matter haloes with a mass of M FOF > 10 10 M h −1 , where M FOF is the mass of the FoF group. Black hole particles then accrete mass with an Eddingtonlimited Bondi accretion rate that is modified if the accreted gas is rotating at a velocity that is significant relative to the sound speed (Rosas-Guevara et al. 2015;Schaye et al. 2015). Black holes that are sufficiently close to each other in position and velocity are allowed to merge, forming a second channel of black hole growth.
Analogous to the implementation of stellar feedback, accreting SMBH particles stochastically inject thermal energy into neighbouring gas particles (Booth & Schaye 2009), with an energy injection ratė whereṁ acc is the gas mass accretion rate on to the SMBH, c is the speed of light, r is the fraction of the accreted rest mass energy that is radiated (set to 0.1), and f is a model parameter that sets the fraction of the radiated energy that couples to the ISM (set to 0.15). The injected thermal energy is stored in the SMBH particle until it is sufficiently large to, on average, heat a single neighbouring gas particle by T = 10 8.5 K, a temperature jump that is an order of magnitude larger than the value used for stellar feedback ( T = 10 7.5 K).

Relating phenomenological feedback modelling to the underlying astrophysics
Having described the salient features of the phenomenological star formation and feedback modelling used in the EAGLE simulations, we briefly discuss here how this relates to the underlying astrophysics of feedback. As presented in Schaye et al. (2015), the conceptual intent for the stellar feedback model in EAGLE is that it represents the combined effects of all stellar feedback processes that are thought to be relevant for galaxy evolution, including supernovae, radiative feedback (e.g. Krumholz & Dekel 2012;Rosdahl et al. 2015), stellar winds (e.g. Gatto et al. 2017), and the effects of cosmic rays seeded by supernovae (e.g. Uhlig et al. 2012;Booth et al. 2013;Girichidis et al. 2016). By imposing a star formation law that reproduces the observed Kennicutt-Schmidt relation, feedback is not required to set the local efficiency of star formation, reducing the need for 'early' stellar feedback that pre-processes the ISM before SNe explode. Furthermore, due to the coarse numerical resolution of EAGLE (compared to the aforementioned numerical studies), and with the equation of state that artificially pressurizes the ISM, the gas phase distribution in the ISM is not expected to be realistic (regardless of whether early stellar feedback is included), which precludes the robust application of much higher resolution calculations that predict (for example) the momentum and energy injection for isolated supernovae as a function of local density and metallicity (e.g. Cioffi, McKee & Bertschinger 1988;Kim & Ostriker 2015;Walch & Naab 2015;Gentry et al. 2017;Gentry, Madau & Krumholz 2020). Similarly, the AGN feedback model used in EAGLE is also a heavily coarse-grained description of the underlying astrophysics. The exact mechanism by which energy is coupled to the surrounding gas is not specified, and the scheme may also need to mimic the outcome of plasma physics in relation to AGN feedback in the intracluster medium, such as the effects of cosmic rays in heating and providing pressure support (e.g. Loewenstein, Zweibel & Begelman 1991;Sijacki et al. 2008;Ruszkowski, Yang & Reynolds 2017). These limitations are important to keep in mind when interpreting results from simulations like EAGLE. The trade-off, however, is that by using simple phenomenological feedback schemes that mitigate immediate radiative losses, we can calibrate a simulation to produce a realistic and representative population of galaxies across a significant dynamic range in stellar mass, and so present a physically viable scenario for how mass and energy fluxes at different scales regulate the growth of galaxies.

Subhalo identification and merger trees
Haloes are first identified from a given simulation output as groups, using an FoF algorithm, with a dimensionless linking length of b = 0.2 (Davis et al. 1985). FoF groups are then split into subhaloes using the SUBFIND algorithm (Springel et al. 2001;Dolag et al. 2009). Each subhalo consists of a set of bound particles (including gas, stars, black holes, and dark matter). For each FoF group, the subhalo containing the particle with the lowest value of the gravitational potential is defined as the central subhalo (and galaxy). Other subhaloes within the FoF group are defined as satellites. The subhalo (and associated galaxy) centre is defined as the position of the particle with the lowest value of the gravitational potential. Finally, for central subhaloes we take an additional step and add/remove particles that are within/outside R 200 , 3 provided the particles are not associated with another subhalo or FoF group. Here, R 200 is the radius enclosing a mean spherical overdensity that is 200 times the critical density of the Universe at a given epoch. Halo masses and virial radii quoted throughout this paper are defined as M 200 and R 200 , respectively, where M 200 is the mass enclosed within R 200 .
We construct merger trees using the algorithm described in appendix A of Jiang et al. (2014). In brief, for each subhalo in a given simulation output (the progenitor in question), the algorithm attempts to identify a single descendant subhalo in the next simulation output. The descendant is selected as the subhalo containing the largest fraction of a set of the progenitor's most-bound particles. Furthermore, if the largest fraction of a set of the most-bound particles of the descendant come from the progenitor in question, the progenitor is identified as the main progenitor of the descendant. In cases where the progenitor in question is not identified as a main progenitor, a number of later simulation outputs are also searched in 3 In practice, this acts to add gas particles within the virial radius that have been raised by feedback to sufficiently high internal plus kinetic energy that they are no longer considered bound to the subhalo by SUBFIND. We need to keep these particles associated with the subhalo in order to ensure that our measurements of halo outflow rates are correct. an attempt to find a descendant for which the progenitor in question is the main progenitor. This procedure accounts for cases where subhaloes temporarily cannot be identified by SUBFIND against the backdrop of a larger subhalo. In post-processing we identify rare cases where the identified main progenitor of a descendant is a clump identified as a subhalo by SUBFIND, but is dominated by star and black hole particles, rather than dark matter particles. In these cases, we find the most massive progenitor of the descendant and set that subhalo as the main progenitor. Put together, this is then the definition of the main progenitor that we use throughout our analysis (in the sense that we measure particles that were present in the ISM/halo of the main progenitor that have since been ejected from the descendant).
We use a number of sets of merger trees constructed with differing numbers of simulation outputs. Most of our results use trees constructed with 200 simulation snipshots, where snipshots are simulation outputs that contain a subset of the information available for each particle from the more sparsely sampled simulation snapshots. The temporal spacing between these 200 snipshots is shown in Appendix A1. In some cases, we use merger trees constructed with different numbers of snipshots or snapshots, either to test the temporal convergence of our method, because processed SUBFIND outputs were not available for a given simulation, or because we required particle information that is only present within the snapshots.

Particle partitioning
Within a given subhalo, we partition the baryonic particles into a discrete number of groups. First, star and black hole particles form two distinct groups. For gas particles, we select particles belonging to the ISM, with the remainder forming a circumgalactic halo component.
Our ISM selection criteria are closely related to the star formation criteria used in the simulation. We define the ISM as the sum of (i) Star-forming gas (i.e. particles with n H > n H and are within 0.5 dex of the temperature floor), irrespective of radius.
The choice to include non-star-forming gas down to n H = 0.01 cm −3 is made primarily to account for dense gas in lowmass haloes with low metallicity, and in effect approximately selects neutral hydrogen out to the imposed radius cut (Rahmati et al. 2013). The effect of this inclusion for our results is to significantly enhance the outflow rates of low-mass galaxies (see Appendix A3), where little star formation and chemical enrichment have occurred. The inclusion also increases the specific angular momentum of the ISM (by effectively selecting more diffuse neutral material in the outskirts of galaxy discs), which we plan to study in the context of inflows/outflows in future work (see also Mitchell et al. 2018b).
We impose a radial cut for the non-star-forming ISM component to exclude dense and low-metallicity infalling and filamentary circumgalactic material (found mostly at high redshift). We do not impose any radial cut for star-forming gas in order to account for stellar feedback that occurs outside of this radius, which is relevant for removing gas from the star-forming gas reservoir of galaxies at high redshift in the simulation (z 2).

Measuring outflow rates
We use a Lagrangian particle tracking method to measure gas outflow rates from galaxies and haloes. We define galaxy-scale outflow rates as the summed mass of particles leaving the ISM per unit time, measured over some finite time interval between two simulation outputs. Halo-scale outflow rates are then defined accordingly for particles leaving the halo virial radius per unit time. In both cases, we apply the additional selection criteria described below to check that the particles are genuinely outflowing. Further details of the rationale, exploration, and testing that was used to arrive at these criteria are described in Appendix A, along with a comparison to simple shell-based outflow rate measurements.
For both galaxy-scale and halo-scale outflows, we require that outflowing particles satisfy and for galaxy-scale outflows, we also require that where V max is the maximum of circular velocity profile of the halo and v rad,1 is the instantaneous radial velocity of the particle at the first simulation output after the particle has left the ISM (output 1).
r 21 t 21 is the time-averaged radial velocity, measured by comparing the particle radius at this output with its radius at a later simulation output (output 2). We choose the time spacing between outputs 1 and 2 to correspond as closely as possible to one quarter of a halo dynamical time. 5 This ensures that our selection criteria are capable of achieving converged answers with respect to the chosen temporal spacing of simulation outputs (see Appendix A1). Further to equations (5) and (6), we also select outflowing particles that have an instantaneous radial velocity greater than V max (at output 1). This catches (rare) cases where particles are feedback accelerated briefly to very high radial velocities but stall 6 before moving a significant distance out into the halo.
Equation (5) is our main criterion for selecting galaxy-scale outflows. It effectively demands that the particles will move outwards by at least one sixteenth of the virial radius within one quarter of a halo dynamical time. Equation (6) is a less stringent secondary criterion that helps to ensure that the particle has already joined the outflow by output 1 (from inspection of particle trajectories, we find that this is only relevant for galaxy-scale outflows).
Particles that leave the ISM/halo that are not selected as outflowing by the aforementioned criteria are added to a list of candidate wind particles that are then propagated down the halo merger tree on subsequent simulation outputs. These particles are re-tested against the same selection criteria at each subsequent simulation output until they either satisfy the criteria or three halo dynamical times have expired (at which point they are removed from the candidate wind list). This procedure ensures that particles that fluctuate over the ISM or virial radius boundary are accounted for in the outflow rate measurements should they be significantly accelerated while just outside the boundary. Including these particles has a negligible effect on outflow rates for lower mass galaxies (M 200 < 10 12 M ), but it does increase the outflow rates of high-mass galaxies appreciably, Figure 1. Mean mass outflow rates from the ISM (top panels) and haloes (bottom panels) for central galaxies, plotted as a function of halo mass. Outflow rates are quantified as a dimensionless mass loading factor (mean outflow rate over mean star formation rate, left-hand panels), and as a mean outflow rate per unit halo mass, scaled by the cosmic baryon fraction, f B ≡ b / m (right-hand panels). Different line colours correspond to different redshift intervals, as labelled, and mean fluxes and star formation rates are computed across all galaxies in each redshift/mass bin. Solid (dashed) lines indicate the halo mass range within which galaxies contain on average more (fewer) than 100 stellar particles. Indicative-power-law scalings for the mass loading factor are shown by the diagonal dashed black lines. and becomes the main contribution to galaxy-scale outflows for halo masses of M 200 > 10 13 M .
Our results are not highly sensitive to the exact values adopted for these selection criteria (as demonstrated in Appendix A3), although it is important to include some cut on time-averaged radial velocity. Fig. 1 presents the main results of this study, showing outflow rates for gas leaving the ISM (top panels) and the halo (bottom panels) of central galaxies. Data are taken from the 100 Mpc reference run, using trees with 200 snipshots. Unless otherwise stated, all subsequent results in this paper are shown for this simulation using these trees. Results are shown here as a function of halo mass; we refer readers interested in the dependence on more readily observable quantities to Section 3.2, where we show outflow rates as functions of stellar mass, star formation rate, and circular velocity. We focus on central galaxies to simplify the interpretation of outflows (which for satellites can also be caused by stripping by gravitational tides or gaseous ram pressure).

R E S U LT S
Following Neistein et al. (2012), the average measurements shown in Fig. 1 (and later figures) are taken by computing the mean of the numerator over the mean of the denominator, including all central galaxies recorded within the quoted redshift range. As demonstrated by Neistein et al. (2012), this approach yields the correct average mass exchange rate, in the sense that taking the time integral over the averaged inflow and outflow rates predicts the correct stellar masses of individual galaxies to within 0.1 dex (because the mean of the time derivative of the mass is equal to the time derivative of the mean of the mass). Taking the mean in this way also helps to average out the discreteness noise that would affect outflow rate measurements of individual galaxies if the numbers of outflowing particles and new stars formed between two simulation outputs are small. We indicate with dashed lines the halo mass range where galaxies contain on average fewer than 100 stellar particles, which we take as an indicator of the range in which galaxies are poorly resolved. We show in Section 4 that this approximately corresponds to the mass scale above which our results are reasonably well converged with respect to numerical resolution.
The left-hand panels of Fig. 1 show average outflow rates normalized by the average star formation rates computed over the same time interval (computed as the total mass of stars formed over the interval, ignoring mass-loss from stellar evolution). This quantity represents a time-averaged dimensionless mass loading factor, η, which can be considered as the efficiency with which outflows are launched from galaxies (top-left) and haloes (bottom-left). Parametric fits to the mass loading factors are provided in Appendix B.
Strong trends with halo mass are visible at both spatial scales, with a local minimum efficiency for outflows found at a halo mass of around M 200 ∼ 10 12 M , approximately independent of redshift. Below this characteristic halo mass, the galaxy-scale wind mass loading scales approximately as M −0.5 200 (the parametric bestfitting value of the exponent is −0.39 − 0.06 z), putting the EAGLE simulations somewhere in between the often considered momentum- Note that these scalings only are only strictly kinetic energy and momentum conserving if the outflow velocity scales linearly with the circular velocity of the system, which we show later is generally not the case for EAGLE. The corresponding mass loading scaling is typically steeper for the halo-scale outflows in the same mass range, with a best-fitting exponent of −1.19 + 0.18 z, matching the energyconserving scaling (∝M −2/3 200 ) by z ≈ 3. Note that the scaling steepens noticeably for the galaxy-scale mass loading in the mass range where more than 20 per cent of the galaxies are not forming stars (indicated by dashed lines). This change in scaling towards very low mass may therefore be related to resolution (and we typically exclude these mass bins from our analysis).
For M 200 > 10 12 M , the mass loading factors start to rise again due to the effects of AGN feedback (we show the explicit comparison with the no-AGN case in Section 3.8). The mass loading factor then declines slightly again for M 200 > 10 13 M for the galaxy-scale outflows, while the mass loading continues to rise monotonically with mass in high-mass haloes for halo-scale outflows for z < 1. Put together, it is clear qualitatively that the scaling of the mass loading factors with halo mass is at least partly responsible for the level of agreement between EAGLE  In the simplistic scenario where outflows alone set the scaling between stellar mass and halo mass, the basic expectation is that M ∝ η −1 M 200 , where η is the mass loading factor (Mitchell et al. 2016). Taking the example of the low-mass regime (where stellar feedback is typically assumed to dominate), empirical constraints indicate that the scaling between stellar mass and halo mass is 200 . This is a stronger dependence compared to what we find in EAGLE for galaxy-scale outflows, but it is consistent (particularly at lower redshifts) with the scaling we find for halo-scale outflows. This implies first that at the spatial scale of galaxies, additional sources of mass scaling must be at play in order to match the observed galaxy stellar mass function. The scaling of the haloscale outflows could in principle be a sufficient explanation (in that they reduce the available reservoir of baryons within the virial radius that can accrete on to the ISM). We defer a more quantitative analysis to a future study where we will present the corresponding picture for gaseous inflows, which is required to fully understand the predicted relationship between stellar mass and halo mass.
The right-hand panels of Fig. 1 show outflow rates without normalizing by the star formation rates, instead normalizing by halo mass to remove the zeroth-order mass scaling to compress the dynamic range. Starting with galaxy-scale outflows (top-right panel), it is interesting to note that the mass scale (M 200 ∼ 10 12 M ) where outflows are least efficient in terms of the mass loading factor is where outflows are most efficient in terms of the mass ejected per unit halo mass. This inversion serves to underline the aforementioned point that the scaling between stellar mass and halo mass is stronger than that between galaxy-scale outflow rate and halo mass, implying that there must be other reasons for the stellar-halo mass scaling. The picture changes markedly when considering instead the halo-scale outflow rates shown in the lower right panel of Fig. 1. The halo-scale outflow rates per unit halo mass are almost independent of halo mass for M 200 ∼ 10 10.5 -10 12.5 M , and for z < 1 even up to 10 14.5 M .
Differing degrees of redshift evolution at a fixed halo mass can be seen in each panel of Fig. 1. The galaxy-scale mass loading factor (top-left) decreases by about 0.5 dex between z = 3 and 0 for haloes of mass M 200 = 10 11 M . We note that the respective positive and negative scalings of energy injected by stellar feedback with gas density and metallicity (equation 3; see also fig. 1 of Crain et al. 2015) could contribute to this redshift evolution, as ISM densities/metallicities increase/decrease, respectively, with redshift at a fixed mass. Interestingly, the redshift dependence is reversed for the halo-scale mass loading factor (bottom-left panel), with the efficiency of halo-scale outflows per unit star formation growing towards low redshift. This presumably reflects an evolution of the properties of circumgalactic gas out to the virial radius. Another possibility is that halo-scale outflows are being driven by energy injected in the past, when star formation rates were higher.
Considering instead the outflow rates normalized by halo mass (right-hand panels) instead of by star formation rate, a trend of outflow rates increasing with increasing redshift is apparent for both galaxy-and halo-scale outflows. This primarily reflects the evolution of galaxy star formation rates at a fixed halo mass, which in turn is related to the slowing of structure formation towards low redshift that occurs in the CDM cosmological model. Indeed, if the outflow rates shown in the right-hand panels are multiplied by the age of the Universe for each redshift bin (in effect removing the redshift scaling of dark matter halo accretion rate), most of the redshift evolution disappears for the galaxy-scale outflows, and almost all of the redshift evolution disappears for the halo-scale outflows.

Comparing outflow rates at galaxy and halo scales
An important feature of the rates shown in Fig. 1 is that in general, substantially more mass is flowing out of the halo virial radius compared to that leaving the ISM. We show this explicitly in Fig. 2. At high redshifts (z > 3), the halo and galaxy-scale outflow rates are roughly equal for halo masses M 200 < 10 12 M (or for M < 10 10 M ). For z < 2, the halo-scale outflow rates evolve to become increasingly elevated over the galaxy-scale rates at lower redshift. The mass dependence becomes stronger at lower redshifts, with halo-scale outflows becoming increasingly elevated over galaxy-scale outflows in both low-mass and high-mass haloes, transitioning around a minimum elevation at M 200 ∼ 10 11.5 M .
All together, the enhanced outflow rates at the halo virial radius will play an important role in the EAGLE simulations, by effectively reducing the reservoir of baryons within the virial radius that can condense down on to the ISM, but without invoking outflow rates at the galaxy scale that are far too high relative to observational The ratio of halo-scale mass loading factor to galaxy-scale mass loading factor, plotted as a function of halo mass (top panel), and stellar mass within a 30 pkpc spherical aperture (bottom panel). Solid (dashed) lines indicate the halo mass range within which galaxies contain on average more (fewer) than 100 stellar particles. In general, substantially more outflowing mass is being removed from the halo than is being removed from the ISM. constraints (see Section 5.3). We explore the origins of the enhancement in the following parts of this section, culminating in the discussion presented in Section 3.7. The question of whether the mass loading factors at the two scales are qualitatively and quantitatively robust with respect to changing numerical resolution is discussed in Section 4. Fig. 3 shows galaxy-scale outflow rates as functions of stellar mass, M , star formation rate,Ṁ , and halo maximum circular velocity, V max , quantities that are more readily observable than halo mass. For outflow rates plotted as a function ofṀ , galaxies are binned according to the mass of stars formed within the last 100 Myr, comparable with the characteristic time-scale of SFR measurements derived from UV luminosities, but to be self-consistent the star formation rate folded into the mass loading factor is always taken from the mass of stars that formed within the same time interval used to measure the outflow rate. The stellar masses and star formation rates plotted along the x-axis are both measured using only star particles within a 30 pkpc spherical aperture. Parametric fits for the mass loading factor as a function of M and V max are given in Appendix B.

Outflow rates as functions of M ,Ṁ , and V max
While trends are similar to those seen in Fig. 1, several notable features do stand out in Fig. 3. While the scaling of galaxy-scale outflow rates plotted as a function of halo mass (upper right panel in Fig. 1) or maximum circular velocity (bottom-right panel in Fig. 3) show a characteristic change in slope around M 200 ∼ 10 12 M or V max ∼ 125 km s −1 , such a change is much less evident in the scaling of outflow rate with stellar mass (top-right panel Fig. 3). This difference reflects in combination the mass scaling of the mass loading factor, the dependence of star formation rate per unit stellar mass on stellar mass (see fig. 5 in Furlong et al. 2015), and the underlying scaling of galaxy stellar mass on halo mass (see fig. 8 in Schaye et al. 2015).
Another feature visible in Fig. 3 is that the negative scaling of the mass loading factor with star formation rate (middle-left) does not flatten or turn over for high star formation rates, unlike for all of the other variables considered. This reflects the strong decrease of galaxy star formation rates per unit stellar mass in massive galaxies (where AGNs power most of the outflow and so change the mass scaling of the mass loading factor; see Section 3.8), such that massive galaxies do not dominate the highest star formation rate bins.

Outflow velocities
While the main focus of this study is on outflow rates, it is also interesting to explore the decomposition of these gas flows as a function of velocity, or gas phase. We defer a detailed analysis to future work, but we do show here the average flux-weighted velocity of outflowing gas in Fig. 4. The median velocities (top panel) exhibit roughly logarithmic scaling with halo mass. Outflowing gas that was ejected from the ISM moves at higher velocities relative to all outflowing gas at a given radius, and exhibits a peak velocity at a characteristic halo mass of 10 12 M at z = 0. This effect is more pronounced for the 90th percentile of the flux-weighted outflow velocity (bottom panel). Except for the scaling of median velocity with halo mass in low-mass haloes (M 200 < 10 12 M ), the scaling of outflow velocity is qualitatively different to the scaling of maximum halo circular velocity with halo mass (shown by the dotted lines). The spread in velocities at a given mass/redshift is large (as can be appreciated by comparing the two percentiles). Outflow velocities at a given halo mass are higher at higher redshifts, with the exception of v 90 around the peak at M 200 ∼ 10 12 M .

Energy and momentum fluxes
While the mass loading factor of galactic winds is one measure of their efficiency, it is also interesting to assess the wind efficiency in terms of energy and radial momentum. Fig. 5 shows measurements of the fluxes of energy (kinetic plus thermal) and momentum, contrasted with the rate of thermal energy injection by feedback processes (Ė inject ). While zero momentum is injected by hand in the simulation, we can define an effective momentum injection rate as p t = 1 t 2 E inject M heated , where E inject of energy is directly injected into M heated of mass over a time interval t. This represents the momentum that the wind would achieve if all of the injected thermal energy were converted to kinetic form, and should be regarded as a rule of thumb rather than as the true momentum that winds are expected to attain.
The top-left panel of Fig. 5 shows the energy flux of outflowing gas close to the galaxy (solid lines), normalized by the kinetic energy that The left-hand panels show the average mass loading factor plotted as a function of different variables, and the right-hand panels show the average outflow rate. Solid (dashed) lines indicate bins within which galaxies contain on average more (fewer) than 100 stellar particles.
would be required to move the entire baryonic content of the halo at the halo circular velocity, V c , assuming that the baryon to dark matter content of the halo matches the universal fraction, f B . At high redshifts, more than sufficient energy is being injected to achieve this within a Gyr, but this is no longer the case at low redshift once the rates of star formation and SMBH accretion have slowed at a fixed halo mass. The upper right panel shows the ratio of the energy flux to the feedback energy injection rate, both close to the galaxy (solid lines) and at the virial radius of the halo (dashed lines). While these measurements are noisier than for the mass loading factor, 7 the . Velocities are measured in spherical shells at a radius 0.1 < r /R vir < 0.2. The median relationship between maximum halo circular velocity, V max , and halo mass is also shown (dotted lines). Outflow velocities are only shown for halo mass bins where more than 80 per cent of the galaxies have non-zero flux within the shell. Median outflow velocities are in general comparable to V max for M 200 < 10 12 M , but saturate (or even decline in some cases) for high-mass haloes. At low (M 200 ≈ 10 11 M ) and high (M 200 > 10 13.5 M ) halo masses, the outflows can carry more energy than is being injected. This serves first to underline that the energy loading factors plotted are upper limits to the efficiency with which the injected energy from feedback is able to power galactic winds. Other sources of energy in outflowing gas include the ultraviolet background (UVB; which could plausibly be responsible for the greater than unity energy loading measured for outflows at the virial radius in low-mass haloes) and gravitational heating (which could plausibly have a larger relative effect in massive haloes, where pressurized hot coronae have developed). Another factor is that the energy/momentum fluxes at the halo virial radius are associated with feedback events that occurred earlier in the history of each galaxy, at which time the star formation and SMBH accretion rates may have been significantly different. We return to this point in Section 3.7.
For intermediate-mass haloes, the energy in outflows close to the galaxy is typically higher than that for outflows close to the virial radius, likely indicating dissipation over the intervening scales. This is less apparent when comparing the momentum flux at the two scales, and by z = 0 the momentum flux is higher at the virial radius than near the galaxy over the entire halo mass range probed (other than the handful of haloes in the highest mass bins). This indicates some level of entrainment of mass at fixed energy, which is consistent with the enhanced mass loading at the virial radius seen in Fig. 2.

Outflows as a function of radius
Entrainment of outflowing mass is shown more directly in Fig. 6, which shows the mass, momentum, and energy fluxes as a function of radius for haloes of mass 12 < log 10 (M 200 /M ) < 12.2 for redshifts 0 < z < 0.3. In this instance, we separate the contribution from gas that has been removed from the ISM (dashed lines), versus gas that has never been in the ISM (dotted lines). Mass flux (top-left panel) is conserved as a function of radius for the former ISM material, but by 0.2 R vir there is a similar mass flux of material that was never in the ISM, and the contribution of this component rises until it dominates the mass flux at the virial radius. A similar picture is seen for the momentum flux (top-right panel).
The total energy flux (solid black line in the bottom-left panel) is approximately constant with radius, with energy seemingly being exchanged from the former ISM component (dashed black line) to gas entrained from the CGM (dotted black line) as outflows propagate outwards. Despite the feedback scheme employed in EAGLE being thermal, the majority of the outflowing energy flux is in kinetic form close the galaxy, but the majority of the energy flux is in thermal form at larger radii. Correspondingly, the mass flux-weighted velocities (bottom-right panel) decline as a function of radius.
Overall, the trends are consistent with a picture whereby gas is entrained on circumgalactic scales, explaining much of the difference between the halo and galaxy-scale outflow rates shown in Fig. 2. A similar picture is seen at lower halo masses at low redshift (not shown), although in that instance the total energy flux actually rises with radius, indicating that another source of energy is involved (possibly the UVB). The picture is again similar at higher halo masses, but in this case the entrainment phenomenon ceases once the outflow reaches half the halo virial radius, thermal energy is more dominant over kinetic energy, and the fractional contribution to the energy flux from outflowing material that has never been in the ISM is higher at the centre. At higher redshifts, the trends are similar but there is systematically less evidence for entrainment, as the mass flux increases much less strongly with radius (as seen also in Fig. 2).

Directionality
The top panel of Fig. 7 shows an example of the angular dependence of galactic outflows in EAGLE. Following the approach used in Nelson et al. (2019), we show mass flux as a function of radius (x-axis) and 'galactocentric' angle, which we take as the angle between each gas particle and the major axis of the galaxy, as viewed in an edgeon projection (defining the mid-plane using the angular momentum vector of the ISM). For a disc galaxy, values of 0 and ±π therefore indicate outflows that are propagating within the plane of the disc, and values of ±π /2 indicate outflows that propagate orthogonally to the disc along the minor axes. The increase of mass flux with radius 2). These can be compared to the input thermal energy injection rate or (pseudo-)input momentum from stellar (dashed lines) and AGN feedback (dotted lines). Because feedback in EAGLE is purely thermal, the (pseudo-)input momentum rate is defined relative to the thermal energy injection rate asṗ = 1 t 2 E inject M heated , where E inject is the energy that is directly injected into M heated of mass over a time interval t (see the main text). Fluxes and injection rates are normalized by the characteristic energy/momenta of the associated haloes. Right-hand panels: fluxes of outflowing gas divided by the corresponding energy/(pseudo-)momentum injection rates, defining effective energy or momentum loading factors. Loading factors are shown for outflowing gas in shells at 0.1 < r/R vir < 0.2 (solid lines), and at 0.9 < r/R vir < 1.0 (dashed lines). In all the panels, data are only shown for mass bins within which galaxies contain on average more than 100 stellar particles. Gas within the ISM is excluded from the flux measurements. Data are taken from the 50 Mpc reference run. Roughly 20 per cent of the energy being injected by feedback is retained in outflows in EAGLE for M 200 ∼ 10 12 M , with this fraction increasing for both higher and lower halo masses.
shows once again the previously discussed entrainment effect. The distribution with angle shows that EAGLE produces a bimodal outflow pattern, aligned with the minor axes, which reflects the relative ease with which outflows can escape the ISM (and propagate through the CGM) in the directions orthogonal to the disc. Note that at r < 0.2 R vir there is also some outflowing flux aligned with the disc, which we interpret as a combination of ISM material (which was not subtracted here) and gas that is settling around the disc after infall.
The bottom panel of Fig. 7 shows as a function of radius the fraction of the virial sphere that is occupied by gas that is on average outflowing with v rad > 0.25V max . The fraction rises from ≈20 per cent at the halo centre up to a peak value close to 40 per cent at r = 0.3R vir , and stays nearly constant out to larger radii. The enhancement of the mass flux with radius is therefore not associated with an increase in the solid angle of the outflow for r > 0.3R vir .

Energy-driven winds and travel-time effects
Having presented information on the mass, momentum and energy fluxes, velocities, and directionality of galactic outflows, we can now put this together to discuss the origin of the enhancement in mass flux with radius (out to the virial radius) seen in EAGLE. We stress at the outset that this is a question that is complicated to address in cosmological simulations because of evolution effects: Galaxies and haloes can grow significantly both in mass and size over time-scales that are comparable to the time-scales for circumgalactic gas flows, the velocity field of circumgalactic gas will also reflect cosmological infall, and the energy and momentum content of circumgalactic gas at different scales will reflect the cumulative injection of feedback energy over a range of time-scales. We can none the less examine some simplified arguments, which we present here.
An obvious mechanism to increase mass flux with radius comes from an 'energy-driven' wind scenario, in which the outflows are overpressurized relative to the ambient ISM or CGM, which generates radial momentum as the hot interior of the outflow does PdV work on the surrounding gas. This is the physical mechanism responsible for increasing the radial momentum of an outflow during the Sedov-Taylor (adiabatic) and pressure-driven snowplow phases of supernova explosions. It has also been discussed within the context of larger scale AGN-driven winds (e.g. King, Zubovas & Power 2011;Faucher-Giguère & Quataert 2012), which have been demonstrated to be capable of driving an increase of mass flux with radius on circumgalactic scales in full cosmological simulations (Costa, Sijacki & Haehnelt 2014).
The top panel of Fig. 8 shows the radial profile of the median thermal pressure, averaging over galaxies with 10 12 < M < 10 12.2 at z = 0.25. We show the estimates of the average thermal pressure of all gas at a given radius (solid/dotted lines, corresponding to mass/volume weighting) as well as the flux-weighted average thermal pressure of outflowing gas with v rad > 0.25V max (dashed line), which we take as a measure of the characteristic thermal pressure within feedback-driven winds. To compute the weighted average of P gas , we weigh by m gas P gas for mass-weighted, (m gas /ρ gas ) P gas for volumeweighted, and (m gas v rad ) P gas for flux-weighted, where m gas is the gas particle mass, ρ gas is the SPH density, v rad is the radial velocity, and P gas is the SPH pressure. We average over spherical shells of width 0.1 R vir , including particles whose centres are inside the shell.
We find that the outflowing gas is overpressurized relative to the total CGM by ≈0.3dex at all radii, which will drive an increase of momentum with time and distance for discrete outflow events as they propagate through the ambient CGM. Referring back to Fig. 6, which shows that the energy flux is roughly constant with radius for this mass/redshift range, it therefore appears that winds are driven across the CGM in an energy-driven configuration. Fig. 6 also shows that the average wind velocity is nearly flat with radius (more precisely, it is slightly declining), which implies that the increase in mass flux must be associated with entrainment of ambient gas, not with an increase in the characteristic velocity of the outflow with radius. The bottom panel of Fig. 7 shows that this entrainment is not associated with an increase in the solid angle occupied by winds as a function of radius. Putting the information from these measurements together implies then that the mass per unit radius in outflows must increase with radius, which is shown explicitly to be the case in the bottom panel of Fig. 8 (see also Schaller et al. 2015b, for a focused analysis of density profiles in EAGLE). With a spherically averaged density profile that is shallower than isothermal [for which, ρ(r) ∝ r −2 ], Figure 7. Directionality (top) and spherical solid-angle covering fraction (bottom) of outflowing gas, plotted as a function of radius, averaging over galaxies with 10 11.75 < M < 10 12.25 at z = 0.25. Top: directionality is defined in terms of the galactocentric angle (see the main text), with values of ±π /2 indicating that gas is flowing orthogonally to the disc (minor axes). The colour scale indicates the radially outflowing mass flux, as labelled, selecting only gas particles with radial velocities v rad > 0.25V max . For r > 0.1 R vir , EAGLE produces a clear bimodal outflow distribution that aligns with the galaxy minor axis. Bottom: the solid angle fraction of the sphere that is covered by solid angle bins within which the net flux-weighted radial velocity satisfies v rad > 0.25V max . For r > 0.2 R vir , about 40 per cent of the virial sphere is occupied by outflowing gas. A colour version is available online. most of the mass in the CGM is in the outer regions of the halo in EAGLE, which helps to explain why the entrainment effect is seen at larger scales.
The middle panel of Fig. 8 shows an estimate of the typical radial acceleration imposed by the radial thermal pressure gradient for gas in winds (dashed red line), and for all gas (solid red line). This is contrasted against the (opposite-sign) gravitational radial acceleration (black line). The CGM is on average undersupported against gravitational infall, and will therefore tend towards a net inflow in the absence of additional sources of inflow/outflow from the ISM and from beyond the virial radius (see Oppenheimer 2018, for a generalized discussion of hydrostatic balance in the EAGLE Figure 8. Radial profiles of thermal pressure (top), radial acceleration (middle), and mass (bottom), averaged over galaxies with 10 12 < M < 10 12.2 at z = 0.25. Top: median thermal pressure profiles for all gas, both mass-weighted (solid line) and volume-weighted (dotted line), as well as the flux-weighted median for outflowing gas with v rad > 0.25V max (dashed line), which we take as a measurement of the characteristic thermal pressure within winds. Winds are on average overpressurized relative to the typical CGM at a given radius, which will drive an increase in momentum as outflows propagate outwards. Middle: median radial acceleration due to the negative thermal pressure gradient (red lines) for all gas (mass-weighted, solid line) and for outflowing gas with v rad > 0.25V max (flux-weighted, dashed line). Also shown is the (inwards) radial acceleration associated with the gradient of the gravitational potential, assuming spherical symmetry. Gas in winds is on average pressure supported against gravitational acceleration for r > 0.4 R vir , but on average gas within the CGM is undersupported at all radii. Bottom: radial mass profiles for total gas (solid), all outflowing gas (dotted), and outflowing gas with v rad > 0.25V max (dashed). Most of the mass in the CGM is in the outer regions of haloes, helping to explain why wind entrainment acts on large spatial scales. A colour version is available online. Figure 9. Top: the mean crossing time for outflowing particles ejected from the ISM to reach the virial radius (recorded at the time that particles reach the virial radius), plotted as a function of halo mass. Bottom: the average change in star formation rate (compared to the main progenitor galaxy) computed over the crossing time-scale shown in the top panel. Since galaxy star formation histories (and AGN activity) in EAGLE on average peak at z ≈ 2, the outflow rate at the virial radius at z = 2 reflects in part the lower star formation/AGN activity in past progenitors, whereas z = 0 outflows at the virial radius reflect the higher past star formation/AGN activity. simulations). Gas within winds is pressure supported against the gravitational acceleration for r > 0.4 R vir , explaining why the radial velocity of outflows (Fig. 6) is nearly flat in this radius range.
Another effect that turns out to be important for understanding the statistical trend of elevated mass fluxes at the virial radius is connected to the finite time taken for energy injected into the ISM to propagate outwards to the virial radius. Comparing the outflow rates at the virial radius to the rates of gas leaving the ISM (Fig. 2), the ratio of the two is of order unity for M 200 ∼ 10 12 M at z ≈ 2 but has increased by ≈0.5 dex by z = 0. The top panel of Fig. 9 shows the mean crossing time for outflowing gas to move from the ISM to the virial radius. This time is not negligible compared to a Hubble time, and outflows at the virial radius will presumably at least partly reflect the energy injection rate at the time outflows were launched from the ISM.
The bottom panel of Fig. 9 then shows how the present star formation rate (at the time the selected particles are leaving the halo) compares to the star formation rate in the past when the particles left the ISM. Due to the shape of star formation histories in EAGLE that peak at z ≈ 2 (see e.g. fig. 9 of Mitchell et al. 2018a), star formation rates (and also AGN activity) were higher in the past for the progenitors of galaxies at z = 0 (black line), but were lower in the past for the progenitors of galaxies at z = 2 (cyan line). While this partially helps to explain the elevated outflow rates at the virial radius at z = 0, the magnitude of the effect is too small to be the main explanation. Time delay effects do, however, present a convincing explanation for the redshift dependence of the ratio of mass loading factors seen in Fig. 2. The offset between the star formation increase/decrease over a crossing time is around 0.3 dex between z = 2 and 0, which is comparable to (and goes in the right direction to explain) the redshift evolution of the mass loading ratio shown in Fig. 2. In addition, the change in star formation rate is more stark for high-mass haloes with M 200 > 10 12 M (due to a longer crossing time), which helps to explain the mass dependence of the mass loading ratio in this mass range.
There are other factors that may contribute to the change in outflow rate with radius in EAGLE, which we now briefly consider. One is that satellites may play a role by injecting energy directly at larger radii. We have checked this explicitly, and find that the energy injection rate is generally negligible compared to the central energy injection rate, and to the energy flux at each radius. A second physical effect that has been discussed recently within the context of stellar feedback-driven outflows is buoyancy (Bower et al. 2017;Keller, Kruijssen & Wadsley 2019), which has also long been considered as an important part of how AGN feedback may operate in the intra-cluster medium (e.g. Churazov et al. 2002;Chandran & Rasera 2007;Pope et al. 2010), albeit generally with additional physics to what is simulated in EAGLE (e.g. cosmic rays and thermal conduction). Since we find in Fig. 8 that outflows in M 200 = 10 12 M mass haloes are overpressured relative to the ambient CGM, we do not expect buoyancy to be the main driver, as buoyancy becomes dynamically important as a mechanism to lift low-entropy gas within a multiphase medium that is locally in pressure equilibrium. This situation may change in higher mass haloes however (M 200 ∼ 10 13 M , not shown), for which winds are still overpressured, but the ambient medium itself is in overall equilibrium with the gravitational potential. Finally, non-feedback related energy sources could in principle act on larger spatial scales to drive elevated mass fluxes at the virial radius. While non-trivial to check, the naive expectation is that gravitation-related motions would peak for gas moving near the halo centre, where the maximum amount of potential energy has been converted into kinetic form. On the other hand, compressive heating associated with gravitational infall (and in particular halo mergers) could overpressurize the CGM and drive large-scale outflows in the same manner as previously described for feedback-driven winds. From looking at individual outflow events in time series (not shown for conciseness), we find that significant outflow at the virial radius is always preceded by an intense but short-lived outflow event at the halo centre, triggered by a period of star formation or AGN activity. This confirms that the large-scale outflows are at least correlated with feedback activity, but on the other hand star formation and feedback will also be correlated with gravitational infall and mergers, so we do not draw any firm conclusions.  zero at z > 2 up to about 40 per cent by z = 0. AGNs provide the majority of energy injection for haloes more massive than 10 12 M at all redshifts recorded.

Impact of AGN feedback
Below z = 5, a strong feature appears at a characteristic halo mass of 10 10 M . This feature arises because of the implementation of SMBH seeding in EAGLE; black hole seeds are placed in FoF groups of that mass. The sudden increase in AGN energy at this specific mass scale is clearly artificial, with the newly formed black hole strongly out of equilibrium with the surrounding ISM. We have checked and verified that this feature has a negligible effect on the median stellar mass as a function of halo mass, by comparing simulations with and without AGN feedback. Fig. 11 compares the outflow rates in simulations with and without AGN feedback. We perform this comparison in terms of mass loading factors to account for the difference in star formation activity between the two simulations at a fixed halo mass. For the galaxy-scale outflows (top panel), AGN feedback is clearly responsible for the upturn in the mass loading factor for haloes with M 200 > 10 12 M . A similar picture emerges for the halo-scale outflows (bottom panel).

N U M E R I C A L C O N V E R G E N C E
As per the results and discussion presented in Schaye et al. (2015), the basic outputs of the EAGLE simulations (e.g. the galaxy stellar mass function; see their fig. 7) are not converged with numerical resolution for a fixed set of model feedback parameters, primarily because the anticipated radiative losses depend on the distribution of ISM densities, which itself changes with numerical resolution. Schaye et al. (2015) argue that this convergence test (dubbed 'strong numerical convergence') is overly stringent for cosmological simulations in which the ISM is unresolved. Because the subgrid parameters of such simulations in any case require calibration, they instead introduce the concept of 'weak numerical convergence', for which the change in radiative losses associated with changing resolution is accounted for by adjusting the efficiency of feedback parameters until agreement with the basic observables used for Figure 11. Impact of AGN feedback on mass loading factors associated with galaxy-scale (top) and halo-scale (bottom) outflows. Solid lines indicate outflow rates for the reference simulation (which includes AGN feedback). Dashed lines indicate the corresponding rates for the no-AGN variant of the reference simulation. Data are taken from the 50 Mpc reference and no-AGN runs, both using trees with 28 snapshots. Data are shown for mass bins within which galaxies contain on average more than 100 stellar particles. AGN feedback starts to appreciably affect outflow rates in haloes with masses M 200 > 10 11.5 M , causing a flattening (or upturn) of the scaling of the mass loading factor with increasing halo mass. calibration is (re)achieved. While clearly less demanding than a conventional ('strong') convergence test, a weak convergence test is still of significant utility, for example to identify the mass scales at which non-convergence of the quenched fraction of galaxies is being driven by sampling issues (e.g. too few star particles), rather than by purely feedback-related issues (Furlong et al. 2015).
It is important then to check if outflow rate scalings change, while still (as far as possible) retaining agreement with the observed galaxy stellar mass function. Since we expect the galaxy stellar mass function to primarily reflect the balance between gaseous inflows and outflows, the naive expectation is that a 'weakly' converged pair of simulations with different resolutions (both calibrated to reproduce the observed galaxy stellar mass function) would produce similar outflow scalings. Fig. 12 compares the mass loading factor (left-hand panels) at ISM (top) and halo (bottom) scales between two EAGLE simulations with equivalent volume (25 3 Mpc 3 ): the first with the reference EAGLE resolution and parameters, and the second with eight times higher mass resolution and recalibrated parameters, which we Bottom-right: energy flux for outflowing gas at r = 0.95 R vir . Dark (light) lines indicate the halo mass range within which galaxies contain on average more (fewer) than 100 stellar particles. Mass loading factors are reasonably well converged at the halo scale for M 200 > 10 11 M , but are not quantitatively converged at the ISM scale for M 200 < 10 12 M , in the stellar feedback regime. Energy fluxes are generally converged at both scales for galaxies with more than 100 stellar particles. refer to as the 'Recal' simulation. Mass loading factors are significantly higher at all redshifts in the higher resolution Recal simulation for M 200 < 10 12 M at the ISM scale, and for M 200 < 10 11 M at the halo scale. (Weak) convergence appears better at higher masses, although statistics are too sparse to make a robust conclusion for group and cluster-mass haloes. Since haloes of M 200 ≈ 10 10.8 M contain on average only 100 stellar particles in the reference model at standard EAGLE resolution, we conclude that there is reasonable weak convergence for mass loading factors at the halo scale for resolved haloes, but not at the ISM scale for M 200 < 10 12 M .
The right-hand panels of Fig. 12 show a comparison of the energy (thermal plus kinetic) of outflowing gas in the two simulations. Outflow energetics are much better converged than mass loading factors at the ISM scale, showing only a significant discrepancy at the halo scale at high redshifts. Given that the Recal model is calibrated against the same observed stellar mass function as the reference model run at lower resolution, this implies that outflow energetics are a better indicator of the efficiency of feedback in regulating galaxy growth. Furthermore, since convergence is better for mass loading factors at the halo scale, we can also infer that galaxy formation in the simulation is being regulated primarily on CGM scales, as a consequence of the work done by energy injected into the CGM by feedback (this interpretation aligns with the analysis of Davies et al. 2019, who find that the CGM mass fraction strongly correlates with the star formation rates in galaxies in EAGLE). This regulation is achieved by shaping inflow rates of gas on to galaxies, which in an upcoming study we will show are higher in the higher resolution Recal simulation (both for recycled and first-time infalling gas), which explains how the simulation produces the same galaxy stellar mass function despite producing different mass loading factors at the ISM scale.
Given the recent focus with cosmological simulations on the question of convergence with numerical resolution in the CGM (for column densities, ionization state, etc.; Hummels et al. 2019;Peeples et al. 2019;Suresh et al. 2019;van de Voort et al. 2019), we briefly mention the possible implications of this for the results presented here. While something that has not been explicitly studied to our knowledge in cosmological simulations, it seems probable that outflows could be affected by CGM resolution, as this will (for example) affect levels of mixing with ambient gas via instabilities. The convergence test we present here is suggestive, in that we find higher inflow and outflow rates in the CGM for the same outflowing energy flux (implying feedback is less effective at disrupting infall at higher resolution), but is inconclusive in that the outflowing mass fluxes also change at the scale of the ISM, before any interaction with the CGM can occur.
To summarize, we find that quantitatively the mass loading scalings in EAGLE are reasonably well converged at the halo scale over the mass range where galaxies are resolved by more than 100 stellar particles at standard resolution (M 200 > 10 11 M ), once feedback parameters are recalibrated against observational constraints. Quantitative convergence is not achieved at the ISM scale for M 200 < 10 12 M , but qualitatively the picture for outflows in EAGLE remains the same at the higher resolution explored: the mass loading factor scales strongly with halo mass, with a minimum value at M 200 ∼ 10 12 M , and outflow rates are elevated at the virial radius compared to at the boundary of the ISM, especially at low redshifts.

L I T E R AT U R E C O M PA R I S O N
Here, we conclude our analysis of outflows by comparing to a range of models, simulations, and observations from the literature, and explore the conclusions that can be drawn from this wider context.

Comparison to semi-analytic models
Semi-analytic models are an established method to study the evolution of galaxies within the full cosmological context (see Baugh 2006;Somerville & Davé 2014, for an overview). Most semi-analytic models assume that stellar feedback drives galactic outflows from the ISM of galaxies, with a mass loading factor that scales negatively with galaxy circular velocity (e.g. Kauffmann, White & Guiderdoni 1993;Cole et al. 2000). This in turn allows the models to achieve a match with the faint end of the galaxy luminosity function (e.g. Benson et al. 2003). 8 Our measurements of outflow rates from EAGLE are (deliberately) suitable for direct comparison to the prescriptions assumed in semi-analytic models, and we show a direct comparison to a subset of recent models from the literature in Fig. 13.
It is immediately apparent from Fig. 13 that there is an enormous dispersion in what is assumed for the mass loading factor from one model to another (up to nearly four orders of magnitude at a given halo mass), despite the fact that all the models shown are calibrated to reproduce the observed distribution of stellar mass. Focusing only on the normalization, the large differences in mass loading factor are driven by two factors. First, each model makes different assumptions regarding the level of dichotomy between outflow rates of gas leaving the ISM (solid lines) versus the halo virial radius (dashed lines). The Henriques et al. (2015), Hirschmann et al. (2016), and Lagos et al. (2018) models (all at least partially adapted from the L-galaxies model of Guo et al. 2011) prescribe the excess energy remaining in galactic winds after they have escaped the ISM, and assume this energy can drive even greater amounts of gas out of the halo. Conversely, the GALFORM and Santa Cruz models assume that the amount of gas ejected from the halo is equivalent (or less than that for the Santa Cruz model) to the amount of gas ejected from the ISM (e.g. Mitchell et al. 2018a;Somerville et al. 2008). Both scenarios are degenerate in terms of stellar mass assembly, in the sense that they both reduce the fraction of baryons that form stars.
The second explanation for the differences in mass loading normalization stems from the assumed efficiency of recycling of Figure 13. A comparison of mass loading factors between EAGLE and a set of semi-analytic models from the literature, plotted as a function of halo mass at z = 0. Semi-analytic models shown include specific implementations of the GALFORM (Mitchell et al. 2018a), L-Galaxies (Henriques et al. 2015), Santa Cruz (Somerville, Popping & Trager 2015), SHARK (Lagos et al. 2018), and GAEA (Hirschmann, De Lucia & Fontanot 2016) models. Outflow rates are plotted for gas being ejected from the ISM (solid lines), and for gas being ejected from the halo virial radius (dashed lines). Each line colour corresponds to a given model, as labelled. For EAGLE, dark (light) lines indicate bins within which galaxies contain on average more (fewer) than 100 stellar particles. All of the semi-analytic models shown (and EAGLE) are tuned to match the local galaxy luminosity function and/or the galaxy stellar mass function. This agreement can apparently be achieved with wildly different scenarios for how much gas outflows from the ISM, and from the halo, emphasizing the deeply degenerate nature of galaxy evolution if only stellar mass constraints are considered. A colour version is available online. ejected wind material. For example, the GALFORM model assumes a very efficient recycling time-scale that is of the order of the halo dynamical time (such that ejected gas returns in only 10 per cent of a Hubble time), whereas the Santa Cruz model assumes that gas returns over a Hubble time. This forces the former model to invoke mass loading factors that are much larger than the latter. Again, these scenarios are degenerate in terms of stellar mass assembly (e.g. Mitchell et al. 2014), at least up until the point that the recycling time-scale becomes so long that galaxy clusters no longer retain the universal baryon fraction (Somerville et al. 2008).
Given this (long-standing) impasse, it is then interesting to consider the picture emerging from modern hydrodynamical simulations. The full simulation picture is shown in Section 5.2, but we choose to show the direct comparison between semi-analytic models and EAGLE here. The outflow rates from EAGLE (blue lines) are qualitatively closer to the scenarios presented by the GAEA (red lines; Hirschmann et al. 2016), SHARK (cyan lines; Lagos et al. 2018), and L-galaxies (black lines; Henriques et al. 2015) models, in that significantly more gas is ejected from halo virial radii than from the ISM. Quantitatively, however, EAGLE differs significantly in both normalization and slope with the L-galaxies model shown. Hirschmann et al. (2016) adopt a mass loading prescription for gas leaving the ISM inspired by the FIRE simulations (Hopkins et al. 2014), as measured by Muratov et al. (2015). Qualitatively, the picture from this model is close to that seen in EAGLE at z = 0, with a relatively low normalization and fairly shallow scaling of the galaxy-scale mass loading factor, combined with a significantly higher normalization for the outflow rates at the halo virial radius. We present a direct comparison with FIRE and other hydrodynamical simulations in the following section.
Finally, we note that the mass loading factors shown for the semianalytic models are for stellar feedback only. The upturn in mass loading factors for high-mass galaxies in EAGLE is caused by AGN feedback. Most semi-analytic models assume that AGN feedback acts only to suppress inflows rather than drive AGN outflows directly, 9 which is qualitatively different from the scenario presented in EAGLE. We note that semi-analytic models where AGNs do eject baryons from haloes have been considered as an explanation for the observed X-ray luminosity of galaxy groups (Bower, McCarthy & Benson 2008;Bower, Benson & Crain 2012).

Comparison to other cosmological simulations
Fig. 14 presents an overview of the mass loading factors from recent cosmological hydrodynamical simulations. Each study shown uses a different method to measure outflow rates, and we have taken care to (as far as is reasonably possible) compare EAGLE to other simulations using equivalent measurements.
The upper panels of Fig. 14 compare EAGLE to the 50 3 Mpc 3 Illustris-TNG (TNG-50) simulation at z = 2, taking measurements from Nelson et al. (2019). Nelson et al. (2019) measure outflow rates in shells at a given physical distance from the halo centre, for gas radially outflowing faster than some minimum radial velocity cut (different line styles in the upper left panel show different cuts). These simple criteria are straightforward to implement, and so we can perform a like-for-like comparison of the simulations at z = 2 (the redshift focused on by Nelson et al. 2019). Taking all outflowing gas with v r > 0 km s −1 at a distance of 10 kpc (solid lines in the top-left panel), EAGLE and TNG-50 display qualitatively similar behaviour for stellar masses M < 10 10.5 M , but are offset in normalization by up to 0.5 dex, with higher mass loading factors in TNG-50 than in EAGLE.
Mass loading for stellar feedback is set by hand at injection for TNG-50 (shown as the dotted red line), with outflows seeded by wind particles that are decoupled from the hydrodynamical scheme until they reach a density below n H ∼ 0.005 cm −3 (Pillepich et al. 2018). In practice, the TNG outflows generally recouple within 10 kpc, and the mass loading factors compared with EAGLE here were measured on scales at which the direct contribution of decoupled particles to the outflow is negligible ). In addition, outflows at these scales in TNG may have to do significant work against the magnetic pressure of circumgalactic gas, which is not accounted for in EAGLE. The TNG mass loading at injection (minus a residual metallicity dependence) is set to scale negatively with circular velocity as V −2 c . Although the measured outflow rate is slightly higher than the injected one, they track each other closely at low mass, where stellar feedback dominates over AGN feedback . No mass loading factor is imposed by hand in EAGLE, but the feedback model is still calibrated against similar observational constraints to those used for TNG, and so it is encouraging (but not surprising) to see that the mass scaling of the mass loading factor is similar between the simulations in the stellar feedback-dominated regime.
At higher stellar masses, Nelson et al. (2019) report a strong upturn in the mass loading factor that is attributed to AGN feedback. A weaker upturn for galaxy-scale outflows at 10 13 M haloes is seen in EAGLE in Fig. 1, but is not visible using the shell-based measurements at 10 kpc, where the mass loading instead flattens at high stellar masses. The upper right panel of Fig. 14 compares shellbased outflows at different radii, and here a clear upturn in the mass loading is visible in EAGLE at a distance of 50 kpc from the halo centre (dotted blue line), similar to that seen in TNG-50 at all radii. This indicates a significant difference in the smaller scale wind launching for AGN feedback between the simulations, with TNG-50 ejecting large amounts of gas from the centre of massive galaxies, while EAGLE launches relatively little gas but with the wind seemingly continuing to load mass as a function of radius, such that the mass loading increases out to the virial radius (dash-dotted blue line).
Comparing the mass loading in the stellar feedback regime in the upper right panel of Fig. 14 reveals further stark differences between the two simulations. While TNG-50 ejects significantly more gas per unit star formation than EAGLE at 10 kpc in low-mass galaxies, the outflows seem to decline strongly as a function of radius in TNG-50. Outflows behave differently in EAGLE, with mass loading that either stays roughly constant with, or grows with, radius. As such, the mass loading factor at 50 kpc is about 0.5 dex higher in EAGLE for galaxies of stellar mass M ∼ 10 9 M . 10 This difference implies that there is likely a large difference in the efficiency of recycling of ejected wind material between EAGLE and TNG-50 (with recycling being a more important source of inflows in TNG-50 than in EAGLE), which presumably affects the observable properties of the CGM as a function of impact parameter from galaxies. Davies et al. (2020) find a very consistent picture by comparing the total baryon content of haloes between EAGLE and Illustris-TNG, which they show is much higher in TNG than in EAGLE at low mass.
The middle-left panel of Fig. 14 compares outflow rates in EAGLE with the FIRE zoom-in simulations (introduced in Hopkins et al. 2014). Relative to EAGLE, the FIRE simulations employ significantly higher mass and spatial resolution (with the improvement scaling negatively with the mass of the targeted haloes), allow a cold ISM phase to form without imposing a temperature floor, and implement a more explicit representation of stellar feedback (separating contributions from radiation, stellar winds, and type II supernova explosions). The FIRE simulations do not include AGN feedback. We show the best-fitting relation to the FIRE simulations at z = 0.25 from Muratov et al. (2015), measured using shells at one quarter of the halo virial radius (dashed magenta line). Mimicking this type of measurement in EAGLE (dashed blue line), the two simulation sets are similar but are offset by a factor of two up until the halo mass scale (M 200 ∼ 10 12 M ) where AGN feedback causes an upturn at high masses in EAGLE. We note that if the comparison is instead performed as a function of halo mass (shown in Appendix C), the two mass loading factors agree almost perfectly between the two simulations over the common mass range between the simulations, which can be explained if the median stellar mass at a fixed halo mass is higher in FIRE than in EAGLE. Anglés-Alcázar et al. (2017) present a complementary measurement to Muratov et al. (2015) using Lagrangian particle tracking to measure particles ejected from the ISM (solid magenta line, taken from the fit presented in Davé et al. 2019), similar to our preferred methodology in this study. These measurements are presented as a cumulative integration over all Figure 14. Comparison of wind mass loading factors and outflow rates between EAGLE and other recent hydrodynamical simulations from the literature. Top-left: compares EAGLE (blue) and Illustris-TNG (red Nelson et al. 2019) at z = 2, showing the median mass loading factor for gas at r = 10 kpc, plotted as a function of stellar mass. Different line styles indicate different minimum radial velocity cuts, as labelled. For Illustris-TNG, we also show the mean mass loading factor applied at injection (dotted red line). Top-Right: mass loading for gas at different distances from the halo centre, as labelled. In this case, gas is selected with a radial velocity of v r > 50 km s −1 . For EAGLE, we only show measurements at 50 kpc for galaxies with M > 10 9 M , below which gas at 50 kpc is outside the halo virial radius (we instead show measurements for a shell at the virial radius with the blue dash-dotted line). Middle-left: compares EAGLE and the FIRE zoom-in simulations (magenta). Note that the FIRE simulations do not include AGN feedback. Dashed blue and magenta lines compare shell-based measurements of the mass loading at r = 0.25R vir , at redshift z = 0.25 (Muratov et al. 2015). For EAGLE, the dashed-dotted blue line shows the same but for a shell at the virial radius. For FIRE, individual galaxies are shown by the magenta points for shells at different radii, as labelled (Muratov et al. 2015(Muratov et al. , 2017. Solid lines show (tracking-based) mass loading factors for gas being ejected from the ISM, time-integrated over the entire history of each galaxy (Anglés-Alcázar et al. 2017). Middle-right: compares EAGLE, NIHAO (orange), and the simulations presented in Christensen et al. (2016) (green) at z = 0, showing (tracking-based) mass loading factors for gas being ejected from the ISM of galaxies (solid lines), plotted as a function of halo circular velocity at the virial radius. For EAGLE and NIHAO, we also show mass loading factors for gas being ejected through the halo virial radius (dashed lines). Note that only EAGLE includes AGN feedback. Bottom: compares EAGLE and the Horizon-AGN simulation (yellow; Beckmann et al. 2017), showing outflow rate as a function of stellar mass. Outflow rates are computed with a shell-based method at r = 0.2 R vir (bottom-left) and at r = 0.95 R vir (bottom-right). A colour version is available online. outflow and star formation events over the entire history of each galaxy shown. We perform an equivalent integration for our particle tracking-based outflow rates in EAGLE, presented as the solid blue line in the middle-left panel of Fig. 14. We note that both the galaxy-and halo-scale outflow selection criteria differ between the two studies, although they are both designed to in principle measure the same thing (the outflow rates of gas being ejected from the ISM/halo by feedback).
As with the comparison to TNG-50, larger differences become apparent when considering the change in the mass loading as a function of radius. Muratov et al. (2017) present measurements of the mass loading in FIRE at the virial radius (magenta crosses), which can be compared to measurements at a quarter of the virial radius (magenta triangles, or the dashed line) from Muratov et al. (2015). In most cases, the mass loading is smaller at larger radii in FIRE, whereas the opposite is true in EAGLE at low redshifts. As with the comparison to TNG-50, this implies that recycling of gas ejected from galaxies is likely much more efficient in FIRE than in EAGLE.
The middle-right panel of Fig. 14 presents a comparison with two additional sets of zoom-in simulations, including the simulations of Christensen et al. (2016, green line), and measurements of the NIHAO simulations presented by Tollet et al. (2019, orange lines), both as a function of stellar mass. Neither of these simulations include AGN feedback. Both these studies utilize particle tracking-based measurements of outflows, which we compare to our particle tracking measurements at z = 0. Tollet et al. (2019) find substantially higher mass loading factors for gas ejected from the ISM (solid orange line) than in EAGLE (solid blue line), with a very steep dependence on mass. They also find that less mass is (on average) ejected from the virial radius (dashed orange line), which is in strong disagreement with the z = 0 measurements from EAGLE (dashed blue line). Christensen et al. (2016) find somewhat lower mass loading factors for gas being ejected from the ISM, and with a slightly steeper mass dependence than in EAGLE. They find that a substantial fraction of this gas is then ejected from the virial radius, but do not present measurements of gas being ejected from haloes that was not previously in the ISM, making it unclear how their simulations compare in terms of outflows at the virial radius.
Finally, the bottom panels of Fig. 14 present a comparison with the Horizon-AGN simulation (Dubois et al. 2014), showing measurements presented in Beckmann et al. (2017). 11 Beckmann et al. (2017) measure outflow rates for two 2 kpc thick shells at 20 per cent (bottom-left) and 95 per cent of the halo virial radius, and we plot their measurements as a function of stellar mass (without any aperture correction). Reproducing these measurements in EAGLE, the comparison shows that outflow rates at a given stellar mass are (in most situations) significantly higher than in Horizon-AGN (for example, by about 0.5 dex at 0.2R vir at z = 0). We note that there are substantial differences between the low redshift galaxy stellar mass function in Horizon-AGN and EAGLE, with Horizon-AGN significantly overpredicting the stellar masses of low-mass galaxies, and EAGLE underpredicting the abundance of galaxies at the knee of the mass function Kaviraj et al. 2017). As such the comparison performed here will be comparing galaxies hosted by dark matter haloes of differing mass.
Also of interest is the comparison between outflow rates at 0.2 versus 0.95R vir . For massive galaxies, both simulations eject similar or greater amounts of gas from haloes than through the inner surface at 0.2R vir . For lower mass galaxies, EAGLE continues to eject similar or greater amounts through the outer surface, whereas Horizon-AGN ejects very little gas through the virial radius compared to the inner surface. This underlines the importance of considering outflowing flux as a function of scale out into the halo.
Taken at face value, the comparisons shown in Fig. 14 indicate that hydrodynamical simulations are seemingly able to reproduce observed stellar masses with different scenarios for gaseous outflows, similar to the situation seen for semi-analytic models in Fig. 13. That said, while we have emphasized the differences it is also important to emphasize that there is qualitative agreement between simulations, in the sense that all predict declining mass loading factors as a function of galaxy mass up to M 200 ∼ 10 12 M , and are in a similar level of qualitative agreement at higher masses if AGN feedback is included. We caution furthermore that some of the differences between the relations shown in this figure will arise from differences in the selection of outflowing particles (this only applies to the Lagrangian measurements), and so the discrepancies could be exaggerated in some cases. In addition, the level of agreement with the observed galaxy stellar mass function is unknown for zoom-in simulations (that must instead rely on comparison to the inferred median relationship between stellar mass and halo mass for central galaxies), and large differences in the stellar mass function could exist between some of the different simulations shown (this is definitely the case for Horizon-AGN).
As for the question of why galactic winds in the EAGLE simulations appear to entrain more circumgalactic gas at larger radii compared to other simulations (at least for those where such a comparison is currently possible), we speculate that is related to the high heating temperatures adopted in the EAGLE feedback model. In reality, energy from feedback is initially injected into a far smaller mass of material compared to the mass that is heated or kicked for the implementations of subgrid feedback models used in all cosmological simulations, such that gas around stars and black holes will (at least locally) achieve much larger velocities and temperatures. The choice made in EAGLE to heat relative few particles to a high temperature was motivated by this realization, and could plausibly lead to outflows escaping the ISM with higher specific energy than in other simulations, allowing the winds to have a greater impact on the ambient CGM. Explicit comparison of the energetics of outflows at different spatial scales between different simulations would show whether or not this is indeed the case.

Comparison to observations
As discussed in the introduction to this work, our priority in this study is to measure the flows of gas leaving galaxies and haloes using the information available from our simulations, independent of observational considerations. Our Lagrangian methodology for selecting outflowing gas does not naturally map on to the way outflowing gas is detected in observations, and in addition we do not explore any phase decomposition of outflowing gas. With these caveats in mind, it is none the less interesting to perform a rudimentary comparison between the outflow rates in EAGLE and a best guess for the outflow rates of real galaxies from observations. We choose to compare to down the barrel observations of seven local galaxies from Chisholm et al. (2017), which use the HST-COS spectrograph to detect multiple ultraviolet (UV) metal ions in absorption against the continuum of the associated galaxy, enabling (via photoionization modelling as a function of velocity) a robust determination of the ionization structure of the outflowing gas, in turn permitting a determination of the outflow rate. These observations are Figure 15. Comparison of mass loading factors between the EAGLE simulations and estimates from 'down-the-barrel' observations of local galaxies from Chisholm et al. (2017), plotted as a function of galaxy stellar mass. Mass loading factors from EAGLE are measured for gas being ejected from the interstellar medium over a redshift range 0 < z < 0.3. Solid (dashed) lines indicate the stellar mass range within which galaxies contain on average more (fewer) than 100 stellar particles. EAGLE is consistent with these observations, although the level of agreement might be fortuitous, given that we are not comparing like-for-like quantities. estimated to probe outflowing gas at small scales with respect to the galaxy; Chisholm et al. (2016) estimate that the detected outflowing gas is within 300 pc from the galaxy along the line of sight.
The EAGLE simulations do not include all of the relevant physics (for example, photoionization from local radiation sources) and do not reach the resolution required to robustly mimic such a selection of gas. Our measurements of how much gas is removed from the ISM by feedback are, however, still conceptually the same quantity as is being inferred from the observations, and so we do none the less present a comparison to the outflow rates of gas being ejected from the ISM in EAGLE, shown in Fig. 15. The observations of Chisholm et al. (2017) probe star-forming galaxies in the stellar mass range where stellar feedback is expected to dominate. They find evidence for an anticorrelation between mass loading factor and stellar mass, with a power-law slope of −1.6 when plotted as a function of circular velocity. Their relation is consistent with our measurements from the EAGLE simulations; we find a best-fitting slope of −1.5 as a function of halo circular velocity, V c ≡ √ GM 200 /R vir , for low-mass haloes (M 200 < 10 12 M ) at z = 0. This agreement is encouraging, and demonstrates that the outflow rates in EAGLE are not implausible given current constraints. At the same time, the level of the quantitative agreement is likely fortuitous to some extent, as we are not comparing like-for-like quantities.

S U M M A RY
We have presented measurements of outflow rates of gas from galaxies and from their associated dark matter haloes, taken from the reference EAGLE hydrodynamical simulation. We find that galactic winds are driven from the ISM in EAGLE with a mass loading factor (η ≡Ṁ out /Ṁ ) that scales approximately as η ∝ M −0.5 200 ∝ V −3/2 c for low-mass galaxies (M 200 < 10 12 M , Fig. 1). For reference, η ∝ M −1 200 would be required to explain the empirically inferred scaling of stellar mass with halo mass for M 200 < 10 12 M using galaxy-scale outflows alone (see the discussion in Section 3), implying that additional sources of mass scaling are required to explain the agreement between EAGLE and the observed galaxy stellar mass function. We do find a scaling close to η ∝ M −1 200 when measuring outflow rates at the virial radius, but a discussion of the complete picture is deferred to a future study where we will present measurements of gaseous inflow rates at different spatial scales. Parametric fits to the mass loading factor as a function of redshift and halo mass (as well as stellar mass and halo maximum circular velocity) are provided in Appendix B.
Similar to the result found in the recent analysis of the TNG-50 simulation of Nelson et al. (2019), we find that AGN feedback causes the scaling of the wind mass loading factor with mass to flatten and then increase for galaxies above a characteristic halo mass of 10 12 M (Fig. 11).
We find that the mass loading factor has a steeper dependence on halo mass when measured at the halo virial radius, and with a much clearer upturn due to AGN feedback at high masses (Fig. 1). We also find typically that significantly more baryons are ejected through the virial radius than out of the ISM, particularly at low and high halo masses, and at low redshift (Fig. 2). We explore a number of simplified explanations for this effect. We find that winds are overpressured relative to the ambient CGM in EAGLE, consistent with an energy-driven scenario in which outflows generate momentum by doing PdV work (Fig. 8). We also find that time delay effects play a role, as outflows at the halo virial radius will reflect energy injection rates in the past, which, for example, were higher for low-redshift galaxies (Fig. 9).
Aside from mass fluxes, we also consider a number of outflow properties in the simulations. We estimate that winds in EAGLE typically retain ≈20 per cent of the energy injected by feedback, modulated by trends with both halo mass and redshift (Fig. 5). Outflow velocities cover a wide range at a given halo mass/redshift, and increase positively with redshift and halo mass up to M 200 ∼ 10 12 M (Fig. 4). Below this mass, the median outflow velocity scales with mass similarly to the halo circular velocity. Outflows exhibit a clear bimodal flow pattern, with strong preferential alignment along the minor axes pointing orthogonal to galaxy discs.
We find that ISM-scale mass fluxes are not fully converged with numerical resolution for M 200 < 10 12 M , despite adjusting feedback parameters to recalibrate the simulation against the observed galaxy stellar mass function (Fig. 12). Convergence is better for outflows at the halo scale, and qualitatively our main conclusions hold at eight times higher mass resolution for both spatial scales. Energy fluxes are well converged at both scales in the recalibrated simulation, indicating that energy fluxes are a better indicator of the impact of feedback on galaxy assembly than ISM-scale mass fluxes (since the galaxy stellar mass functions of the two simulations agree despite discrepant ISM-scale mass fluxes).
Comparing to other cosmological hydrodynamical simulations (Fig. 14), we demonstrate that while substantial quantitative differences are found for gas being driven from the ISM (up to 0.5 dex), most simulations show qualitatively similar trends, although for M 200 > 10 12 M this is only the case if AGN feedback is included. The largest uncertainty in the current picture for outflows comes from the dichotomy between outflow rates measured at different spatial scales. For example, we show that the EAGLE and Illustris-TNG simulations present completely different scenarios for gas outflows at 50 kpc from galaxies versus outflows at 10 kpc. At z = 2 and at a galaxy stellar mass of 10 9 M , outflow rates are an order of magnitude higher at 10 kpc than at 50 kpc in Illustris-TNG, whereas there is little difference in flux between these spatial scales in the EAGLE simulation. At a high mass (M ∼ 10 11 M ), outflows in TNG stay approximately constant with radius from 10 to 50 kpc, whereas outflow rates increase with radius in EAGLE by nearly an order of magnitude over the same range. EAGLE therefore presents an ejective (but not ballistic) scenario for galactic winds driven by stellar feedback, where comparatively few baryons are removed from the ISM but are driven out to relatively large distances while sweeping up circumgalactic gas. Illustris-TNG instead presents a comparatively more fountain-like scenario, where more baryons are removed from the ISM by supernovae but are not driven as far, and so (presumably) can be recycled more efficiently. Davies et al. (2019) have shown that this is reflected in the total baryon content of lowmass (M 200 < 10 12 M ) haloes between the two simulations, with much higher baryon fractions for TNG than for EAGLE. The FIRE zoom-in simulations are similar to TNG in the sense that they report lower outflow rates at the virial radius than at one quarter of the virial radius (although FIRE agrees well with EAGLE at r = 0.25 R vir if the comparison is performed at a fixed halo mass).
The differences between simulations closely echo the picture encapsulated by semi-analytic galaxy formation models (Fig. 13), where acceptable matches to the observed galaxy luminosity function can be achieved using a very wide range in mass loading factor, with high outflow rates from the ISM (but not from the halo) being degenerate with high outflow rates through the halo virial radius (but not from the ISM), and both scenarios also being degenerate with the time-scale for ejected gas to return. Measurements of the distribution of metals both within and outside galaxies presumably represent a means to move beyond this impasse, as well as observational estimates of outflow rates that span a range of spatial scales. With some of the clearest differences between simulations seen for low-mass galaxies (M ∼ 10 9 M ), the regime where AGN feedback is not predicted to play an important role, observations that probe metals in the vicinity of dwarf galaxies may represent a particularly promising avenue to distinguish between ejective and fountain-like scenarios (e.g. Burchett et al. 2016;Johnson et al. 2017).

A1 Temporal convergence
An important caveat of Lagrangian flux measurements is that any mass element that crosses the chosen surface more than once (over the finite time interval adopted) will lead to an underestimate of the flux. This is particular pertinent for measuring fluxes at the interface to the ISM, where the time-scales for gas to be accreted Figure A2. The impact of changing the time-averaged radial velocity cut used to select outflowing galaxy-scale wind particles (top/middle), and of changing the ISM definition (bottom). Adjusting the radial velocity cut has negligible effects on the halo-scale outflows (middle panel). A factor of two change relative to our fiducial velocity cut of 0.25 V max changes the galaxyscale outflow rates by about 0.1 dex, although larger differences are seen in high-mass haloes. The lower panel shows the impact of changing our fiducial ISM definition to a selection of star-forming gas only. The main effect of including non-star-forming gas in our ISM definition is to enhance the outflow rates in low-mass galaxies (where metallicities are low and less gas can pass the metallicity-dependent star formation threshold). and then ejected by feedback can be short. The top panel of Fig. A1 shows the temporal convergence properties for measurements of gas ejected from the ISM, using a 25 3 Mpc 3 simulation with a higher frequency of simulation outputs (1000 in total, compared to 200 for our fiducial simulation). We find that the measurements start to be reasonably well converged once 200 simulation outputs are used (cyan lines), apart from at low redshift, where outflow rates are underestimated by ≈0.2 dex with respect to the measurements made using 1000 snapshots (solid yellow line). There are only 200 processed simulation outputs available for the larger 100 Mpc reference simulation, and so we use this set of outputs (and associated merger trees) for our fiducial analysis in this study. The temporal spacing of these 200 outputs is shown by the bottom panel of Fig. A1. Fig. A1 shows that temporal convergence issues tend to affect outflow rates with a fairly constant fractional offset as a function of halo mass at a given redshift. The main effect on our results is that the offset between our measurements of galaxy-and halo-scale outflow rates (which are much better converged due to the longer associated time-scales) will increase spuriously at low redshift, partly explaining the trends seen in Fig. 2.

A2 Lagrangian outflow rates
We present here a more detailed explanation of how we arrived at the selection criteria described in Section 2.6. These criteria were chosen to find a reasonable balance between completeness and purity, with temporal convergence (as described above) as another consideration. Given that the aim of this study is to measure the flux of gas being evacuated from the ISM or from the halo, we take genuine outflowing particles to be those that leave a given component, and then proceed to move a significant distance outwards in radius.
For a number of example galaxies (spanning a wide range in mass and redshift), we compute the maximum change of radius, r max , for particles that leave the ISM (or halo virial radius), up until the time that the particle ceases to be outflowing, or otherwise rejoins the ISM (or halo). The distribution of r max is typically characterized by a peak of particles that do not move significantly outwards, and then a long extended tail of particles that move over a wide range of radii. For at least a subset of the total mass and redshift range, the fraction of particles leaving the ISM (or halo) that do not move a significant distance is substantial. From detailed inspection of individual particle trajectories, these are often particles that have recently been accreted on to the ISM but are still in the process of settling into the disc, and so fluctuate across the ISM boundary a number of times.
We find that adopting cuts in instantaneous velocity or energy (as with the criteria adopted by, for example, Hopkins et al. 2012;Christensen et al. 2016) yields relatively poor completeness/purity in terms of the radial distance then travelled by outflowing particles. At the same time, computing r max from the full future radial trajectory of all particles from a large simulation would be prohibitively expensive. We compromise in this by computing the radial displacement reached by particles after one quarter of a halo dynamical time has passed since they left the ISM, which we find to be an excellent proxy for the r max computed using many simulation outputs. This thus motivates our choice of equation (5) in Section 2.6.

A3 Impact of radial velocity cuts and ISM definition
Our results are not highly sensitive to the choice of (time-averaged) radial velocity cut in equation (5). This is demonstrated in Fig. A2, which shows that galaxy-scale outflow rates change by small amounts when varying the cut (although there is a more significant impact Figure A3. A comparison of our fiducial Lagrangian measurement of wind mass loading factors with simple shell-based measurements, performed at z = 2 (left) and at z = 0.25 (right). Lagrangian outflow rates (red) are shown for gas leaving the ISM (solid), and the halo virial radius (dashed). Eulerian measurements (blue) are computed from the instantaneous radial momentum, summed over particles within shells of width 0.1R vir , including only particles outflowing faster than some minimum (instantaneous) radial velocity. Top-row: a comparison with shells at different radii, selecting all outflowing gas. Middle-row: a comparison with shells at one quarter of the halo virial radius, selecting gas with different radial velocity cuts. Bottom-row: a comparison with shells at the halo virial radius, selecting gas with different radial velocity cuts. A colour version is available online. on outflow rates in high-mass haloes). Note that removing the cut completely would have a much larger impact (gas leaving the ISM can often move inwards over quarter of a halo dynamical time).
The bottom panel of Fig. A2 demonstrates the impact of including/excluding non-star-forming gas from our ISM criteria. Including this material slightly enhances the galaxy-scale outflow rates at all mass scales, but the main effect is to substantially enhance the outflow rates at a low halo mass. Given that this is the regime where galaxies often have not formed a single star particle over the entire redshift interval (meaning results are likely not well converged at a low mass), the impact on our results is modest. Fig. A3 presents a comparison of our Lagrangian measurement of wind mass loading factors with a simple Eulerian measurement performed by summing the radial momentum of outflowing particles within spherical shells. For shells placed at the halo virial radius, it is evident that our Lagrangian criteria are equivalent to selecting all outflowing gas with v rad > 0 km s −1 , reflecting the looseness of our Lagrangian selection criteria for halo-scale outflows.

A4 Comparison of Lagrangian and Eulerian fluxes
Comparing our Lagrangian galaxy-scale outflows to shell-based measurements at one quarter of the halo virial radius, it is clear that the Lagrangian measurements are always lower. At z = 2 (top-left), there is some evidence for entrainment seen in the shell measurements at different radii, with outflow rate increasing with radius by about 0.3 dex in both low-and high-mass haloes.
At z = 0, there is an increased entrainment effect seen in the shell measurements, with outflow rates higher by about 0.5 dex at the virial radius when compared to one quarter of the virial radius. Our Lagrangian measurements for gas leaving the ISM are again lower. A lack of temporal convergence at z = 0 has a small systematic contribution to this effect (see Fig. A1).

A P P E N D I X B : PA R A M E T R I C F I T S TO G A L A X Y-S C A L E A N D H A L O -S C A L E M A S S L OA D I N G FAC TO R S
To facilitate comparisons with other studies, we provide parametric fits to the mass loading factors for both galaxy-and halo-scale outflows (as shown in Fig. 1). We only fit to data from bins where more than 80 per cent of the galaxies have formed at least one star particle, integrated over the redshift bins indicated in Fig. 1. We find that a reasonable fit to the mass loading factors as a function of halo mass is given by the five-parameter function, where two of the parameters are fitted as a linear function of redshift, z. For halo-scale outflows, we find that a reasonable fit is given by where in this case log 10 (N) is fitted as a function of expansion factor, a.
The parameter α sets the low-mass power-law slope of the mass loading factor as a function of halo mass (primarily related to stellar feedback), and β sets the power-law slope of the upturn at higher masses (primarily related to AGN feedback). M 1 sets the transition halo mass scale between these regimes, 12 and N sets the overall normalization. δ and M cut are responsible for the third (flatter or negative) power-law slope that becomes evident in group/cluster mass haloes. Both of these fits (galaxy and halo scale) provide a reasonable description of the data shown in Fig. 1 (within at least ≈0.1 dex) up until z = 4.
We also provide parametric fits to the galaxy-scale mass loading factor as a function of galaxy stellar mass and halo maximum circular velocity. We again adopt the five parameter functional form given by equation (B1), switching the dependent variable from halo mass to maximum circular velocity or stellar mass.
As a function of halo maximum circular velocity, V max , we find log 10 (N ) = −0.31 + 0.14 z In both cases, the fits provide a reasonable description of the data up to z = 3.

A P P E N D I X C : C O M PA R I S O N TO T H E F I R E S I M U L AT I O N S A S A F U N C T I O N O F H A L O M A S S
Further to the comparison between simulations shown in Fig. 14 of Section 5.2, Fig. C1 shows a comparison between the mass loading factors in FIRE and EAGLE at low redshifts, measured at r = 0.25 R vir , and plotted in this case as a function of halo mass, rather than as a function of stellar mass. Compared at a given halo mass, EAGLE and the best-fitting relation from FIRE are in remarkably good agreement over the common mass range. The level of agreement is significantly better than when the simulations are compared as a function of galaxy stellar mass (as shown in Fig. 14), where the Figure C1. A comparison of mass loading factors between EAGLE and the FIRE simulations, in this case plotted as a function of halo mass. Mass loading factors are measured in shells with 0.2 < r/R vir < 0.3, selecting all outflowing gas. The magenta line shows the best-fitting relation from Muratov et al. (2015) plotted at z = 0.25. The blue line shows the average mass loading factor from EAGLE, also plotted at z = 0.25. Plotted as a function of halo mass, EAGLE and FIRE are in excellent agreement over the common mass range at r = 0.25R vir . A colour version is available online.
best-fitting mass loading factor relation in FIRE is about 0.3 dex higher than the average from EAGLE at M ∼ 10 9 M . This implies that there is a systematic difference in the median stellar mass versus halo mass relation between the two sets of simulations. Individual galaxies in FIRE at M H ∼ 10 12 M fall below the plotted best-fitting relation at low redshifts, and are observed to be relatively quiescent in terms of outflow activity, with residual outflowing flux attributed to non-feedback sources (Muratov et al. 2015). This is the mass scale where AGN feedback (which is not implemented in the FIRE simulations) starts to play a significant role in EAGLE, causing the upturn of the mass loading factor at higher halo masses. This paper has been typeset from a T E X/L A T E X file prepared by the author.