# Semi-Analytic Galaxy Evolution (SAGE): Model Calibration and Basic Results

Darren J. Croton, Adam R. H. Stevens, Chiara Tonini, Thibault Garel, Maksym Bernyk, Antonio Bibiano, Luke Hodkinson, Simon J. Mutch, Gregory B. Poole, Genevieve M. Shattow

► To cite this version:

Darren J. Croton, Adam R. H. Stevens, Chiara Tonini, Thibault Garel, Maksym Bernyk, et al. Semi-Analytic Galaxy Evolution (SAGE): Model Calibration and Basic Results. The Astrophysical Journal Supplement Series, 2016, 222, 10.3847/0067-0049/222/2/22. insu-03710587

HAL Id: insu-03710587

https://insu.hal.science/insu-03710587v1

Submitted on 21 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

SEMIANALYTIC GALAXY EVOLUTION (SAGE): MODEL CALIBRATION AND BASIC RESULTS

DARREN J. CROTON<fnref n="1" />, ADAM R. H. STEVENS<fnref n="1" />, CHIARA TONINI <fnref n="1" />,<fnref n="2" />, THIBAULT GAREL <fnref n="1" />,<fnref n="3" />,<fnref n="4" />, MAKSYM BERNYK<fnref n="1" />, ANTONIO BIBIANO<fnref n="1" />, LUKE HODKINSON<fnref n="1" />, SIMON J. MUTCH <fnref n="1" />,<fnref n="2" />, GREGORY B. POOLE <fnref n="1" />,<fnref n="2" />, AND GENEVIEVE M. SHATTOW<fnref n="1" />

Received 2015 October 6; accepted 2016 January 17; published 2016 February 18

ABSTRACT

This paper describes a new publicly available codebase for modeling galaxy formation in a cosmological context, the "Semi-Analytic Galaxy Evolution" model, or SAGE for short.<fnref n="5" /> SAGE is a significant update to the 2006 model of Croton et al. and has been rebuilt to be modular and customizable. The model will run on any N-body simulation whose trees are organized in a supported format and contain a minimum set of basic halo properties. In this work, we present the baryonic prescriptions implemented in SAGE to describe the formation and evolution of galaxies, and their calibration for three N-body simulations: Millennium, Bolshoi, and GiggleZ. Updated physics include the following: gas accretion, ejection due to feedback, and reincorporation via the galactic fountain; a new gas cooling-radio mode active galactic nucleus (AGN) heating cycle; AGN feedback in the quasar mode; a new treatment of gas in satellite galaxies; and galaxy mergers, disruption, and the build-up of intra-cluster stars. Throughout, we show the results of a common default parameterization on each simulation, with a focus on the local galaxy population.

Key words: galaxies: active – galaxies: evolution – galaxies: halos – methods: numerical

# 1. INTRODUCTION

Developing a complete theory of galaxy evolution is a formidable task. Without the ability to construct real universes in a laboratory, we are left to test ideas through conducting supercomputer simulations and comparing their results against what we observe. Arguably the most thorough way of doing this currently is through cosmological hydrodynamic simulations (e.g., Carlberg et al. 1990; Dubois et al. 2014b; Vogelsberger et al. 2014; Khandai et al. 2015; Schaye et al. 2015), where the physics of baryons and the dark sector are self-consistently considered. However, one can successfully reproduce the large-scale structure and formation sites of galaxies inside halos with pure N-body simulations (e.g., Davis et al. 1985; Springel et al. 2005; Kim et al. 2011; Klypin et al. 2011; Skillman et al. 2014; Poole et al. 2015), into which galaxies can be later added in post-processing. Such simplicity can be significantly advantageous.

By foregoing a simulation of costly hydrodynamic processes, simulators can invest their computational resources in increasing the number of particles. Structures on both smaller and larger scales are then resolved through improved particle mass resolution and by simulating larger volumes, respectively. Larger volumes also lead to less biased sampling and more halos overall, which in turn allow for better statistical significance. For example, the main run of the hydrodynamic EAGLE simulations (Schaye et al. 2015) used 1504³ dark matter particles, but required approximately the same number of floating point operations as the main run of the pure N-body Dark Sky simulations (Skillman et al. 2014) with 10,240³ particles. The addition of hydrodynamics comes at the cost of approximately two orders of magnitude in particle number.

Semi-analytic models of galaxy evolution take advantage of the relative computational efficiency of N-body simulations by

adding the bound baryons to a simulation as a post-processing step. By using information about the gravitationally bound halos, such as their mass, size, spin, substructure, and merger history, the properties of galaxies hosted within these structures can be inferred through differential equations that describe the relevant physics macroscopically. The runtime of such a model is orders of magnitude less than the N-body simulation itself. Hence, after that initial investment one can explore large regions of the parameter space that underlie the baryonic physics by running the model many times.

In the original implementation of the semi-analytic method (White & Frenk 1991), halos were drawn statistically from the Press–Schechter formalism for gravitational collapse (Press & Schechter 1974; later extended by Bond et al. 1991; Bower 1991), and the analytic baryonic physics was based on the milestone theory of White & Rees (1978). Kauffmann et al. (1993) introduced the use of merger trees for a model, where individual systems could then be tracked across time statistically within Press–Schechter theory. Merger trees produced from an N-body simulation should always be more representative of the real universe, and Kauffmann et al. (1999) were the first to utilize this technique. Now, it is standard practice for semi-analytic models to be coupled to an N-body simulation (e.g., Hatton et al. 2003; Bower et al. 2006; Croton et al. 2006; Benson 2012; Henriques et al. 2015). For a more in-depth look at the history of semi-analytic models, see Baugh (2006).

In this paper, we present our new semi-analytic model, SAGE (Semi-Analytic Galaxy Evolution), which updates the work of Croton et al. (2006, hereafter C06). The new model revamps many prescriptions for the treatment of baryons, including the suppression of cooling within halos from active galactic nucleus (AGN) feedback, reincorporation of ejected gas, and the stripping of gas from satellite systems. While semi-analytic models have historically been designed and calibrated for a single N-body simulation, SAGE is designed to be run on any

<footnote n="1">¹Centre for Astrophysics & Supercomputing, Swinburne University of Technology, P.O. Box 218, Hawthorn, Victoria 3122, Australia</footnote>

<footnote n="2">²School of Physics, University of Melbourne, Parkville, Victoria 3010, Australia</footnote>

<footnote n="3">³Centre de Recherche Astrophysique de Lyon, Université de Lyon, Université Lyon 1, CNRS, Observatoire de Lyon,
9 avenue Charles André, F-69561 Saint-Genis Laval Cedex, France</footnote>

<footnote n="4">⁴ Australian Research Council Super Science Fellow.</footnote>

<footnote n="5">⁵ https://github.com/darrencroton/sage</footnote>
simulation, so long as the merger trees are provided in an
appropriate format. We show SAGE's performance on three
cosmological N-body simulations; namely the Millennium
(Springel et al. 2005), Bolshoi (Klypin et al. 2011), and
GiggleZ (Poole et al. 2015) simulations, all of which follow the
standard ΛCDM cosmological paradigm. That said, there is
nothing to stop SAGE being run on simulations that explore
alternate gravity models, or different dark matter candidates,
for example.

Along with those from other semi-analytic models, SAGE galaxy catalogs built on various N-body simulations are available through the Theoretical Astrophysical Observatory (Bernyk et al. 2016).<fnref n="6" /> The codebase of SAGE is also publicly available (see footnote 5), allowing the community to build further models, or modify those described here. The repository includes an IPYTHON notebook for conducting simple analysis with SAGE output, specifically showing how we produced some of the figures presented in this paper.<fnref n="7" />

We lay out the sections of this paper in the following
manner. The *N*-body simulations used as input for SAGE are
summarized in Section 2. We provide an overview of SAGE in
Section 3, covering which components of the model have been
upgraded from that of C06. Sections 4–12 then describe the
physics of the model in more detail, covering the following: gas
infall in halos (4); the role of reionization (5); cooling of gas
from the hot halo (6); the consideration of cold gas, star
formation, and metal enrichment in galactic disks (7); the role
of supernova feedback (8); the growth of black holes and their
associated AGN feedback (9); dealing with mergers and intra-
cluster stars (10); disk instabilities (11); and starbursts (12). We
finish with some discussion and concluding remarks in
Section 13.

Throughout, we present all results assuming *h* = 0.73, based
on the cosmology used for the Millennium simulation, the
primary simulation used for calibrating the model. Where
relevant, we also use a Chabrier (2003) initial mass function to
produce stellar masses.

## 2. N-BODY SIMULATIONS

In this paper, we present the performance of SAGE on three cosmological simulations. Each simulation not only varies in terms of its cosmological parameters, having used the best available measurements at their various times of being run, but also in terms of the codes and pipelines used to generate the final data products. Even for simulations with identical initial conditions and cosmological parameters, the use of different codes for either running or post-processing cosmological simulations can lead to non-trivial differences in their results. Knebe et al. (2011, 2013) assessed how and why the choice of (sub)halo finder can affect results, while complimentary studies investigated how the choice of merger tree code can change the derived structure formation histories and the consequences this has for semi-analytic models (Srisawat et al. 2013; Lee et al. 2014, respectively). Understanding these non-trivial differences is key to understanding the theoretical and numerical uncertainties associated with semi-analytic models, although we do not attempt to provide a detailed analysis of

such effects here. We describe each simulation below and
summarize their properties in Table 1.

2.1. *The Millennium Simulation*

The Millennium simulation (Springel et al. 2005) signifi-
cantly upped the ante in cosmological simulations, boasting
unrivalled detail for its time. It has since been the focus of a
plethora of scientific studies.<fnref n="8" /> Eleven years on from its
completion, the simulation remains a benchmark, and continues
to be used for science. Of the simulations used in this paper, it
remains the most balanced in terms of (cosmologically
representative) size and resolution (cf. Table 1).

Millennium was run using the popular GADGET-2 code (Springel 2005). The cosmological parameters followed those from a combined analysis of WMAP1 (Wilkinson Microwave Anisotropy Probe, first year) data (Spergel et al. 2003) and the Two-Degree-Field Galaxy Redshift Survey (Colless et al. 2001). Arguably, the biggest weakness of Millennium is its dated cosmological parameters, which now differ significantly from the best-fitting values which are more precisely measured (for the latest observational results, see Planck Collaboration et al. 2015).

Structure and subhalos were identified in the Millennium
simulation with the SUBFIND algorithm (Springel et al. 2001).
Parent halos are found through a friends-of-friends procedure,
while subhalos are defined as having at least 20 gravitationally
bound particles by the halo finder. The merger trees which feed
SAGE were constructed with the L-HALOTREE code (described in
the supplementary information of Springel et al. 2005).

## 2.2. *The Bolshoi Simulation*

