SEMI-ANALYTIC GALAXY EVOLUTION (SAGE): MODEL CALIBRATION AND BASIC RESULTS

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 underly 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 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).⁶ 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.⁷

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.

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.

Table 1.  Details of the N-body Simulations Used in the Analysis of sage in This Paper

Simulation	Npart	Mparth	lboxh	ΩM	σ8	Code	Subhalo Finder	Tree Constructor
		(${M}_{\odot }$)	(cMpc)
Millennium	21603	8.60 × 108	500	0.250	0.900	gadget-2	subfind	l-halotree
Bolshoi	20483	1.35 × 108	250	0.270	0.820	art	rockstar	consistent-trees
GiggleZ-MR	5203	9.50 × 108	125	0.273	0.812	gadget-2	subfind	G. B. Poole et al. (2016, in preparation)

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, Ω_M (for which the equivalent for dark energy is Ω_Λ = 1 − Ω_M); the redshift-zero extrapolation of the root-mean-square linear mass fluctuation within a sphere of radius 8 h⁻¹ Mpc, σ₈; the code the simulation was run with; the code used to identify subhalos; and the code used to build the merger trees.

Download table as:  ASCIITypeset image

2.1. The Millennium Simulation

2.2. The Bolshoi Simulation

2.3. The GiggleZ Simulation Suite

2.4. Halo Merger Tree Structure and Required Properties

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.

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).

Table 2.  Fiducial sage Parameters used Throughout This Text, and Also Compared to Those From C06

Parameter	Description	Value	C06 Value	Fixed	Section(s)
${f}_{b}^{{\rm{(cosmic)}}}$	(Cosmic) baryon fraction	0.17, 0.13	0.17	No	4, 5
z0	Redshift when H ii regions overlap	8.0	8.0	Yes	5
zr	Redshift when the intergalactic medium is fully reionized	7.0	7.0	Yes	5
αSF	Star formation efficiency	0.05	0.07	No	7
Y	Yield of metals from new stars	0.025	0.03	No	7
${ \mathcal R }$	Instantaneous recycling fraction	0.43	0.30	Yes	7, 8
disk	Mass-loading factor due to supernovae	3.0	3.5	No	8
halo	Efficiency of supernovae to unbind gas from the hot halo	0.3	0.35	No	8
kreinc	Sets velocity scale for gas reincorporation	0.15	N/A	Yes	8
κR	Radio mode feedback efficiency	0.08	1.0	No	9.1
κQ	Quasar mode feedback efficiency	0.005	N/A	No	9.2
fBH	Rate of black hole growth during quasar mode	0.015	0.03	No	9.2
ffriction	Threshold subhalo-to-baryonic mass for satellite disruption or merging	1.0	N/A	Yes	10
fmajor	Threshold mass ratio for merger to be major	0.3	0.3	Yes	10

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.

Download table as:  ASCIITypeset image

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, Φ. 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.

