https://orcid.org/0000-0003-2842-9434
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és-Alcázar et al. 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 |${\sim } 10^3\,\,{\rm km}\, {\rm s}^{-1}$| at high Eddington ratios, increasing to |${\sim } 8000\,\,{\rm km}\, {\rm s}^{-1}$| at Eddington ratios below 2 per cent, with a constant momentum input of 20L/c. simba further includes an on-the-fly dust production, growth, and destruction model. Our simba run with |$(100h^{-1}\, {\rm Mpc})^3$| and 10243 gas elements reproduces numerous observables, including galaxy stellar mass functions at z = 0−6, the stellar mass–star formation rate main sequence, H i and H2 fractions, the mass–metallicity relation at z ≈ 0, 2, star-forming 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 too-large sizes for low-mass quenched galaxies. We show that simba’s jet feedback is primarily responsible for quenching massive galaxies.
1 INTRODUCTION
The formation and evolution of galaxies are governed by a wide range of physical processes spanning from sub-parsec to giga-parsec scales. These include the growth of large-scale 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 ever-advancing 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 regard to feedback processes that are currently among its least well-understood aspects.
Modern galaxy formation models are generally built on the premise that feedback at the very lowest masses is dominated by photoionization 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 discs around supermassive central black holes (see Somerville & Davé 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 with 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.
Cosmological-scale simulations that model galaxy growth and feedback dynamically within evolving large-scale structure are an increasingly valuable tool for testing and constraining galaxy formation physics, owing both to rapidly advancing computational power and to 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 (Genel et al. 2014; Vogelsberger 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.
The mufasa simulation project (Davé, Thompson & Hopkins 2016) has added to the pantheon of such simulations, employing several novel approaches that distinguish it from others. First, it utilizes 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 parametrizations 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.
An aspect where mufasa was physically less well motivated than other state-of-the-art simulations was in its depiction of black hole growth and active galactic nucleus (AGN) feedback. mufasa did not directly grow black holes and utilize the accretion energy for feedback to quench galaxies. Instead, following Gabor & Davé (2015), mufasa implemented a heuristic ‘quenching feedback’ in which diffuse gas within massive haloes 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-fitting analytic equilibrium model of Mitra, Davé & Finlator (2015) and Mitra et al. (2017), and evolved slowly upwards with redshift. Despite its simplicity, this prescription displayed impressive successes at reproducing a red sequence and massive central galaxy properties in excellent accord with observations (Davé, Rafieferantsoa & Thompson 2017b), albeit with some non-trivial discrepancies such as overquenching satellites particularly in the outskirts of large haloes (Rafieferantsoa & Davé 2018). This demonstrated that a model based primarily on starvation of the central galaxy via ‘radio mode’ feedback (Bower et al. 2006; Croton 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.
This 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 formation-driven 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és-Alcázar, Özel & Davé 2013; Anglés-Alcázar et al. 2015, 2017a). In this model, black hole growth is regulated by the ability for gas in the inner disc to lose angular momentum via disc 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 on to the black hole accretion disc with properties of the inner galactic disc. They showed that even at ∼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.
Anglés-Alcázar et al. (2013) explored the torque-limited accretion model via the post-processing of zoom simulations, and Anglés-Alcázar 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 that represents the fraction of material entering the accretion disc that accretes on to the black hole; a plausible choice of ∼10 per cent provided a good match to data, insensitive to the choice of the (uncertain) black hole seed mass. Anglés-Alcázar et al. (2017a) extended these previous works to self-consistently incorporate torque-limited 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és-Alcázar 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 use 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.
The second part of simba’s new black hole model involves a novel sub-grid prescription for AGNs feedback. AGN feedback connects flows coming off the black hole accretion disc to energy release on scales of tens or hundreds of kpc. To model this transfer of energy from small to large scales, simba utilizes 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 ionized outflows with velocities of |${\sim } 1000\,\,{\rm km}\, {\rm s}^{-1}$| or more (Sturm et al. 2011; Greene, Zakamska & Smith 2012; Maiolino et al. 2012; Liu et al. 2013; Perna et al. 2017a), and jets at velocities up to |${\sim } 10^4\,\,{\rm km}\, {\rm s}^{-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_{\rm Edd}\lesssim$|few per cent), while lower-velocity outflows typically arise in systems with higher fEdd (Best & Heckman 2012; Heckman & Best 2014). Extreme systems such as bright quasars often show both types of outflows. simba’s black hole outflows are parametrized 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−1 in low-mass, fast-accreting black holes, up to many thousands of |$\,\,{\rm km}\, {\rm s}^{-1}$| for slower-accreting black holes.
A key unique feature is that simba’s kinetic feedback is purely bipolar, with the ejection vector given by the angular momentum of the inner disc. This direction is relatively stable over galaxy dynamical time-scales. To be conservative in minimizing 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 randomize the kinetic feedback’s direction on short time-scales (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).
The 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 self-regulation 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 sphericalizes at large distances, sufficient in massive haloes 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 feedback that is more closely connected to observations of AGN winds and jets.
In addition to kinetic AGN feedback, simba also includes X-ray feedback input into surrounding gas. The importance of this feedback channel has been emphasized 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 cosmological-volume simulation to include such X-ray feedback.
Another 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 2 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.
In this paper, we describe the simulation methodology in simba, which, besides the new black hole model, also makes various other minor improvements to mufasa (Section 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é et al. 2016, 2017a; Davé et al. 2017b; Rafieferantsoa & Davé 2018), paying particular attention to black hole and massive galaxy properties that represent the most direct test of simba’s new AGN feedback model (Section 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 (Section 4), showing that the high-velocity jet feedback is crucial for producing a quenched massive galaxy population. Finally, we summarize our results in Section 5.
2 SIMULATION METHODOLOGY
2.1 Code and input physics
The simba simulations utilize much of the framework of the mufasa simulations as described in Davé 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.
simba utilizes a forked version of the gizmo cosmological gravity plus hydrodynamics solver (Hopkins 2015, 2017), in its MFM version. This code, based on 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 discs, 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.
Radiative cooling and photoionization 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 grackle-2.1 used in mufasa, offering two advantages: First, the adiabatic and radiative terms are evolved together during the cooling sub-time-step, unlike the previous operator-split approach where first the system was evolved adiabatically over the full time-step 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 self-shielding self-consistently during the simulation run, based on the Rahmati et al. (2013) prescription in which the ionizing background strength is attenuated depending (primarily) on gas density. A spatially uniform ionizing background is assumed as specified by Haardt & Madau (2012), modified to account for self-shielding (A. Emerick, private communication). These changes do not make a noticeable difference to the resulting galaxy population, but they do improve the accuracy of the baryonic thermal evolution that 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é et al. 2017a).
As in mufasa, we use an H2-based star formation rate, where the H2 fraction is computed based on the sub-grid model of Krumholz & Gnedin (2011) based on the metallicity and local column density, with minor modifications as described in Davé et al. (2016) to account for variations in numerical resolution. The star formation rate is given by the H2 density divided by the dynamical time: SFR=ϵ*ρH2/tdyn, where we use ϵ* = 0.02 (Kennicutt 1998). The chemical enrichment model tracks 11 elements (H, He, C, N, O, Ne, Mg, Si, S, Ca, and 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é (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 Section 2.5. Type Ia SNe and AGB wind heating are also included as in mufasa, along with ISM pressurization at a minimum level as required to resolve the Jeans mass in star-forming gas as described in Davé et al. (2016).
Mass loading factor η versus stellar mass M* from a suite of FIRE simulations analysed via particle tracking in Anglés-Alcázar et al. (2017b). The points show values measured at various redshifts, while the orange line is the best-fitting relation. The grey solid line shows the Muratov et al. (2015) scaling and the dashed and dotted grey lines show the mass loading factors used in the Illustris and Illustris-TNG simulations (η(Mhalo) fitting functions from Pillepich et al. 2018, converted to η(M*) using the M*–Mhalo relation of Moster, Naab & White 2013).
simba uses an on-the-fly approximate friends-of-friends (FOF) finder applied to stars and dense gas as described in Davé et al. (2016) in order to compute the galaxy properties such as M*, and as in mufasa obtains vcirc using a scaling based on the baryonic Tully–Fisher relation. Besides algorithmic improvements to improve parallel performance, the only 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 two-mode accretion model. The first mode closely follows the torque-limited accretion model presented in Anglés-Alcázar 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és-Alcázar et al., in preparation) 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
The normalization factor ϵT ≡ ϵm × α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 disc), where αT = 5 is the original normalization of |$\dot{M}_{\rm Torque}$| in Hopkins & Quataert (2011) and we set ϵm = 0.1 to match the normalization of the local MBH–M⋆ relation as in Anglés-Alcázar et al. (2017a). One can view αT as corresponding to an efficiency of transport of material from the inner galactic disc on to the black hole accretion disc, and ϵm as the efficiency of transport from the accretion disc on to the black hole, for which 10 per cent is a canonical value. However, αT is itself fairly uncertain, and in the end the meaningful subgrid parameter is only the combination ϵT.
2.2.2 Bondi accretion from hot gas
Hot gas can also accrete on to 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 > 105 K.
2.2.3 Numerical implementation
Numerically, we follow Springel, Di Matteo & Hernquist (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 fm of the mass of gas particles within R0 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 per cent 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 for 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_{\rm Edd}\equiv \dot{M}_{\rm BH}/\dot{M}_{\rm Edd} \gtrsim$| few per cent), in which AGNs are seen to drive multiphase winds at velocities of |${\sim } 1000 \,\,{\rm km}\, {\rm s}^{-1}$| that include molecular and warm ionized gas (e.g. Sturm et al. 2011; Maiolino et al. 2012; Perna et al. 2017a); and a ‘jet mode’ at low Eddington ratios, where AGNs mostly drive hot gas in collimated jets at high velocities (of order |${\sim } 10^4\,\,{\rm km}\, {\rm s}^{-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’ versus ‘low excitation’ radio galaxies) above and below roughly a per cent 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 large-scale 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
Based on observations of AGN outflows and the inferred momentum and energy input (e.g. Fiore et al. 2017; Ishibashi, Fabian & Maiolino 2018), we set the amount of material ejected in the AGN winds in order to obtain a momentum input of |$\dot{P}_{\rm out} = 20\, L/c$|, where |$L = \eta \, \dot{M}_{\rm BH} \, c^2$| is the bolometric luminosity of the AGN, η = 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_{\rm BH} = 10^9\, \, \mathrm{M}_\odot$| in the high-fEdd mode injects outflows with |$v_{\rm w,EL} = 1000\, \,\,{\rm km}\, {\rm s}^{-1}$|, mass loading |$\dot{M}_{\rm out,EL} / \dot{M}_{\rm BH} \approx 600$|, and kinetic energy efficiency |$\dot{E}_{\rm kin,EL} \approx 0.03\, L$|, while in the low-fEdd mode at full jet speed reaches |$v_{\rm w,jet}=8000\,\,{\rm km}\, {\rm s}^{-1}$|, |$\dot{M}_{\rm out,jet} / \dot{M}_{\rm BH} \approx 75$|, and |$\dot{E}_{\rm kin,jet} \approx 0.3\, L$|.
All outflows are ejected in a purely bipolar fashion. That is, we eject gas elements in a direction ±parallel to the angular momentum vector of the inner disc 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 disc 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 (|$\lesssim1$| kpc). We note that the outflow direction can precess owing to variations in the inner disc, but is in practice typically stable over tens to hundreds of Myr. Hence any effect of ‘sphericalizing’ 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.
Since 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_{\rm vir}=9.52\times 10^7 (M_{\rm halo}/10^{15}\, \mathrm{M}_\odot)^{1/3}$| K (Voit 2005). This choice is motivated by observations showing that jets contain mostly synchrotron-emitting plasma, and eventually thermalize their energy into surrounding hot gas at around Tvir (Fabian 2012). The extra energy input required for this is typically less than a few per cent of the jet kinetic energy, so it does not figure significantly into the overall energy budget.
We apply a short hydrodynamic and radiative cooling decoupling time of 10−4tH to the outflowing wind gas elements, where tH 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 for from the assumption of constant momentum input of 20L/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) time-step limiter in order to ensure proper interactions of the high-speed winds and their surrounding gas as they recouple.
Our model has similarities to the two-mode thermal and kinetic AGN feedback model employed in Illustris-TNG (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 fEdd rather than kinetic feedback. This may owe to the fact that torque-limited accretion does not require self-regulation, while Bondi accretion does (Anglés-Alcázar 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 low-fEdd, Illustris-TNG randomizes the direction of the jets at each time-step, rather than always ejecting jets perpendicular to the inner disc. Our approach seems more physically motivated, since jets are not known to dramatically precess on time-scales of ∼Myr (though they may occasionally do so over hundreds of Myr). (iv) Illustris-TNG uses, at maximum, 200 per cent of the AGN bolometric luminosity (assuming a 10 per cent 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 disc, 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 fgas < 0.2, where fgas = Mgas/M* as computed by our galaxy finder, and we only include X-ray heating in galaxies with full-velocity jets.
The 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, pressurized ISM as we assume in simba would cool quickly and not be physically well motivated, 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 neighbours the X-ray energy input exceeds this value, then no further X-ray heating is done for that black hole at that time-step. 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 Section 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és-Alcázar et al. 2017a). If a galaxy reaches a stellar mass M* > γBH × Mseed 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 Mseed = 104 M⊙/h and γBH = 3 × 105, which places black holes in galaxies with M* ≳ 109.5 M⊙.
This 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és-Alcázar 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, Volonteri & Dubois 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 |$\dot{M}_{\rm Bondi} \propto M_{\rm BH}^{2}$|). 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 Mseed and stellar mass threshold (Anglés-Alcázar et al. 2015).
We 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 × R0, where R0 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 ‘well-behaved’ 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 R0 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
s imba includes a dust physics module to track the lifecycle 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 μm. We ignore the active dust cooling, which will be applied in future work.
Dust 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_{i,d}^j$| refers to the dust mass of the ith element (C, O, Mg, Si, S, Ca, Fe) produced by the jth stellar process (Type II SNe or AGB stars), whereas |$m_{i,{\rm ej}}^j$| refers to the mass of ejecta from the jth process.
We take a fixed dust condensation efficiency |$\delta ^{\rm AGB}_{i,\rm 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 |$\delta ^{\rm SN II}_{i,\rm dust}=0.15$| to match the low-metallicity end of the observed relation between dust-to-gas mass ratios (DGRs) and gas-phase metallicities (Rémy-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, Torrey & Vogelsberger (2016) and Popping, Somerville & Galametz (2017), where Type Ia SNe are assumed to have the same condensation efficiency as Type II SNe.
We additionally destroy dust, as well as molecular hydrogen, completely in hot winds and during star formation (Section 2.1) and in any gas that is impacted by AGN X-ray heating or jets (Section 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, Ménard & Corrales 2015).
2.6 Runs and analysis
The primary simba runs have 10243 dark matter particles and 10243 gas elements. We are running four volumes: |$100h^{-1}\, {\rm Mpc}$| down to z = 0, |$50h^{-1}\, {\rm Mpc}$| to z = 1, |$25h^{-1}\, {\rm Mpc}$| to z = 2, and |$12.5h^{-1}\, {\rm Mpc}$| to z = 5. All runs have identical input physics, begin at z = 249, and assume a Planck Collaboration XIII (2016) concordant cosmology of Ωm = 0.3, |$\Omega _\Lambda =0.7$|, Ωb = 0.048, |$H_0=68\,\,{\rm km}\, {\rm s}^{-1}\,\,{\rm Mpc}^{-1}$|, σ8 = 0.82, and ns = 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 |$100h^{-1}\, {\rm Mpc}$| run, as the other runs are at various stages of completion.
The simba simulation suite.
| Name | |$L_{\rm box}^a$| | |$\epsilon _{\rm min}^b$| | |$z_{\rm end}^c$| | |$m^d_{\rm gas}$| | |$m^e_{\rm DM}$| | |$M^f_{\rm *,min}$| |
|---|---|---|---|---|---|---|
| m100n1024 | 100 | 0.5 | 0 | 1.82 × 107 | 9.6 × 107 | 5.8 × 108 |
| m50n1024 | 50 | 0.25 | 1 | 2.28 × 106 | 1.2 × 107 | 7.3 × 107 |
| m25n1024 | 25 | 0.125 | 2 | 2.85 × 105 | 1.5 × 106 | 9.1 × 106 |
| m12.5n1024 | 12.5 | 0.0625 | 5 | 3.56 × 104 | 1.88 × 105 | 1.14 × 106 |
| Name | |$L_{\rm box}^a$| | |$\epsilon _{\rm min}^b$| | |$z_{\rm end}^c$| | |$m^d_{\rm gas}$| | |$m^e_{\rm DM}$| | |$M^f_{\rm *,min}$| |
|---|---|---|---|---|---|---|
| m100n1024 | 100 | 0.5 | 0 | 1.82 × 107 | 9.6 × 107 | 5.8 × 108 |
| m50n1024 | 50 | 0.25 | 1 | 2.28 × 106 | 1.2 × 107 | 7.3 × 107 |
| m25n1024 | 25 | 0.125 | 2 | 2.85 × 105 | 1.5 × 106 | 9.1 × 106 |
| m12.5n1024 | 12.5 | 0.0625 | 5 | 3.56 × 104 | 1.88 × 105 | 1.14 × 106 |
Notes.aBox length in comoving |$h^{-1}\, {\rm Mpc}$|.
bMinimum gravitational softening length in comoving |$h^{-1}\, {\rm kpc}$|.
cEnding redshift (all begin at z = 249).
dInitial gas element mass resolution in |$\, \mathrm{M}_\odot$|.
eDark matter particle mass resolution in |$\, \mathrm{M}_\odot$|.
fMinimum stellar mass of a resolved galaxy in |$\, \mathrm{M}_\odot$|.
The simba simulation suite.
| Name | |$L_{\rm box}^a$| | |$\epsilon _{\rm min}^b$| | |$z_{\rm end}^c$| | |$m^d_{\rm gas}$| | |$m^e_{\rm DM}$| | |$M^f_{\rm *,min}$| |
|---|---|---|---|---|---|---|
| m100n1024 | 100 | 0.5 | 0 | 1.82 × 107 | 9.6 × 107 | 5.8 × 108 |
| m50n1024 | 50 | 0.25 | 1 | 2.28 × 106 | 1.2 × 107 | 7.3 × 107 |
| m25n1024 | 25 | 0.125 | 2 | 2.85 × 105 | 1.5 × 106 | 9.1 × 106 |
| m12.5n1024 | 12.5 | 0.0625 | 5 | 3.56 × 104 | 1.88 × 105 | 1.14 × 106 |
| Name | |$L_{\rm box}^a$| | |$\epsilon _{\rm min}^b$| | |$z_{\rm end}^c$| | |$m^d_{\rm gas}$| | |$m^e_{\rm DM}$| | |$M^f_{\rm *,min}$| |
|---|---|---|---|---|---|---|
| m100n1024 | 100 | 0.5 | 0 | 1.82 × 107 | 9.6 × 107 | 5.8 × 108 |
| m50n1024 | 50 | 0.25 | 1 | 2.28 × 106 | 1.2 × 107 | 7.3 × 107 |
| m25n1024 | 25 | 0.125 | 2 | 2.85 × 105 | 1.5 × 106 | 9.1 × 106 |
| m12.5n1024 | 12.5 | 0.0625 | 5 | 3.56 × 104 | 1.88 × 105 | 1.14 × 106 |
Notes.aBox length in comoving |$h^{-1}\, {\rm Mpc}$|.
bMinimum gravitational softening length in comoving |$h^{-1}\, {\rm kpc}$|.
cEnding redshift (all begin at z = 249).
dInitial gas element mass resolution in |$\, \mathrm{M}_\odot$|.
eDark matter particle mass resolution in |$\, \mathrm{M}_\odot$|.
fMinimum stellar mass of a resolved galaxy in |$\, \mathrm{M}_\odot$|.
We will also explore parameter space and compare to our previous mufasa simulations using |$50h^{-1}\, {\rm Mpc}$| runs with 2 × 5123 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 ‘No-jet’ 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.
To 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 interparticle 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 nH >0.13 H atoms cm−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.
Because 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 gravitationally 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 minimized.
Halos 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 interparticle spacing. We do not identify or consider sub-haloes in this work.
Galaxies and haloes 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 haloes so that any other properties can be computed via user-written python scripts.
All results shown here are obtained from the caesar catalogues 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 ≈250 GB in size, and the c aesar catalogues are typically ∼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.
To begin, we show in Fig. 2 a projected temperature map from the (50 Mpc h−1)3, 5123simba simulation. The slice shown is 10 Mpc h−1 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-hand panel), a significant volume of the IGM has been heated to high temperatures 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 neighbouring 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 ionization 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.
Temperature map projected through a random 10 Mpc h−1 slice from a 50 Mpc h−1simba 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 haloes, with cooler dense filamentary structures penetrating the hot gas.
Fig. 3 shows some examples of individual galaxies. We choose a Milky Way-sized disc galaxy at z = 0, with |$M_*\approx 4.7\times 10^{10}\, \mathrm{M}_\odot$| and SFR|$=1.3\ \, \mathrm{M}_\odot$| yr−1, and show the face-on (upper row) and edge-on (lower row) views, in both H2 surface density (left) and stellar mass surface density (right). The z = 2 galaxy shown in the bottom four panels has essentially the same M*, but with SFR|$=45\ \, \mathrm{M}_\odot$| yr−1 that is typical of a main sequence galaxy at Cosmic Noon. The z = 0 disc 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\times 10^7\, \mathrm{M}_\odot$| black hole accreting at |$0.005 \, \mathrm{M}_\odot$| yr−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 disc 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.
Examples of the molecular gas (left) and stellar (right) surface density distributions in star-forming disc galaxies with |$M_*\approx 4.7\times 10^{10}\, \mathrm{M}_\odot$| 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 disc in both H i and H2, while the high-z galaxy is clumpier and thicker.
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 characterizes the efficiency by which haloes 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* (e.g. Behroozi, Wechsler & Conroy 2013; Moster et al. 2013), and current models attribute this to self-regulation by star formation-driven feedback below L* and quenching of galaxies due to AGN feedback above L* (Somerville & Davé 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.
Fig. 4 shows the GSMF at z = 0.1, 1, 2, 3, 4, 6 from simba (green lines). Observational data are 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−1.8 + 0.3zGyr−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.
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 eight simulations sub-octants. Red and blue dashed lines show the mass functions of central galaxies below and above sSFR=10−1.8 + 0.3zGyr−1, respectively. Cyan-dotted line shows the results from EAGLE for comparison.
simba produces generally good agreement with the observed GSMFs at all redshifts, overall comparably well to EAGLE. There is excellent agreement at z ≥ 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 that can miss massive dusty galaxies.
At 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 underpredicts the GSMF around M⋆ by a factor of up to 2; 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.
The 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_*\gtrsim10^{11}\, \mathrm{M}_\odot$| objects. There are also increased observational uncertainties at the massive end. For instance, it is a matter of debate as to how many 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, Vikhlinin & Meshcheryakov 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^{-1}\, {\rm Mpc}$| volume, by z = 0, the largest halo has a virial mass of |$M_{\rm halo}=1.16\times 10^{15}\, \mathrm{M}_\odot$|, which is larger than expected by about 50 per cent 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 discrepancies reflect a significant failing of simba’s input physics.
Examining the SF versus Q samples, we see that massive quenched galaxies begin to appear in significant numbers at |$z\gtrsim2$|. By z = 1 they outnumber the SF galaxies among the most massive galaxies with |$M_*\gtrsim10^{11}\, \mathrm{M}_\odot$|, and by z = 0, they dominate at |$M_*\gtrsim2\times 10^{10}\, \mathrm{M}_\odot$|. The quenched population grows quickly at low redshifts, and the number of massive star-forming galaxies drops quickly since z ∼ 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 artefact of the friends-of-friends halo finder.
In 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−M*) relation. Unlike the GSMF that is often used as a primary constraint on models, SFR−M* is not, making it more of a true prediction of models. The SFR−M* 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.
During Cosmic Noon, it has long been seen that cosmological models tend to underpredict the main sequence amplitude (Daddi et al. 2007; Davé 2008; Narayanan et al. 2012; Somerville & Davé 2015; Sparre et al. 2015), typically by a factor of 2–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-ionizing stellar populations at high-z, or having a more top-heavy IMF.
Fig. 5 shows the specific SFR−M* 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 star-forming galaxies only (simba-SF), defined as before by sSFR>10−1.8 + 0.3zGyr−1 (dotted horizontal yellow line). The error bars show the 1σ 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 MBH/M*; points at |$M_*\lesssim10^{10}\, \mathrm{M}_\odot$| 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 low-z are shown from the GALEX-SDSS-WISE Legacy Catalog Salim et al. (GSWLC; 2016); Salim, Boquien & Lee (GSWLC; 2018), shown as grey hexbins, and the running median to the star-forming galaxies with the same criterion as above are shown as the black dashed line, along with error bars showing the 1σ spread around the median. At high-z, we show the median sSFR−M* relation measured for 2 < z < 2.5 galaxies by Whitaker et al. (2014). Finally, results from EAGLE are shown as the magenta dotted line.
Star formation rate–stellar mass relation at z ≈ 0 (left) and z ≈ 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σ spread around the running median value, typically 0.3–0.4 dex. At high-z we show the black dashed line as the best-fitting 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 high-z sSFR determinations, though in small galaxies it appears to overpredict the SFR at z ∼ 0 and underpredict at z ∼ 2.
At z = 0 (left-hand panel), simba nicely reproduces the observed GSWLC main sequence slope and amplitude at |$M_*\gtrsim10^{10}M_{\odot }$|. 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 magenta-coloured points with higher black hole mass for their M*. Indeed, there is a very strong trend that the galaxies that are quenching are specifically the ones with a high MBH/M* ratio; massive galaxies left on the main sequence at z ∼ 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−M* for star-forming galaxies over the entire mass range plotted is −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–0.4 dex, with a mild tendency to drop with M*; this is very comparable to that seen in the GSWLC data.
At z = 2.2 (right-hand panel), the simba main sequence generally tracks the observed one from Whitaker et al. (2014), but is low in amplitude by ≈×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) point out that more sophisticated SED fitting applied to the latest data sets can lead to a systematic increase in the inferred M* while lowering the SFR that results in a combined ≈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 ∼ 2 observations at long last.
Finally, we show in Fig. 6 histograms of the specific SFR, broken up into mass bins of |$10^9\lt M_*\lt 10^{10}\, \mathrm{M}_\odot$|, |$10^{10}\lt M_*\lt 10^{11}\, \mathrm{M}_\odot$|, and |$M_*\gt 10^{11}\, \mathrm{M}_\odot$|. 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≤10−2.5 Gyr−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 quenched. There is impressive agreement in the lowest sSFR bin, showing that simba well reproduces the quenched fractions at various M*. However, the low-mass galaxies in simba have somewhat too high sSFR values, reflecting the same excess as seen in Fig. 5.
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 ∼ 0.1 observations from GSWLC-D2. All galaxies with sSFR<10−2.5 Gyr−1 are placed in the lowest bin. There is good agreement, particularly in the quenched fractions in more massive galaxies, though Simba produces somewhat too high sSFRs at low-M*.
In summary, simba generally reproduces the main sequence of star-forming galaxies as seen at z ≈ 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*\gtrsim10^{10}\, \mathrm{M}_\odot$| is encouraging because it suggests that 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* 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 sub-M⋆ 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 homogenized database for the SFRD that can be more robustly compared.
Fig. 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_*\gt 5.8\times 10^8 \, \mathrm{M}_\odot$|) makes a negligible difference.
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 best-fitting line).
The overall shape of the predicted SFRD versus time is in good agreement with observations. simba matches the present-day SFRD very well, and generally reproduces the order-of-magnitude rise in SFRD towards the peak at z ∼ 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-fitting line from Madau & Dickinson (2014). The peak SFRD at z ∼ 2 is also slightly lower than observed, following the trend shown in Fig. 5 that simba has slightly lower main sequence than observed. Despite these minor differences, the overall shape and amplitude are in very good agreement with observations, comparable 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 (H2) 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 H2 gas fractions versus M*. Recent millimetre and radio observation data have greatly expanded our knowledge of gas contents for low-z galaxies, with constraints at higher z promising continued rapid advancement in the near future.
Fig. 8 shows the scaling relations for H2 (top) and H i (bottom) mass fractions, versus M*. The points show individual galaxies at z = 0 colour-coded by their sSFR deviation from the main sequence at that M* (Δ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 H2 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.
Molecular (top) and neutral (bottom) gas fractions MH2/M* and MH i/M* as a function of M*. The points show z = 0 values from simba colour-coded by the deviation in sSFR from the star-forming main sequence – bluer points have higher-than 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 fH2 from xCOLDGASS (Saintonge et al. 2017) are shown in the top panel, and observations of fH i 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.
Overall, 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-H2 conversion factor that is poorly determined particularly in the low-mass (|$M_*\lt 10^{10}\, \mathrm{M}_\odot$|) regime (e.g. Narayanan et al. 2012; Bolatto, Wolfire & Leroy 2013), as well as theoretical uncertainties in the approximate way that self-shielding is applied, galaxy gas contents can be considered to be a remarkably successful prediction of simba.
We note that our massive galaxies have a non-trivial amount of cold gas, despite their very low sSFR. This was not the case in mufasa (Davé 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 per cent 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 ρ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 are 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.
Finally, the colours of the simba points in Fig. 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* have a higher SFR. This is unsurprising for H2 in our models, given that star formation is tied to the molecular content. It is somewhat more surprising to see this for the H i fraction, but such a correlation was also seen fairly strongly in mufasa (Davé et al. 2017a), though not as strong as for H2 (Rafieferantsoa, Davé & Naab 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 Section 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.
Metals 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é 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.
Fig. 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 per cent (Asplund et al. 2009). Points show central galaxies colour-coded by their deviation from the main sequence, as in Fig. 8; black points are satellite galaxies. A running median for star-forming galaxies (sSFR>10 − 1.8 + 0.3zGyr−1) is shown as the dashed green line. Fits to observations at z = 0 are shown from strong emission line fitting (black dashed; Tremonti et al. 2004), stacked measurement of direct metallicities (black dot-dashed; Andrews & Martini 2013), and individual semidirect metallicities (black dotted; Yates et al. 2019). Observations at z = 2.3 are shown from the MOSDEF survey (Sanders et al. 2015).
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 low-z, best-fits to observations from Tremonti et al. (2004), Andrews & Martini (2013), and Yates et al. (2019) are shown as the black lines. At z ≈ 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 low-z.
simba 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.
mufasa 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 SN ii yields (Davé et al. 2017a). In simba, metals are locked into dust, 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 preparation).
Since 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 ∼ 2 → 0 is approximately correct in simba.
simba also shows a clear second parameter dependence on the specific SFR such that at a given M*, 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ópez et al. 2010). The existence of the FMR remains somewhat controversial (Salim et al. 2014; Sánchez et al. 2017; Cresci, Mannucci & Curti 2018), but a careful analysis of the MOSDEF data has revealed such a trend at z ∼ 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, Gallazzi & van den Bosch 2012).
Fig. 10 shows the mass-weighted stellar metallicity as a function of M* at z = 0. Points show centrals colour-coded 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).
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) and 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.
simba 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 s imba produces a low-mass stellar MZR that is somewhat too high compared to observations.
In 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 normalization 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 Section 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 gas 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.
To conduct a fair comparison to observed sizes, we compute projected galaxy half-light radii in the Rband. 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, and z axes. Fig. 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 (simba-Q) and star-forming (simba-SF) populations (magenta and cyan lines), where we divide the two populations at sSFR=10−1.8 Gyr−1 as before. For comparison we show observations from SDSS (Zhang & Yang 2019). The SDSS sample has been subdivided into red and blue galaxies, albeit with a criterion based on photometry, not sSFR.
R-band 2-D projected half-light radii of simba galaxies at z = 0, as a function of M*. Points are colour-coded by sSFR. We fit median relations to the red and blue sub-samples, delineated by 10−1.8 Gyr−1, as the cyan and magenta dashed lines, respectively. Observational relations are shown from Zhang & Yang (2019) 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.
Star-forming galaxy sizes in simba show an amplitude and scaling with M* that agrees quite well with the observed slope, which is encouraging. There is a suggestion that low-mass galaxies are too small, but this occurs in the mass range 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.
In 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 ∼× 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 disc 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 haloes. Such energy input can concurrently have a strong impact on the amount and thermal state of hot gas within those haloes. In particular, it can evacuate gas even from fairly sizeable haloes, 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 intragroup and intracluster gas. Such observations thus provide an important constraint on the AGN feedback model.
The 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 haloes compared to observations (Genel et al. 2014), which provided one motivation for the new AGN feedback model in Illustris-TNG (Weinberger et al. 2018). None the less, 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.
Fig. 12 shows the baryon fractions as a function of halo mass (M500) from simba. M500 is computed as the radius enclosing 500 times the critical density, centred on the most bound halo particle. Black points show the total baryon fraction, red points show hot gas (T > 105K) fractions, and blue points show cold gas (T < 105K). 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 Ωb = 0.048, so a halo at unity has its cosmic share of baryons.
Gas fractions as a function of M500 at z = 0. Black points show the total baryon fraction within the halo, and red and blue points show the gas fractions subdivided at 105 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 M500, which has been difficult for previous simulations to achieve without tuning.
At |$10^{12}\, \mathrm{M}_\odot$| haloes have about 40 per cent of the cosmic baryon fraction within R500, with a large scatter. This dips to |${\sim } 30{{\ \rm per\ cent}}$| at |$10^{12.5-13}\, \mathrm{M}_\odot$|, 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 per cent for all haloes 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 haloes, 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.
The 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 105K. Comparing to observations, we see that simba’s haloes 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 preparation), we will examine quantitatively where these evacuated baryons end up.
3.8 Black hole mass versus 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 ϵm = 10 per cent, while most other cosmological simulations (based on Bondi accretion) tune the AGN feedback efficiency. In previous works, Anglés-Alcázar et al. (2015) and Anglés-Alcázar et al. (2017a) showed that the MBH−M* 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.
Fig. 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 reverberation-mapped 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 forming (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.
MBH−M* 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) show the relationship more appropriate for spiral star-forming systems at lower M*. simba broadly reproduces these observed relations in its appropriate galaxy populations.
simba black holes generally lie in the range of observations. Although we tuned ϵ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 sub-samples. 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*, 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.
At the low-mass end, consistent with Anglés-Alcázar et al. (2015), torque-limited accretion grows black holes very quickly once the galaxy stellar mass exceeds |$3\times 10^9\, \mathrm{M}_\odot$|, which is where we choose to seed the black holes in simba. Hence the rapid rise is not directly physical but a numerical artefact 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és-Alcázar 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.
Overall, 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 DGR.
Fig. 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, Eales & Edmunds (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 DMF is dominated by star-forming galaxies (blue dashed line), which is as observed (Beeston et al. 2018).
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.
The z = 2 DMF is compared to observations from Dunne et al. (2003), and shows a deficit of ∼× 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), owing in part to tuning of each model, but our z = 2 DMF is significantly higher than theirs.
Fig. 15 shows the z = 0 DGR as a function of gas-phase metal abundance. Points are colour-coded by sSFR. simba is in good agreement with the data shown (crosses Rémy-Ruyer et al. 2014), showing a slope of increasing DGR with metallicity as observed. In massive quenched systems, the DGR drops quickly.
Dust-to-gas mass 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émy-Ruyer et al. (2014). simba reproduces the observed trend and amplitude in DGRs.
Finally, Fig. 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_*\gtrsim10^{9.5}\, \mathrm{M}_\odot$|, the fraction is typically 30–40 per cent, 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 × 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. We will examine galaxy dust content and evolution in significantly more detail in forthcoming work (Li et al., in preparation), but these preliminary comparisons suggest that simba’s dust tracking model yields plausible galaxy dust contents.
Metal mass fraction locked in dust, as a function of M*, 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.
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 |$50h^{-1}\, {\rm Mpc}$|, 2 × 5123 simulations run to z = 0, with simba input physics and parameters except for the AGN feedback variations.
Fig. 17 shows the GSMF for the full s imba physics run, a No-X run turning off only X-ray feedback, and a No-jet run turning off both jet and X-ray feedback, at z = 2, 1, 0. Observations are overplotted as described in Fig. 4. We do not show z ≥ 3 results because these variants’ GSMFs are indistinguishable there.
Stellar mass function evolution from z = 2 → 0 in |$50h^{-1}\, {\rm Mpc}$|, 5123 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 (No-Xray; red dashed). For comparison we show the main |$100h^{-1}\, {\rm Mpc}$|simba run (green dashed) reproduced from Fig. 4, as well as selected observations as indicated. Turning on just the jet feedback (No-X-ray) results in a substantial truncation of the GSMF that does not occur without jets (No-jet).
Looking at the z = 0 panel, without jets, the GSMF is strongly overproduced at high masses. Turning on jets results 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 ∼ 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_*=10^{12}\, \mathrm{M}_\odot$| 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.
A more detailed view of how AGN feedback impacts both stellar and black hole growth can be obtained by examining the MBH−M* relation in these variants, shown in Fig. 18, with galaxies colour-coded by specific SFR as in Fig. 13. For clarity we show only central galaxies.
MBH−M* 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 Fig. 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.
Comparing the No-Jet version (top panel) to the original simba in Fig. 13 highlights several key points. As expected, the sub-M⋆ objects show little difference in the trends, in either MBH or sSFR. However, for massive galaxies, the full simba run shows significantly lower sSFR and somewhat higher MBH, particularly for the most massive galaxies. This demonstrates more explicitly that the AGN jet feedback is crucial for quenching galaxies.
Interestingly, the No-jet run still shows a few quenched galaxies at high MBH 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š et al. 2005; Gabor & Davé 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).
So 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* (Thomas et al., in preparation) 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 MBH (equation 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.
The 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* values in this simulation are substantially higher, which moves galaxies leftwards in the MBH−M* diagram; the black hole masses themselves are not substantially different. The relative roles of Bondi versus 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és-Alcázar et al., in preparation).
The 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 X-ray 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−1 box run with 10243 dark matter and 10243 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 hole, 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.
Here are our main findings:
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 substantially in numbers at |$z\lesssim2$|, and by z = 0 they dominate at |$M_*\gt 2\times 10^{10}\, \mathrm{M}_\odot$|.
simba’s star-forming main sequence is in good agreement with observations at z = 0, and low at z ≈ 2 by only a factor of 2, which is explainable via observational systematics. Predicted quenched fractions at z = 0 as a function of M* are in good agreement with observations. Galaxies at a given M* that quench first in simba have preferentially larger black holes.
simba gas fractions, both neutral and molecular, show a dropping trend with M* 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.
The gas–phase and stellar mass–metallicity relations generally agree with observations at z = 0 and z ∼ 2. The MZR evolves upwards by a factor of ∼× 3 for our smallest systems, but at |$M_*\gtrsim10^{11}\, \mathrm{M}_\odot$| there is little evolution.
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.
There is substantial evacuation of baryons from haloes 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.
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* have lower black hole masses.
simba predicts a z = 0 dust mass function and DGRs in agreement with observations, with more massive star-forming galaxies having higher ratios. At a given metallicity, galaxies that have higher SFR have higher DGRs. Roughly one-third of metals are locked in dust.
These 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 suggestions that s imba overproduces the stellar metallicities as well as the specific SFRs in low-mass present-day star-forming systems. Finally, in many low-z scaling relations (e.g. gas and metallicities versus M*) there is an abrupt break in the typical properties of galaxies above and below |$M_*\approx 2\times 10^{10}\, \mathrm{M}_\odot$|, 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 s imba (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 e-Infrastructure.