Bolshoi (Klypin et al. 2011) was run using the adaptive-
mesh-refinement code ART (Adaptive Refinement Tree, Kravt-
sov et al. 1997). The chosen cosmological parameters were
very close to those of the WMAP7 data (Jarosik et al. 2011),
while maintaining consistency with WMAP5 (Dunkley
et al. 2009; also see Komatsu et al. 2009). When compared
with WMAP1, the data from these WMAP releases describe a
universe with a greater average matter density, that presently
expands more slowly, with smaller mass fluctuations. While
smaller in box size, Bolshoi complements the results from
Millennium due to its higher mass resolution, allowing us to
probe the low-mass end of the mass function in more detail.

Subhalos in Bolshoi were identified with the **ROCKSTAR** algorithm (Behroozi et al. 2013a),<fnref n="9" /> which builds a hierarchy of friends-of-friends subgroups and utilizes one temporal and six phase-space quantities to determine which particles constitute those subhalos. Merger trees were subsequently constructed using **CONSISTENT-TREES** (Behroozi et al. 2013b).<fnref n="10" />

## 2.3. *The GiggleZ Simulation Suite*

The Gigaparsec WiggleZ simulation suite (GiggleZ, Poole
et al. 2015) was run as a theoretical counterpart to the WiggleZ
Dark Energy Survey (Drinkwater et al. 2010). Each simulation
was performed using GADGET-2, with cosmological parameters
based on data from WMAP5, baryonic acoustic oscillations, and

<footnote n="6">⁶ https://tao.asvo.org.au</footnote>

<footnote n="7">⁷ This can be viewed directly on or downloaded from GitHub at the above link. Once downloaded, the notebook will function without needing to run SAGE first (it will fetch output).</footnote>

<footnote n="8">⁸ See http://www.mpa-garching.mpg.de/millennium/ for an up-to-date list.</footnote>

<footnote n="9">⁹ Robust Overdensity Calculation using K-Space Topologically Adaptive Refinement, available at https://bitbucket.org/gfcstanford/rockstar.</footnote>

<footnote n="10">¹⁰ https://bitbucket.org/pbehroozi/consistent-trees.</footnote>
<table><caption>**Table 1**
Details of the *N*-body Simulations Used in the Analysis of SAGE in This Paper</caption><thead><tr><th>Simulation</th><th>N<sub>part</sub></th><th>M<sub>part</sub>h<br>(M<sub>☉</sub>)</th><th>l<sub>box</sub>h<br>(cMpc)</th><th>Ω<sub>M</sub></th><th>σ<sub>8</sub></th><th>Code</th><th>Subhalo Finder</th><th>Tree Constructor</th></tr></thead><tbody><tr><td>Millennium</td><td>2160<sup>3</sup></td><td>8.60 × 10<sup>8</sup></td><td>500</td><td>0.250</td><td>0.900</td><td>GADGET-2</td><td>SUBFIND</td><td>L-HALOTREE</td></tr><tr><td>Bolshoi</td><td>2048<sup>3</sup></td><td>1.35 × 10<sup>8</sup></td><td>250</td><td>0.270</td><td>0.820</td><td>ART</td><td>ROCKSTAR</td><td>CONSISTENT-TREES</td></tr><tr><td>GiggleZ-MR</td><td>520<sup>3</sup></td><td>9.50 × 10<sup>8</sup></td><td>125</td><td>0.273</td><td>0.812</td><td>GADGET-2</td><td>SUBFIND</td><td>G. B. Poole et al. (2016, in preparation)</td></tr></tbody></table>

**Note.** The columns provide particle number, $N_{part}$; particle mass, $M_{part}$; periodic box length, $l_{box}$, in comoving units; contribution of matter to the average universal energy density at the present epoch, $\Omega_M$ (for which the equivalent for dark energy is $\Omega_\Lambda = 1 - \Omega_M$); the redshift-zero extrapolation of the root-mean-square linear mass fluctuation within a sphere of radius $8 h^{-1} \text{ Mpc}$, $\sigma_8$; the code the simulation was run with; the code used to identify subhalos; and the code used to build the merger trees.

supernovae (Komatsu et al. 2009). While the main box of GiggleZ boasts shear size ($1h^{-1}$ comoving Gpc on a side), due to its lower mass resolution we instead assess the complimentary GiggleZ-MR simulation for this paper. GiggleZ-MR was run in a smaller box of side length $125 h^{-1}$ comoving Mpc but with a particle mass much closer to Millennium.

GiggleZ subhalos were identified with SUBFIND. Trees were built following the method in Poole et al. (2016, in preparation). This approach repairs pathological defects in merger trees introduced by the halo-finding process (e.g., over linking, or the disappearance of halos during pericentric passages) through a process of forward and backward matching which scans both ways over multiple snapshots.

## 2.4. Halo Merger Tree Structure and Required Properties

SAGE is modular enough that it should run on any halo merger tree that is structured in a supported format and contains a minimum set of information per halo. The initial public release of SAGE requires halo measurements for the following:

1. $M_{vir}$, the halo virial mass, and “Len,” the number of particles in the (sub)halo;

2. $V_{vir}$ and $R_{vir}$, measured using the standard virial relations;

3. $V_{max}$, a fit to the maximum circular velocity of the halo; and

4. the cartesian position, velocity, and spin vector of each halo.

Additional input properties may be required depending on the version of SAGE being used. This can be checked in `coreSimulation.h`, where the tree file structure is defined, and also directly in each relevant core or science module.

In terms of tree structure, SAGE assumes the halo trees are organized in depth-first order, as described in the supplementary information of Springel et al. (2005, see, in particular, their Figure 5). This is the same output format used for the trees of the Millennium simulation (Section 2.1), which were constructed with the `GADGET-2/L-HALOTREE` codebase. `L-HALOTREE` produces additional output properties required by SAGE that enable the membership and history of each halo to be determined:

1. “FirstProgenitor,” “NextProgenitor,” “FirstHaloIn-FOFgroup,” and “NextHaloInFOFgroup.”

Given the modular nature of SAGE, other tree formats, such as those produced by `ROCKSTAR/CONSISTENT-TREES`, may be supported in the future.

## 3. AN OVERVIEW OF SAGE

The new SAGE model of galaxy formation is an updated version of that described in C06. One of the most important aspects of this update is that the code has been cleaned, significantly optimized, generalized to run on a wide variety of simulations, and is ready to be used by the wider community. Key changes to the physics that distinguish SAGE from the C06 model are summarized below. A more detailed description of each model component is then given in the subsequent sections, where we also explore their consequences on the galaxy population and its evolution.

1. **Gas cooling and AGN heating:** Cooling and heating of halo gas is now more directly coupled in the new SAGE model. We introduce a heating radius from radio mode AGN feedback, interior to which gas will never cool. This radius can only move outward, thus retaining the memory of previous heating episodes. The cooling rate is then calculated using only the hot gas between the cooling and heating radii.

2. **Quasar mode feedback:** We have added feedback from the quasar mode, the dominant growth channel of the black hole, triggered by mergers or disk instabilities. This feedback is most effective at removing disk (and sometimes halo) gas at high redshift in gas rich galaxies.

3. **Ejected gas reincorporation:** Gas ejected from the halo is now reincorporated according to the dynamical time of the halo modulated by $V_{vir}/V_{crit}$, with $V_{crit}$ a parameter set at the galaxy group halo mass scale. Previously, reincorporation was dependent on the dynamical time alone.

4. **Satellite galaxies:** Hot gas is no longer instantaneously removed from the satellite/subhalo system upon infall, but stripped in proportion to the subhalo dark matter mass stripping rate. Satellites are treated in the same way as central galaxies for the longest time possible (e.g., we allow cooling in subhalos).

5. **Mergers and intra-cluster stars:** At the moment of infall, an expected satellite merger time is calculated. The satellite is then tracked until its baryon-to-subhalo mass falls below a critical threshold (taken at ~1). At this point the current survival time is compared to the expected merger time. If the subhalo has survived longer than expected we say it is more resistant to disruption and the satellite is merged with the central in the usual way. Otherwise the satellite is disrupted and its stars are added to a new intra-“cluster” mass component. As a consequence, SAGE no longer produces satellite galaxies lacking a subhalo, the so-called orphan population.
<table><caption>**Table 2**
Fiducial SAGE Parameters used Throughout This Text, and Also Compared to Those From C06</caption><thead><tr><th>Parameter</th><th>Description</th><th>Value</th><th>C06 Value</th><th>Fixed</th><th>Section(s)</th></tr></thead><tbody><tr><td>f<sub>b</sub><sup>(cosmic)</sup></td><td>(Cosmic) baryon fraction</td><td>0.17, 0.13</td><td>0.17</td><td>No</td><td>4, 5</td></tr><tr><td>z<sub>0</sub></td><td>Redshift when H &infin; regions overlap</td><td>8.0</td><td>8.0</td><td>Yes</td><td>5</td></tr><tr><td>z<sub>r</sub></td><td>Redshift when the intergalactic medium is fully reionized</td><td>7.0</td><td>7.0</td><td>Yes</td><td>5</td></tr><tr><td>&alpha;<sub>SF</sub></td><td>Star formation efficiency</td><td>0.05</td><td>0.07</td><td>No</td><td>7</td></tr><tr><td>Y</td><td>Yield of metals from new stars</td><td>0.025</td><td>0.03</td><td>No</td><td>7</td></tr><tr><td>R</td><td>Instantaneous recycling fraction</td><td>0.43</td><td>0.30</td><td>Yes</td><td>7, 8</td></tr><tr><td>&epsilon;<sub>disk</sub></td><td>Mass-loading factor due to supernovae</td><td>3.0</td><td>3.5</td><td>No</td><td>8</td></tr><tr><td>&epsilon;<sub>halo</sub></td><td>Efficiency of supernovae to unbind gas from the hot halo</td><td>0.3</td><td>0.35</td><td>No</td><td>8</td></tr><tr><td>k<sub>reinc</sub></td><td>Sets velocity scale for gas reincorporation</td><td>0.15</td><td>N/A</td><td>Yes</td><td>8</td></tr><tr><td>&kappa;<sub>R</sub></td><td>Radio mode feedback efficiency</td><td>0.08</td><td>1.0</td><td>No</td><td>9.1</td></tr><tr><td>&kappa;<sub>Q</sub></td><td>Quasar mode feedback efficiency</td><td>0.005</td><td>N/A</td><td>No</td><td>9.2</td></tr><tr><td>f<sub>BH</sub></td><td>Rate of black hole growth during quasar mode</td><td>0.015</td><td>0.03</td><td>No</td><td>9.2</td></tr><tr><td>f<sub>friction</sub></td><td>Threshold subhalo-to-baryonic mass for satellite disruption or merging</td><td>1.0</td><td>N/A</td><td>Yes</td><td>10</td></tr><tr><td>f<sub>major</sub></td><td>Threshold mass ratio for merger to be major</td><td>0.3</td><td>0.3</td><td>Yes</td><td>10</td></tr></tbody></table>

**Note.** The fifth column indicates whether the value was kept fixed or used in the calibration. The last column indicates the section where the parameter is discussed. The cosmic baryon fraction used for Millennium and GiggleZ-MR is 0.17, while 0.13 was used for Bolshoi. All other parameters are common across all three simulations.

### 3.1. Model Calibration