Zoom In Zoom Out Reset image size

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.

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.

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 time-step, the baryons in each halo grow by ${f}_{{\rm{b}}}{M}_{{\rm{vir}}}-{m}_{{\rm{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 Milky Way-sized halos, and cold gas dominating at lower mass scales.

Zoom In Zoom Out Reset image size

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 ΛCDM halo mass function.¹¹

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:

$$\begin{eqnarray}&&{f}_{{\rm{b}}}^{{\rm{halo}}}(z,{M}_{{\rm{vir}}})=\displaystyle \frac{{f}_{{\rm{b}}}^{{\rm{cosmic}}}}{{(1+0.26{M}_{{\rm{F}}}(z)/{M}_{{\rm{vir}}})}^{3}}.\end{eqnarray} \tag{ 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₀, 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₀ = 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 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).

Zoom In Zoom Out Reset image size

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,

$$\begin{eqnarray}&&{T}_{{\rm{vir}}}=\displaystyle \frac{1}{2}\;\displaystyle \frac{\bar{\mu }{m}_{p}}{k}\;{V}_{{\rm{vir}}}^{2}=35.9{\left(\displaystyle \frac{{V}_{\mathrm{vir}}}{\mathrm{km}{{\rm{s}}}^{-1}}\right)}^{2}\;{\rm{K}},\end{eqnarray} \tag{ 2 }$$

where V_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

$$\begin{eqnarray}&&{\rho }_{g}(r)=\displaystyle \frac{{m}_{{\rm{hot}}}}{4\pi {R}_{{\rm{vir}}}{r}^{2}},\end{eqnarray} \tag{ 3 }$$

where m_hot is the total hot gas mass in the halo which extends out to the virial radius R_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_cool, as the radius at which the gas cooling time is equal to the dynamical time of the system, t_cool = R_vir/V_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 ρ_g(r) and temperature T_vir, the cooling time can be approximated by the ratio of its specific thermal energy to the cooling rate per unit volume,

$$\begin{eqnarray}&&{t}_{{\rm{cool}}}=\displaystyle \frac{3}{2}\displaystyle \frac{\bar{\mu }\;{m}_{p}\;k\;{T}_{{\rm{vir}}}}{{\rho }_{g}(r){\rm{\Lambda }}({T}_{{\rm{vir}}},Z)}.\end{eqnarray} \tag{ 4 }$$

Here k is the Boltzmann constant and Λ (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_cool above, allows us to find r_cool for each system. We then solve the continuity equation

$$\begin{eqnarray}&&{\dot{m}}_{{\rm{cool}}}=4\pi \;{\rho }_{g}({r}_{{\rm{cool}}}){r}_{{\rm{cool}}}^{2}\;{\dot{r}}_{{\rm{cool}}}\end{eqnarray} \tag{ 5 }$$

to determine the instantaneous cooling rate as

$$\begin{eqnarray}&&{\dot{m}}_{{\rm{cool}}}=\displaystyle \frac{1}{2}\left(\displaystyle \frac{{r}_{{\rm{cool}}}}{{R}_{{\rm{vir}}}}\right)\left(\displaystyle \frac{{m}_{{\rm{hot}}}}{{t}_{{\rm{cool}}}}\right).\end{eqnarray} \tag{ 6 }$$

Equation (6) is valid as long as r_cool < R_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_vir/V_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.

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

$$\begin{eqnarray}&&{m}_{{\rm{crit}}}=3.8\times {10}^{9}\left(\displaystyle \frac{{V}_{{\rm{vir}}}}{200\;\mathrm{km}\;{{\rm{s}}}^{-1}}\right)\left(\displaystyle \frac{{r}_{{\rm{disk}}}}{10\;{\rm{kpc}}}\right){M}_{\odot }.\end{eqnarray} \tag{ 7 }$$

Here, the disk radius is defined as

$$\begin{eqnarray}&&{r}_{{\rm{disk}}}=\displaystyle \frac{3}{\sqrt{2}}\lambda {R}_{{\rm{vir}}},\end{eqnarray} \tag{ 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, λ (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):

$$\begin{eqnarray}&&{\dot{m}}_{*}={\alpha }_{{\rm{SF}}}\;\displaystyle \frac{({m}_{{\rm{cold}}}-{m}_{{\rm{crit}}})}{{t}_{\mathrm{dyn},\mathrm{disk}}},\end{eqnarray} \tag{ 9 }$$

where m_cold is the total mass of cold gas and α_SF is the star formation efficiency. In other words, a fraction α_SF of gas above the threshold is converted into stars in a disk dynamical time ${t}_{\mathrm{dyn},{\rm{disk}}}={r}_{{\rm{disk}}}/{V}_{{\rm{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 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 ∼ 2 for Bolshoi, z ∼ 2.5 for GiggleZ, and z ∼ 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 ≳ 2, which can be attributed to its larger σ₈ 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.

Zoom In Zoom Out Reset image size

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 ${\mathrm{log}}_{10}(({m}_{*}+{m}_{{\rm{cold}}})/{M}_{\odot })$ = $(3.94\pm 0.07){\mathrm{log}}_{10}({V}_{\mathrm{max}}/\mathrm{km}\;{{\rm{s}}}^{-1})$ + $(1.79\pm 0.26)+2{\mathrm{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_max. We also note that the contours for Bolshoi show more galaxies at low V_max on account of the simulation's higher mass resolution.

Zoom In Zoom Out Reset image size

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 ${ \mathcal R }$ of the mass of newly 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 ${ \mathcal 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}_{*,\mathrm{Chabrier}}=(5/6){m}_{*,\mathrm{Kroupa}}$. Our choice for the value of the yield is almost entirely driven by this constraint.

Zoom In Zoom Out Reset image size

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.

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}}_{{\rm{reheated}}}$, is proportional to the rate at which new stars are forming:

$$\begin{eqnarray}&&{\dot{m}}_{{\rm{reheated}}}={\epsilon }_{{\rm{disk}}}\;{\dot{m}}_{*}.\end{eqnarray} \tag{ 10 }$$

The proportionality constant, _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

$$\begin{eqnarray}&&{\dot{E}}_{{\rm{SN}}}=\displaystyle \frac{1}{2}\;{\epsilon }_{{\rm{halo}}}\;{\dot{m}}_{*}\;{V}_{{\rm{SN}}}^{2},\end{eqnarray} \tag{ 11 }$$

where $\frac{1}{2}{V}_{{\rm{SN}}}^{2}$ is the mean energy in supernova ejecta per unit mass of stars formed, V_SN = 630 km s⁻¹ taking the commonly accepted value, and _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

$$\begin{eqnarray}&&{\dot{E}}_{{\rm{hot}}}=\displaystyle \frac{1}{2}\;{\dot{m}}_{{\rm{reheated}}}{V}_{{\rm{vir}}}^{2}.\end{eqnarray} \tag{ 12 }$$

Equations (11) and (12) allow us to determine the excess energy in the hot gas after reheating: ${\dot{E}}_{{\rm{excess}}}={\dot{E}}_{{\rm{SN}}}-{\dot{E}}_{{\rm{hot}}}$. When ${\dot{E}}_{{\rm{excess}}}\;\lt 0$ supernovae has failed to transfer enough energy with the reheated disk gas to unbind any hot halo gas. However ${\dot{E}}_{{\rm{excess}}}\;\gt 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:

$$\begin{eqnarray}&&{\dot{m}}_{{\rm{ejected}}}=\displaystyle \frac{{\dot{E}}_{{\rm{excess}}}}{{E}_{{\rm{hot}}}}{m}_{{\rm{hot}}}=\ \left({\epsilon }_{{\rm{halo}}}\displaystyle \frac{{V}_{{\rm{SN}}}^{2}}{{V}_{{\rm{vir}}}^{2}}-{\epsilon }_{{\rm{disk}}}\right){\dot{m}}_{*}.\end{eqnarray} \tag{ 13 }$$

Here, ${E}_{{\rm{hot}}}=\frac{1}{2}\;{m}_{{\rm{hot}}}{V}_{{\rm{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}}_{{\rm{excess}}}\lt 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⁻¹, 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

$$\begin{eqnarray}&&{\dot{m}}_{{\rm{reinc}}}=\left(\displaystyle \frac{{V}_{{\rm{vir}}}}{{V}_{{\rm{crit}}}}-1\right)\;\displaystyle \frac{{m}_{{\rm{ejected}}}}{{t}_{{\rm{dyn}}}},\end{eqnarray} \tag{ 14 }$$

where ${t}_{{\rm{dyn}}}={R}_{{\rm{vir}}}/{V}_{{\rm{vir}}}$, and assuming V_vir > V_crit (${\dot{m}}_{{\rm{reinc}}}=0$ otherwise). Here, ${V}_{{\rm{crit}}}={k}_{{\rm{reinc}}}{V}_{{\rm{esc}}}$, where ${V}_{{\rm{esc}}}={V}_{{\rm{SN}}}/\sqrt{2}$ is the critical halo virial velocity, above which the supernova wind velocity is sufficient for the gas to escape, and k_reinc parameterizes how efficient this idealized process actually is (cf. Table 2). In effect, V_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.

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):

$$\begin{eqnarray}&&{\dot{m}}_{{\rm{Bondi}}}=2.5\pi {G}^{2}\;\displaystyle \frac{{m}_{{\rm{BH}}}^{2}\;{\rho }_{0}}{{c}_{{\rm{s}}}^{3}}.\end{eqnarray} \tag{ 15 }$$

Here, the sound speed, c_s, is approximated by the virial velocity of the parent halo, V_vir, while ρ₀ is the density of hot gas around the black hole. To find ρ₀, C06 equated the sound travel time across a shell of diameter twice the Bondi radius (${r}_{{\rm{Bondi}}}\equiv 2{{Gm}}_{{\rm{BH}}}/{c}_{{\rm{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:

$$\begin{eqnarray}&&{\dot{m}}_{{\rm{BH,R}}}={\kappa }_{{\rm{R}}}\ \displaystyle \frac{15}{16}\pi G\ \bar{\mu }{m}_{p}\ \displaystyle \frac{{kT}}{{\rm{\Lambda }}}\ {m}_{{\rm{BH}}}.\end{eqnarray} \tag{ 16 }$$

In a departure from C06, we insert a "radio mode efficiency" parameter, κ_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}_{{\rm{BH,R}}}=\eta \;{\dot{m}}_{{\rm{BH,R}}}\;{c}^{2}$, where η = 0.1 is the standard¹² 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}_{{\rm{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

$$\begin{eqnarray}&&{\dot{m}}_{{\rm{heat}}}=\displaystyle \frac{{L}_{{\rm{BH,R}}}}{\frac{1}{2}{V}_{{\rm{vir}}}^{2}}.\end{eqnarray} \tag{ 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_heat, interior to which the gas retains the memory of its past heating. This gas will never cool thereafter. To determine r_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_heat would lose if it were to cool onto the galaxy disk from Equation (6) alone. This is given by

$$\begin{eqnarray}&&{r}_{{\rm{heat}}}=\displaystyle \frac{{\dot{E}}_{{\rm{heat}}}}{{\dot{E}}_{{\rm{cool}}}}\;{r}_{{\rm{cool}}},\end{eqnarray} \tag{ 18 }$$

where

$$\begin{eqnarray}&&{\dot{E}}_{{\rm{heat}}}=\displaystyle \frac{1}{2}{\dot{m}}_{{\rm{heat}}}{V}_{{\rm{vir}}}^{2}\ \ {\rm{and}}\ \ {\dot{E}}_{{\rm{cool}}}=\displaystyle \frac{1}{2}{\dot{m}}_{{\rm{cool}}}{V}_{{\rm{vir}}}^{2},\end{eqnarray} \tag{ 19 }$$

and ${\dot{m}}_{{\rm{cool}}}$ and ${\dot{m}}_{{\rm{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

$$\begin{eqnarray}&&{\dot{m}}_{{\rm{cool}}}^{\prime }=\left(1-\displaystyle \frac{{r}_{{\rm{heat}}}}{{r}_{{\rm{cool}}}}\right){\dot{m}}_{{\rm{cool}}}.\end{eqnarray} \tag{ 20 }$$

This model assumes that only the gas between r_heat and r_cool can cool, and when ${r}_{{\rm{heat}}}\gt {r}_{{\rm{cool}}}$, cooling is effectively quenched (i.e., ${\dot{m}}_{{\rm{cool}}}^{\prime }=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.¹³ 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₅₀₀), 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.

Zoom In Zoom Out Reset image size

The upper-limit nature of the C06 model is highlighted in Figure 7 by the tight band of points between 10^0.4 < T_vir/keV < 10^0.6 relative to the observations. In this model, the radio feedback required to offset the excessive cooling demanded κ_R = 1, with correspondingly higher absolute heating rates. With less cooling to counteract, our nominal value for κ_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_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_heat moving back inward places a lower limit on the cooling rates of halos. The addition of r_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

$$\begin{eqnarray}&&{\rm{\Delta }}{m}_{{\rm{BH,Q}}}=\displaystyle \frac{{f}_{{\rm{BH}}}^{\prime }\ {m}_{{\rm{cold}}}}{1+{(280\mathrm{km}{{\rm{s}}}^{-1}/{V}_{{\rm{vir}}})}^{2}},\end{eqnarray} \tag{ 21 }$$

where

$$\begin{eqnarray}&&{f}_{{\rm{BH}}}^{\prime }={f}_{{\rm{BH}}}\ ({m}_{{\rm{sat}}}/{m}_{{\rm{central}}})\end{eqnarray} \tag{ 22 }$$

is an accretion efficiency parameter with constant f_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 implementation, which is similar to that used in C06. For the sake of instability-driven accretion here, we modify Equation (21) by taking ${f}_{{\rm{BH}}}^{\prime }={f}_{{\rm{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}_{{\rm{BH,Q}}}=\eta \;{\dot{m}}_{{\rm{BH,Q}}}\;{c}^{2}$, where η = 0.1 as before. This is used to calculate the total energy contained in the quasar wind,

$$\begin{eqnarray}&&{E}_{{\rm{BH,Q}}}={\kappa }_{{\rm{Q}}}\;\displaystyle \frac{1}{2}\;\eta \;{\rm{\Delta }}{m}_{{\rm{BH,Q}}}\;{c}^{2},\end{eqnarray} \tag{ 23 }$$

where κ_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:

$$\begin{eqnarray}&&{E}_{{\rm{cold}}}=\displaystyle \frac{1}{2}\;{m}_{{\rm{cold}}}\;{V}_{{\rm{vir}}}^{2}\ \ {\rm{and}}\ \ {E}_{{\rm{hot}}}=\displaystyle \frac{1}{2}\;{m}_{{\rm{hot}}}\;{V}_{{\rm{vir}}}^{2}.\end{eqnarray} \tag{ 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.

Zoom In Zoom Out Reset image size

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.

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:

$$\begin{eqnarray}&&{t}_{{\rm{friction}}}=1.17\displaystyle \frac{{V}_{{\rm{vir}}}{R}_{{\rm{vir}}}^{2}}{{{Gm}}_{{\rm{sat}}}\mathrm{ln}{\rm{\Lambda }}}.\end{eqnarray} \tag{ 25 }$$

Here, m_sat is the total mass of the subhalo/satellite system (dark matter plus baryonic), and the Coulomb logarithm $\mathrm{ln}{\rm{\Lambda }}$ is approximated by $\mathrm{ln}(1+{M}_{{\rm{vir}}}/{m}_{{\rm{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 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_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_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¹⁴ galaxies from each N-body simulation, determined as those with specific star formation rates < 10⁻¹¹ yr⁻¹, 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,¹⁵ 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).

Zoom In Zoom Out Reset image size

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.

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_disk; radius, r_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:

$$\begin{eqnarray}&&\displaystyle \frac{{V}_{c}}{{(G{m}_{\mathrm{disk}}/{r}_{\mathrm{disk}})}^{1/2}}\geqslant 1.\end{eqnarray} \tag{ 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).

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 given by

  $$\begin{eqnarray}&&{e}_{{\rm{burst}}}={\beta }_{{\rm{burst}}}{({m}_{{\rm{sat}}}/{m}_{{\rm{central}}})}^{{\alpha }_{{\rm{burst}}}}.\end{eqnarray} \tag{ 27 }$$

  The two parameters above are fixed at α_burst = 0.7 and β_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_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).

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 supermassive 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 ¹⁶ (Benson 2012), as well as many others (see the Astrophysics Source Code Library¹⁷ 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 build 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-IDEX-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. Gottloeber 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.