We introduce the SIMBA simulations the next generation of the MUFASA cosmological galaxy formation simulations run with GIZMOs meshless finite mass hydrodynamics. SIMBA includes updates to MUFASAs subresolution star formation and feedback prescriptions and introduces black hole growth via the torquelimited accretion model of AnglsAlczar et al. from cold gas and Bondi accretion from hot gas along with black hole feedback via kinetic bipolar outflows and Xray energy. Ejection velocities are taken to be   103 km s1 at high Eddington ratios increasing to   8000 km s1 at Eddington ratios below 2 per cent with a constant momentum input of 20Lc. SIMBA further includes an onthefly dust production growth and destruction model. Our SIMBA run with 100h1 Mpc3 and 1024SUP3SUP gas elements reproduces numerous observables including galaxy stellar mass functions at z  06 the stellar massstar formation rate main sequence H I and HSUB2SUB fractions the massmetallicity relation at z  0 2 starforming galaxy sizes hot gas fractions in massive haloes and z  0 galaxy dust properties. However SIMBA also yields an insufficiently sharp truncation of the z  0 mass function and toolarge sizes for lowmass quenched galaxies. We show that SIMBAs jet feedback is primarily responsible for quenching massive galaxies.
['galaxies: evolution', 'galaxies: formation', 'Astrophysics - Astrophysics of Galaxies', 'Astrophysics - Cosmology and Nongalactic Astrophysics']
SIMBA Cosmological simulations with black hole growth and feedback
{'Simba : Cosmological Simulations with Black Hole Growth and Feedback': "Romeel Dav'e 1 ; 2 ; 3 , Daniel Angl'es-Alc'azar 4 , Desika Narayanan 5 ; 6 ; 7 , Qi Li 5 , Mika H. Rafieferantsoa 2 ; 3 , & Sarah Appleby 1  \n- 1 Institute for Astronomy, Royal Observatory, Univ. of Edinburgh, Edinburgh EH9 3HJ, UK\n- 2 University of the Western Cape, Bellville, Cape Town 7535, South Africa\n- 3 South African Astronomical Observatories, Observatory, Cape Town 7925, South Africa\n- 4 Center for Computational Astrophysics, Flatiron Institute, 162 Fifth Avenue, New York, NY 10010, USA\n- 5 Department of Astronomy, University of Florida, 211 Bryant Space Sciences Center, Gainesville, FL, USA\n- 6 University of Florida Informatics Institute, 432 Newell Drive, CISE Bldg E251, Gainesville, FL, USA\n- 7 Cosmic Dawn Center at the Niels Bohr Institute, University of Copenhagen and DTU-Space, Technical University of Denmark  \nAccepted XXX. Received YYY; in original form ZZZ", 'ABSTRACT': "We introduce the Simba simulations, the next generation of the Mufasa cosmological galaxy formation simulations run with Gizmo 's meshless finite mass hydrodynamics. Simba includes updates to Mufasa 's sub-resolution star formation and feedback prescriptions, and introduces black hole growth via the torque-limited accretion model of Angl'es-Alc'azar et al. (2017a) from cold gas and Bondi accretion from hot gas, along with black hole feedback via kinetic bipolar outflows and X-ray energy. Ejection velocities are taken to be GLYPH<24> 10 3 kms GLYPH<0> 1 at high Eddington ratios, increasing to GLYPH<24> 8000 kms GLYPH<0> 1 at Eddington ratios below 2%, with a constant momentum input of 20 L GLYPH<157> c . Simba further includes an on-the-fly dust production, growth, and destruction model. Our Simba run with ' 100 h GLYPH<0> 1 Mpc ' 3 and 1024 3 gas elements reproduces numerous observables, including galaxy stellar mass functions at z = 0 GLYPH<0> 6 , the stellar mass-star formation rate main sequence, H i and H 2 fractions, the mass-metallicity relation at z GLYPH<25> 0 ; 2 , star-forming galaxy sizes, hot gas fractions in massive halos, and z = 0 galaxy dust properties. However, Simba also yields an insufficiently sharp truncation of the z = 0 mass function, and too-large sizes for low-mass quenched galaxies. We show that Simba 's jet feedback is primarily responsible for quenching massive galaxies.  \nKey words: galaxies: formation, galaxies: evolution, methods: N-body simulations, galaxies: mass function", '1 INTRODUCTION': "The formation and evolution of galaxies is governed by a wide range of physical processes spanning from sub-parsec to giga-parsec scales. These include the growth of largescale structure, cooling and heating of astrophysical plasmas, and the formation of stars and central black holes along with associated energy return processes collectively known as feedback. As observations of galaxies improve at an everadvancing pace, a key goal of modern astrophysics is to understand how such observations can be used to constrain the balance of underlying physical processes driving galaxy evolution, particularly in regards to feedback processes that are currently among its least well-understood aspects.  \nModern galaxy formation models are generally built on  \nthe premise that feedback at the very lowest masses is dominated by photoionisation from the metagalactic ultraviolet background, feedback at scales above this but below L ? is driven primarily by energy and momentum from young stars and supernovae, and feedback in massive galaxies predominantly owes to energetic release from accretion disks around supermassive central black holes (see Somerville & Dav'e 2015, for a review). While this framework has been broadly successful at reproducing many key observed characteristics of the galaxy population at a range of cosmic epochs, the physical understanding of most of these small-scale feedback processes remains coarse and heuristic (see Naab & Ostriker 2017, for a review). Galaxy formation models have now progressed to a point where numerous models can reproduce similar core sets of observations, but they often do so  \nwith significantly different underlying physical models and assumptions. To discriminate among these physical drivers, it becomes important to develop modeling methodologies that are as well-motivated and realistic as possible, and to test such models against the widest possible suite of observations quantifying the stellar, gas, metal, and black hole properties of galaxies along with their surrounding gas.  \nCosmological-scale simulations that model galaxy growth and feedback dynamically within evolving largescale structure are an increasingly valuable tool for testing and constraining galaxy formation physics, owing both to rapidly advancing computational power and the commensurate ability to concurrently model a large range of scales and physical processes. State of the art models now simultaneously predict the co-evolution of stars, interstellar media, black holes, and circum-galactic gas, enabling a holistic approach towards testing the input physics against observations across a wide range of scales, environments, and wavelengths. Modern cosmological-scale simulations such as Illustris (Vogelsberger et al. 2014; Genel et al. 2014), Magneticum (Hirschmann et al. 2014), Horizon-AGN (Dubois et al. 2014; Volonteri et al. 2016; Kaviraj et al. 2017), Eagle (Schaye et al. 2015), MassiveBlack (Khandai et al. 2015), Blue Tides (Feng et al. 2016), Romulus (Tremmel et al. 2017), and Illustris-TNG (Springel et al. 2018) have implemented ever-improving sub-grid models aimed at more successfully reproducing the stellar, gaseous, and black hole contents of galaxies in bulk, while numerous associated hierarchically-situated zoom simulations using more detailed input physics can examine the internal structural and dynamical properties of galaxies with increasing fidelity.  \nThe Mufasa simulation project (Dav'e et al. 2016) has added to the pantheon of such simulations, employing several novel approaches that distinguish it from others. First, it utilises meshless finite mass (MFM) hydrodynamics as implemented in the Gizmo code (Hopkins 2015, 2017), which offers important accuracy advantages over Smoothed Particle Hydrodynamics (SPH), and owing to the mass conserving nature of its gas elements greater ease of analysis as compared to adaptive mesh refinement ( Ramses ) and moving mesh ( Arepo ) codes. Second, rather than employing simple parameterisations or a cooling shutoff to drive galactic outflows, the kinetic mass outflow rate is taken directly from very high-resolution simulations from the Feedback in Realistic Environments (FIRE; Hopkins et al. 2014, 2018; Muratov et al. 2015) simulations, providing a synergy between ISM-resolving simulations of individual galaxies and cosmological-scale simulations of galaxy populations.  \nAn aspect where Mufasa was physically less well motivated than other state of the art simulations was in its depiction of black hole growth and AGN feedback. Mufasa did not directly grow black holes and utilise the accretion energy for feedback to quench galaxies. Instead, following Gabor & Dav'e (2015), Mufasa implemented a heuristic 'quenching feedback' in which diffuse gas within massive halos was prevented from cooling. This energy was envisioned to be putatively supplied by an AGN, but Mufasa did not explicitly model black hole accretion and the interaction of its energy release with surrounding gas. The halo mass scale above which quenching feedback was applied was taken from the best-fit analytic equilibrium model of Mitra et al. (2015, 2017), and evolved slowly upwards with redshift. Despite its  \nsimplicity, this prescription displayed impressive successes at reproducing a red sequence and massive central galaxy properties in excellent accord with observations (Dav'e et al. 2017b), albeit with some non-trivial discrepancies such as over-quenching satellites particularly in the outskirts of large halos (Rafieferantsoa & Dav'e 2018). This demonstrated that a model based primarily on starvation of the central galaxy via 'radio mode' feedback (Croton et al. 2006; Bower et al. 2006) is able to quench the galaxy population in a hydrodynamic simulation in broad agreement with observations. While this model represented an interesting numerical experiment, it would clearly be valuable to implement a more physically-motivated black hole growth and feedback model that retains and perhaps even extends the successes of Mufasa 's more heuristic model.  \nThis is the primary goal of the Simba project. As the descendent of Mufasa , Simba marries two lines of investigation to achieve this. First, it builds on the successful Mufasa model, including its representation of star formationdriven feedback and other modern features. To this, Simba adds a novel and promising model for black hole growth: Torque-limited accretion (Hopkins & Quataert 2011; Angl'esAlc'azar et al. 2013, 2015, 2017a). In this model, black hole growth is regulated by the ability for gas in the inner disk to lose angular momentum via disk instabilities. Hopkins & Quataert (2011) developed an analytic formalism for this, tested and calibrated using sub-pc scale numerical simulations, which yielded a formula that connects the infall rate of material onto the black hole accretion disk with properties of the inner galactic disk. They showed that even at GLYPH<24> 1 kpc resolution typical of cosmological simulations, their gravitational torque accretion formula provides a significantly better match to the measured accretion rate in their high-resolution simulations than employing the canonical Bondi accretion formula used in all other current cosmological black hole growth simulations.  \nAngl'es-Alc'azar et al. (2013) explored the torque-limited accretion model via the post-processing of zoom simulations, and Angl'es-Alc'azar et al. (2015) extended this approach to cosmological simulations. Their most significant result was that, unlike Bondi accretion, torque-limited accretion does not require the black hole to self-regulate its own growth. In particular, torque-limited accretion naturally results in black holes growing along the observed galaxy-black hole scaling relations, even without any black hole feedback. There is one free parameter in the model which represents the fraction of material entering the accretion disk that accretes onto the black hole; a plausible choice of GLYPH<24> 10 % provided a good match to data, insensitive to the choice of the (uncertain) black hole seed mass. Angl'es-Alc'azar et al. (2017a) extended these previous works to self-consistently incorporate torquelimited accretion into Gizmo , along with bipolar black hole winds, and demonstrated that the results obtained without feedback were reproduced in this case - in particular, the inclusion of feedback self-consistently confirmed the primary result obtained in the post-processed case that black hole-galaxy scaling relations arise naturally without the black hole self-regulating its own growth. Simba builds on this work to employ the torque-limited black hole accretion model of Angl'es-Alc'azar et al. (2017a) when accreting from cold or star-forming gas, in order to self-consistently grow black holes within galaxies during the simulation run. The  \nuse of torque-limited black hole growth is unique among current cosmological simulations. Simba also includes Bondi accretion, but only from hot gas when present close to the black hole since it is the physically appropriate model in that case.  \nThe second part of Simba 's new black hole model involves a novel sub-grid prescription for active galactic nuclei (AGN) feedback. AGN feedback connects flows coming off the black hole accretion disk to energy release on scales of tens or hundreds of kpc. To model this transfer of energy from small to large scales, Simba utilises kinetic outflows with outflow parameters based on observed AGN feedback. While there is still no well-defined theoretical consensus on the generation of black hole outflows and jets, recent observational progress has been rapid, showing that AGN can drive molecular and ionised outflows with velocities of GLYPH<24> 1000 kms GLYPH<0> 1 or more (Sturm et al. 2011; Greene et al. 2012; Maiolino et al. 2012; Liu et al. 2013; Perna et al. 2017a), and jets at velocities up to GLYPH<24> 10 4 kms GLYPH<0> 1 and more (Fabian 2012). Generally, high-velocity jets are observed to arise from early-type galaxies hosting massive black holes with low accretion rates relative to its Eddington rate ( f Edd < GLYPH<24> few percent), while lower-velocity outflows typically arise in systems with higher f Edd (Best & Heckman 2012; Heckman & Best 2014). Extreme systems such as bright quasars often show both types of outflows. Simba 's black hole outflows are parameterised to broadly follow such observed trends. Simba employs kinetic feedback (i.e. gas element kicks) for both feedback modes, with the kick velocity ranging from many hundreds of km/s in low-mass, fast-accreting black holes, up to many thousands of kms GLYPH<0> 1 for slower-accreting black holes.  \nA key unique feature is that Simba 's kinetic feedback is purely bipolar, with the ejection vector given by the angular momentum of the inner disk. This direction is relatively stable over galaxy dynamical timescales. To be conservative in minimising black hole feedback impact on the galaxy interstellar medium, we employ an opening angle of zero. This is in contrast to other simulations that successfully reproduce massive galaxy properties using Bondi accretion, which employ either spherical thermal input (e.g. EAGLE; Schaye et al. 2015) or randomise the kinetic feedback's direction on short timescales (e.g. Illustris-TNG; Weinberger et al. 2017). Horizon-AGN employed Bondi accretion with a two-mode feedback scheme (Dubois et al. 2012), but still used spherical thermal feedback during the high-Eddington growth phase, as later also done in Illustris-TNG. More detailed isolated elliptical simulations have also highlighted the importance of radiative mechanical feedback (Gan et al. 2014) that can reproduce observations of AGNs in ellipticals such as their duty cycle (Gan et al. 2019).  \nThe reason Simba is able to be successful with a genuinely bipolar model likely traces back to its accretion model: Torque-limited accretion does not require selfregulation of black hole growth, whereas Bondi accretion requires quasi-spherical energy distribution close to the black hole in order to self-regulate its own growth. In Simba 's case, the energy input sphericalises at large distances, sufficient in massive halos to quench inflow into the galaxy by keeping the halo gas hot. In this way, Simba 's accretion and feedback models work together to build a sub-grid description of AGN  \nfeedback that is more closely connected to observations of AGN winds and jets.  \nIn addition to kinetic AGN feedback, Simba also includes X-ray feedback input into surrounding gas. The importance of this feedback channel has been emphasised in zoom simulations by Choi et al. (2012), showing that it can potentially drive the quenching of massive galaxies. We adapt this model to operate under the lower-resolution conditions present in Simba 's cosmological-scales runs, and show that it plays a minor but non-trivial role in quenching the most massive galaxies. Simba is the first cosmologicalvolume simulation to include such X-ray feedback.  \nAnother novel aspect of Simba is that it includes a model for on-the-fly dust production and destruction, broadly following McKinnon et al. (2017)'s implementation into Arepo where the dust is passively advected with the gas. We include dust production from Type II supernovae (SNe) and Asymptotic Giant Branch (AGB) stars, and further growth via condensation from metals, while destruction can occur from sputtering, consumption by star formation, or SNe shocks. The fraction of metals locked into dust can be substantial, leading to significant changes in the predicted mass-metallicity relations. In Mufasa , we found it necessary to reduce the SN yields arbitrarily by a factor of two in order to match the observed gas-phase mass-metallicity relation, but in Simba we can reproduce this as well or better without such arbitrary factors, since a substantial fraction of the metals ends up locked in dust.  \nIn this paper, we describe the simulation methodology in Simba which, besides the new black hole model, also makes various other minor improvements to Mufasa ( § 2). We then present a range of observational comparisons to predicted stellar, gas, and metal properties, analogous to a sampling of results presented for Mufasa in a series of recent papers (Dav'e et al. 2016, 2017a,b; Rafieferantsoa & Dav'e 2018), paying particular attention to black hole and massive galaxy properties which represent the most direct test of Simba 's new AGN feedback model ( § 3). We show that Simba reproduces many observations comparably well or better than Mufasa , but now with a more realistic and self-consistent description of black hole growth and feedback. We then examine variants of Simba 's AGN feedback model in order to isolate the impact of its various components ( § 4), showing that the high-velocity jet feedback is crucial for producing a quenched massive galaxy population. Finally, we summarize our results in § 5.", '2.1 Code and input physics': "The Simba simulations utilise much of the framework of the Mufasa simulations as described in Dav'e et al. (2016), but there are a number of updates and additions based on recent theoretical and observations results, in addition to the major change of modeling black hole growth and feedback as well as dust. Here we recap the main features of the model, and then describe in more detail the new aspects of Simba .  \nSimba utilises a forked version of the Gizmo cosmological gravity plus hydrodynamics solver (Hopkins 2015, 2017), in its Meshless Finite Mass (MFM) version. This code, based  \non Gadget-3 (Springel 2005), evolves dark matter and gas elements together including gravity and pressure forces, handling shocks via a Riemann solver with no artificial viscosity. It performs very well in numerous standard hydrodynamics tests including strong shocks, rotating disks, and cold blobs moving through hot media (Hopkins 2015). It also preserves the mass within each fluid element during the evolution, thereby enabling detailed tracking of gas flows. It thus marries the advantages of a particle-based code such as adaptivity in space and time, with the hydrodynamics accuracy of a Riemann solved-based mesh code, without the imposition of a Cartesian mesh that can cause difficulties with Galilean invariance and rotating shear flows.  \nRadiative cooling and photoionisation heating are modeled using the Grackle-3.1 library (Smith et al. 2017), including metal cooling and non-equilibrium evolution of primordial elements. This is an updated version of Grackle2.1 used in Mufasa , offering two advantages: First, the adiabatic and radiative terms are evolved together during the cooling sub-timestep, unlike the previous operator-split approach where first the system was evolved adiabatically over the full timestep and then cooling was applied; this results in more accurate and stable thermal evolution particularly in the stiff regions of the cooling curve. Second, it includes selfshielding self-consistently during the simulation run, based on the Rahmati et al. (2013) prescription in which the ionising background strength is attenuated depending (primarily) on gas density. A spatially-uniform ionising background is assumed as specified by Haardt & Madau (2012), modified to account for self-shielding (A. Emerick, priv. comm.). These changes do not make a noticeable difference to the resulting galaxy population, but they do improve the accuracy of the baryonic thermal evolution which may be particularly important within circum-galactic gas. Furthermore, computing the neutral hydrogen content of gas particles is now done self-consistently within the code rather than via a post-processed application of self-shielding (Dav'e et al. 2017a).  \nAs in Mufasa , we use an H 2 -based star formation rate, where the H 2 fraction is computed based on the subgrid model of Krumholz & Gnedin (2011) based on the metallicity and local column density, with minor modifications as described in Dav'e et al. (2016) to account for variations in numerical resolution. The star formation rate is given by the H 2 density divided by the dynamical time: SFR = GLYPH<15> GLYPH<3> GLYPH<26> H 2 GLYPH<157> t dyn , where we use GLYPH<15> GLYPH<3> = 0 : 02 (Kennicutt 1998). The chemical enrichment model tracks eleven elements (H,He,C,N,O,Ne,Mg,Si,S,Ca,Fe) during the simulation, with enrichment tracked from Type II supernovae (SNe), Type Ia SNe, and Asymptotic Giant Branch (AGB) stars. The yield tables employed are the same as in Mufasa , namely Nomoto et al. (2006) for SNII yields, Iwamoto et al. (1999) for SNIa yields, and AGB star enrichment following Oppenheimer & Dav'e (2006). However, we no longer apply an arbitrary reduction of yields by a factor of 2 that was previously needed to match the mass-metallicity relation, and instead lock individual metals into dust; we detail the dust model implementation in § 2.5. Type Ia SNe and AGB wind heating are also included as in Mufasa , along with ISM pressurisation at a minimum level as required to resolve the Jeans mass in star-forming gas as described in Dav'e et al. (2016).  \nThe model for star formation-driven galactic winds closely follows that in Mufasa ; we continue to use decoupled two-phase winds, with 30% of wind particles ejected 'hot' i.e. with a temperature set by the supernova energy minus the wind kinetic energy. However, we make a significant update to the mass loading factor scaling with stellar mass. Mufasa used the scalings taken from Muratov et al. (2015), who computed the outflow rates based on mass advection across a boundary at one-quarter of the virial radius in the FIRE zoom simulations. Angl'es-Alc'azar et al. (2017b) used similar FIRE simulations, but instead tracked individual particles in order to quantify the mass outflow rates out of the star-forming region, thus providing a more direct measurement of the amount of material leaving the ISM. This yields a cumulative mass loading factor versus stellar mass as shown in Figure 1, which is well fit by a broken power law at M 0 = 5 : 2 GLYPH<2> 10 9 M GLYPH<12> :  \nGLYPH<17> ' M GLYPH<3> ' / 8 > > > < > > > 9 GLYPH<16> M GLYPH<3> M 0 GLYPH<17> GLYPH<0> 0 : 317 ; if M GLYPH<3> < M 0 9 GLYPH<16> M GLYPH<3> M 0 GLYPH<17> GLYPH<0> 0 : 761 ; if M GLYPH<3> > M 0 (1)  \n:  \nand is independent of redshift. This has a similar slope to Muratov et al. (2015) below M 0 but with roughly double the amplitude, and is much steeper above M 0 . It is also similar to the assumed mass loading factor in IllustrisTNG(Pillepich et al. 2018). Similar to TNG, Simba employs a flat GLYPH<17> ' M GLYPH<3> ' below an M GLYPH<3> corresponding to 16 star particles ( M GLYPH<3> GLYPH<20> 2 : 9 GLYPH<2> 10 8 M GLYPH<12> for the 100 h GLYPH<0> 1 Mpc run), otherwise poorly-resolved forming galaxies are unable to grow owing to excessive feedback. As in Mufasa , we further apply a reduction in GLYPH<17> at high redshifts, in order to allow for early galaxy growth in poorly resolved situations. Mufasa was found to underproduce z > 6 galaxies, and hence we strengthen the suppression factor at z > 3 from ' a GLYPH<157> 0 : 25 ' to ' a GLYPH<157> 0 : 25 ' f a , where a is the expansion factor. We tune the value of f a based on the resolution of the simulation, since the origin of the lack of early galaxy formation owes to poor resolution. Testing has shown that we obtain converged results that match z > GLYPH<24> 6 observations (shown later) if we use f a = 2 at our largest ( 100 h GLYPH<0> 1 Mpc ) volume's resolution, f a = 1 : 5 at 8 GLYPH<2> higher mass resolution, and so on. Fortunately, because galaxy growth is very rapid at high redshifts, this choice makes little difference to galaxy predictions at z > GLYPH<24> 3 over most of cosmic time. Note that the FIRE simulations do not make strong predictions for GLYPH<17> ' M GLYPH<3> ' at z > GLYPH<24> 3 owing to the limited dynamic range covered by their zooms at early epochs, so this choice is not in obvious conflict with using FIRE scalings at lower redshifts. However, the ad hoc nature of this correction means that results for galaxy stellar growth at high redshifts from Simba should be considered as tuned rather than predictive.  \nA new feature in Simba is metal-loaded winds. When a wind particle is launched, it extracts some metals from nearby particles to represent the local enrichment by the supernovae driving the wind. The metallicity added to the wind particle is given by  \ndZ = f SNII y SNII ' Z 'GLYPH<157> MAX ' GLYPH<17>; 1 ' (2)  \nwhere f SNII = 0 : 18 is the stellar mass fraction lost to supernova (assumed to be instantaneous in Simba ), y SNII ' Z ' is the metal-dependent Type II SN yield for each species, and GLYPH<17> is the mass loading factor. This amount is subtracted from nearby gas in a kernel-weighted manner. If there is not  \nFigure 1. Mass loading factor GLYPH<17> versus stellar mass M GLYPH<3> from a suite of FIRE simulations analysed via particle tracking in Angl'es-Alc'azar et al. (2017b). The points show values measured at various redshifts, while the orange line is the best-fit relation. The gray solid line shows the Muratov et al. (2015) scaling and the dashed and dotted gray lines show the mass loading factors used in the Illustris and Illustris-TNG simulations ( GLYPH<17> ' M halo ' fitting functions from Pillepich et al. 2018, converted to GLYPH<17> ' M GLYPH<3> ' using the M GLYPH<3> -M halo relation of Moster et al. 2013).  \n<!-- image -->  \n∗  \n/circledot  \nenough metals nearby (as can happen early on), then dZ is reduced appropriately. In all circumstances, the total metal mass is conserved. The metal loading factor (i.e. the ejected metallicity relative to surrounding gas) can be a factor of two or larger when the ISM has a metallicity GLYPH<28> y SNII , but is more typically around 10-20%, in broad agreement with metal loading estimates from zoom simulations (Muratov et al. 2017; Christensen et al. 2018).  \nSimba continues to use the wind velocity scalings from Muratov et al. (2015) as in Mufasa , since the scaling follows the expected quasi-linear scaling of wind speed with escape velocity as observed (e.g. Martin 2005), and also because a full analysis of the velocity distribution of outflowing gas is not available from FIRE. In Mufasa , the amplitude was taken as a tunable parameter, and set to 2 v circ at v circ = 200 kms GLYPH<0> 1 . Owing to the increase in the mass loading factor from low-mass galaxies, we find that a somewhat lower value of the wind velocity is required to compensate for this, so in Simba we reduce the normalisation to 1.6:  \nv w = 1 : 6 GLYPH<16> v circ 200 kms GLYPH<0> 1 GLYPH<17> 0 : 12 v circ + GLYPH<1> v ' 0 : 25 R vir ' (3)  \nwhere GLYPH<1> v ' 0 : 25 R vir ' is the velocity corresponding to the potential difference between the launch point and one-quarter of the virial radius (see Dav'e et al. 2016). A related new aspect in Simba is that we limit the wind kinetic energy to the available supernova energy by attenuating the wind speed when needed, although this only has noticeable effect in small galaxies at very early epochs.  \nSimba uses an on-the-fly approximate friends-of-friends (FOF) finder applied to stars and dense gas as described in Dav'e et al. (2016) in order to compute the galaxy properties such as M GLYPH<3> , and as in Mufasa obtains v circ using a scaling based on the baryonic Tully-Fisher relation. Besides algorithmic improvements to improve parallel performance, the  \nonly change to this is that the FOF finder now also groups black holes into galaxies.", '2.2 Black hole growth': "The most significant change in Simba relative to Mufasa is that black holes are seeded and grown during the simulation, and the accretion energy is used to drive feedback that serves to quench galaxies. In this section we describe Simba 's twomode accretion model. The first mode closely follows the torque-limited accretion model presented in Angl'es-Alc'azar et al. (2017a), and we refer the reader there for full details. The second mode uses Bondi accretion, but solely from the hot gas component. In future work (Angl'es-Alc'azar et al., in prep.) we will show that this contribution from Bondi accretion is sub-dominant for all but the highest mass black holes.", '2.2.1 Torque-limited accretion from cold gas': "We model the gas inflow rate GLYPH<219> M Torque driven by disk gravitational instabilities from galactic scales down to the accretion disk surrounding the central black hole following Hopkins & Quataert (2011):  \nGLYPH<219> M Torque GLYPH<25> GLYPH<15> T f 5 GLYPH<157> 2 d GLYPH<2> GLYPH<18> M BH 10 8 M GLYPH<12> GLYPH<19> 1 GLYPH<157> 6 GLYPH<18> M enc ' R 0 ' 10 9 M GLYPH<12> GLYPH<19> GLYPH<2> GLYPH<18> R 0 100 pc GLYPH<19> GLYPH<0> 3 GLYPH<157> 2 GLYPH<18> 1 + f 0 f gas GLYPH<19> GLYPH<0> 1 M GLYPH<12> yr GLYPH<0> 1 ; (4)  \nwhere f d is the disk mass fraction (including both stars and gas), M enc ' R 0 ' is the total gas+stellar mass, f gas is the gas mass fraction in the disk component, f 0 GLYPH<25> 0 : 31 f 2 d ' M d ' R 0 'GLYPH<157> 10 9 M GLYPH<12> ' GLYPH<0> 1 GLYPH<157> 3 , and all quantities are evaluated within a distance R 0 of each black hole enclosing the nearest 256 gas elements, with an upper limit R 0 GLYPH<20> 2 h GLYPH<0> 1 kpc (comoving) imposed throughout the simulation. Evaluating this equation for GLYPH<219> M Torque requires separating the spheroidal and disk components within R 0 , which we do by means of the kinematic decomposition implemented in Angl'esAlc'azar et al. (2017a); see the Appendix of Angl'es-Alc'azar et al. (2015) for further tests of this. Unlike previous work, we evaluate torque-limited accretion only for the cold ( T < 10 5 K) gas within the black hole kernel, since it relies on instabilities in a cold gaseous disk to drive mass inflow. We consider all ISM gas to be in the cold phase, where ISM gas is taken to be gas that has been artificially pressurised in order to resolve the Jeans mass as described in Dav'e et al. (2016). In our 100 h GLYPH<0> 1 Mpc run this corresponds to gas above a hydrogen number density of n H > 0 : 13 cm GLYPH<0> 3 and within 0.5 dex of the pressurized ISM temperature floor.  \nThe normalization factor GLYPH<15> T GLYPH<17> GLYPH<15> m GLYPH<2> GLYPH<11> T encapsulates processes that affect the radial transport of gas on unresolved scales (e.g. nuclear star formation and stellar feedback and mass loss in winds from the accretion disk), where GLYPH<11> T = 5 is the original normalization of GLYPH<219> M Torque in Hopkins & Quataert (2011) and we set GLYPH<15> m = 0 : 1 to match the normalization of the local M BH -M ? relation as in Angl'es-Alc'azar et al. (2017a). One can view GLYPH<11> T as corresponding to an efficiency of transport of material from the inner galactic disk onto the black hole accretion disk, and GLYPH<15> m as the efficiency of transport  \nfrom the accretion disk onto the black hole, for which 10% is a canonical value. However, GLYPH<11> T is itself fairly uncertain, and in the end the meaningful subgrid parameter is only the combination GLYPH<15> T .", '2.2.2 Bondi accretion from hot gas': "Hot gas can also accrete onto black holes, but in this case Bondi accretion is a more appropriate physical mechanism since the hot gas is typically more spherically distributed. Thus we also account for Bondi accretion, but only from non-ISM gas with a temperature of T > 10 5 K.  \nBondi-Hoyle-Lyttleton accretion is computed via the standard Bondi (1952) formula:  \nGLYPH<219> M Bondi = GLYPH<15> m 4 GLYPH<25> G 2 M 2 BH GLYPH<26> ' v 2 + c 2 s ' 3 GLYPH<157> 2 (5)  \nwhere GLYPH<26> is the mean density of hot ( T > 10 5 K) gas computed within the black hole accretion kernel, c s is the kernelaveraged sound speed of the hot gas, and v is the kernelaveraged velocity of the hot gas relative to the black hole. In practice, we neglect the gas relative velocity and set v = 0 , since the dynamics of the black hole particle is controlled by the repositioning algorithm ( § 2.4). We do not include any boost factor, since hot gas is likely distributed quite smoothly over the size of the black hole kernel. For consistency with the gravitational torque rate, we suppress Bondi accretion by the same efficiency GLYPH<15> m = 0 : 1 .", '2.2.3 Numerical implementation': "The total accretion rate for a given black hole is then  \nGLYPH<219> M BH = ' 1 GLYPH<0> GLYPH<17> ' GLYPH<2> ' GLYPH<219> M Torque + GLYPH<219> M Bondi ' ; (6)  \nwhere we adopt a constant radiative efficiency GLYPH<17> = 0 : 1 (e.g. Yu & Tremaine 2002). We apply an overall limit to the accretion rate, based on the given black hole's Eddington accretion rate. For torque-limited accretion, we apply a limit of 3 times the Eddington rate, based on the idea that non-spherical accretion can potentially exceed Eddington as occasionally observed particularly at higher redshifts (Martınez-Aldama et al. 2018) and consistent with recent accretion disk simulations (e.g. Jiang et al. 2014). For Bondi accretion, we apply a strict Eddington limit, since this is intended to represent quasi-spherical accretion from hot gas where the Eddington limit is directly applicable.  \nNumerically, we follow Springel et al. (2005) and track separately the physical black hole mass, which grows continuously according to equation 6, and the dynamical black hole particle mass, which grows by stochastically accreting a fraction f m of the mass of gas particles within R 0 with a probability that statistically satisfies mass conservation and the desired mass outflow rate in AGN winds (see below). A time step limiter is imposed on black hole particles such that black holes do not grow by more than 0.1 % of their current mass in a single time step.", '2.3 Black hole feedback': "We incorporate a kinetic subgrid model for black hole feedback, along with X-ray energy feedback. The motivation  \nfor the kinetic feedback model comes from the observed dichotomy in black hole growth modes that is reflected in their outflow characteristics (e.g. Heckman & Best 2014): A 'radiative mode'at high Eddington ratios ( f Edd GLYPH<17> GLYPH<219> M BH GLYPH<157> GLYPH<219> M Edd > GLYPH<24> few percent), in which AGN are seen to drive multi-phase winds at velocities of GLYPH<24> 1000 kms GLYPH<0> 1 that include molecular and warm ionised gas (e.g. Sturm et al. 2011; Maiolino et al. 2012; Perna et al. 2017a); and a 'jet mode' at low Eddington ratios, where AGN mostly drive hot gas in collimated jets at high velocities (of order GLYPH<24> 10 4 kms GLYPH<0> 1 ) that in some circumstances are seen to inflate super-virial temperature bubbles in surrounding hot gas (e.g. McNamara & Nulsen 2007; Fabian 2012). This dichotomy also appears in radio jet activity ('high excitation'vs.'low excitation'radio galaxies) above and below roughly a percent of the Eddington rate (Best & Heckman 2012), with the former tending to be found in lower-mass, bluer host galaxies and the latter in massive early-types. Simba 's AGN feedback model is designed to directly mimic the energy injection into largescale surrounding gas from these two modes, using purely bipolar feedback with velocities and temperatures taken as much as possible from AGN outflow observations. We also include X-ray heating from black holes broadly following the model introduced by Choi et al. (2012), modified to operate at the lower resolution of Simba 's cosmological simulations. We now describe these subgrid models in more detail.", '2.3.1 Kinetic feedback': "For the highf Edd mode outflows, we choose an outflow velocity based on ionised gas linewidth observations of X-ray detected AGN from SDSS by Perna et al. (see Fig. 8 of 2017a), which we parameterise in terms of the black hole mass M BH (in M GLYPH<12> ) as  \nv w ; EL = 500 + 500 ' log M BH GLYPH<0> 6 'GLYPH<157> 3 kms GLYPH<0> 1 ; (7)  \nwe refer to these as radiative AGN winds. The gas is ejected without modifying its temperature, meaning that it is ejected at the ISM temperature given by our ISM pressurisation model. This is generally consistent with observations suggesting typical electron temperature of GLYPH<24> 10 4 K for the ionised gas outflows (Perna et al. 2017b). This is broadly similar to the AGN feedback model implemented into Gizmo by Angl'es-Alc'azar et al. (2017a), except here with a variable outflow velocity.  \nIf the Eddington ratio is below f Edd < 0 : 2 , we begin to transition to a jet mode with a velocity that becomes increasingly strong as f Edd drops, as follows:  \nv w ; jet = v w ; EL + 7000 log ' 0 : 2 GLYPH<157> f Edd ' kms GLYPH<0> 1 ; (8)  \nwith a cap to the velocity increase at 7000 kms GLYPH<0> 1 . In this way, full speed jets are achieved only once f Edd < GLYPH<24> 0 : 02 . To trigger jet mode, we also include an additional criterion requiring M BH > M BH ; lim , motivated by observations showing that radio jets only arise in galaxies with velocity dispersions corresponding to black holes with M BH > GLYPH<24> 10 8 M GLYPH<12> (Bariˇsi'c et al. 2017). To be conservative we choose M BH ; lim lower than this, namely M BH ; lim = 10 7 : 5 M GLYPH<12> . Physically, this mass limit is implemented in order to prevent small black holes that temporarily have low accretion rates from driving highpowered jets.  \nBased on observations of AGN outflows and the inferred momentum and energy input (e.g. Fiore et al. 2017; Ishibashi et al. 2018), we set the amount of material ejected in the AGN winds in order to obtain a momentum input of GLYPH<219> P out = 20 L GLYPH<157> c , where L = GLYPH<17> GLYPH<219> M BH c 2 is the bolometric luminosity of the AGN, GLYPH<17> = 0 : 1 , and c is the speed of light. This value is kept constant for both modes, resulting in the mass loading factor in AGN winds scaling inversely with the outflow velocity. For our parameter choices, a black hole with M BH = 10 9 M GLYPH<12> in the highf Edd mode injects outflows with v w ; EL = 1000 kms GLYPH<0> 1 , mass loading GLYPH<219> M out ; EL GLYPH<157> GLYPH<219> M BH GLYPH<25> 600 , and kinetic energy efficiency GLYPH<219> E kin ; EL GLYPH<25> 0 : 03 L , while in the lowf Edd mode at full jet speed reaches v w ; jet = 8000 kms GLYPH<0> 1 , GLYPH<219> M out ; jet GLYPH<157> GLYPH<219> M BH GLYPH<25> 75 , and GLYPH<219> E kin ; jet GLYPH<25> 0 : 3 L .  \nParticles are selected to be ejected randomly from within the black hole accretion kernel, with probability  \np j = 1 GLYPH<0> f m f m GLYPH<2> w j m j GLYPH<2> GLYPH<219> M BH GLYPH<1> t ; (9)  \nwhere w j is a kernel weight ( GLYPH<6> j w j = 1 ) and f m is the fraction of mass accreted by the black hole. The desired mass loading factor relative to the black hole accretion rate ( GLYPH<219> M out GLYPH<157> GLYPH<219> M BH = ' 1 GLYPH<0> f m 'GLYPH<157> f m ) is achieved by setting f m such that:  \nGLYPH<219> P out L GLYPH<157> c = 20 = v w GLYPH<17> c GLYPH<18> 1 GLYPH<0> f m f m GLYPH<19> : (10)  \nAll outflows are ejected in a purely bipolar fashion. That is, we eject gas elements in a direction GLYPH<6> parallel to the angular momentum vector of the inner disk that we use to compute the black hole accretion (typically, the 256 nearest gas particle neighbours to the black hole). We assume zero opening angle for all winds; this is probably conservative for the radiative mode winds, as the opening angles are likely to be wider, but for the jet mode it is a good approximation to observed highly collimated jets. Even in the case of initially spherical radiative winds, it is likely that there is substantial collimation from the inner galaxy disk on scales that we cannot resolve in our cosmological runs, so the assumption of collimated winds is likely to be closer to correct. Since the wind particles are launched from their current location, this results in a collimated outflow with a small but finite extent ( < GLYPH<24> 1 kpc). We note that the outflow direction can precess owing to variations in the inner disk, but is in practice typically stable over tens to hundreds of Myr. Hence any effect of 'sphericalising' the jet energy input on super-galactic scales is done self-consistently via the hydrodynamic interactions of the outflows with ambient gas at larger scales.  \nSince jets are observed to carry very hot gas, we raise the temperature of jet mode (only) outflows to the virial temperature of the halo, specifically T vir = 9 : 52 GLYPH<2> 10 7 ' M halo GLYPH<157> 10 15 M GLYPH<12> ' 1 GLYPH<157> 3 K (Voit 2005). This choice is motivated by observations showing that jets contain mostly synchrotron-emitting plasma, and eventually thermalise their energy into surrounding hot gas at around T vir (Fabian 2012). The extra energy input required for this is typically less than a few percent of the jet kinetic energy, so it does not figure significantly into the overall energy budget.  \nWe apply a short hydrodynamic and radiative cooling decoupling time of 10 GLYPH<0> 4 t H to the outflowing wind gas elements, where t H is the Hubble time at launch. This is in order to avoid further entrainment within the unresolved ISM close to the black hole, since the mass loading is accounted  \nfor from the assumption of constant momentum input of 20 L GLYPH<157> c . This also avoids some numerical inaccuracies from high Mach number shocks in very dense gas. We note that for the jet mode, this can result in a decoupled distance of up to tens of kpc at the present epoch. Hence the jet energy begins to be deposited at a distance comparable to the extent of observed radio lobes. Gizmo employs a Durier & Dalla Vecchia (2012) timestep limiter in order to ensure proper interactions of the high-speed winds and their surrounding gas as they recouple.  \nOur model has similarities to the two-mode thermal and kinetic AGN feedback model employed in IllustrisTNG (Weinberger et al. 2017). The main differences are as follows: (i) Illustris-TNG uses Bondi accretion rather than torque-limited accretion for cold, rotationally supported gas. (ii) Illustris-TNG uses spherical thermal feedback at high f Edd rather than kinetic feedback. This may owe to the fact that torque-limited accretion does not require selfregulation, while Bondi accretion does (Angl'es-Alc'azar et al. 2013), and hence Simba can employ non-spherical feedback during the growth phase and yield black holes consistent with observed scaling relations. (iii) At lowf Edd , IllustrisTNG randomises the direction of the jets at each timestep, rather than always ejecting jets perpendicular to the inner disk. Our approach seems more physically motivated, since jets are not known to dramatically precess on timescales of GLYPH<24> Myr (though they may occasionally do so over hundreds of Myr). (iv) Illustris-TNG uses, at maximum, 200% of the AGN bolometric luminosity (assuming a 10% radiative efficiency), whereas for our model, the maximum is approximately a third of the bolometric luminosity. Despite these and other minor differences, the use of two-mode AGN feedback as in Illustris-TNG and Simba seems to be a reasonably successful approach in state of the art AGN feedback models.", '2.3.2 X-ray feedback': "We include energy input into surrounding gas from X-rays off the accretion disk, as motivated and discussed in Choi et al. (2012). Specifically, we compute the volume heating rate owing to X-rays following equation 12 of Choi et al. (2012), assuming (as they did) a radiative efficiency of 0.1. We only apply this heating when jet mode is activated, as the lower velocity winds typically arise in more gaseous blue galaxies for which radiative losses would be severe (Best & Heckman 2012). To be more explicit, we assume that more gas-rich galaxies are able to absorb and radiate away the X-ray energy, so we implement a gas fraction threshold such that we only apply X-ray heating if f gas < 0 : 2 , where f gas = M gas GLYPH<157> M GLYPH<3> as computed by our galaxy finder, and we only include X-ray heating in galaxies with full-velocity jets.  \nThe X-ray heating is applied to gas within the black hole accretion kernel, scaled inversely with the square of the distance between the gas elements and the black hole, including Plummer softening based on the gas's smoothing length in order to mitigate large energy deposition in gas close to the black hole. For non-ISM gas, we directly increase the gas's temperature according to the heating flux at the gas's position. For ISM gas, because depositing such heat into a low-resolution, pressurised ISM as we assume in Simba would cool quickly and not be physically well moti-  \nvated, we instead take half of the X-ray energy and apply it a radial outwards kick; the remainder is added as heat. We further limit the total energy input in both kinetic and thermal forms to the overall available heating energy; if while looping over BH neighbors the X-ray energy input exceeds this value, then no further X-ray heating is done for that black hole at that timestep. The X-ray heating has a fairly minimal effect on the galaxy mass function, but it provides an important additional energy input to more fully quench massive galaxies, as we discuss in § 4.", '2.4 Black hole seeding and dynamics': "We use the on-the-fly FOF algorithm to seed black holes in galaxies dynamically during the simulation (e.g. Di Matteo et al. 2008; Angl'es-Alc'azar et al. 2017a). If a galaxy reaches a stellar mass M GLYPH<3> > GLYPH<13> BH GLYPH<2> M seed and it does not already contain a black hole particle, then the star particle closest to the center of mass of the galaxy is converted into a black hole particle. For our fiducial simulations, we employ M seed = 10 4 M GLYPH<12> GLYPH<157> h and GLYPH<13> BH = 3 GLYPH<2> 10 5 , which places black holes in galaxies with M GLYPH<3> glyph[greaterorsimilar] 10 9 : 5 M GLYPH<12> .  \nThis somewhat high stellar mass threshold for black hole seeding is motivated by recent simulations from the FIRE project, showing that stellar feedback strongly suppresses black hole growth in low mass galaxies by evacuating the nuclear gas reservoir on < 100 pc scales (Angl'es-Alc'azar et al. 2017c). A qualitatively similar effect was also found in EAGLE (Bower et al. 2017; McAlpine et al. 2018) and Ramses -based simulations (Dubois et al. 2015; Habouzit et al. 2017), though their use of Bondi accretion may inhibit the growth of low mass black holes even in the absence of resolved stellar feedback (owing to the strong dependence GLYPH<219> M Bondi / M 2 BH ). Owing to poorer cosmological resolution as well as Simba 's decoupled kinetic winds that explicitly avoids interaction of star formation feedback with ISM gas, Simba does not reproduce this effect self-consistently. Hence we simply seed black holes in the regime where they are expected to grow more efficiently. We note that our results are insensitive to the exact choice of M seed and stellar mass threshold (Angl'es-Alc'azar et al. 2015).  \nWe assume that dynamical friction is efficient enough to maintain black holes near the host galaxy's center. At every time step, black hole particles are repositioned to the location of the potential minimum within the FOF host group, if it is found within a distance < 4 GLYPH<2> R 0 , where R 0 is the size of the black hole kernel used to compute the accretion rate. The black hole particle velocity is then set to the center of mass velocity of the FOF group. While current cosmological large volume simulations cannot self-consistently model the dynamics of black holes within galaxies, this algorithm is sufficient to capture the mass growth and feedback of 'wellbehaved'central black holes (see Tremmel et al. 2017, for an attempt to include sub-grid dynamic friction for black holes in cosmological simulations). Any two black holes located within R 0 are allowed to merge instantaneously if their relative velocity is lower than three times their mutual escape velocity.", '2.5 Dust production, growth and destruction': 'Simba includes a dust physics module to track the lifecy-  \ncle of cosmic dust. In this implementation, dust is passively advected following the gas particles. This treatment is essentially accurate, as gas drag is usually able to decelerate grains on very short time scales especially when the radiative pressure is weak, so the drift cannot be resolved in our simulations. We additionally assume all dust grains have the same physical properties with a fixed radius a = 0 : 1 GLYPH<22> m . We ignore the active dust cooling, which will be applied in future work.  \nDust is produced by condensation of metals from ejecta of SNe and AGB stars. We follow the prescription described in the work of Dwek (1998), with updated condensation efficiencies based on recent studies. In the following, m j i ; d refers to the dust mass of the i th element (C, O, Mg, Si, S, Ca, Fe) produced by the j th stellar process (Type II SNe or AGB stars), whereas m j i ; ej refers to the mass of ejecta from the j th process.  \nThe mass of dust produced by AGB stars with a carbonto-oxygen mass ratio C/O > 1 is expressed as  \nm AGB i ; d = ( GLYPH<14> AGB C \' m AGB C ; ej GLYPH<0> 0 : 75 m AGB O ; ej \' ; i = C 0 ; otherwise ; (11)  \nwhere GLYPH<14> AGB i is the condensation efficiency of element i for AGB stars. The mass of dust produced by AGB stars with a carbon-to-oxygen mass ratio C/O < 1 is expressed as  \nm AGB i ; d = 8 > > > > < > > > > : 0 ; i = C 16 " i = Mg ; Si ; S ; Ca ; Fe GLYPH<14> AGB i m AGB i ; ej ; i = O GLYPH<14> AGB i m AGB i ; ej ; otherwise ; (12)  \nwhere GLYPH<22> i is the mass of element i in atomic mass units. The mass of dust produced by Type II SNe is described as  \nm SNII i ; d = 8 > > < > > : 16 " i = Mg ; Si ; S ; Ca ; Fe GLYPH<14> SNII i m SNII i ; ej ; i = O GLYPH<14> SNII i m SNII i ; ej ; otherwise ; (13)  \nwhere GLYPH<27> SNII i is the condensation efficiency of element i for Type II SNe.  \nWe take a fixed dust condensation efficiency GLYPH<14> AGB i ; dust = 0 : 2 based on the theoretical models of Ferrarotti & Gail (2006). Guided by computations of Bianchi & Schneider (2007), we choose the dust condensation efficiency of Type II SNe GLYPH<14> SNII i ; dust = 0 : 15 to match the low-metallicity end of the observed relation between dust-to-gas mass ratios (DGR) and gas-phase metallicities (R\'emy-Ruyer et al. 2014). We omit the condensation of Type Ia SNe ejecta, as recent work suggests that Type Ia SNe are not significant sources of dust production (see Nozawa et al. 2011; Dwek 2016; Gioannini et al. 2017). This is different from McKinnon et al. (2016) and Popping et al. (2017) where Type Ia SNe are assumed to have the same condensation efficiency as Type II SNe.  \nOnce dust grains are produced, they can grow by accreting gas-phase metals. Derived by Dwek (1998), the growth rate of grain radius can be expressed as:  \nGLYPH<18> d M dust d t GLYPH<19> grow = GLYPH<18> 1 GLYPH<0> M dust M metal GLYPH<19> GLYPH<18> M dust GLYPH<28> accr GLYPH<19> ; (14)  \nwhere M metal is the total mass of dust and local gas-phase metals. Following Hirashita (2000) and Asano et al. (2013)  \nwhich assume the accretion is a two-body collisional process, the accretion time scale GLYPH<28> accr is  \nGLYPH<28> accr = GLYPH<28> ref GLYPH<18> GLYPH<26> ref GLYPH<26> g GLYPH<19> GLYPH<18> T ref T g GLYPH<19> GLYPH<18> Z GLYPH<12> Z g GLYPH<19> : (15)  \nwhere GLYPH<26> g , T g and Z g are the local gas density, temperature and metallicity, respectively. GLYPH<26> ref , T ref and Z ref are the reference values correspondingly. We take GLYPH<26> ref = 100 H atoms cm GLYPH<0> 3 , T ref = 20 K and GLYPH<28> ref = 10 Myr. Inclusion of the multiplier \' Z GLYPH<12> GLYPH<157> Z g \' , unlike McKinnon et al. (2017), is integral to reproduce the observed relation between the dust to gas ratio and gas-phase metal abundance ( § 3.9).  \nDust grains can be eroded by colliding with thermally excited gas especially in hot halos. A number of works have calculated the thermal sputtering rate in detail (e.g. Barlow 1978; Draine & Salpeter 1979; Tielens et al. 1994). In this work, we adopt an analytic approximation of the growth rate of grain radii of Tsai & Mathews (1995) (also adopted by McKinnon et al. 2017 and Popping et al. 2017) described as  \nGLYPH<18> d a d t GLYPH<19> sp = GLYPH<0> a GLYPH<28> sp ; (16)  \nwhere the sputtering time scale  \nGLYPH<28> sp = a GLYPH<12> GLYPH<12> GLYPH<12> GLYPH<12> d a d t GLYPH<12> GLYPH<12> GLYPH<12> GLYPH<12> GLYPH<0> 1 GLYPH<24> \' 0 : 17Gyr \' GLYPH<18> a 0 : 1 GLYPH<22> m GLYPH<19> GLYPH<18> 10 GLYPH<0> 27 g cm GLYPH<0> 3 GLYPH<26> g GLYPH<19> GLYPH<20> GLYPH<18> T 0 T g GLYPH<19> ! + 1 GLYPH<21> ; (17)  \nwhere ! = 2 : 5 controls the low-temperature scaling of the sputtering rate and T 0 = 2 GLYPH<2> 10 6 K is the temperature above which the sputtering rate flattens. The corresponding dust mass changes as  \nGLYPH<18> d M dust d t GLYPH<19> sp = GLYPH<0> M dust GLYPH<28> sp GLYPH<157> 3 (18)  \nBecause SN blast waves are not resolved in our simulations, we implement an additional dust destruction mechanism by SN shocks which enhance inertia and the thermal sputtering of dust grains (Dwek & Scalo 1980; Seab & Shull 1983; McKee et al. 1987; McKee 1989). We follow the prescription outlined by McKinnon et al. (2016) in this work. The growth rate of the dust particle mass due to SN destruction is  \nGLYPH<18> d M dust d t GLYPH<19> de = GLYPH<0> M dust GLYPH<28> de ; (19)  \nwhere the characteristic time scale GLYPH<28> de is  \nGLYPH<28> de = M g GLYPH<15> GLYPH<13> M s ; (20)  \nwhere M g is the local gas mass, GLYPH<15> = 0 : 3 is the efficiency with which grains are destroyed in SN shocks (McKee 1989), GLYPH<13> is the local SN II rate, and M s is the mass of local gas shocked to at least 100 km/s. Considering that our simulations are unable to resolve multi-phase ISM, we apply the Sedov-Taylor solution to a homogeneous medium of n H = 0 : 13 H atoms cm GLYPH<0> 3 (the minimum SF threshold density of our simulations) and obtain  \nM s = 6800 E SNII ; 51 GLYPH<18> v s 100 km s GLYPH<0> 1 GLYPH<19> ; (21)  \nwhere E SNII ; 51 is the energy released by a SN II in units of 10 5 1 erg, and v s GLYPH<24> 100 km s GLYPH<0> 1 is the shock wave speed.  \nWe additionally destroy dust, as well as molecular hydrogen, completely in hot winds and during star formation ( § 2.1) and in any gas that is impacted by AGN X-ray heating or jets ( § 2.3). This is done instantaneously, with all dust mass and metals being returned to the gaseous phase. Note that we do not do this for cold star-forming winds or AGN winds in the high-Eddington mode, so these outflows carry molecular gas and dust out of the galaxy. We leave for future work an investigation into whether this reproduces observations of AGN-driven molecular outflows (e.g. Sturm et al. 2011) and circum-galactic dust (e.g. Peek et al. 2015).', '2.6 Runs and analysis': "The primary Simba runs have 1024 3 dark matter particles and 1024 3 gas elements. We are running four volumes: 100 h GLYPH<0> 1 Mpc down to z = 0 , 50 h GLYPH<0> 1 Mpc to z = 1 , 25 h GLYPH<0> 1 Mpc to z = 2 , and 12 : 5 h GLYPH<0> 1 Mpc to z = 5 . All runs have identical input physics, begin at z = 249 , and assume a Planck Collaboration et al. (2016) concordant cosmology of GLYPH<10> m = 0 : 3 , GLYPH<10> GLYPH<3> = 0 : 7 , GLYPH<10> b = 0 : 048 , H 0 = 68 kms GLYPH<0> 1 Mpc GLYPH<0> 1 , GLYPH<27> 8 = 0 : 82 , and n s = 0 : 97 . Other parameters such as the minimum gravitational softening length and mass resolutions are listed in Table 1. In this paper we will only present results from the main 100 h GLYPH<0> 1 Mpc run, as the other runs are at various stages of completion.  \nWe will also explore parameter space and compare to our previous Mufasa simulations using 50 h GLYPH<0> 1 Mpc runs with 2 GLYPH<2> 512 3 particles that match Mufasa 's size. We run a full physics Simba simulation at this resolution, and in order to directly assess the impact of our new quenching feedback modules, namely jet and X-ray feedback, we also run a 'Nojet' simulation where these modules are turned off, and a 'No-Xray' run where jets are kept on but X-ray feedback is turned off. All other input physics in these runs, including stellar feedback and radiative mode black hole feedback, remains identical to that in Simba .  \nTo analyse the simulation outputs, we employ a suite of tools as described below. First, galaxies are identified using a friends-of-friends galaxy finder, assuming a spatial linking length of 0.0056 times the mean inter-particle spacing (equivalent to twice the minimum softening length). In our tests, this gives very similar results to the more comprehensive Spline Kernel Interpolative Denmax (SKID) galaxy finder. Galaxy finding is applied to all stars and black holes plus all gas elements with a density above the minimum SF threshold density of n H > 0.13 H atoms cm GLYPH<0> 3 ; this captures all the stars and molecular gas in galaxies. Black holes are assigned to the galaxy to which they are most gravitationally bound; large galaxies can have many black holes. We take the central black hole to be the most massive black hole in the galaxy, and use this when we discuss black hole masses. In most cases, the other black holes are very small and add no significant black hole mass compared to the central one.  \nBecause significant amounts of neutral hydrogen can lie in an extended configuration beyond the star-forming region of galaxies, we assign H i to galaxies in a separate step. To do this, we consider all gas elements with H i fractions above 0.001, and assign them to the galaxy to which they are most  \n<!-- image -->  \nFigure 2. Temperature map projected through a random 10 Mpc/h slice from a 50 Mpc/h Simba volume, at z = 2 (left) and z = 0 (right). At z = 2 , warm-hot gas traces large-scale filaments, with energetic bipolar outflows owing to jets evident from the nodes where the most massive galaxies and black holes reside. At z = 0 , high-speed AGN outflows have shocked the IGM gas throughout much of this volume to well beyond the virial radii of halos, with cooler dense filamentary structures penetrating the hot gas.  \n<!-- image -->  \nTable 1. The Simba simulation suite.  \n- b Minimum gravitational softening length in comoving h GLYPH<0> 1 kpc .\n- c Ending redshift (all begin at z = 249 ).\n- d Initial gas element mass resolution in M GLYPH<12> .\n- e Dark matter particle mass resolution in M GLYPH<12>  \n.  \n- f Minimum stellar mass of a resolved galaxy in M GLYPH<12> .  \ngravitationally bound, i.e. its kinetic energy relative to the galaxy's center of mass velocity minus the potential energy from the galaxy at the gas element's location is minimised.  \nHalos are identified on the fly during the simulation run using a 3-D friends-of-friends algorithm within Gizmo , which is identical to the one in Gadget-3 written by V. Springel. The linking length is taken to be 0.2 times the mean inter-particle spacing. We do not identify or consider sub-halos in this work.  \nGalaxies and halos are cross-matched in post-processing using the yt -based package Caesar , which outputs a single hdf5 catalogue containing all galaxy and halo information with many key properties pre-computed, as well as particle lists of individual galaxies and halos so that any other properties can be computed via user-written python scripts.  \nAll results shown here are obtained from the Caesar catalogs generated from simulation snapshots at specified redshifts. We output 151 snapshots to z = 0 , 105 to z = 1 , and 78 to z = 2 . Each snapshot is GLYPH<25> 250 GB in size, and the Caesar catalogues are typically GLYPH<24> 15 GB each.", '3 RESULTS': "In this section we provide a comprehensive suite of predictions for Simba for a range of key global galaxy properties. The purpose is to ascertain how well Simba reproduces observed galaxy stellar and gas properties that have historically provided stringent constraints on feedback models in previous simulations, and thereby demonstrate the suitability of Simba as a platform to study detailed galaxy evolution.  \nTo begin, we show in Figure 2 a projected temperature map from the (50 Mpc/h) 3 , 512 3 Simba simulation. The slice shown is 10 Mpc/h thick, and is arbitrarily chosen to contain representative structures in the volume. At z = 2 , the familiar Cosmic Web is evident as traced out by warmer gas arising from mild shock heating on filamentary structures. Closer inspection reveals the earliest AGN jets becoming active, with characteristic bipolar outflows that are typically perpendicular to the large-scale filaments. But these large early black holes are sparse, and most of the IGM is unaffected by feedback. In contrast, by z = 0 (right panel), a significant volume of the IGM has been heated to high  \nFigure 4 shows the GSMF at z = 0 : 1 ; 1 ; 2 ; 3 ; 4 ; 6 from Simba (green lines). Observational data is shown at z = 0 from Bernardi et al. (2017). At z = 1 ; 2 ; 3 we show observations from Tomczak et al. (2014) combining CANDELS and zFOURGE data, while at z = 4 ; 6 we show observations based on CANDELS from Song et al. (2016). We also show the GSMF of central galaxies only, subdivided into star-forming (SF) and quenched (Q) samples at a specific SFR = 10 GLYPH<0> 1 : 8 + 0 : 3 z Gyr GLYPH<0> 1 . Error bars are shown from jacknife re-sampling over eight simulation sub-octants. Finally, we show the results from the EAGLE simulation as the dotted cyan line at selected redshifts.  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \n<!-- image -->  \nFigure 3. Examples of the molecular gas (left) and stellar (right) surface density distributions in star-forming disk galaxies with M GLYPH<3> GLYPH<25> 4 : 7 GLYPH<2> 10 10 M GLYPH<12> at z = 0 (top four panels) and z = 2 (bottom four panels), showing face-on and edge-on views. At z = 0 there is a thin, well-ordered disk in both H i and H 2 , while the highz galaxy is clumpier and thicker.  \n<!-- image -->  \ntemperatures from AGN feedback, and the hot bubbles encroach upon regions untouched by AGN feedback containing the canonical warm filaments. These bubbles are reasonably spherical since they arise from clustered massive galaxies, each one ejecting jets that are relatively stable in direction but overlap quickly with neighboring outflows. In some cases individual bipolar jets and the resulting bow shocks can still be picked out. Such a dramatic impact on the IGM may have significant consequences for the ionisation state of diffuse neutral hydrogen and the statistics of Lyman alpha forest absorbers (Kollmeier et al. 2014), as well as the diffuse IGM pressure measurable via the Sunyaev-Zel'dovich effect (Lim et al. 2018); we will explore these in future work. In this paper, we focus on the demographics of the galaxy population predicted by Simba .  \nFigure 3 shows some examples of individual galaxies. We choose a Milky Way-sized disk galaxy at z = 0 , with M GLYPH<3> GLYPH<25> 4 : 7 GLYPH<2> 10 10 M GLYPH<12> and SFR = 1 : 3 M GLYPH<12> yr GLYPH<0> 1 , and show the face-on (upper row) and edge-on (lower row) views, in both H 2 surface density (left) and stellar mass surface density  \n<!-- image -->  \n(right). The z = 2 galaxy shown in the bottom four panels has essentially the same M GLYPH<3> , but with SFR = 45 M GLYPH<12> yr GLYPH<0> 1 that is typical of a main sequence galaxy at Cosmic Noon. The z = 0 disk is a grand design spiral, with a thin cold gas distribution. There is a small central hole in cold gas that owes to the AGN feedback from its 6 GLYPH<2> 10 7 M GLYPH<12> black hole accreting at 0 : 005 M GLYPH<12> yr GLYPH<0> 1 . The stellar distribution does not show the spiral structure owing to the relatively low resolution of Simba , compared to zooms or higher-resolution simulations such as Illustris-TNG and EAGLE. The z = 2 system shows more prominent star forming clumps and a thicker gas distribution, and is overall more compact (note the scale bar). While Simba 's numerical resolution smooths out many of the detailed internal features, this shows that it still produces galaxies that have features broadly like star-forming disk galaxies in the real Universe. We do not show more massive quenched examples, but as expected they tend to be elliptical in their stellar morphology, with little cold gas.", '3.1 Galaxy stellar mass functions': "Since galaxies are a collection of stars, the most basic property of a galaxy is its stellar mass. Given that the concordance cosmological model strongly constraints the halo mass function, the galaxy stellar mass function (GSMF) thus characterises the efficiency by which halos convert their baryons into stars. It is well established that (under the abundance-matching ansatz) the stellar-to-halo mass ratio drops quickly to low and high masses away from the peak at L GLYPH<3> (e.g. Moster et al. 2013; Behroozi et al. 2013), and current models attribute this to self-regulation by star formationdriven feedback below L GLYPH<3> and quenching of galaxies due to AGN feedback above L GLYPH<3> (Somerville & Dav'e 2015). Since the GSMF is reasonably well measured over much of cosmic time (albeit with non-trivial systematic uncertainties; Mobasher et al. 2015), it represents a stringent test for the key feedback modules of a galaxy formation model. Indeed, simulations these days including Simba tend to use the z = 0 GSMF as a primary constraint to tune feedback models.  \nSimba produces generally good agreement with the observed GSMFs at all redshifts, overall comparably well to EAGLE. There is excellent agreement at z GLYPH<21> 3 , especially given the systematic uncertainties in stellar mass determinations at higher redshifts (Mobasher et al. 2015). At z = 2 , there starts to be a slight excess at the massive end in Simba . This may owe to insufficient quenching of the most massive galaxies, or may represent an underestimate of the observed GSMF owing to selection effects in the rest-optical surveys used for the GSMF determinations which can miss massive dusty galaxies.  \nFigure 4. Stellar mass function evolution from z = 6 ! 0 , compared to observations as indicated in the legends. Green band shows the results from all Simba galaxies, with the spread computed from jackknife resampling the 8 simulations sub-octants. Red and blue dashed lines show the mass functions of central galaxies below and above sSFR = 10 GLYPH<0> 1 : 8 + 0 : 3 z Gyr GLYPH<0> 1 , respectively. Cyan dotted line shows the results from EAGLE for comparison.  \n<!-- image -->  \nAt lower redshifts, there is a clear truncation at the massive end, but a mild overproduction of the most massive galaxies remains all the way to z = 0 . Like EAGLE, Simba under-predicts the GSMF around M ? by a factor of up to two; this was an advantage of Mufasa that is unfortunately not retained in Simba . This highlights that it continues to be a challenge to achieve such a sharp turndown in the GSMF using a physically-motivated AGN feedback model.  \nThe overproduction of the most massive galaxies could owe to a number of effects. First off, there are numerical uncertainties in quantifying the most massive systems, because they tend to have large extended envelopes of stars and many satellites that, owing to poor resolution, can be overmerged into the central object either during the dynamical evolution or during the post-processing galaxy finding stage. These tend to artificially boost the mass in the simulated massive galaxies. One way to mitigate this is to compare the stellar mass within fixed apertures to data, which Schaye et al. (2015) showed using EAGLE can significantly reduce the mass of M GLYPH<3> > GLYPH<24> 10 11 M GLYPH<12> objects. There are also in-  \ncreased observational uncertainties at the massive end. For instance, it is a matter of debate as to how much of the surrounding stars should be classified as part of the central galaxy and how much should be intracluster light; this can strongly impact the stellar mass (Kravtsov et al. 2018). There is also the issue of the stellar initial mass function (IMF) - stellar population (Conroy et al. 2013) and dynamical (Cappellari et al. 2013) studies suggest that the most massive galaxies have bottom-heavy IMFs relative to Milky Way-like galaxies, which can result in the stellar mass being underestimated by a factor of 2 or more for the most massive systems. Finally, there is an issue particular to this Simba run - it turns out, in the 100 h GLYPH<0> 1 Mpc volume, by z = 0 , the largest halo has a virial mass of M halo = 1 : 16 GLYPH<2> 10 15 M GLYPH<12> , which is larger than expected by about 50% for its volume; this may contribute to the excess of the very most massive galaxies. Hence although at face value there is some disagreement at the massive end in comparing Simba with recent observations, more work must be done to determine whether these  \ndiscrepancies reflect a significant failing of Simba 's input physics.  \nExamining the SF vs. Q samples, we see that massive quenched galaxies begin to appear in significant numbers at z > GLYPH<24> 2 . By z = 1 they outnumber the SF galaxies among the most massive galaxies with M GLYPH<3> > GLYPH<24> 10 11 M GLYPH<12> , and by z = 0 , they dominate at M GLYPH<3> > GLYPH<24> 2 GLYPH<2> 10 10 M GLYPH<12> . The quenched population grows quickly at low redshifts, and the number of massive star-forming galaxies drops quickly since z GLYPH<24> 1 , in broad agreement with observations (e.g. Bell et al. 2004). There are a few very small quenched centrals, but this is likely an artifact of the friends-of-friends halo finder.  \nIn summary, to within systematic uncertainties, Simba produces a GSMF that is in quite good agreement with observations across most of cosmic time, with the possible exeption of the z = 0 massive end. Simba passes this primary check at a level comparable to that seen for Mufasa , EAGLE, and Illustris-TNG (Pillepich et al. 2018). In no small part, this owes to these various models tuning feedback parameters to match such data, but even the fact that such a tuning is now possible is a recent and important step forward for cosmological hydrodynamic simulations. It does mean that the growth of galaxies' stellar component over time is no longer a strong discriminant between current galaxy formation models. Instead, for this we must rely on the many other predicted observables that are not used to tune the models. We now examine some of these predictions for Simba .", '3.2 Star formation rate-stellar mass relation': "Another key barometer of galaxy formation models is the star formation rate-stellar mass (SFR GLYPH<0> M GLYPH<3> ) relation. Unlike the GSMF that is often used as a primary constraint on models, SFR GLYPH<0> M GLYPH<3> is not, making it more of a true prediction of models. The SFR GLYPH<0> M GLYPH<3> relation consists of a star-forming 'main sequence' of galaxies, and a population of quenched galaxies falling below the main sequence that dominates at high masses at later epochs. Getting the balance of these populations in accord with observations over cosmic time, as well as predicting their growth rates, has traditionally been difficult to reproduce in cosmological simulations.  \nDuring Cosmic Noon, it has long been seen that cosmological models tend to underpredict the main sequence amplitude (Daddi et al. 2007; Dav'e 2008; Narayanan et al. 2012; Sparre et al. 2015; Somerville & Dav'e 2015), typically by a factor of 2 GLYPH<0> 3 . Fixing this requires rather substantially changing the star formation histories, not just the overall SFRs, since a multiplicative constant on the SFR will tend to move galaxies along the relation rather than increase its amplitude. There are also potential observationally-oriented systematics that may be overestimating the SFR owing to one or more of many possible factors, such as galaxies being dominated by harder-ionising stellar populations at highz , or having a more top-heavy initial mass function.  \nFigure 5 shows the specific SFR GLYPH<0> M GLYPH<3> relation at z = 0 : 1 (left) and z = 2 : 3 (right) for Simba galaxies. The SFRs are computed as instantaneous SFRs from the gas elements, which corresponds well to the SFR computed from young star particles when averaged over several tens of Myr. The running median (green curve) includes starforming galaxies only ( Simba -SF), defined as before by  \nsSFR > 10 GLYPH<0> 1 : 8 + 0 : 3 z Gyr GLYPH<0> 1 (dotted horizontal yellow line). The error bars show the 1 GLYPH<27> spread around the median value in each bin. Points are colour-coded by the ratio of black hole to stellar mass, with magenta points having higher M BH GLYPH<157> M GLYPH<3> ; points at M GLYPH<3> < GLYPH<24> 10 10 M GLYPH<12> in cyan have no or very small growing black holes, as we will discuss later. Galaxies with very low or zero SFR are plotted near the bottom for visibility. Observations at lowz are shown from the GALEX-SDSSWISE Legacy CatalogSalim et al. (GSWLC; 2016, 2018), shown as grey hexbins, and the running median to the starforming galaxies with the same criterion as above is shown as the black dashed line, along with error bars showing the 1 GLYPH<27> spread around the median. At highz , we show the median sSFR GLYPH<0> M GLYPH<3> relation measured for 2 < z < 2 : 5 galaxies by Whitaker et al. (2014). Finally, results from EAGLE are shown as the magenta dotted line.  \nAt z = 0 (left panel), Simba nicely reproduces the observed GSWLC main sequence slope and amplitude at M GLYPH<3> > GLYPH<24> 10 10 M GLYPH<12> . Below this mass, Simba shows noticeably higher SFRs. This mass corresponds to the onset of massive black holes, as shown by the growing number of magentacoloured points with higher black hole mass for their M GLYPH<3> . Indeed, there is a very strong trend that the galaxies that are quenching are specifically the ones with a high M BH GLYPH<157> M GLYPH<3> ratio; massive galaxies left on the main sequence at z GLYPH<24> 0 in Simba are only those that for some reason have not grown their black hole as rapidly. A similar trend is seen in EAGLE (Matthee & Schaye 2019), which arises owing to a spread in halo formation times (Davies et al. 2019). We will investigate the detailed reasons for this dichotomy in Simba in future work, but for now we note the tight connection between quenching and black holes already appearing in Simba , which will be a recurring theme throughout this paper. The average slope of sSFR GLYPH<0> M GLYPH<3> for star-forming galaxies over the entire mass range plotted is GLYPH<0> 0 : 27 , which is in reasonable agreement with observations (e.g. Noeske et al. 2007; Speagle et al. 2014). The scatter around the main sequence in Simba is 0 : 3 GLYPH<0> 0 : 4 dex, with a mild tendency to drop with M GLYPH<3> ; this is very comparable to that seen in the GSWLC data.  \nAt z = 2 : 2 (right panel), the Simba main sequence generally tracks the observed one from Whitaker et al. (2014), but is low in amplitude by GLYPH<25> GLYPH<2> 2 . This continues the trend in models that the main sequence at Cosmic Noon remains too low, though not quite as strongly as in some previous models. However, Leja et al. (2018) points out that more sophisticated SED fitting applied to the latest datasets can lead to a systematic increase in the inferred M GLYPH<3> while lowering the SFR that results in a combined GLYPH<25> 0 : 3 dex lower sSFR compared to previous determinations. If confirmed, then at face value this would bring Simba 's (and other models') simulated main sequence into agreement with z GLYPH<24> 2 observations at long last.  \nFinally, we show in Figure 6 histograms of the specific SFR, broken up into mass bins of 10 9 < M GLYPH<3> < 10 10 M GLYPH<12> , 10 10 < M GLYPH<3> < 10 11 M GLYPH<12> , and M GLYPH<3> > 10 11 M GLYPH<12> . Solid lines show the results from Simba at z = 0 : 1 , while dashed lines show identically selected galaxies from the GSWLC-D2 catalog (Salim et al. 2018). Galaxies with sSFR GLYPH<20> 10 GLYPH<0> 2 : 5 Gyr GLYPH<0> 1 have been placed in the lowest sSFR bin. There is an overall bimodal distribution, with low-mass galaxies being predominantly star-forming, while massive galaxies are almost uniformly  \nFigure 5. Star formation rate-stellar mass relation at z GLYPH<25> 0 (left) and z GLYPH<25> 2 (right). Points show Simba galaxies, colour-coded by their black hole to stellar mass ratio. The thick green line shows the running median to star-forming galaxies (i.e. above the horizontal dotted yellow line). Observations at z = 0 from GSWLC-X2 are shown as the grey hexbins, with the black dashed line showing the median to galaxies using the same sSFR cut as shown for Simba . The errorbars show the 1 GLYPH<27> spread around the running median value, typically 0 : 3 GLYPH<0> 0 : 4 dex. At highz we show the black dashed line as the best-fit relation for 2 < z < 2 : 5 galaxies from Whitaker et al. (2014). Results from EAGLE are shown as the magneta dotted line for comparison. Simba reproduces the star-forming main sequence at both redshift reasonably well, especially accounting for systematics in highz sSFR determinations, though in small galaxies it appears to overpredict the SFR at z GLYPH<24> 0 and underpredict at z GLYPH<24> 2  \n<!-- image -->  \nFigure 6. Histogram of sSFR in three bins of stellar mass. Solid lines show the results for Simba at z = 0 : 1 , while dotted lines show z GLYPH<24> 0 : 1 observations from GSWLC-D2. All galaxies with sSFR < 10 GLYPH<0> 2 : 5 Gyr GLYPH<0> 1 are placed in the lowestmost bin. There is good agreement, particularly in the quenched fractions in more massive galaxies, though Simba produces somewhat too high sSFRs at lowM GLYPH<3> .  \n<!-- image -->  \nquenched. There is impressive agreement in the lowest sSFR bin, showing that Simba well reproduces the quenched fractions at various M GLYPH<3> . However, the low-mass galaxies in Simba have somewhat too high sSFR values, reflecting the same excess as seen in Figure 5.  \nIn summary, Simba generally reproduces the main sequence of star-forming galaxies as seen at z GLYPH<25> 0 ; 2 , to within current systematic uncertainties. Potential disagreements from data lie mostly at lower masses, where the observations are less certain and more subject to selection effects. The success at M GLYPH<3> > GLYPH<24> 10 10 M GLYPH<12> is encouraging because it suggests  \nthat the balance of quenching and quenched galaxies in this transition mass near M ? is being reproduced roughly correctly in Simba . There is also a strong trend that quenched galaxies at a given M GLYPH<3> tend to have larger fractional black hole masses, which is an interesting prediction that can be tested in future samples of black hole mass measurements in subM ? galaxies.", '3.3 Global star formation rate evolution': 'The evolution of the cosmic SFR density (SFRD) has long been a key test for cosmological galaxy formation models. While proper model comparisons to data can be challenging owing to the variety of different selection effects used to measure SFR over time, the recent compilation by Madau & Dickinson (2014) has provided a homogenised database for the SFRD that can be more robustly compared.  \nFigure 7 shows the comparison of the cosmic SFRD as a function of cosmic age in Simba versus the Madau & Dickinson (2014) compilation. The Simba SFRD values include all the star formation in the volume at each epoch, but we have checked that including only star formation in resolved galaxies ( M GLYPH<3> > 5 : 8 GLYPH<2> 10 8 M GLYPH<12> ) makes a negligible difference.  \nThe overall shape of the predicted SFRD versus time is in good agreement with observations. Simba matches the preset-day SFRD very well, and generally reproduces the order-of-magnitude rise in SFRD towards the peak at z GLYPH<24> 2 . There is a slight tendency for Simba to form more stars globally at earlier epochs, with the peak shifted very slightly towards higher redshift compared to the best-fit line from Madau & Dickinson (2014). The peak SFRD at z GLYPH<24> 2 is also slightly lower than observed, following the trend shown in Figure 5 that Simba has slightly lower main sequence than observed. Despite these minor differences, the overall shape and amplitude is in very good agreement with observations,  \nFigure 7. Star formation rate density evolution versus age of the Universe in Simba (curve), compared to the observational compilation of Madau & Dickinson (2014) (black points and bestfit line).  \n<!-- image -->  \ncomparable to the agreement seen versus other recent simulations such as Illustris-TNG (Pillepich et al. 2018).', '3.4 Neutral and molecular gas fractions': "Simba tracks the neutral (H i ) and molecular (H 2 ) hydrogen separately during its evolution, via sub-grid prescriptions to account for molecular gas production and destruction, and approximate self-shielding that results in neutral gas. Thus Simba lends itself to testing against a complementary set of constraints: the scaling relations of H i and H 2 gas fractions versus M GLYPH<3> . Recent millimetre and radio observation data has greatly expanded our knowledge of gas contents for lowz galaxies, with constraints at higher z promising continued rapid advancement in the near future.  \nFigure 8 shows the scaling relations for H 2 (top) and H i (bottom) mass fractions, versus M GLYPH<3> . The points show individual galaxies at z = 0 colour-coded by their sSFR deviation from the main sequence at that M GLYPH<3> ( GLYPH<1> sSFR), where the main sequence is defined by fitting a running median to the main sequence. The running mean of the gas fractions is shown as the cyan dashed line. Observations are shown from the mass-selected GASS H i survey (Catinella et al. 2012) and its follow-up COLDGASS survey (Saintonge et al. 2017) that obtained H 2 masses from CO emission measurements. We further show the mean predicted trends at z = 1 (green dashed) and z = 2 (magenta dashed), to illustrate how these quantities evolve.  \nOverall, Simba does an excellent job of reproducing the trends in both molecular and neutral gas fractions with stellar mass. There is a hint that the amplitude of both is low by about 0.1 dex, but given observational uncertainties such as the CO-to-H 2 conversion factor which is poorly determined particularly in the low-mass ( M GLYPH<3> < 10 10 M GLYPH<12> ) regime (e.g. Narayanan et al. 2012; Bolatto et al. 2013), as well as theoretical uncertainties in the approximate way that selfshielding is applied, galaxy gas contents can be considered to be a remarkably successful prediction of Simba .  \nWe note that our massive galaxies have a non-trivial  \nFigure 8. Molecular (top) and neutral (bottom) gas fractions M H 2 GLYPH<157> M GLYPH<3> and M HI GLYPH<157> M GLYPH<3> as a function of M GLYPH<3> . The points show z = 0 values from Simba colour-coded by the deviation in sSFR from the star-forming main sequence - bluer points have higherthan typical SFR, redder have lower. A running median at z = 0 is shown as the cyan dashed line. For comparison we show the running medians at z = 1 ; 2 (green, magneta lines). Observations of f H 2 from xCOLDGASS (Saintonge et al. 2017) are shown in the top panel, and observations of f HI from GASS (Catinella et al. 2012) are shown in the bottom panel. Simba predicts gas fraction scalings in good agreement with data, and predicts a small but significant amount of gas even in the most massive quenched systems.  \n<!-- image -->  \namount of cold gas, despite their very low sSFR. This was not the case in Mufasa (Dav'e et al. 2017a), where the most massive galaxies were devoid of essentially any cold gas. Recent observations seem to suggest that, perhaps surprisingly, many massive quenched galaxies contain substantial cold gas fractions of up to a percent or more (e.g. Young et al. 2014), which is consistent with Simba 's predictions. The efficiency of star formation from the molecular gas is, however, low. Simba qualitatively reproduces this trend, possibly because the cold gas generally sits in a more extended configuration where the densities are not as high. Since the star formation rate is proportional to GLYPH<26> 1 : 5 , this means that even if the gas has high molecular content, its low density will curtail star formation relative to the same gas being in a compact configuration. The origin and fate of this cold dense gas is unclear; it could be a transient phase brought in by satellites, or else a stable phase maintained in a more diffuse configuration owing to the presence of hot gas and AGN feedback. We will examine the exact nature of cold gas in quenched galaxies in future work.  \nFinally, the colours of the Simba points in Figure 8 indicate the deviation of a given galaxy from the main sequence. There is a clear trend that galaxies that are more gas-rich at a given M GLYPH<3> have a higher SFR. This is unsurprising for H 2 in our models, given that star formation is tied to the molecular content. It is somewhat more surprising to see  \nthis for the H i fraction, but such a correlation was also seen fairly strongly in Mufasa (Dav'e et al. 2017a), though not as strong as for H 2 (Rafieferantsoa et al. 2018). Such trends have also been noted in observations (Bothwell et al. 2013) and in EAGLE (Lagos et al. 2016).", '3.5 Gas-phase and stellar mass-metallicity relations': "Simba tracks the production and distribution of various heavy elements, through several nucleosynthetic channels. Produced metals can be carried out from galaxies via outflows, which in Simba are typically mildly enriched compared to the mean ISM metallicity (see § 2.1). Additionally, Simba locks individual metals into dust, removing them from the gas phase. Hence predictions for the relationship between galaxy stellar mass and metallicity, which is observed to be among the tightest relations known for galaxies, tests how numerous aspects of Simba work together to establish galaxy metallicities.  \nMetals can be associated with gas, stars, or dust. Measurements of the gas-phase metallicity reflect a balance between relatively pristine inflow and ejection of enriched material via outflows, and thus provide a direct constraint on the mass outflow rate in gas-rich (star-forming) galaxies (Finlator & Dav'e 2008). The stellar metallicity is measured from stellar atmospheric absorption lines that reflect the accumulated metals from both gas and dust that ended up locked in the stars. The inclusion of dust production and destruction model can in principle therefore decouple the stellar and gas phase metallicities. Here we present predictions for the gas-phase and stellar metallicity scaling relations from Simba .  \nFigure 9 shows the gas-phase mass-metallicity relation (gMZR) at z = 0 (left) and z = 2 : 3 (right). The gas-phase metallicity is computed as the SFR-weighted oxygen abundance in all galaxy gas particles, normalized to the solar value of 1.34% (Asplund et al. 2009). Points show central galaxies colour-coded by their deviation from the main sequence, as in Figure 8; black points are satellite galaxies. A running median for star-forming galaxies (sSFR > 10 GLYPH<0> 1 : 8 + 0 : 3 z Gyr GLYPH<0> 1 ) is shown as the dashed green line. Fits to observations at z = 0 are shown from strong emission line fitting (Tremonti et al. 2004, black dashed), stacked measurement of direct metallicities (Andrews & Martini 2013, black dot-dashed), and individual semi-direct metallicities (Yates et al. 2019, black dotted). Observations at z = 2 : 3 are shown from the MOSDEF survey (Sanders et al. 2015).  \nSimba predicts a gas phase mass-metallicity relation that agrees quite well with observations, lying generally in agreement with the range of current observational determinations. The metallicities may be slightly too high at the highest masses, but this turns out to be strongly dependent on the assumed cut for star-forming galaxies; a more stringent cut would lower the massive end fit, and highlights the sensitivity of MZR predictions there to precise selection effects.  \nMufasa produced a gMZR that was slightly too steep, in addition to having an amplitude that roughly agreed with data only because of an arbitrary halving of the SNII yields (Dav'e et al. 2017a). In Simba , metals are locked into  \ndust, increasingly so at higher metallicities and hence larger masses. This likely leads to a suppression of the gas-phase metallicity in massive galaxies and thus a flatter gMZR, as well as a lower amplitude. We will examine the impact of dust on the metal content of galaxies in more detail in future work (Li et al., in prep.).  \nSince z = 2 : 3 , Simba produces more metal evolution at the low-mass end, in general agreement with observations suggesting that the most metal-rich galaxies are in place at early epochs (e.g. Zahid et al. 2014). Star-forming galaxy metallicities are in good agreement with the MOSDEF data, suggesting that the amount of metal evolution from z GLYPH<24> 2 ! 0 is approximately correct in Simba .  \nSimba also shows a clear second parameter dependence on the specific SFR such that at a given M GLYPH<3> , galaxies with lower sSFR tend to have higher metallicities. This has been noted observationally as the Fundamental Metallicity Relation (FMR; Mannucci et al. 2010; Lara-L'opez et al. 2010). The existence of the FMR remains somewhat controversial (Salim et al. 2014; S'anchez et al. 2017; Cresci et al. 2018), but a careful analysis of the MOSDEF data has revealed such a trend at z GLYPH<24> 2 (Sanders et al. 2018). The trend is quite obvious at z = 0 with the bluest galaxies clearly having the lowest metallicities, but is not quite so evident at z = 2 : 3 except at the massive end where a population of quenched galaxies has appeared. Finally, the small black points showing the satellites tend to lie above the mean relation, in qualitative agreement with data (Pasquali et al. 2012).  \nFigure 10 shows the mass-weighted stellar metallicity as a function of M GLYPH<3> at z = 0 . Points show centrals colourcoded by sSFR, with satellites in black. The yellow dashed line shows a running median. Observations are shown from Gallazzi et al. (2005) and Panter et al. (2008).  \nSimba nicely reproduces the stellar MZR for massive, quenched galaxies. At lower masses, the star-forming population dominates, and these tend to have a stellar MZR that is typically slightly higher, with larger scatter, than expected from an extrapolation of the massive galaxy relation. This owes to the fact that these galaxies have continued to form stars after their more massive counterparts have quenched. However, no such feature is evident in the observations, and hence Simba produces a low-mass stellar MZR that is somewhat too high compared to observations.  \nIn summary, Simba does a reasonable job reproducing observed galaxy metallicities, both stellar and gas phase, and the evolution out to Cosmic Noon. The fact that no arbitrary normalisation was required as in Mufasa is a step forward, suggesting that Simba is locking metals into dust in a realistic manner; we will explore this more in § 3.9.", '3.6 Galaxy photometric sizes': 'Modern cosmological simulations typically have sufficent resolution to resolve the size of galaxies, even if the detailed structure of the ISM remains unresolved. Illustris highlighted the ability for simulations to produce galaxies populating the full range of the Hubble sequence (Vogelsberger et al. 2014). For EAGLE, galaxy sizes provided a key constraint on their star formation feedback implementation (Schaye et al. 2015), namely that they employ a steeper dependence of the star formation rate on the density in dense  \nFigure 9. Gas-phase mass-metallicity relation at z = 0 (left) and z = 2 : 3 (right) from Simba . Points are colour-coded by deviation in sSFR from the star-forming main sequence. Running median values are shown as the green lines. At lowz , best-fits to observations from Tremonti et al. (2004); Andrews & Martini (2013); Yates et al. (2019) are shown as the black lines. At z GLYPH<25> 2 : 3 , observations are shown from the MOSDEF survey (Sanders et al. 2015). Simba reproduces the gas-phase MZR well at both redshifts, with a noticeable second-parameter dependence on SFR particularly at lowz .  \n<!-- image -->  \nFigure 11. R -band 2-D projected half-light radii of Simba galaxies at z = 0 , as a function of M GLYPH<3> . Points are colour-coded by sSFR. We fit median relations to the red and blue sub-samples, delineated by 10 GLYPH<0> 1 : 8 Gyr GLYPH<0> 1 , as the cyan and magenta dashed lines lines, respectively. Observational relations are shown from Zhang & Yang (2017) from SDSS, split into red and blue galaxies. Simba broadly reproduces the sizes of star-forming galaxies, but fails to show the observed trend that quiescent galaxies have a much steeper slope and are much more compact at low masses.  \n<!-- image -->  \nFigure 10. Stellar mass-stellar metallicity relation at z = 0 from Simba , with a running median shown as the dashed yellow line. Points are colour-coded by specific SFR. Observations are shown from Gallazzi et al. (2005); Panter et al. (2008). Simba reproduces the stellar metallicities of galaxies fairly well, although it appears that low-mass star-forming galaxies tend to have somewhat higher metallicities than typically observed.  \n<!-- image -->  \ngas in order to prevent galaxies from being overly compact. Simba did not use sizes to tune the feedback model, so instead they provide a test of it.  \nTo conduct a fair comparison to observed sizes, we compute projected galaxy halflight radii in the R -band. We obtain R -band luminosities from the Bruzual & Charlot (2003) models interpolated to the age and metallicity of each star particle. The radius is determined for each galaxy by averaging 2-D half-light projections along the x ; y ; z axes. Figure 11 shows the galaxy half-light sizes at z = 0 from Simba versus stellar mass. We colour-code the points by sSFR as before, and fit separate running medians for quenched  \n( Simba -Q) and star-forming ( Simba -SF) populations (magenta and cyan lines), where we divide the two populations at sSFR = 10 GLYPH<0> 1 : 8 Gyr GLYPH<0> 1 as before. For comparison we show observations from SDSS (Zhang & Yang 2017). The SDSS sample has been subdivided into red and blue galaxies, albeit with a criterion based on photometry, not sSFR.  \nStar-forming galaxy sizes in Simba show an amplitude and scaling with M GLYPH<3> that agrees quite well with the observed slope, which is encouraging. There is a suggestion that lowmass galaxies are too small, but this occurs in the mass  \nFigure 12. Gas fractions as a function of M 500 at z = 0 . Black points show the total baryon fraction within the halo, and red and blue points show the gas fractions subdivided at 10 5 K. Purple hexbins show an observational compilation of hot gas fractions from McCarthy et al. (2017), to be compared to the red points. Simba does a reasonable job of reproducing the observed trend of hot gas fraction with M 500 , which has been difficult for previous simulations to achieve without tuning.  \n<!-- image -->  \nrange where the number of particles is below a few hundred, and given that the sizes are light-weighted, stochasticity can give rise to smaller-than-expected sizes. We note that a stellar mass-weighted size does not show this drop-off at low masses. But for well-resolved star-forming galaxies, the sizes are in quite good agreement with data. This is an important success that did not require any specific tuning of the feedback model.  \nIn contrast to the star-forming systems, Simba shows quenched galaxy sizes that are quite discrepant with observations. Massive galaxy sizes are in reasonable agreement with data, but the lowest-mass quiescent galaxies are up to GLYPH<24> GLYPH<2> 3 larger than the comparable sample in SDSS, showing that the size-mass trend for passive galaxies is incorrect in Simba ; indeed, in Simba the low-mass passive galaxies are actually larger than the star-forming ones, which is opposite to the observed trend. There are a number of potential reasons for this. The large number of stellar orbits in older quiescent galaxies tends to puff out the distribution numerically. The discrepancy could also represent a failing in physics, if for instance low-mass galaxies are preferentially quenched via some rapid mode such as merging, violent disk instability, or stripping that is simultaneously associated with compactification (e.g. Tacchella et al. 2016). Alternatively, it could be a failing of the feedback physics associated with quenching low-mass galaxies. We can test this issue directly with higher resolution runs once they complete. For now, we note that the sizes of small quenched galaxies in Simba is a clear failing of the current model, in contrast to its success at reproducing the sizes of star-forming galaxies.', '3.7 Halo gas fractions': "In Simba , AGN feedback provides the primary energy input that serves to quench massive galaxies in large halos. Such energy input can concurrently have a strong impact on the amount and thermal state of hot gas within those halos.  \nIn particular, it can evacuate gas even from fairly sizeable halos, somewhat by entraining gas in jets but mostly by depositing heat that results in the gas becoming unbound to the halo. These processes result in halo gas fractions that deviate strongly from the mean cosmic baryon fraction, a departure that can be measured in real systems via X-ray emission from intra-group and intra-cluster gas. Such observations thus provide an important constraint on the AGN feedback model.  \nThe hot gas fraction as a function of halo mass has been a challenging constraint for modern cosmological AGN feedback models to reproduce (McCarthy et al. 2017). In Illustris, it was found that the AGN feedback mechanism overevacuated hot gas from group-sized halos compared to observations (Genel et al. 2014), which provided one motivation for the new AGN feedback model in Illustris-TNG (Weinberger et al. 2018). Nonetheless, while closer than Illustris, TNG somewhat overpredicts the observed hot gas fractions (Barnes et al. 2018). EAGLE likewise overpredicts the hot gas fractions, while the Bahamas simulation suite was able to match this with mild tuning (McCarthy et al. 2017). Here we examine this constraint for Simba .  \nFigure 12 shows the baryon fractions as a function of halo mass ( M 500 ) from Simba . M 500 is computed as the radius enclosing 500 times the critical density, centered on the most bound halo particle. Black points show the total baryon fraction, red points show hot gas ( T > 10 5 K) fractions, and blue show cold gas ( T < 10 5 K). Colour-coordinated lines show the running median values. Note that the black points include the stellar (and black hole) contribution, which is not explicitly shown. A compilation of observations from McCarthy et al. (2017) is shown as the purple hexbins. All fractions have been scaled to GLYPH<10> b = 0 : 048 , so a halo at unity has its cosmic share of baryons.  \nAt 10 12 M GLYPH<12> halos have about 40% of the cosmic baryon fraction within R 500 , with a large scatter. This dips to GLYPH<24> 30% at 10 12 : 5 GLYPH<0> 13 M GLYPH<12> , before rising again to large halo masses. The dip owes to jet AGN feedback, which we have checked by comparing to the No-jet test run, which shows baryon fractions around 90% for all halos over this mass range (and a much flatter trend of hot gas fraction versus halo mass). This shows that the energy input required to quench galaxies can cause substantial evacuation of group-sized halos, as has been noted in e.g. Illustris (Genel et al. 2014). This strong evacuation has important implications for using these systems as probes of cosmology, which we will probe in future work.  \nThe hot baryon fraction is shown in red, which can be compared to the observations shown in purple. In massive systems, the total baryons are dominated by this hot phase. Most of this hot gas is near the virial temperature, so the results are insensitive to the exact value of the cut at 10 5 K. Comparing to observations, we see that Simba 's halos have hot baryon fractions that are well in agreement with data in the overlapping mass range, in both amplitude and scatter. We note that there was no tuning done to obtain this agreement. The halo hot baryon fraction is thus a non-trivial success of Simba 's AGN feedback model, and shows that Simba evacuates halo baryons in a manner that is concordant with observations. In Borrow et al. (in prep.) we will examine quantitatively where these evacuated baryons end up.  \nFigure 13. M BH GLYPH<0> M GLYPH<3> relation at z = 0 in Simba . Points show galaxies, with centrals colour-coded by specific SFR, and satellites as the grey points. Observations are shown from Kormendy & Ho (2013) for comparison with bulge-dominated (redder) systems, while Bentz & Manne-Nicholas (2018) shows the relationship more appropriate for spiral star-forming systems at lower M GLYPH<3> . Simba broadly reproduces these observed relations in its appropriate galaxy populations.  \n<!-- image -->", '3.8 Black hole mass vs stellar mass': "The canonical relation that highlights the connection between galaxies and their central supermassive black holes is the relationship between the black hole mass and the galaxy bulge mass or stellar velocity dispersion (e.g. Magorrian et al. 1998; Kormendy & Ho 2013; McConnell & Ma 2013; Graham 2016; Bentz & Manne-Nicholas 2018). Modern galaxy formation models that track black holes typically have free parameter(s) that are tuned to match these relations; in Simba , this is set by the accretion efficiency tuned to GLYPH<15> m = 10 %, while most other cosmological simulations (based on Bondi accretion) tune the AGN feedback efficiency. In previous works, Angl'es-Alc'azar et al. (2015) and Angl'es-Alc'azar et al. (2017a) showed that the M BH GLYPH<0> M GLYPH<3> relation emerged naturally from the torque-limited accretion model, without or with AGN feedback, respectively. But these studies were done via post-processing or without star formation feedback. Here we examine whether the full physics model in Simba likewise reproduces the relationship between black hole mass and galaxy properties.  \nFigure 13 shows the black hole mass-stellar mass relation at z = 0 for Simba galaxies. Central galaxies are shown colour-coded by specific SFR, while satellite galaxies are indicated by grey points. The relationship for galaxy bulges is shown from Kormendy & Ho (2013) as the magenta dashed line; this is an appropriate comparison sample for bulge-dominated galaxies, which are expected to be the quiescent systems with redder points. Meanwhile, Bentz & Manne-Nicholas (2018) assembled a sample of reverberationmapped galaxies, and found the steeper relation shown as the blue dotted line. In their case, the lower-mass systems are predominantly late-type galaxies, while the most massive systems are early-type. Hence in the region plotted, the Bentz & Manne-Nicholas (2018) sample is probably best compared to later-type systems, and therefore star form-  \ning (bluer points). We note that all observational relations show a large scatter, typically at least 0.3 dex, which is not represented on this plot.  \nSimba black holes generally lie in the range of observations. Although we tuned GLYPH<15> m = 0 : 1 in order to obtain the correct amplitude of the relation, the slope of the relation is not tunable, particularly in terms of different galaxy subsamples. Hence the agreement of the quiescent galaxy slope with the bulge-dominated galaxy black holes, and likewise the general agreement of the star-forming galaxies with the lower black hole masses at a given M GLYPH<3> , is in good agreement with observations. We note that there is some disagreement on whether the late-type galaxies have a steeper slope or the same slope but offset to lower black hole masses (e.g. Graham 2016; Savorgnan et al. 2016; Bentz & Manne-Nicholas 2018), but Simba predictions are broadly compatible with either scenario.  \nAt the low-mass end, consistent with Angl'es-Alc'azar et al. (2015), torque-limited accretion grows black holes very quickly once the galaxy stellar mass exceeds 3 GLYPH<2> 10 9 M GLYPH<12> , which is where we choose to seed the black holes in Simba . Hence the rapid rise is not directly physical but a numerical artifact of our seeding prescription, though it is intended to mimic the physical effect of black hole growth suppression due to early star formation seen in e.g. FIRE (Angl'es-Alc'azar et al. 2017c). Also, we note that we attempt to keep black holes fixed to the centre of the galaxy potential well, but in dense regions this does not always work owing to the shallow potential wells in poorly resolved galaxies, so black holes can move between galaxies and thus merge. We continue to test for approaches for better handling this given the poor cosmological resolution.  \nOverall, Simba predicts a relationship between black hole and stellar masses in agreement with observations. In upcoming work we will examine black hole scaling relations in more detail, but for now, the good agreement corroborates the idea that black holes in Simba grow in accord with observations, and thus the feedback energy released by black holes and used by Simba to quench galaxies is plausible.", '3.9 Dust Properties': "Simba includes a model to form and destroy dust from metals within the ISM of galaxies during the simulation run. As a basic check on the production of dust, here we examine two measurables tracking dust in galaxies: The dust mass function, and the dust-to-gas ratio.  \nFigure 14 shows the z = 0 (bottom panel) and z = 2 (top) dust mass function (DMF) from Simba (green line), versus two z = 0 observational determinations from Dunne et al. (2011) and Clemens et al. (2013), and a z = 2 determination from Dunne et al. (2003). At z = 0 , Simba agrees well with Dunne et al. (2011), but not Clemens et al. (2013). The difference between the two can be traced to their assumption of the dust mass opacity coefficient used to infer the dust mass from far-IR data; Clemens et al. (2013) showed that under the same assumption of this quantity, the two results agree. Hence given current uncertainties in inferring dust masses, it is probably premature to use the DMF as a strong constraint on models. But Simba 's DMF is at least within the ballpark of currently observed values, with good agreement in the overall DMF shape. Unsurprisingly, the  \nFigure 14. Dust mass function from Simba at z = 0 , shown as the green shaded region. DMF is split into star-forming and quenched samples, shown as blue and red dashed lines, respectively. Observations are shown from Dunne et al. (2011) and Clemens et al. (2013); differences owe to assumptions regarding inferring dust mass from far-IR flux. Simba reproduces the observed shape of the DMF, and is in the range of observed amplitudes, modulo systematic uncertainties.  \n<!-- image -->  \nDMF is dominated by star-forming galaxies (blue dashed line), which is as observed (Beeston et al. 2018).  \nThe z = 2 DMFis compared to observations from Dunne et al. (2003), and shows a deficit of GLYPH<24> GLYPH<2> 3 in the number density of galaxies at a given dust mass. We note that the observational DMF by Dunne et al. (2003) is from surveys of sub-mm sources with large beam sizes, which could result in multiple objects being blended within one beam therefore overestimating their dust masses. If one regards the Clemens et al. (2013) results at z = 0 to be more accurate, as confirmed by Beeston et al. (2018), then the shortfall in the predicted DMF is very similar at both redshifts. This suggests that the evolution in dust masses in Simba is viable, but the overall dust production is short, or else destruction is too efficient. It may be possible to remedy this with differing choices of dust parameters; we are exploring this. We note that our z = 0 DMF agrees well with the predictions from cosmological simulations of McKinnon et al. (2017),  \nFigure 15. Dust to gas ratio as a function of gas-phase metallicity at z = 0 in Simba galaxies, colour-coded by specific SFR. Observations are shown as crosses from R'emy-Ruyer et al. (2014). Simba reproduces the observed trend and amplitude in dust-togas ratios.  \n<!-- image -->  \nFigure 16. Metal mass fraction locked in dust, as a function of M GLYPH<3> , for star-forming galaxies in Simba . Plot is colour-coded by the mean sSFR within each hexbin. Except for the smallest galaxies, typically one-third of the metals are locked in dust.  \n<!-- image -->  \nowing in part to tuning of each model, but our z = 2 DMF is significantly higher than theirs.  \nFigure 15 shows the z = 0 dust-to-gas ratio (DGR) as a function of gas-phase metal abundance. Points are colourcoded by sSFR. Simba is in good agreement with the data shown (R'emy-Ruyer et al. 2014, crosses), showing a slope of increasing DGR with metallicity as observed. In massive quenched systems, the DGR drops quickly.  \nFinally, Figure 16 shows the fraction of metals locked into dust at z = 0 for star-forming galaxies. The green line shows the running median. For galaxies with M GLYPH<3> > GLYPH<24> 10 9 : 5 M GLYPH<12> , the fraction is typically 30-40%, but drops significantly to lower masses. This mostly explains why the mass-metallicity relation agrees with observations in Simba without the ad hoc reduction of the yields by GLYPH<2> 2 as in Mufasa . Low-sSFR galaxies also have fewer metals locked into dust, as AGN feedback returns metals locked in dust into the gas phase.  \nFigure 17. Stellar mass function evolution from z = 2 ! 0 in 50 h GLYPH<0> 1 Mpc , 512 3 test runs with different AGN feedback variants: Original Simba (green solid), Jet and X-ray feedback both turned off (No-jet; blue dashed), only X-ray feedback turned off (NoXray; red dashed). For comparison we show the main 100 h GLYPH<0> 1 Mpc Simba run (green dashed) reproduced from Figure 4, as well as selected observations as indicated. Turning on just the jet feedback (No-Xray) results in a substantial truncation of the GSMF that does not occur without jets (No-jet).  \n<!-- image -->  \nWe will examine galaxy dust content and evolution in significantly more detail in forthcoming work (Li et al., in prep.), but these preliminary comparisons suggest that Simba 's dust tracking model yields plausible galaxy dust contents.", '4 AGN FEEDBACK VARIATIONS': "AGN feedback is believed to be responsible for quenching galaxies. Simba includes three different forms of AGN feedback: Radiative mode AGN winds, AGN jets, and X-ray feedback. In this section we examine the importance of these various modules in producing a quenched galaxy population, by running simulations with AGN jets and X-ray feedback off, and with only X-ray feedback off. We always include radiative AGN winds. For these tests we use 50 h GLYPH<0> 1 Mpc , 2 GLYPH<2> 512 3 simulations run to z = 0 , with Simba input physics and parameters except for the AGN feedback variations.  \nLooking at the z = 0 panel, without jets, the GSMF is strongly overproduced at high masses. Turning on jets  \nFigure 17 shows the GSMF for the full Simba physics run, a No-X run turning off only X-ray feedback, and a Nojet run turning off both jet and X-ray feedback, at z = 2 ; 1 ; 0 . Observations are overplotted as described in Figure 4. We do not show z GLYPH<21> 3 results because these variants' GSMFs are indistinguishable there.  \n<!-- image -->  \nFigure 18. M BH GLYPH<0> M GLYPH<3> relations in test simulations with jet and X-ray black hole feedback turned off (No-jet, top panel), and jets on but X-ray feedback turned off (No-Xray, bottom). By comparing to Figure 13, the jet feedback is seen to enact most of the quenching, but the X-ray feedback is important for fully quenching the most massive galaxies.  \n<!-- image -->  \nresults in much better agreement with full Simba . This is even true when including X-ray feedback, which makes only a small change to the GSMF. The redshift evolution shows that the impact of jets is fairly minor at z GLYPH<24> 2 in terms of the GSMF, only impacting the very largest few galaxies. The importance of jets in truncating the GSMF grows steadily with time, to where without jets the number density of M GLYPH<3> = 10 12 M GLYPH<12> galaxies would be an order of magnitude higher, in strong disagreement with data. These results clearly show that the main driver in truncating the massive end of the GSMF in Simba is AGN jet feedback.  \nA more detailed view of how AGN feedback impacts both stellar and black hole growth can be obtained by examining the M BH GLYPH<0> M GLYPH<3> relation in these variants, shown in Figure 18, with galaxies colour-coded by specific SFR as in Figure 13. For clarity we show only central galaxies.  \nComparing the No-Jet version (top panel) to the original Simba in Figure 13 highlights several key points. As expected, the subM ? objects show little difference in the trends, in either M BH or sSFR. However, for massive galaxies, the full Simba run shows significantly lower sSFR and somewhat higher M BH , particularly for the most massive  \ngalaxies. This demonstrates more explicitly that the AGN jet feedback is crucial for quenching galaxies.  \nInterestingly, the No-jet run still shows a few quenched galaxies at high M BH around M ? . These are clearly correlated with the presence of a massive black hole, and would not occur in a model with no AGN feedback at all. This arises from the fact that we still have radiative AGN winds in our No-jet run. These become effective around M ? because it is at the corresponding halo mass that a significant hot gaseous halo begins to form (Kereˇs et al. 2005; Gabor & Dav'e 2012). The AGN energy can then be deposited into the hot gas, providing a mechanism for quenching the galaxy by shutting off the fuel supply (Dekel et al. 2009).  \nSo why do radiative winds cease to be effective at higher masses? An examination of the energetics shows the reason. In Simba 's AGN feedback model the momentum input is assumed to be constant, which means that the AGN feedback energy scales as the wind velocity. Since Simba 's black hole accretion rates are a quite weak function of M GLYPH<3> (Thomas et al., in prep) while the number of hot halo baryons is growing, this means that the energy injected per hot halo baryon is dropping with the halo mass. The logarithmic increase in velocity with M BH (eq. 7) is too slow to compensate for this, so one quickly ends up in a situation where the energy injection is insufficient to keep the hot halo baryons near the virial temperature. What is required is a strong increase in the outflow velocity, and hence energy input, in this halo mass regime. This is why high-velocity AGN jets are crucial for quenching massive galaxies.  \nThe black hole masses in the No-jet run are also appear to be significantly lower. However, this can primarily be explained by the effect that M GLYPH<3> values in this simulation are substantially higher, which moves galaxies leftwards in the M BH GLYPH<0> M GLYPH<3> diagram; the black hole masses themselves are not substantially different. The relative roles of Bondi vs. torque-limited accretion in growing black holes across the full mass range over cosmic time will be examined more fully in a forthcoming paper (Ang'es-Alc'azar et al. in prep.).  \nThe No-Xray run is fairly similar to the full Simba run, but there is a slight if noticeable increase in the sSFR in the massive galaxies. This is not so much as to contribute significant mass growth, hence the GSMF is only modestly affected, but it is higher than typical observed values for massive red and dead ellipticals. This suggests that the Xray feedback is important for fully quenching massive galaxies in accord with observations, even if it does not play a leading role in regulating mass growth.", '5 SUMMARY': "We have introduced the new Simba suite of cosmological galaxy formation simulations, and explored predictions from a 100 Mpc/h box run with 1024 3 dark matter and 1024 3 gas elements. The most novel aspect of Simba is its implementation of black hole growth via torque-limited accretion, and two-mode black hole feedback via bipolar kinetic outflows. Simba further includes numerous updates over its predecessor simulation Mufasa , including a dust production and destruction model. In this paper we present comparisons to a range of different observational probes measuring the stellar mass, star formation rate, neutral and molecular gas, black  \nhole, and dust properties in Simba . We show that, in all cases, Simba produces galaxies that are in quite reasonable agreement with observations. While our feedback parameters were generally chosen to follow observations or expectation from high-resolution simulations, some of these observations were used to further tune these parameters. However, many of them were not, and these represent model predictions that demonstrate the viability of Simba as a platform for studying cosmological-scale galaxy evolution.  \nHere are our main findings:  \n- GLYPH<15> Simba produces a stellar mass function evolution that is in very good agreement with data across all masses at all cosmic epochs, although it may overproduce the massive end slightly by z = 0 . Quenched galaxies grow substantial in numbers at z < GLYPH<24> 2 , and by z = 0 they dominate at M GLYPH<3> > 2 GLYPH<2> 10 10 M GLYPH<12> .\n- GLYPH<15> Simba 's star forming main sequence is in good agreement with observations at z = 0 , and low at z GLYPH<25> 2 by only a factor of two which is explainable via observational systematics. Predicted quenched fractions at z = 0 as a function of M GLYPH<3> are in good agreement with observations. Galaxies at a given M GLYPH<3> that quench first in Simba have preferentially larger black holes.\n- GLYPH<15> Simba gas fractions, both neutral and molecular, show a dropping trend with M GLYPH<3> that is in good agreement with observations. Gas fractions evolve downwards with time, but even at z = 0 , massive quenched galaxies still typically have some cold gas.\n- GLYPH<15> The gas-phase and stellar mass-metallicity relations generally agree with observations at z = 0 and z GLYPH<24> 2 . The MZR evolves upwards by a factor of GLYPH<24> GLYPH<2> 3 for our smallest systems, but at M GLYPH<3> > GLYPH<24> 10 11 M GLYPH<12> there is little evolution.\n- GLYPH<15> Galaxy photometric projected sizes are in good agreement with observations for star-forming systems, but are too large for quenched systems particularly at lower masses.\n- GLYPH<15> There is substantial evacuation of baryons from halos at group scales, with Local Group-sized objects retaining typically only a third of their cosmic baryon fraction. Hot halo gas fractions show a rising trend with halo mass in good agreement with data.\n- GLYPH<15> The black hole mass-stellar mass relation shows a population of quenched galaxies that agrees well with observations of bulge-dominated systems, while star-forming galaxies at a given M GLYPH<3> have lower black hole masses.\n- GLYPH<15> Simba predicts a z = 0 dust mass function and dust-togas ratios in agreement with observations, with more massive star-forming galaxies having higher ratios. At a given metallicity, galaxies that have higher SFR have higher dustto-gas ratios. Roughly one-third of metals are locked in dust.  \nThese results demonstrate that Simba broadly reproduces a wide range of observed galaxy properties including their stellar content, star formation rates, gas content, metallicities, sizes, black hole properties, and dust properties. Clearly there are many other more detailed constraints that would be tested, and in subsequent work we aim to examine these in more detail. It is important to note that Simba also displays some aspects that are in conflict with observations. It fails to produce as sharp a knee in the z = 0 stellar mass function as is observed. It produces low-mass quenched galaxy sizes that are larger than for star-forming systems, opposite to the observed trend in SDSS. There are  \nsuggestions that Simba overproduces the stellar metallicities as well as the specific SFRs in low-mass present-day star-forming systems. Finally, in many lowz scaling relations (e.g. gas and metallicities vs. M GLYPH<3> ) there is an abrupt break in the typical properties of galaxies above and below M GLYPH<3> GLYPH<25> 2 GLYPH<2> 10 10 M GLYPH<12> , which qualitatively agrees observations but quantitatively may be too sharp. This rapid transition contrasts with the too-gradual turn-down in the GSMF, suggesting a tension in how Simba (and similar cosmological models) quench massive galaxies at low redshifts. Despite these minor issues, Simba provides a state of the art platform for studying galaxy evolution across a range of masses, environments, and cosmic epochs, and promises to yield numerous new insights into the physical drivers of galaxy evolution over cosmic time.", 'ACKNOWLEDGEMENTS': 'The authors acknowledge helpful discussions with Josh Borrow, Weiguang Cui, Shuiyao Huang, Katarina Kraljic, and Neal Katz. We thank Philip Hopkins for making Gizmo public, and providing our group with early access. We thank Robert Thompson for developing Caesar , and the yt team for development and support of yt . RD acknowledges support from the Wolfson Research Merit Award program of the U.K. Royal Society. DAA acknowledges support by a Flatiron Fellowship. The Flatiron Institute is supported by the Simons Foundation. DN and QL were supported in part by NSF Award AST-1715206 and HST Theory Award 15043.0001. This work used the DiRAC@Durham facility managed by the Institute for Computational Cosmology on behalf of the STFC DiRAC HPC Facility. The equipment was funded by BEIS capital funding via STFC capital grants ST/P002293/1, ST/R002371/1 and ST/S002502/1, Durham University and STFC operations grant ST/R000832/1. DiRAC is part of the National eInfrastructure.', 'REFERENCES': 'Yates R. M., Schady P., Chen T.-W., Schweyer T., Wiseman P., 2019, arXiv e-prints, Young L. M., et al., 2014, MNRAS, 444, 3408 Yu Q., Tremaine S., 2002, MNRAS, 335, 965 Zahid H. J., et al., 2014, ApJ, 792, 75 Zhang Y., Yang X., 2017, arXiv e-prints,  \nThis paper has been typeset from a T E X/L A T E X file prepared by the author.'}