A summary of the primary SAGE parameters, along with the choices that define our fiducial galaxy model, are given in Table 2 for quick reference. These parameters were manually selected to simultaneously perform well across all three simulation sets, with a slightly higher emphasis on Millennium. All figures and results in this paper are calculated using a model with these parameter choices. Note that, due to the order-of-magnitude higher resolution of Bolshoi, we found it necessary to lower its baryon fraction parameter (only) from 0.17 to 0.13 to obtain comparable results to Millennium and GiggleZ-MR (see Section 4).

The primary constraint for our parameter choices was the ability to reproduce the $z=0$ stellar mass function: that is, the number density of galaxies as a function of their stellar mass, $\Phi$. We show in Figure 1 that SAGE produces an excellent stellar mass function for each simulation with our fiducial parameter set, tightly matching the observational uncertainty range presented by Baldry et al. (2008). To conservatively account for the different simulation resolution limits, a minimum stellar mass equal to the median of galaxies in Millennium (sub)halos made of 50 particles is adopted. We also compare the SAGE mass functions to the C06 model, but note that C06 used the luminosity function as its primary constraint instead.

We further used a set of secondary constraints for setting our fiducial parameter values. These constraints include the star formation rate density history (Somerville et al. 2001, our Figure 4), the Baryonic Tully-Fisher relation (Stark et al. 2009, our Figure 5), the mass-metallicity relation of galaxies (Tremonti et al. 2004, our Figure 6), and the black hole-bulge mass relation (Scott et al. 2013, our Figure 8). We detail which parameters were allowed to be varied in the model calibration in order to meet these constraints in Table 2. All of the free parameters and several of the fixed parameters have different values to C06.

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Stellar mass functions at $z=0$ for SAGE galaxies for each of the N-body simulations compared to observational data from Baldry et al. (2008) and the published C06 model, taking $h=0.73$.</figcaption></figure>

For the remainder of the paper, we walk through the different baryonic reservoirs in the model and the physics that describe how mass and energy move between them.

### 4. GAS INFALL IN HALOS

Perturbations in the primordial density field lead to gravitational instability that drives mass to continually collapse into “halo”-like structures. This process is dominated by dark matter, being the main mass component of the universe, with the baryons tending to follow rather than lead the dynamical evolution. In a cosmological N-body simulation, such halo growth is typically well measured by a halo finder. Hence, the total baryonic content of each halo, and how it changes with time, can be inferred.

In the C06 model, as in SAGE, the total baryonic content of each dark matter halo at (nearly) all times is maintained at the universal fraction, $f_b$. This requires that for each simulation
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Fraction of baryonic content to total mass for Millennium halos at $z=0$ from SAGE, binned by virial mass. Each curve represents the mean fraction of the labeled component.</figcaption></figure>

time-step, the baryons in each halo grow by $f_b M_{vir} - m_b$, where $m_b$ is the total mass of baryons present in the previous time-step. This mass difference is added to the hot gas reservoir of the system and assumed to be pristine. If halo mass were to decrease over a time-step, then ejected gas (see Section 8) or, secondarily, hot gas is removed from the system (cold gas and stars, located deep in the potential well of the halo, are always unaffected here). However, in low-mass halos and at high redshift, little to none of the baryons are expected to be hot (Birnboim & Dekel 2003; Kereš et al. 2005); a decrease in halo mass can then lead to a temporary increase in the baryon fraction, i.e., above universal.

However, unlike C06, SAGE considers the baryon fraction to be a free parameter during model calibration, which sets the baseline level of baryons within the virial radius of each halo. For the Millennium and GiggleZ-MR simulations, the original value of 0.17 provides the baryonic mass required to produce a well-calibrated model. However, we find that the Bolshoi simulation, with its order-of-magnitude higher mass resolution, needs a lower baseline fraction of baryons (assuming the same remaining model parameters), $f_b = 0.13$, to obtain comparable results.

Investigating this further, it is interesting to see how the baryon fraction and simulation resolution limit together influence which halos get gas at each epoch and when this gas turns into stars. Since halos are identified earlier in Bolshoi, and given the power-law steepness of the halo mass function, a considerable fraction of baryons accumulate in low-mass halos but are delayed in their opportunity to turn into stars by the reionization prescription (see next section). Our choice of a lower baryon fraction compensates for this accumulation by lessening the significance of such baryons once they are later able to contribute to galaxy evolution. A simple back-of-the-envelope calculation reveals that a model with $f_b = 0.17$ on Millennium results in approximately the same baryonic mass per unit volume as a model with $f_b = 0.13$ on Bolshoi.

In Figure 2 we show the baryon fraction contained within Millennium halos of a given mass at $z=0$, broken down by mass component: galaxy stars, “intra-cluster” stars, hot gas, cold gas, and ejected gas (these will each be defined in the subsequent sections). One can see the dominance of hot gas in group and cluster mass systems, with stellar mass peaking in

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Effect of reionization in SAGE on the $z=0$ stellar mass function for Bolshoi. Because of the simulation’s relatively high resolution, reionization is needed to get the low-mass end of the mass function right, while also affecting the high-mass end.</figcaption></figure>

Milky Way-sized halos, and cold gas dominating at lower mass scales.

## 5. REIONIZATION

At high redshift, the baryonic content of low-mass halos is most likely suppressed as a result of photoionization heating of the intergalactic medium (IGM) by strong feedback from the first stars. This heating acts to lower the concentration of baryons in shallow potentials (Efstathiou 1992). Observationally, this can be seen in the low abundance of local dwarf galaxies relative to the prediction that results from the $\Lambda$CDM halo mass function.<fnref n="11" />

Gnedin (2000) showed using high-resolution SLH-P³M (softened Lagrangian hydrodynamics) simulations that the effect of photoionization heating can be modeled by defining a filtering mass, $M_F$, below which the baryonic fraction in the halo is reduced relative to the universal value:

$$f_b^{\text{halo}}(z, M_{\text{vir}}) = \frac{f_b^{\text{cosmic}}}{(1 + 0.26 M_F(z)/M_{\text{vir}})^3} \quad (1)$$

Notice that the filtering mass is a function of time. In their simulations, the change was sharpest around the epoch of reionization.

Kravtsov et al. (2004) fitted an analytic model to approximate the behavior of the Gnedin (2000) filtering mass. They defined two parameters that characterize this transition: $z_0$, which marks the redshift where the first H II regions overlap, and $z_r$, which marks the time when the IGM is fully reionized. The best-fit values to the Gnedin simulation are $z_0 = 8$ and $z_r = 7$. We adopt these in our model, leaving it identical to that used in C06. See Appendix B of Kravtsov et al. (2004) for the expression of $M_F(z)$.

In Figure 3, we show how reionization can have an important and positive effect on both the low-mass and high-mass ends of the stellar mass function at $z=0$. We show this for Bolshoi, using all galaxies more massive than the median

<footnote n="11">¹¹ This has been referred to in the literature as the “missing satellites problem” (Klypin et al. 1999; Moore et al. 1999), for which baryonic physical processes, including reionization heating, are capable of solving (e.g., Zolotov et al. 2012; Brooks et al. 2013).</footnote>
stellar mass in (sub)halos made of 50 particles, which keeps us above the simulation resolution limit. Because reionization suppresses gas infall (and therefore star formation) at high redshift, satellite galaxies at $z=0$ have less stars when reionization is turned on. Instead this gas cools at later epochs, being brought into more massive halos by the satellites, and contributes directly to the higher-mass central galaxy. As such, these systems form more stars, and the high-mass end of the mass function kicks out. Higher-resolution simulations have more subhalos, especially at high redshift, in which more gas can cool, and the strength of this effect is resolution-dependent (hence the effect is stronger for Bolshoi).

## 6. THE HOT GAS HALO

The physics of infalling gas is complicated and a detailed accounting is beyond the scope of the semi-analytic methodology. Its evolution can be approximated at a level of accuracy consistent with the rest of the model by assuming that new gas added is shock heated to the virial temperature of the system as it loses its gravitational potential energy. This results in the formation of a quasistatic hot halo of baryons. However, at early times and in low-mass systems, the shock heating may be unstable, or the shock may not form at all. In this case, the gas rapidly radiates away any gained energy and instead is channeled to the center as a cold stream, the so-called “cold accretion” or “rapid cooling” mode. See Section 3.2.1 of C06 for further discussion.

In SAGE, we follow the hot-halo model described in C06. That model, based on the original work of White & Frenk (1991), assumes that halo gas initially settles into an isothermal density profile at the virial temperature of the system,

$$T_{\text{vir}} = \frac{1}{2} \frac{\bar{\mu} m_p}{k} V_{\text{vir}}^2 = 35.9 \left( \frac{V_{\text{vir}}}{\text{km s}^{-1}} \right)^2 \text{K}, \quad (2)$$

where $V_{\text{vir}}$ is the halo virial velocity and $\bar{\mu} m_p$ is the mean particle mass, with $\bar{\mu} = 0.59$ nominally. The density profile of the gas is thus

$$\rho_g(r) = \frac{m_{\text{hot}}}{4\pi R_{\text{vir}} r^2}, \quad (3)$$

where $m_{\text{hot}}$ is the total hot gas mass in the halo which extends out to the virial radius $R_{\text{vir}}$.

To calculate the rate at which gas cools out of this distribution, the similarity solutions of Bertschinger (1989) are used. We define the cooling radius, $r_{\text{cool}}$, as the radius at which the gas cooling time is equal to the dynamical time of the system, $t_{\text{cool}} = R_{\text{vir}}/V_{\text{vir}}$ (although other timescales can be used, for example the age of the system—see Lu et al. 2011). Bertschinger (1989) showed that the cooling mass flux across this radius is proportional to the mass deposition rate at the center of the halo, with the proportionality close to one.

For a parcel of gas with local density $\rho_g(r)$ and temperature $T_{\text{vir}}$, the cooling time can be approximated by the ratio of its specific thermal energy to the cooling rate per unit volume,

$$t_{\text{cool}} = \frac{3}{2} \frac{\bar{\mu} m_p k T_{\text{vir}}}{\rho_g(r) \Lambda(T_{\text{vir}}, Z)}. \quad (4)$$

Here $k$ is the Boltzmann constant and $\Lambda(T, Z)$ is the cooling function, dependent on both the gas temperature and metallicity $Z$ (metal production is described below in Section 7).

Combining Equations (3) and (4), and assuming $t_{\text{cool}}$ above, allows us to find $r_{\text{cool}}$ for each system. We then solve the continuity equation

$$\dot{m}_{\text{cool}} = 4\pi \rho_g(r_{\text{cool}}) r_{\text{cool}}^2 \dot{r}_{\text{cool}} \quad (5)$$

to determine the instantaneous cooling rate as

$$\dot{m}_{\text{cool}} = \frac{1}{2} \left( \frac{r_{\text{cool}}}{R_{\text{vir}}} \right) \left( \frac{m_{\text{hot}}}{t_{\text{cool}}} \right). \quad (6)$$

Equation (6) is valid as long as $r_{\text{cool}} < R_{\text{vir}}$, which marks when the system is in the hot-halo regime, as discussed above. Otherwise, infalling gas is in the cold-accretion regime, and we assume it deposits at the center of the halo on a free-fall timescale, $R_{\text{vir}}/V_{\text{vir}}$. In the presence of radio mode AGN heating, this cooling rate is modified with the addition of an inner heating radius. We define our new hybrid cooling/heating model further in Section 9.1.

## 7. THE DISK: COLD GAS, STAR FORMATION, AND METAL ENRICHMENT

Cooling gas that settles in the center of a halo is assumed to conserve angular momentum and spin up to form a rotationally supported disk. It is from this disk of cold gas that stars form. In the SAGE model, star formation proceeds as in C06. A more sophisticated accounting of the formation of new stars, including the tracking of atomic and molecular hydrogen to fuel star formation, will be introduced in future work.

As discussed in C06, we assume that only cold gas above a critical surface density threshold can form stars. Following the work of Kauffmann (1996), and assuming the gas is evenly distributed across the disk (note that assuming an exponential disk profile simply results in a minor adjustment to the normalization), we find a critical mass given by

$$m_{\text{crit}} = 3.8 \times 10^9 \left( \frac{V_{\text{vir}}}{200 \text{ km s}^{-1}} \right) \left( \frac{r_{\text{disk}}}{10 \text{ kpc}} \right) M_{\odot}. \quad (7)$$

Here, the disk radius is defined as

$$r_{\text{disk}} = \frac{3}{\sqrt{2}} \lambda R_{\text{vir}}, \quad (8)$$

which is thrice the disk scale length proposed by Mo et al. (1998), based on the properties of the Milky Way (van den Bergh 2000). The spin parameter of the dark halo, $\lambda$ (Bullock et al. 2001), is taken directly from the N-body simulation.

A star formation rate can now been calculated from a Kennicutt–Schmidt-type relation (Kennicutt 1998):

$$\dot{m}_{*} = \alpha_{\text{SF}} \frac{(m_{\text{cold}} - m_{\text{crit}})}{t_{\text{dyn,disk}}}, \quad (9)$$

where $m_{\text{cold}}$ is the total mass of cold gas and $\alpha_{\text{SF}}$ is the star formation efficiency. In other words, a fraction $\alpha_{\text{SF}}$ of gas above the threshold is converted into stars in a disk dynamical time $t_{\text{dyn,disk}} = r_{\text{disk}}/V_{\text{vir}}$.

This model of star formation, combined with the feedback processes described in the below sections, produces a galaxy population whose combined star formation rate density evolution is consistent with the observed universe, as shown in Figure 4. The supporting observational data were compiled by Somerville et al. (2001), which we used as a constraint for
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Average star formation rate density history produced by SAGE for each cosmological simulation compared to observations and the C06 model. Observational data were originally compiled and corrected by (Somerville et al. 2001, see their Table A2 for a complete list of references).</figcaption></figure>

SAGE. Observational constraints on star formation rate density are non-trivial to obtain, naturally leading to large uncertainties (for discussion, see Somerville et al. 2001; Springel & Hernquist 2003). For example, in SAGE the time of peak star formation shows dependence on the *N*-body simulation: $z \sim 2$ for Bolshoi, $z \sim 2.5$ for GiggleZ, and $z \sim 3$ for Millennium, yet they all agree with the broadness of the real data. Millennium exhibits a systematically higher star formation rate than the other simulations for $z \gtrsim 2$, which can be attributed to its larger $\sigma_8$ value. Despite a largely identical star formation law to C06, variations in other parts of the model in SAGE (e.g., AGN feedback) affect the star formation histories of galaxies noticeably. As seen in Figure 4, the C06 model predicted star formation rates that were too high on average at high redshift.

A well-known constraint on the properties of nearby disk galaxies is the correlation between luminosity (a proxy for stellar mass) and rotation velocity (Tully & Fisher 1977). An even tighter correlation is found if one considers the contribution of cold gas to the mass of the systems: the so-called Baryonic Tully–Fisher relation (McGaugh et al. 2000). For each of the *N*-body simulations, SAGE successfully produces Sb and Sc Hubble-type galaxies (proxied by a bulge-to-total ratio cut between 0.1 and 0.5) that closely match this relation, as shown in Figure 5. The observational backdrop shows the fitted relation from Stark et al. (2009) with their random uncertainties; i.e., we plot $\log_{10}((m_* + m_{\text{cold}})/M_{\odot}) = (3.94 \pm 0.07)\log_{10}(V_{\text{max}}/\text{km s}^{-1}) + (1.79 \pm 0.26) + 2\log_{10}(0.75/0.73)$, where the last term accounts for the assumed value of $h$. We note that we have approximated the flat rotation velocity (as used by Stark et al. 2009) of SAGE galaxies by their halo’s maximum value, $V_{\text{max}}$. We also note that the contours for Bolshoi show more galaxies at low $V_{\text{max}}$ on account of the simulation’s higher mass resolution.

As stellar populations evolve, they enrich the interstellar medium with elements heavier than hydrogen and helium. Such metal enrichment is modeled in SAGE by assuming a yield $Y$ of metals are returned for each solar mass of stars formed. We deposit these metals into the cold disk of the galaxy, but more complicated models may wish to explore direct metal enrichment of the hot halo gas and beyond (e.g., Shattow et al. 2015). Furthermore, a fraction $R$ of the mass of newly

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Baryonic Tully–Fisher relation for SAGE galaxies at $z = 0$. 2500 representative galaxies (randomly sampled within the axes) with bulge-to-total ratios between 0.1 and 0.5 are plotted from Millennium, while contours encapsulating 68% of galaxies of the same cut are shown for the other *N*-body simulations. The maximum Keplerian velocity for the subhalos, $V_{\text{max}}$, represents the rotational velocity of galaxies. Appropriately, only central galaxies were considered. The thick strip represents the scatter in the observed relation (Stark et al. 2009).</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 6.** Stellar mass–gas metallicity relationship for $z = 0$ central galaxies in SAGE. 2500 representative galaxies are plotted in the axes from Millennium, while vertical lines cover the 16th–84th percentile range for bins of stellar mass for the other *N*-body simulations. The shaded region compares the 16th–84th percentile range from observational data published by Tremonti et al. (2004), binned identically to the Bolshoi and GiggleZ-MR data.</figcaption></figure>

formed stars is recycled immediately back to the cold gas disk: the so-called “instantaneous recycling approximation” (see Cole et al. 2000). Our treatment of metals and recycling are identical to that described in C06 (and De Lucia et al. 2004 before), aside from differences in the assumed $Y$ and $R$ parameters, which are guided by a change to the Chabrier (2003) initial mass function.

Through the measurement of optical nebular emission lines supplemented by galaxy photometry, Lequeux et al. (1979) were among the first to show evidence for a correlation between total mass and gas metallicity for galaxies. Furthering this, Tremonti et al. (2004, and see references therein) showed there is a near quadratic relationship in logarithmic space for the oxygen abundance in galaxies, calculated based on the model of Charlot & Longhetti (2001), as a function of their stellar mass. In Figure 6, we show that the redshift-zero population of SAGE
galaxies is representative of this observed relation for each N-body simulation by comparing against the 16th–84th percentile range of Tremonti et al. (2004, cf. their Table 3), adjusting for a conversion from a Kroupa (2001) to Chabrier (2003) initial mass function: $m_{*,Chabrier} = (5/6)m_{*,Kroupa}$. Our choice for the value of the yield is almost entirely driven by this constraint.

The scatter seen in Figure 6 for the simulations is larger than the observational data on average (up to ~75%), although Bolshoi in particular shows strikingly good agreement. The scatter for Millennium and GiggleZ-MR are almost identical across the full mass range, which can be seen in more detail by the individual Millennium points.

## 8. SUPERNOVA FEEDBACK AND THE GALACTIC FOUNTAIN

For each episode of new star formation, very massive stars have lifetimes much shorter than the typical time resolution of the simulations on which SAGE is built. Such stars end as supernovae. Supernovae play an important role in the life-cycle of a galaxy, injecting metals (as discussed above), mass and energy into the surrounding interstellar medium. Our modeling of supernova feedback follows that used in C06.

First, we assume that supernova winds remove cold gas from the disk, which in turn acts to suppresses star formation. The rate at which disk gas is driven from the galaxy, $\dot{m}_{\text{reheated}}$, is proportional to the rate at which new stars are forming:

$$ \dot{m}_{\text{reheated}} = \epsilon_{\text{disk}} \dot{m}_{*} \quad (10) $$

The proportionality constant, $\epsilon_{\text{disk}}$, is typically referred to as a mass-loading factor.

More generally, the energy released by supernovae during the star formation episode can be approximated as

$$ \dot{E}_{\text{SN}} = \frac{1}{2} \epsilon_{\text{halo}} \dot{m}_{*} V_{\text{SN}}^2, \quad (11) $$

where $\frac{1}{2}V_{\text{SN}}^2$ is the mean energy in supernova ejecta per unit mass of stars formed, $V_{\text{SN}} = 630 \text{ km s}^{-1}$ taking the commonly accepted value, and $\epsilon_{\text{halo}}$ quantifies the efficiency with which this energy is able to reheat disk gas. If this reheated disk gas were moved to the hot halo without changing its specific energy, the total thermal energy of the hot gas would change by

$$ \dot{E}_{\text{hot}} = \frac{1}{2} \dot{m}_{\text{reheated}} V_{\text{vir}}^2 \quad (12) $$

Equations (11) and (12) allow us to determine the excess energy in the hot gas after reheating: $\dot{E}_{\text{excess}} = \dot{E}_{\text{SN}} - \dot{E}_{\text{hot}}$. When $\dot{E}_{\text{excess}} < 0$ supernovae has failed to transfer enough energy with the reheated disk gas to unbind any hot halo gas. However $\dot{E}_{\text{excess}} > 0$ signals that there is enough energy in the system to unbind some of the hot halo and maintain virial balance. Ejected hot gas is accounted for in the model by placing it in an external reservoir:

$$ \dot{m}_{\text{ejected}} = \frac{\dot{E}_{\text{excess}}}{E_{\text{hot}}} m_{\text{hot}} = \left( \epsilon_{\text{halo}} \frac{V_{\text{SN}}^2}{V_{\text{vir}}^2} - \epsilon_{\text{disk}} \right) \dot{m}_{*} \quad (13) $$

Here, $E_{\text{hot}} = \frac{1}{2} m_{\text{hot}} V_{\text{vir}}^2$ is the total thermal energy of the hot gas. No gas is ejected when the right-hand side of Equation (13) is less than zero; this signals that $\dot{E}_{\text{excess}} < 0$, as mentioned above.

Equation (13) has many desirable properties. In low-mass halos with shallow potentials, the feedback can be very destructive, with the possibility that the entire disk and halo gas be expelled from the system if the feedback (i.e., star formation) is strong enough. Conversely, in halos with virial velocities above 200 km s$^{-1}$, the potential well is sufficiently deep that no amount of feedback will remove hot gas. Such halos develop and maintain very stable hot atmospheres throughout their lives (barring mergers or other cataclysmic events).

### 8.1. Reincorporation of Ejected Gas

In a dynamically evolving universe, gas that is ejected may not stay ejected forever. In SAGE, we adopt a modified version of that used in C06 to determine ejected-gas reincorporation. Previously, C06 assumed that a fixed fraction of the ejected material returned to the hot halo over a halo dynamical time, and that this was true for all halos. As addressed by Mutch et al. (2013), this model was poorly constrained, where alterations to the reincorporation parameter could significantly alter the number of low-mass systems. In fitting SAGE for multiple simulations, we found that a better match to the data can be obtained when we allow the reincorporation rate to increase for the more massive halos, and limit it to zero for the very lowest-mass halos. We do this by assuming the mass of ejected gas reincorporated per time-step is

$$ \dot{m}_{\text{reinc}} = \left( \frac{V_{\text{vir}}}{V_{\text{crit}}} - 1 \right) \frac{m_{\text{ejected}}}{t_{\text{dyn}}}, \quad (14) $$

where $t_{\text{dyn}} = R_{\text{vir}}/V_{\text{vir}}$, and assuming $V_{\text{vir}} > V_{\text{crit}}$ ($\dot{m}_{\text{reinc}} = 0$ otherwise). Here, $V_{\text{crit}} = k_{\text{reinc}} V_{\text{esc}}$, where $V_{\text{esc}} = V_{\text{SN}}/\sqrt{2}$ is the critical halo virial velocity, above which the supernova wind velocity is sufficient for the gas to escape, and $k_{\text{reinc}}$ parameterizes how efficient this idealized process actually is (cf. Table 2). In effect, $V_{\text{crit}}$ sets the velocity scale below which no gas can reincorporate back into a hot halo. Conversely, far above this scale, ejected gas reincorporates in a fraction of the dynamical time, meaning it quickly becomes available for future cooling and hence star formation.

It has been suggested by Henriques et al. (2013) that using a virial-mass dependence rather than virial-velocity one for the reincorporation rate could be conducive to reproducing observations across multiple epochs, in particular the galaxy stellar mass function. Our initial exploration of this variation has not been quite as positive however (nor was it for White et al. 2015), and hence we do not adopt it here. In truth, this is an area for improvement for both SAGE and other popular models. How to rectify the mass evolution of galaxies is an ongoing area of research deserving of far more investigation than the scope of this paper.

Note that the ejected gas component of our model need not be physically removed from the system. It merely marks gas that is unable to cool into the disk to form stars for a period of time. Such evacuated gas and metals may still play a part in a galaxy's evolution, simply one at a later epoch in its history.

## 9. SUPERMASSIVE BLACK HOLES AND THEIR FEEDBACK

Feedback from active supermassive black holes is now known to play a critical role regulating the life-cycle of many
types of galaxies. C06 distilled a number of popular ideas about the interplay of AGNs and galaxies into a simple picture of how this feedback shapes the properties of the low-redshift galaxy population (see also Bower et al. 2006; Cattaneo et al. 2006). Kitzbichler & White (2007) extended their comparison to higher redshift, specifically looking at galaxy masses, luminosities and number counts. Other authors have since adopted similar models and developed the basic framework to further explore how AGNs evolve (e.g., Marulli et al. 2008; Somerville et al. 2008; Guo et al. 2011).

C06 broke AGNs into two general classes, loosely dubbed the “quasar mode” and “radio mode.” These two classes can be distinguished by their triggering mechanism, lifetime, and accretion rate. We maintain this two-mode approach in SAGE but have updated both, for which we describe each in turn below.

We note that we do not consider the AGN luminosity function in the calibration or results of the model. This would require additional levels of modeling beyond the scope of this paper. We do, however, plan to include more complex AGN physics to explore more observables in later developments of the model.

## 9.1. The Radio Mode

Radio mode feedback was introduced into semi-analytic models to solve the cooling flow problem, where the over-accretion of cooling gas onto the central galaxy led to massive galaxy properties that were inconsistent with observations (too massive, too blue, too disky). The first implementations were simple and framed either phenomenologically, which attempted to infer a black hole accretion rate (and hence feedback) based on the local black hole and gas properties (C06), or in terms of a sharp cutoff in cooling when the halo or galaxy evolved past a chosen critical state (Bower et al. 2006), e.g., a mass threshold.

The model we employ in SAGE is an enhancement of the Bondi–Hoyle accretion model described in Section 5.2 of C06, which is also used by Somerville et al. (2008). In this model, hot gas accretes onto the central black hole at a rate approximated using the Bondi–Hoyle formula (Bondi 1952):

$$ \dot{m}_{\text{Bondi}} = 2.5\pi G^2 \frac{m_{\text{BH}}^2 \rho_0}{c_s^3}. \qquad (15) $$

Here, the sound speed, $c_s$, is approximated by the virial velocity of the parent halo, $V_{vir}$, while $\rho_0$ is the density of hot gas around the black hole. To find $\rho_0$, C06 equated the sound travel time across a shell of diameter twice the Bondi radius ($r_{\text{Bondi}} \equiv 2Gm_{\text{BH}}/c_s^2$) to the local cooling time (Equation (4)), the so-called “maximal cooling flow” model of Nulsen & Fabian (2000). Inserting this density back into Equation (15) allows us to write the accretion rate as a function of local temperature and black hole mass alone:

$$ \dot{m}_{\text{BH,R}} = \kappa_R \frac{15}{16} \pi G \bar{\mu} m_p \frac{kT}{\Lambda} m_{\text{BH}}. \qquad (16) $$

In a departure from C06, we insert a “radio mode efficiency” parameter, $\kappa_R$, to the right-hand side of Equation (16). While we have added this term by hand, it allows us to counteract the approximations used in the derivation of Equation (16) and to modulate the strength of black hole accretion (and subsequently

radio mode feedback) within SAGE. As such, the “best” value need not be 1. In the present work, we employ a default value of 0.08 (see below).

The accretion rate given by Equation (16) enables us to estimate the luminosity of the black hole in the radio mode, $L_{\text{BH,R}} = \eta \dot{m}_{\text{BH,R}} c^2$, where $\eta = 0.1$ is the standard$^{12}$ efficiency with which inertial mass is liberated upon approaching the event horizon, and c the speed of light. We assume this luminosity acts as a source of heating that offsets the energy losses from the cooling gas. If enough energy from the central AGN is injected into the cooling flow, it can be turned off entirely, leading to longer-term quenching, which was the focus of the work by C06. If we define the specific energy of gas in the hot halo as $\frac{1}{2}V_{vir}^2$, and heating as adding this amount of energy per unit mass to offset cooling, then the heating rate from radio mode feedback can be written as

$$ \dot{m}_{\text{heat}} = \frac{L_{\text{BH,R}}}{\frac{1}{2}V_{\text{vir}}^2}. \qquad (17) $$

### 9.1.1. A More Self-consistent Treatment of the Cooling–Heating Cycle

In C06, the heating rate given by Equation (17) was subtracted off the cooling rate given by Equation (6) to determine an effective AGN-adjusted cooling rate. However, despite its ability to reproduce many key local galaxy properties, there remained important limitations. Significantly, in the C06 model (and many similar radio mode models), cooling and heating are decoupled, meaning the latter only offsets the former after both have been independently calculated. In the real universe, one would expect episodes of AGN activity to have a lasting (or, at least, extended) effect on the surrounding gas, modifying its temperature and density profile in such a way that would alter the later cooling. A coupled cooling–heating model is clearly a desirable refinement if one wants to make realistic predictions for, e.g., the X-ray halo cooling luminosity or AGN jet power. Unfortunately, there is no natural way to achieve such cooling–heating coupling within the paradigm of current semi-analytic models.

To achieve a coarsely equivalent behavior, we implement a simple idea in SAGE to create an updated AGN model. We assume that cooling gas is heated by radio mode feedback out to a particular radius, called $r_{\text{heat}}$, interior to which the gas retains the memory of its past heating. This gas will never cool thereafter. To determine $r_{\text{heat}}$ we find the radius at which the energy deposited into the gas due to the radio mode equals the energy the halo gas interior to $r_{\text{heat}}$ would lose if it were to cool onto the galaxy disk from Equation (6) alone. This is given by

$$ r_{\text{heat}} = \frac{\dot{E}_{\text{heat}}}{\dot{E}_{\text{cool}}} r_{\text{cool}}, \qquad (18) $$

<footnote>$^{12}$ “Standard” in terms of what has been adopted in the literature (e.g., Di Matteo et al. 2005; Sijacki et al. 2007; Booth & Schaye 2009; Dubois et al. 2014a, 2014b; Henriques et al. 2015; Schaye et al. 2015), a value which falls between the efficiency expected for a non-spinning and maximally spinning black hole (Bardeen et al. 1972; also see Figure 1 of Maio et al. 2013). Many authors cite Shakura & Sunyaev (1973) as the justification for $\eta = 0.1$, despite their using $\eta = 0.06$. Recent hydrodynamic works have started to use $\eta = 0.2$ instead (e.g., Hirschmann et al. 2014; Sijacki et al. 2015), motivated by observations of more-luminous AGNs (Yu & Tremaine 2002; Davis & Laor 2011, anti-respectively).</footnote>
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Cooling luminosity of hot gas in halos as a function of virial temperature at $z = 0$. Here we compare what the C06 model for cooling would have returned against SAGE for Millennium halos. 2500 randomly points are plotted within the axes for both the SAGE and C06 cooling methods. Squares and starred points with error bars compare observational data of galaxy clusters from Peres et al. (1998), measured with the High Resolution Imager (HRI) and Position Sensitive Proportional Counter (PSPC) at ROSAT, respectively. Triangular points with errors show observational data of galaxy groups (Ponman et al. 1996; Bharadwaj et al. 2015). Stacked X-ray observations surrounding locally brightest galaxies from Anderson et al. (2015) are given by the horizontal marks with errors.</figcaption></figure>

where

$$ \dot{E}_{\text{heat}} = \frac{1}{2} \dot{m}_{\text{heat}} V_{\text{vir}}^2 \quad \text{and} \quad \dot{E}_{\text{cool}} = \frac{1}{2} \dot{m}_{\text{cool}} V_{\text{vir}}^2, \qquad (19) $$

and $\dot{m}_{\text{cool}}$ and $\dot{m}_{\text{heat}}$ are determined from Equations (6) and (17) above. The cooling rate given in Section 6 is then modified by the presence of this heating radius. The new cooling rate, in the presence of current and past radio mode episodes, now becomes

$$ \dot{m}'_{\text{cool}} = \left(1 - \frac{r_{\text{heat}}}{r_{\text{cool}}}\right) \dot{m}_{\text{cool}}. \qquad (20) $$

This model assumes that only the gas between $r_{\text{heat}}$ and $r_{\text{cool}}$ can cool, and when $r_{\text{heat}} > r_{\text{cool}}$, cooling is effectively quenched (i.e., $\dot{m}'_{\text{cool}} = 0$). To retain the memory of past heating episodes, the heating radius is never allowed to move inwards, only outwards.

This new model behaves almost identically to that used in C06, but with the added benefit that a more realistic cooling rate can be extracted from the cooling-heating cycle rather than simply an upper limit. To demonstrate this point, we show the cooling rates predicted by SAGE (i.e., Equation (20)) against those from the C06 model (i.e., Equation (6) less Equation (17)) for the same halos in Figure 7 at $z = 0$. The cooling rates are shown against halo virial temperature (Equation (2)). For comparison, also plotted are various X-ray observations of hot gas halos surrounding galaxy clusters and groups (Ponman et al. 1996 (P96); Peres et al. 1998 (P98); Anderson et al. 2015 (A15); Bharadwaj et al. 2015 (B15)), as marked.<fnref n="13" /> Note that P98 measure the cooling luminosity (X-ray luminosity inside

the cooling radius), whereas P96 and B15 measure the bolometric luminosity (X-ray luminosity inside $R_{500}$), for which the cooling luminosity of such systems will be lower. Furthermore, the observations statistically favor the brightest systems at any given mass scale, while the model draws more broadly from the wider galaxy and halo population. The A15 data are stacked observations of X-ray emission around “locally brightest galaxies.” These data should, in principle, be more representative of an average halo (suffering optical selection bias rather than X-ray selection bias), and therefore more directly comparable to SAGE. We have plotted the maximum of their measurement and bootstrapping errors for each binned point.

The upper-limit nature of the C06 model is highlighted in Figure 7 by the tight band of points between $10^{0.4} < T_{\text{vir}}/\text{keV} < 10^{0.6}$ relative to the observations. In this model, the radio feedback required to offset the excessive cooling demanded $\kappa_R = 1$, with correspondingly higher absolute heating rates. With less cooling to counteract, our nominal value for $\kappa_R$ has been reduced to 0.08 for SAGE. In other words, our required heating rates are significantly lower than before, the results of which can be more meaningfully compared with observations.

In truth, if an AGN has been off for an extended period of time, then one might expect the hot gas in the halo to return to its previous state and cool at closer to the maximal rate. Ideally, $r_{\text{heat}}$ would shrink during this time, which is most likely why our cooling rates lie below the observations. In this sense, our prevention of $r_{\text{heat}}$ moving back inward places a lower limit on the cooling rates of halos. The addition of $r_{\text{heat}}$ in SAGE is still a step in the right direction though, and goes some way to include often missed effects, such as the observed entropy floors in many cluster hot gas profiles. We leave a more thorough treatment of the cooling-heating cycle for a future study.

## 9.2. The Quasar Mode

In most simulations of galaxy formation, quasars are triggered by mergers or from some form of instability in the disk. The key requirement is to funnel gas into the galactic center on a very short timescale, which results in high black hole accretion rates and rapid black hole growth. Hence, the quasar mode is the dominant mode through which a black hole gets its mass.

Galaxy mergers are described below in Section 10. Of relevance here, to model the effect of mergers on black hole growth we follow the work of Kauffmann & Haehnelt (2000), as well as the enhancements in C06, and assume that mergers trigger black hole growth according to the phenomenological relation

$$ \Delta m_{\text{BH,Q}} = \frac{f'_{\text{BH}} m_{\text{cold}}}{1 + (280 \text{ km s}^{-1}/V_{\text{vir}})^2}, \qquad (21) $$

where

$$ f'_{\text{BH}} = f_{\text{BH}} (m_{\text{sat}}/m_{\text{central}}) \qquad (22) $$

is an accretion efficiency parameter with constant $f_{\text{BH}}$, which controls the fraction of cold gas accreted by a black hole and is modulated by the satellite-to-central galaxy merger mass ratio.

Disk instabilities can similarly lead to rapid black hole growth. In Section 11 we describe our instability

<footnote n="13">¹³ We have used a subset of the observational data in the literature to cover the plotted range of virial temperatures. Many other data exist (see, e.g., Figure 7 of Anderson et al. 2015 and references therein) and they are all consistent with the conclusions we draw here.</footnote>
implementation, which is similar to that used in C06. For the sake of instability-driven accretion here, we modify Equation (21) by taking $f'_{BH} = f_{BH}$ and substitute the mass of unstable cold gas for $m_{cold}$.

Although C06 used the quasar mode to grow black holes, they did not include quasar feedback when accounting for the evolution of the surrounding baryons. In SAGE, quasar mode feedback is included. This mode has little effect on the local galaxy population but can have a significant impact on early universe galaxy formation. An exploration of quasar mode feedback is left for a future paper; here we simply describe its implementation.

In the absence of a detailed understanding of how quasar accretion and feedback operates (but see Lynden-Bell 1969; Novikov & Thorne 1973; Costa et al. 2014; Stevens 2015, and references therein), we adopt a simple phenomenological model that is consistent with the quasar mode narrative. When a merger or disk instability occurs and the black hole has undergone some form of rapid accretion, we assume a quasar wind follows with luminosity $L_{BH,Q} = \eta \dot{m}_{BH,Q} c^2$, where $\eta = 0.1$ as before. This is used to calculate the total energy contained in the quasar wind,

$$E_{BH,Q} = \kappa_Q \frac{1}{2} \eta \Delta m_{BH,Q} c^2, \quad (23)$$

where $\kappa_Q$ parametrizes the efficiency with which the wind influences the surrounding gas as it escapes the galaxy and halo. Next, we calculate the total thermal energies in both the cold disk gas and hot halo gas:

$$E_{\text{cold}} = \frac{1}{2} m_{\text{cold}} V_{\text{vir}}^2 \quad \text{and} \quad E_{\text{hot}} = \frac{1}{2} m_{\text{hot}} V_{\text{vir}}^2. \quad (24)$$

Simply put, if the total energy in the quasar wind (Equation (23)) exceeds the total energy in the cold disk gas we blow out the cold gas and associated metals to the ejected gas reservoir. If the quasar energy is greater than the combined total energy in the cold gas and hot halo gas, the quasar wind instead ejects both the cold and hot gas (and metals) from the halo. This is an “all or nothing” approach that is ripe for development.

### 9.2.1. Black Hole Population

While the radio mode regulates cooling in SAGE, the quasar mode is the dominant channel for black hole growth. In Figure 8, we show that our treatment of black hole evolution is in general agreement with the observed black hole-bulge mass relation. We compare SAGE galaxies to the observed sample published in Scott et al. (2013, S13), which considers Sèrsic and core-Sèrsic galaxies (the latter with typically more-massive bulges) separately (also see Graham & Scott 2015). Once again, observational statistics favor large bulge masses, while numbers are naturally greater for low-bulge-mass systems from the theory perspective. The region of overlap spans over 1 dex in width though, and shows clear agreement.

We note that we do not consider either the quasar or radio AGN luminosity functions in the current work. This would require additional layers of modeling that are beyond the scope of this paper. We do, however, plan to include more complex AGN physics in later developments of SAGE that will allow us to explore more observables.

<figure><img src="image_10.png" /><figcaption>Figure 8. black hole–bulge mass relation for SAGE galaxies. 2500 representative galaxies are plotted for Millennium within the axes, while contours encapsulating 68% of systems containing black holes and bulges are shown for the other N-body simulations. Error bars compare observational data from Scott et al. (2013), where cyan represents Sèrsic galaxies, and purple represents core-Sèrsic.</figcaption></figure>

## 10. GALAXY MERGERS AND INTRA-CLUSTER STARS

The new SAGE model treats satellite galaxies somewhat differently to the previous C06 model. In C06, satellites had their hot halo instantly stripped upon infall, and their orbits were followed using the host subhalo position until the subhalo dark matter mass stripped below the resolution limit of the simulation. Upon losing the subhalo, an analytic merger time was calculated assuming the dynamical friction model of Binney & Tremaine (1987) and the properties of the subhalo in the time-step before it was lost:

$$t_{\text{friction}} = 1.17 \frac{V_{\text{vir}} R_{\text{vir}}^2}{G m_{\text{sat}} \ln \Lambda}. \quad (25)$$

Here, $m_{sat}$ is the total mass of the subhalo/satellite system (dark matter plus baryonic), and the Coulomb logarithm $\ln \Lambda$ is approximated by $\ln(1 + M_{vir}/m_{sat})$. The “orphan” (i.e., subhalo-less) galaxy was then allowed to survive until this clock ran out, after which it was merged with the central galaxy.

A few important consequences resulted from this satellite galaxy model. In particular, it was realized that the colors of satellite galaxies were too red, as noted by many authors (C06; Weinmann et al. 2006; Guo et al. 2011). This was because the instantaneous stripping of hot gas was imposing an artificial quenching mechanism on satellites, leading to premature suppression of star formation (also see Font et al. 2008). Furthermore, in the C06 model, all satellites were assumed to merge with the central galaxy once their dynamical friction clock reached zero, whereas we know from observations that many satellites are shredded to pieces well before merging, and instead instead become part of an intra-group or intra-cluster stellar halo. While scatter in a mass-dependent dynamical friction formula is inevitable, hydrodynamic simulations have since suggested there is a better generic fit than Equation (25) (Jiang et al. 2008, 2010).

In SAGE, we allow for the evolution of the satellite population by treating them more like central galaxies. Hot-halo stripping now happens in proportion to the dark matter subhalo stripping, rather than instantaneously. Any hot gas present in the subhalo
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Fraction of quiescent (i.e., low specific star formation rate) satellites at $z=0$ in SAGE as a function of stellar mass. Compared are observed galaxies at $z<0.05$ from SDSS along with the C06 model. Too many quiescent galaxies are consistently seen in the models for $m_* \lesssim 10^{10} M_\odot$.</figcaption></figure>

is allowed to cool onto the satellite in the usual way. Upon infall, a merger time is calculated for the satellite, using the same dynamical friction formula as above (Equation (25)). We take this as the *average* merger time expected for systems of similar properties. We then follow the satellite with time and measure the ratio of subhalo-to-baryonic mass. When this ratio falls below a critical threshold, $f_{\text{friction}}$ (typically 1: cf. Table 2), we compare its current survival time with the average time determined at infall. If the subhalo has survived longer than average, then we say that the subhalo/satellite system was more bound than average and merge it with the central in the standard way. On the other hand, if the subhalo/satellite mass ratio has fallen below the threshold sooner than average, then we argue that the system was instead loosely bound and more susceptible to disruption. In this case, we add the satellite stars to a new intra-cluster stars component, and any remaining gas goes to the parent hot halo. This omission of orphan galaxies is one of SAGE's notable points of difference to other semi-analytic models currently in the literature.

Once the occurrence of a galaxy–galaxy merger has been identified, we check the satellite-to-central baryonic mass ratio. If the ratio is above a threshold $f_{\text{major}}$–in SAGE set at a default of 0.3—we say the merger is “major.” In a major merger the disks of both galaxies are destroyed and all stars are combined to form a spheroid. Otherwise the merger is “minor,” and only the satellite stars are added to the central galaxy bulge. Furthermore, any cold gas present in either system can lead to a starburst, as described below in Section 12.

SAGE still overproduces the fraction of quiescent satellite galaxies, a problem shared with other semi-analytic models (e.g., Guo et al. 2011). In fact, this is true for galaxies at lower masses ($m_* \lesssim 10^{10} M_\odot$) in general. In Figure 9, we show the fraction of quiescent$^{14}$ galaxies from each N-body simulation, determined as those with specific star formation rates $< 10^{-11}$ yr$^{-1}$, as a function of stellar mass. To compare to real galaxies, we calculate a quiescent fraction (by the same definition as the model galaxies) from Sloan Digital Sky Survey (SDSS; York et al. 2000) galaxies, specifically Data Release 7 (Abazajian

et al. 2009). Stellar masses and star formation rates come from the MPA-JHU catalog,$^{15}$ where star formation rates are based on Brinchmann et al. (2004). We bin these data in stellar mass of width 0.1 dex and display Poisson errors for the quiescent fraction. We only include galaxies at $z < 0.05$. Also compared against these data is the quiescent fraction for C06 galaxies. While SAGE has improved the quiescent fraction at masses around $10^{10} M_\odot$, the problem at lower masses, i.e., in satellites, persists (as does an overproduction of star-forming galaxies at high masses). We note that because Millennium was our primary simulation for constraining the model, and GiggleZ-MR is of similar resolution and used the same halo finder, these simulations produce more promising results than Bolshoi. Due to a number of effects, including Bolshoi being more affected by reionization, baryons cool in Bolshoi halos at a later time, therefore systematically raising the specific star formation rates of galaxies at $z=0$ (also seen in Figure 4).

We acknowledge that removing orphans may have consequences for galaxy clustering (Saro et al. 2008; Guo et al. 2011), which we intend to address in a follow-up study. As it stands, however, most cosmological simulations have sufficiently high subhalo resolution to reproduce the observed clustering of the galaxy population through the halo occupation distribution model (see Berlind & Weinberg 2002), which also lacks an orphan population. An orphan-less semi-analytic model is significantly more modular and transportable between differing simulations, which is a desirable property in the SAGE codebase.

## 11. DISK INSTABILITIES

Galaxies also transform through instabilities that occur within a galaxy disk itself. Again, we follow C06: after each episode of star formation, we determine a stability criterion (Mo et al. 1998) using the disk mass, $m_{\text{disk}}$; radius, $r_{\text{disk}}$; and circular velocity, $V_c$, approximated by the maximum circular velocity of the halo. For the disk to be stable the following inequality must be met:

$$ \frac{V_c}{(G m_{\text{disk}}/r_{\text{disk}})^{1/2}} \ge 1. \qquad (26) $$

When the left side of Equation (26) is less than unity, we transfer (in proportion) enough stellar and cold gas mass to the bulge to return the disk to stability. Unstable cold gas can both grow the central black hole (as described in Section 9.2) and lead to the formation of new stars in a starburst (as described in the next section).

## 12. STARBURSTS

Starbursts mark the rapid formation of new stars, triggered as a result of a specific event. This is opposed to the more typical quiescent star formation that goes on in the galactic disk and described in Section 7. In SAGE, both galaxy mergers and disk instabilities act as triggers for starbursts.

**Mergers:** When a galaxy merger has been identified, the resulting starburst occurs in proportion to the total sum of the cold gas in the two merging galaxies. SAGE treats starbursts using the implementation of Somerville et al. (2001), where the fraction of cold gas converted to stars is

<footnote>$^{14}$ It is worth noting that a red fraction, defined by a cut in color (typically $g-r$), is not equivalent to a quiescent fraction, defined by a cut in specific star formation rate.</footnote>

<footnote>$^{15}$ http://wwwmpa.mpa-garching.mpg.de/SDSS/DR7/</footnote>
given by

$$e_{\text{burst}} = \beta_{\text{burst}} (m_{\text{sat}}/m_{\text{central}})^{\alpha_{\text{burst}}}. \quad (27)$$

The two parameters above are fixed at $\alpha_{\text{burst}} = 0.7$ and $\beta_{\text{burst}} = 0.56$, which provides a good fit to the numerical results of Cox et al. (2004) and also Mihos & Hernquist (1994, 1996) for merger mass ratios ranging from 1:10 to 1:1. For minor mergers, the new stars are added to the galactic bulge and the stability of the disk is subsequently checked. As mentioned in Section 10, for major mergers, all stars go to the spheroid, including those newly formed.

Instabilities: For instability-driven starbursts, $e_{\text{burst}}$ is taken as the fraction of cold disk gas that is unstable (Section 11), minus any gas that is accreted onto the central black hole (Section 9.2). All newly formed stars from a disk instability burst are then added to the bulge (adding them to the disk would simply leave the disk unstable).

Our starburst implementation for mergers is in contrast to the recent semi-analytic work of Padilla et al. (2014); rather than applying Equation (27), those authors instead check for an instability to drive a starburst after merging, by virtue of evolving the size of disks rather than assuming Equation (8).

## 13. DISCUSSION AND SUMMARY

We have presented our new, publicly available semi-analytic model of galaxy formation and evolution, SAGE. The model is based on that of Croton et al. (2006), but has updated many of the transfer processes between baryonic reservoirs, including gas infall, cooling, heating, and reincorporation. It further includes previously omitted quasar mode feedback for super-massive black holes and an intra-cluster star component for central galaxies. In addition, SAGE has a unique take on satellite galaxies, whereby they are assumed to be disrupted or merge when they lose their dark matter subhalo, rather than surviving as orphans.

The codebase that describes SAGE is modular and well suited to run on the halo merger trees of any cosmologically representative N-body simulation. In this paper we have shown the performance of SAGE on each of the Millennium, Bolshoi, and GiggleZ simulations using a common set of default parameters. With these, SAGE can successfully reproduce many observational properties of the local galaxy population simultaneously. Our calibrations include the redshift-zero stellar mass function, baryonic Tully-Fisher relation, stellar mass-gas metallicity relationship, black hole-bulge mass relation, and the average star formation rate density history. We find only minor differences between simulations in the predicted scaling relations. This insensitivity to halo finding and merger tree construction reflects the robust nature of the included physical prescriptions.

However, no model of the galaxy population can ever be considered complete. The physics that governs galaxy evolution can always be examined at a higher level of complexity, and in the current era of survey science our instruments continue to produce increasingly vast and rich data sets. Within these data, the properties of galaxies are being measured in increasingly refined detail. To stay current, SAGE is ripe for further development in a number of key areas.

1. More detailed modeling of gas outside of and within the halo and galaxy, including radio predictions for the large-scale distribution of atomic hydrogen, new baryonic

reservoirs such as the circumgalactic medium and warm IGM, and the different phases of gas in the galactic disk, specifically its neutral and molecular hydrogen content.

2. Improved modeling of the various stellar growth channels of disks and bulges in galaxies, their size and evolution, and how this modifies their broader predicted observed properties.

3. An improved understanding of gas and star formation in satellite galaxies; in particular, the formation and abundance of low mass galaxies at high redshift, and the discrepancies found by many models with the quiescent fraction of satellites at low redshift.

4. Expanding the model to produce predictions for the AGN population, such as radio jet luminosities and sizes, and the abundance of radio AGNs as a function of host galaxy mass and redshift.

Orthogonal to its use to explore questions of science, the philosophy behind SAGE is one of transparency and reproducibility of scientific results. By making SAGE an open and community project, we are hoping to (a) widen the accessibility of such models to more astronomers, especially students; (b) enable wider development of the science modules, which will hopefully be useful to astronomers with more specific interests; and (c) increase the scrutiny of how such models are produced, their uncertainty and limitations, and their use in the literature. In this sense, SAGE is following a similar path to that already established by scientific codes such as GADGET-2 (Springel 2005) and GALACTICUS<fnref n="16" /> (Benson 2012), as well as many others (see the Astrophysics Source Code Library<fnref n="17" /> for further examples).

The default catalogs produced by SAGE and used here are publicly available for download at Swinburne University's Theoretical Astrophysical Observatory (TAO, Bernyk et al. 2016). With TAO, users can additionally add a wide range of observational filters to produce apparent and absolute magnitudes, build custom light cones to mimic popular surveys, and create custom images of a mock galaxy population. These data sets form a solid basis for comparisons with survey data, although users may wish to locally re-run them using SAGE and tweak the parameters to further refine the model. When new simulations of significance become publicly available, and new models are built on top of them, we expect to use TAO to distribute these as well.

Having an array of simulations and models on-hand, and being able to generate new ones as needed, will enable astronomers to explore the theoretical uncertainty between different mock data sets. This is especially important for models which claim to follow a similar underlying physical narrative, but with different technical implementation. Expanding the pallet of theoretical predictions available to observers will add valuable context when using such models to compare with and interpret observational results.

D.C. acknowledges receipt of a QEII Fellowship from the Australian Research Council. T.G. acknowledges support from an Australian Research Council Super Science Fellowship, and is grateful to the LABEX Lyon Institute of Origins (ANR-10-LABX-0066) at the Université de Lyon for its financial support within the program "Investissements d'Avenir" (ANR-11-

<footnote n="16">¹⁶ https://sites.google.com/site/galacticusmodel/</footnote>

<footnote n="17">¹⁷ http://ascl.net/</footnote>
INDEX-0007) of the French government, operated by the National Research Agency (ANR).

The Millennium Simulation was carried out by the Virgo Supercomputing Consortium at the Computing Centre of the Max Plank Society in Garching. It is publicly available at http://www.mpa-garching.mpg.de/Millennium/. The Bolshoi Simulation was carried out by A. Klypin, J. Primack and S. Gottlieber at the NASA Ames Research Centre. The simulation and data products can be found at http://astronomy.nmsu.edu/aklypin/Bolshoi/. The GiggleZ simulations were conducted by G. Poole on the Green supercomputer at Swinburne University with subsequent analysis conducted with Swinburne's gSTAR cluster.

Some data used in this work were generated using Swinburne University's Theoretical Astrophysical Observatory (TAO). TAO is part of the Australian All-Sky Virtual Observatory (ASVO) and is freely accessible at https://tao.asvo.org.au.

ARHS thanks Toby Brown for sourcing the relevant SDSS data used in Figure 9 and discussion on the quiescent fraction. Funding for SDSS-III has been provided by the Alfred P. Sloan Foundation, the Participating Institutions, the National Science Foundation, and the U.S. Department of Energy Office of Science. The SDSS-III web site is http://www.sdss3.org/.

Finally, the authors would like to thank the anonymous referee, whose careful reading of this paper resulted in many valuable improvements.

## REFERENCES

Abazajian, K. N., Adelman-McCarthy, J. K., Agüeros, M. A., et al. 2009, ApJ, 182, 543

Anderson, M. E., Gaspari, M., White, S. D. M., Wang, W., & Dai, X. 2015, MNRAS, 449, 3806

Baldry, I. K., Glazebrook, K., & Driver, S. P. 2008, MNRAS, 388, 945

Bardeen, J. M., Press, W. H., & Teukolsky, S. A. 1972, ApJ, 178, 347

Baugh, C. M. 2006, RPh, 69, 3101

Behroozi, P. S., Wechsler, R. H., & Wu, H.-Y. 2013a, ApJ, 762, 109

Behroozi, P. S., Wechsler, R. H., Wu, H.-Y., et al. 2013b, ApJ, 763, 18

Benson, A. J. 2012, NewA, 17, 175

Berlind, A. A., & Weinberg, D. H. 2002, ApJ, 575, 587

Bernyk, M., Croton, D. J., Tonini, C., et al. 2016, arXiv:1403.5270

Bertschinger, E. 1989, ApJ, 340, 666

Bharadwaj, V., Reiprich, T. H., Lovisari, L., & Eckmiller, H. J. 2015, A&A, 573, A75

Binney, J., & Tremaine, S. 1987, Galactic Dynamics (Princeton, NJ: Princeton Univ. Press)

Birnboim, Y., & Dekel, A. 2003, MNRAS, 345, 349

Bond, J. R., Cole, S., Efstathiou, G., & Kaiser, N. 1991, ApJ, 379, 440

Bondi, H. 1952, MNRAS, 112, 195

Booth, C. M., & Schaye, J. 2009, MNRAS, 398, 53

Bower, R. G. 1991, MNRAS, 248, 332

Bower, R. G., Benson, A. J., Malbon, R., et al. 2006, MNRAS, 370, 645

Brinchmann, J., Charlot, S., White, S. D. M., et al. 2004, MNRAS, 351, 1151

Brooks, A. M., Kuhlen, M., Zolotov, A., & Hooper, D. 2013, ApJ, 765, 22

Bullock, J. S., Dekel, A., Kolatt, T. S., et al. 2001, ApJ, 555, 240

Carlberg, R. G., Couchman, H. M. P., & Thomas, P. A. 1990, ApJL, 352, L29

Cattaneo, A., Dekel, A., Devriendt, J., Guiderdoni, B., & Blaizot, J. 2006, MNRAS, 370, 1651

Chabrier, G. 2003, PASP, 115, 763

Charlot, S., & Longhetti, M. 2001, MNRAS, 323, 887

Cole, S., Lacey, C. G., Baugh, C. M., & Frenk, C. S. 2000, MNRAS, 319, 168

Colless, M., Dalton, G., Maddox, S., et al. 2001, MNRAS, 328, 1039

Costa, T., Sijacki, D., & Haehnelt, M. G. 2014, MNRAS, 444, 2355

Cox, T. J., Primack, J., Jonsson, P., & Somerville, R. S. 2004, ApJL, 607, L87

Croton, D. J., Springel, V., White, S. D. M., et al. 2006, MNRAS, 365, 11

Davis, M., Efstathiou, G., Frenk, C. S., & White, S. D. M. 1985, ApJ, 292, 371

Davis, S. W., & Laor, A. 2011, ApJ, 728, 98

De Lucia, G., Kauffmann, G., & White, S. D. M. 2004, MNRAS, 349, 1101

Di Matteo, T., Springel, V., & Hernquist, L. 2005, Natur, 433, 604

Drinkwater, M. J., Jurek, R. J., Blake, C., et al. 2010, MNRAS, 401, 1429

Dubois, Y., Volonteri, M., & Silk, J. 2014a, MNRAS, 440, 1590

Dubois, Y., Pichon, C., Welker, C., et al. 2014b, MNRAS, 444, 1453

Dunkley, J., Komatsu, E., Nolta, M. R., et al. 2009, ApJS, 180, 306

Efstathiou, G. 1992, MNRAS, 256, 43P

Font, A. S., Bower, R. G., McCarthy, I. G., et al. 2008, MNRAS, 389, 1619

Gnedin, N. Y. 2000, ApJ, 542, 535

Graham, A. W., & Scott, N. 2015, ApJ, 798, 54

Guo, Q., White, S., Boylan-Kolchin, M., et al. 2011, MNRAS, 413, 101

Hatton, S., Devriendt, J. E. G., Ninin, S., et al. 2003, MNRAS, 343, 75

Henriques, B. M. B., White, S. D. M., Thomas, P. A., et al. 2013, MNRAS, 431, 3373

Henriques, B. M. B., White, S. D. M., Thomas, P. A., et al. 2015, MNRAS, 451, 2663

Hirschmann, M., Dolag, K., Saro, A., et al. 2014, MNRAS, 442, 2304

Jarosik, N., Bennett, C. L., Dunkley, J., et al. 2011, ApJS, 192, 14

Jiang, C. Y., Jing, Y. P., Faltenbacher, A., Lin, W. P., & Li, C. 2008, ApJ, 675, 1095

Jiang, C. Y., Jing, Y. P., & Lin, W. P. 2010, A&A, 510, A60

Kauffmann, G. 1996, MNRAS, 281, 475

Kauffmann, G., Colberg, J. M., Diaferio, A., & White, S. D. M. 1999, MNRAS, 303, 188

Kauffmann, G., & Haehnelt, M. 2000, MNRAS, 311, 576

Kauffmann, G., White, S. D. M., & Guiderdoni, B. 1993, MNRAS, 264, 201

Kennicutt, R. C. 1998, ApJ, 498, 541

Kereš, D., Katz, N., Weinberg, D. H., & Davé, R. 2005, MNRAS, 363, 2

Khandai, N., Di Matteo, T., Croft, R., et al. 2015, MNRAS, 450, 1349

Kim, J., Park, C., Rossi, G., Lee, S. M., & Gott, J. R., III 2011, JKAS, 44, 217

Kitzbichler, M. G., & White, S. D. M. 2007, MNRAS, 376, 2

Klypin, A., Kravtsov, A. V., Valenzuela, O., & Prada, F. 1999, ApJ, 522, 82

Klypin, A. A., Trujillo-Gomez, S., & Primack, J. 2011, ApJ, 740, 102

Knebe, A., Knollmann, S. R., Muldrew, S. I., et al. 2011, MNRAS, 415, 2293

Knebe, A., Pearce, F. R., Lux, H., et al. 2013, MNRAS, 435, 1618

Komatsu, E., Dunkley, J., Nolta, M. R., et al. 2009, ApJS, 180, 330

Kravtsov, A. V., Gnedin, O. Y., & Klypin, A. A. 2004, ApJ, 609, 482

Kravtsov, A. V., Klypin, A. A., & Khokhlov, A. M. 1997, ApJS, 111, 73

Kroupa, P. 2001, MNRAS, 322, 231

Lee, J., Yi, S. K., Elahi, P. J., et al. 2014, MNRAS, 445, 4197

Lequeux, J., Peimbert, M., Rayo, J. F., Serrano, A., & Torres-Peimbert, S. 1979, A&A, 80, 155

Lu, Y., Kereš, D., Katz, N., et al. 2011, MNRAS, 416, 660

Lynden-Bell, D. 1969, Natur, 223, 690

Maio, U., Dotti, M., Petkova, M., Perego, A., & Volonteri, M. 2013, ApJ, 767, 37

Marulli, F., Bonoli, S., Branchini, E., Moscardini, L., & Springel, V. 2008, MNRAS, 385, 1846

McGaugh, S. S., Schombert, J. M., Bothun, G. D., & de Blok, W. J. G. 2000, ApJL, 533, L99

Mihos, J. C., & Hernquist, L. 1994, ApJL, 425, L13

Mihos, J. C., & Hernquist, L. 1996, ApJ, 464, 641

Mo, H. J., Mao, S., & White, S. D. M. 1998, MNRAS, 295, 319

Moore, B., Ghigna, S., Governato, F., et al. 1999, ApJL, 524, L19

Mutch, S. J., Poole, G. B., & Croton, D. J. 2013, MNRAS, 428, 2001

Novikov, I. D., & Thorne, K. S. 1973, in Black Holes (Les Astres Occlus), ed. C. Dewitt, & B. S. Dewitt (Paris: Gordon & Breach)

Nulsen, P. E. J., & Fabian, A. C. 2000, MNRAS, 311, 346

Padilla, N. D., Salazar-Albornoz, S., Contreras, S., Cora, S. A., & Ruiz, A. N. 2014, MNRAS, 443, 2801

Peres, C. B., Fabian, A. C., Edge, A. C., et al. 1998, MNRAS, 298, 416

Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2015, preprint (arXiv:1502.01589)

Ponman, T. J., Bourner, P. D. J., Ebeling, H., & Böhringer, H. 1996, MNRAS, 283, 690

Poole, G. B., Blake, C., Marín, F. A., et al. 2015, MNRAS, 449, 1454

Press, W. H., & Schechter, P. 1974, ApJ, 187, 425

Saro, A., De Lucia, G., Dolag, K., & Borgani, S. 2008, MNRAS, 391, 565

Schaye, J., Crain, R. A., Bower, R. G., et al. 2015, MNRAS, 446, 521

Scott, N., Graham, A. W., & Schombert, J. 2013, ApJ, 768, 76

Shakura, N. I., & Sunyaev, R. A. 1973, A&A, 24, 337

Shattow, G. M., Croton, D. J., & Bibiano, A. 2015, MNRAS, 450, 2306

Sijacki, D., Springel, V., Di Matteo, T., & Hernquist, L. 2007, MNRAS, 380, 877

Sijacki, D., Vogelsberger, M., Genel, S., et al. 2015, MNRAS, 452, 575

Skillman, S. W., Warren, M. S., Turk, M. J., et al. 2014, arXiv:1407.2600
Somerville, R. S., Hopkins, P. F., Cox, T. J., Robertson, B. E., & Hernquist, L. 2008, MNRAS, 391, 481

Somerville, R. S., Primack, J. R., & Faber, S. M. 2001, MNRAS, 320, 504

Spergel, D. N., Verde, L., Peiris, H. V., et al. 2003, ApJS, 148, 175

Springel, V. 2005, MNRAS, 364, 1105

Springel, V., & Hernquist, L. 2003, MNRAS, 339, 289

Springel, V., White, S. D. M., Jenkins, A., et al. 2005, Natur, 435, 629

Springel, V., White, S. D. M., Tormen, G., & Kauffmann, G. 2001, MNRAS, 328, 726

Srisawat, C., Knebe, A., Pearce, F. R., et al. 2013, MNRAS, 436, 150

Stark, D. V., McGaugh, S. S., & Swaters, R. A. 2009, AJ, 138, 392

Stevens, A. R. H. 2015, PASA, 32, 30

Tremonti, C. A., Heckman, T. M., Kauffmann, G., et al. 2004, ApJ, 613, 898

Tully, R. B., & Fisher, J. R. 1977, A&A, 54, 661

van den Bergh, S. 2000, The Galaxies of the Local Group (Cambridge: Cambridge Univ. Press)

Vogelsberger, M., Genel, S., Springel, V., et al. 2014, Natur, 509, 177

Weinmann, S. M., van den Bosch, F. C., Yang, X., et al. 2006, MNRAS, 372, 1161

White, C. E., Somerville, R. S., & Ferguson, H. C. 2015, ApJ, 799, 201

White, S. D. M., & Frenk, C. S. 1991, ApJ, 379, 52

White, S. D. M., & Rees, M. J. 1978, MNRAS, 183, 341

York, D. G., Adelman, J., Anderson, J. E., et al. 2000, AJ, 120, 1579

Yu, Q., & Tremaine, S. 2002, MNRAS, 335, 965

Zolotov, A., Brooks, A. M., Willman, B., et al. 2012, ApJ, 761, 71