THE AVERAGE STAR FORMATION HISTORIES OF GALAXIES IN DARK MATTER HALOS FROM z = 0–8

Constraining the buildup of stellar mass in galaxies provides fundamental constraints on galaxy formation models. An ever growing number of galaxy observations have been made, covering a time period from 500 Myr after the big bang (e.g., Zheng et al. 2012) to the present day. A model of galaxy formation that can match observations over this entire stretch of time would represent a significant aid to our understanding. So far, matching the evolution of the stellar masses and star formation rates (SFRs) of galaxies over this entire epoch with has proved difficult (see, e.g., Lu et al. 2012; Borgani & Kravtsov 2011; Weinmann et al. 2012). Despite challenges of reproducing these observations with detailed models of galaxy formation, significant progress has been made in recent years with empirical models that connect the evolution of galaxy properties to the evolution of dark matter halos.

Over the past decade, a range of studies have associated galaxies with dark matter halos at a given epoch, using a variety of techniques, including Halo Occupation Distribution modeling (e.g., Berlind & Weinberg 2002; Bullock et al. 2002) the Conditional Luminosity Function modeling (e.g., Yang et al. 2003), and variants of the abundance matching technique (e.g., Colín et al. 1999; Kravtsov & Klypin 1999; Neyrinck et al. 2004; Kravtsov et al. 2004; Vale & Ostriker 2004, 2006; Tasitsiomi et al. 2004; Conroy et al. 2006; Shankar et al. 2006; Berrier et al. 2006; Marín et al. 2008; Guo et al. 2010; Moster et al. 2010; Behroozi et al. 2010). The simplest abundance matching models assign the most massive observed galaxies in rank order to the largest halos in an equal simulation volume. With only slight modifications (e.g., using the peak mass for satellite halos), this technique accurately reproduces the redshift—and scale-dependent clustering of galaxies (e.g., Conroy et al. 2006; Reddick et al. 2012).

Because the cosmological model provides a prediction for the buildup of dark matter halos, one can combine knowledge of the galaxy-halo connection at different epochs with knowledge of the mass accretion and merger histories of halos to constrain the buildup of stellar mass in galaxies over cosmic time. This was first done in a comprehensive way by Conroy & Wechsler (2009), who provided an empirical constraint on the star formation histories (SFHs) and stellar mass growth in galaxies from z = 2 to the present. This approach has been explored in several studies using variants of the techniques (Zheng et al. 2007b; White et al. 2007; Firmani & Avila-Reese 2010; Leitner 2012; Béthermin et al. 2012; Wang et al. 2013; Moster et al. 2013), which together provide important constraints on the basic picture for the buildup of stars in galaxies and their connection to dark matter halos.

These studies represented important advances, albeit with shortcomings. For example, most of these studies only modeled galaxy evolution from z ∼ 2 to the present due to perceived conflicts between the integrated cosmic star formation and the cosmic stellar mass density at z > 1 (Wilkins et al. 2008; Hopkins & Beacom 2006). In addition, because passive stellar mass loss (e.g., supernovae of massive stars) depends on galaxies' SFHs, these studies had to make assumptions about the ages of stars already formed at z > 2. Finally, these studies paid limited attention to the vast array of observational uncertainties as well as modeling uncertainties affecting their derived constraints; indeed, many results are presented without error bars (Conroy & Wechsler 2009; Firmani & Avila-Reese 2010; Béthermin et al. 2012; Wang et al. 2013).

Here we present a new technique and a new compilation of observations, aimed squarely at resolving these issues. On the observational side, new observations of galaxies at very high redshifts (e.g., Bouwens et al. 2011, 2012; Bradley et al. 2012; McLure et al. 2011) have made it possible to place constraints on galaxy formation all the way to z ∼ 8. At lower redshifts, we present the first constraints from stellar mass functions (SMFs) based on the PRIMUS survey, as well as from a new analysis at z = 0 based on Sloan Digital Sky Survey (SDSS) and GALEX data (Moustakas et al. 2013). These are obtained from a consistent methodology from z = 0 to z = 1, and thus alleviate many concerns about matching inconsistently identified galaxies at different epochs. Using these combined data sets, we find that observations of the evolution of galaxy SFRs and stellar masses can now be reconciled (see also Bernardi et al. 2010 and Moster et al. 2013).

Our new method constrains the galaxy–halo relation using observed galaxy SMFs as well as specific star formation rates (SSFRs) and cosmic star formation rates (CSFRs). As with previous studies, we match observed galaxies to halos, but the additional information on SFRs allows us to break degeneracies and directly constrain the buildup of the intracluster light (ICL), as well as the amount of stars that can transfer from satellites to the central galaxy during mergers. We account for a number of statistical and systematic effects, including uncertainties from stellar population synthesis models, dust models, SFH models, the faint-end slope of the SMF, observational completeness, scatter between stellar mass and halo mass.

Given a parameterization for the intrinsic stellar mass–halo mass (SMHM) relationship as well as for these observational effects, we use a Markov Chain Monte Carlo (MCMC) method to find the allowed posterior distribution. This is important not only for addressing concerns that our results may be biased by limited observational constraints at high redshifts, but it also gives us the power to determine which observations would best improve the resulting constraints on the relationship between stellar mass, star formation, and halo mass. The direct results of this analysis are constraints on the distribution of galaxy stellar masses as a function of halo mass and redshift. From these, many other constraints relevant to galaxy formation are derived, including the average SFR as a function of halo mass and redshift, the average SFH in galaxies at a given epoch as a function of galaxy stellar mass or host halo mass, the instantaneous baryon conversion efficiency of galaxies as a function of mass and redshift, the buildup of the ICL, and the evolution of the stellar mass to halo mass ratio for progenitors of today's galaxies, all including full uncertainties and covering a redshift range from z = 0 to z = 8.

We provide a broad overview of our methodology in Section 2, including our parameterization of the SMHM relation and the relevant uncertainties, with additional details in Appendices A and B. We discuss the observational data sets relevant to our method in Section 3, with special attention to new measurements of CSFRs. We discuss the simulations that we use in Section 4, along with recalibrations of the halo mass function (Appendix G), halo mass accretion rates (Appendix H), and subhalo merger/disruption rates (Appendix I). We present our main results in Section 5, with discussion in Section 6 and a summary of our conclusions in Section 7.

Throughout this work, we assume a Chabrier (2003) initial mass function (IMF), the stellar population synthesis model of Bruzual & Charlot (2003), and the dust model employed in Blanton & Roweis (2007). We convert data sets from other papers to these models as necessary. We additionally assume a flat, ΛCDM cosmology with parameters Ω_M = 0.27, $\Omega _\Lambda = 0.73$, h = 0.7, n_s = 0.95, and σ₈ = 0.82.

2.1. Overview

Much previous work has gone into determining the relation between halo mass and stellar mass as a function of redshift (e.g., Moster et al. 2010, 2013; Yang et al. 2009a, 2012; Leauthaud et al. 2012; Behroozi et al. 2010; Guo et al. 2010; Conroy & Wechsler 2009; Wang & Jing 2010; Zheng et al. 2007b; Mandelbaum et al. 2006; Gavazzi et al. 2007; Hansen et al. 2009; Lin & Mohr 2004). The redshift evolution of this relation is due to three contributing physical effects for stars in galaxies (new star formation, merging satellite galaxies, and stellar mass loss) and one effect for halos (continuing mass accretion). Merging galaxies and mass accretion are well-constrained by dark matter simulations, and stellar mass loss is well-constrained with the assumption of an IMF for stars. The most uncertain of these effects, in relation to dark matter halos, are new star formation and stars ejected during galaxy mergers. With observational constraints on the redshift evolution of the stellar mass and the SFR as well as computed constraints on satellite mass loss, we may effectively constrain both of these factors.

Schematically, the SFR for a given galaxy over a given timestep is constrained by

$$\begin{eqnarray} {\rm SFR} \cdot \Delta t & = & \hbox{Expected Stellar Mass Now}\nonumber \\ & - & \hbox{Remaining Stellar Mass from Previous Timestep}\nonumber \\ & - & (\hbox{Stellar Mass from Mergers})\nonumber \\ &&\times (1 - \hbox{Fraction Ejected}). \end{eqnarray} \tag{ 1 }$$

In this way, we can combine a specific assignment of stellar mass to halos with the evolution of the halo mass function over time, to derive SFRs for halo trajectories. Our approach is to flexibly parameterize the possibilities for the SMHM relation, M_*(M_h, z), as well as the uncertainties affecting the remaining terms in Equation (1). Then, using a MCMC method, we can determine the allowed parameter space for M_*(M_h, z) by comparing the implied SFRs and stellar mass abundances to observations across a wide range of redshifts.

Parameterizing M_*(M_h, z) separately for individual halos is beyond the scope of this work (although it will be explored in future papers). Here, we parameterize the median value of M_* for halos at a given mass and redshift (i.e., a M_*(M_h, z)) as well as the scatter in stellar mass at fixed halo mass as a function of redshift. By taking a summation of Equation (1) over halos in a specific mass and redshift bin, we can derive the total amount of stars formed in that bin. Thus, dividing by the number of halos in that bin gives the average SFR, which is similar to approaches taken in previous works (Conroy & Wechsler 2009; Leitner 2012; Wang et al. 2013; Moster et al. 2013). This approach is accurate except in terms of small second-order effects that come from different stellar mass histories having a range of stellar mass loss rates. This aspect is discussed in detail in Appendix A. Further discussion of how individual galaxy SFHs can differ from the average is presented in Section 5.8.

In this section, we present our method for parameterizing the SMHM relation (Section 2.2.1), our parameterization of the uncertainties affecting observational data (Section 2.2.2), our model for stellar mass loss (Section 2.3), and our model for ejection of stars into the ICL (Section 2.4). To streamline the presentation of this section, we present the methodology for incorporating the effects of stellar mass accreted in mergers (which involves calibrating merger rates from simulations to z = 8) and for calculating observables (such as the SMF and SSFRs) to Appendices B and C. We present a summary of the methodology and details of the MCMC method in Section 2.5.

2.2. Determining the Stellar Mass–Halo Mass Relation to z = 8

2.2.1. Intrinsic Relation

As in Behroozi et al. (2010), we parameterize the SMHM relation at a given epoch. The most commonly used function is a double power law (Yang et al. 2012; Moster et al. 2013), which has a characteristic halo mass (M₁) and stellar mass (M_{*, 0}), a low-mass slope (α), and a high-mass slope (β). Because we include data down to low stellar masses (M_* ∼ 10^7.25 M_☉) and have tight observational constraints from the SDSS, we find that a double power-law cannot accurately fit the unique shape of the SMF. As shown in Appendix D, a double power-law results in an SMF off by as much as 0.1 dex at z = 0. For this reason, we choose a different, five-parameter form as the fitting function for the SMHM relation. The chosen form retains the low-mass power-law slope (α) as well as a characteristic stellar (M_{*, 0}) and halo mass (M₁). The high-mass behavior is trickier to fit. We find the best match with a subpower (superlogarithm) function with index γ:

$$\begin{equation} \log _{10}(M_\ast (M_h \rightarrow \infty)) \propto \left(\log _{10}\left(\frac{M_h}{M_1}\right)\right)^\gamma. \end{equation} \tag{ 2 }$$

If γ is 1, this equation becomes an ordinary power law. As γ approaches 0, this equation approaches a logarithm. In between, this equation grows more slowly than any power law, but faster than any logarithm. While this equation can become multiply valued for M_h < M₁, it is straightforward to mitigate this behavior as well as allowing a smooth connection to a power-law for M_h < M₁. The specific parameterization we adopt is⁵

$$\begin{eqnarray} \log _{10}(M_\ast (M_h)) & = &\! \log _{10}(\epsilon M_1) + f\left(\log _{10}\left(\frac{M_h}{M_1}\right)\right) - f(0) \nonumber \\ f(x) & = &\! -\log _{10}(10^{\alpha x} + 1) + \delta \frac{(\log _{10}(1+\exp (x)))^\gamma }{1+\exp (10^{-x})}. \nonumber\\ \end{eqnarray} \tag{ 3 }$$

This is a power law with slope −α for M_h ≪ M₁ and a subpower law with index γ for M_h ≫ M₁. The characteristic stellar mass to halo mass ratio is at the characteristic halo mass M₁. The maximum errors of this fit to the SMHM relation are about 0.025 dex, or roughly four times better than a double power-law fit, as discussed in Appendix D.

We use the virial mass (as defined in Bryan & Norman 1998) to define the halo mass of central galaxies, and for satellites, we use the peak progenitor virial mass (M_peak). As this paper was being prepared, Reddick et al. (2012) found that the peak v_max for halos is an even better proxy for stellar mass, at least at z = 0.05, in agreement with physical expectations that the depth of the halo potential well before it is impacted by stripping should be most correlated with the galaxy mass. The use of M_peak is at present easier both conceptually (to consider mass ratios between galaxies and halos) and operationally (because it allows the use of several previously calibrated relationships in the literature), so we stick with this choice in the present work for simplicity. The use of M_peak may result in slight underestimates of the clustering and the satellite fraction (see discussion in Reddick et al. 2012 for a comparison of these quantities in the local universe). However, Section 5.6 demonstrates that the differences between our results and the best-fitting SMHM of Reddick et al. (2012) are small, and we expect that they are within the current systematic errors.

The redshift scaling of the relation poses a unique challenge. At z = 0, the parameter fits must be flexible enough to match the tight constraints from SDSS observations. At higher redshifts, the fit must be flexible enough to allow exploration of possible SFHs allowed by higher uncertainties, but not so flexible that it allows significant overfitting of the systematic and statistical errors in the data. Moreover, the scaling must be physical: obviously, negative SFRs cannot be allowed, nor can stellar mass to halo mass ratios that exceed the halo's baryon fraction. These physical constraints help ensure basic sanity for the functional form where no observational data exist (for redshifts z > 8, as well as for faint galaxies).

By necessity, the error bars we obtain will be sensitive to the choice of the redshift fit. This is not only so for the variation allowed in the stellar mass histories, but it is also true for the variation allowed in the nuisance parameters for systematic errors. Because the latter are often poorly constrained by available data, the flexibility of the redshift fit also influences how much of the available nuisance parameter space can be explored. For this reason, we have tried a large number of different redshift scalings, including cubic spline interpolation with many control points (which results in overfitting) and many different choices of scaling parameters with a or z—many of which have unphysical behavior at very high redshifts. We have also tried fitting directly to abundance matching results (Appendix E), but no robust trends with redshift emerged.

The choice that we settled on has several conventions. For many variables in Equation (3), we have three parameters: one each for the low-redshift value, one for the scaling at intermediate redshifts (0.5 ≲ z ≲ 2), and one for scaling at high redshifts (z ≳ 2). We also apply an exponential shutoff to the redshift scaling toward lower redshift to isolate the highly constrained low-redshift parameters from higher-redshift ones. Our final parameterization is the following:

$$\begin{eqnarray} \nu (a) & = & \exp (-4 a^2)\nonumber \\ \log _{10}(M_1) & = & M_{1,0} + (M_{1,a}(a-1) + M_{1,z}z)\nu \nonumber \\ \log _{10}(\epsilon) & = & \epsilon _0 + (\epsilon _a (a-1) + \epsilon _z z)\nu + \epsilon _{a,2} (a-1)\nonumber \\ \alpha & = & \alpha _0 + (\alpha _a (a-1))\nu \nonumber \\ \delta & = & \delta _0 + (\delta _a (a-1) + \delta _{z} z)\nu \nonumber \\ \gamma & = & \gamma _0 + (\gamma _a (a-1) + \gamma _{z} z)\nu. \end{eqnarray} \tag{ 4 }$$

Additionally, there will be scatter in stellar mass at fixed halo mass (ξ). While current studies have indicated that this is approximately 0.16–0.2 dex at z = 0 (More et al. 2009; Yang et al. 2009b; Reddick et al. 2012) with no evidence for a mass trend at least down to 10¹² M_☉ in halo mass (Reddick et al. 2012), little is known about the redshift evolution of this scatter. For this work, we parameterize the possible evolution of scatter with redshift via a two-parameter scaling:

$$\begin{equation} \xi = \xi _0 + \xi _a (a-1), \end{equation} \tag{ 5 }$$

and we take the prior on ξ₀ to be 0.20 ± 0.03 dex, as determined by Reddick et al. (2012) using a combination of constraints from the correlation function and conditional SMF.

2.2.2. Observational Systematics

A full discussion of the systematic uncertainties affecting SMFs may be found in Behroozi et al. (2010). These include uncertainties from the IMF, the stellar population synthesis model, the dust model, the SFH model, sample variance, Eddington bias, redshift errors, and magnification bias. Most of these effects result in a constant systematic offset (μ) in stellar masses:

$$\begin{equation} \log _{10}\left(\frac{M_{\ast,\mathrm{meas}}}{M_{\ast,\mathrm{true}}}\right) = \mu. \end{equation} \tag{ 6 }$$

However, some effects can have different behavior as a function of stellar mass, largely on account of assumptions for passive galaxies (e.g., declining SFRs and low dust fractions) not being true for active galaxies. We account for these effects by introducing an additional parameter (κ) to capture offsets in the stellar masses of active galaxies:

$$\begin{equation} \log _{10}\left(\frac{M_{\ast,\mathrm{meas,active}}}{M_{\ast,\mathrm{true,active}}}\right) = \mu + \kappa. \end{equation} \tag{ 7 }$$

For the fraction of galaxies that are quiescent, we use a fitting formula for recent measurements from Brammer et al. (2011), corrected to the stellar mass estimates we use in this paper:

$$\begin{equation} f_\mathrm{passive}(M_{\ast,\mathrm{meas}},z) = \left[\left(\frac{M_{\ast,\mathrm{meas}}}{10^{10.2+0.5z}\,M_{\odot }}\right)^{-1.3}+1\right]^{-1}. \end{equation} \tag{ 8 }$$

At low redshifts, most galaxies below 10^10.2 M_☉ in stellar mass are active; beyond z = 3, nearly all galaxies are considered active.

As in Behroozi et al. (2010), we do not model uncertainties in the IMF in this work. The most significant remaining uncertainties come from the stellar population synthesis (SPS) model (∼0.1 dex), the dust model (∼0.1 dex), and the SFH (∼0.2 dex). By definition, passive galaxies have had no recent star formation, so that the SFH uncertainties are substantially smaller. For that reason, we take the priors on μ to be a log-normal distribution centered at zero with width 0.14 dex (combined SPS and dust model errors), but the priors on κ to be of width 0.24 dex (all three sources combined). The functional form of how μ and κ evolve with redshift is unknown. We therefore adopt the following fiducial formulae:

$$\begin{equation} \mu = \mu _0 + (a-1)\mu _a \end{equation} \tag{ 9 }$$

$$\begin{equation} \kappa = \kappa _0 + (a-1)\kappa _a. \end{equation} \tag{ 10 }$$

At higher redshifts (z ∼ 2.3), Muzzin et al. (2009) suggest that uncertainties from SPS modeling become more important, on the order of 0.2 dex instead of 0.1 dex. We thus set the widths of the priors on μ_a and κ_a to be 0.22 dex and 0.30 dex, respectively, which take these increased uncertainties into account.

Because stellar populations cannot be fully constrained with limited photometric information, stellar mass estimates for individual galaxies have intrinsic scatter relative to the true galaxy stellar mass. This causes an Eddington bias (Eddington 1913, 1940) in the SMF: on the high-mass end of the SMF, many low-mass galaxies are upscattered, but there is a limited number of higher-mass galaxies that can be downscattered. This results in a net increase in the observed numbers of high stellar-mass galaxies. This effect is best estimated for each redshift range and modeling technique individually; while few authors have corrected for this error, many at least provide an estimate of its effect. As in Behroozi et al. (2010), we find that these estimates depend most significantly on the redshift. We thus model the distribution in the observed stellar mass estimates (compared to the true stellar masses) as a log-normal Gaussian with mean 0 and a redshift-dependent standard deviation given by

$$\begin{equation} \sigma (z) = \sigma _0 + \sigma _z z. \end{equation} \tag{ 11 }$$

Following Conroy et al. (2009) and Behroozi et al. (2010), we take σ₀ = 0.07 and we take the prior on σ_z to be σ_z = 0.04 ± 0.015, consistent with estimates from the literature (Conroy et al. 2009; Kajisawa et al. 2009; Pérez-González et al. 2008; Marchesini et al. 2009, 2010; Stark et al. 2009; Caputi et al. 2011; Lee et al. 2011).

Finally, there are some concerns that, due to the Lyman break techniques used to detect high-redshift galaxies, not all high-redshift galaxies may be detected. If galaxies have bursty star formation, or if some fraction of star forming galaxies have extreme quantities of dust, then not all galaxies will be detected in Lyman break surveys (see also Stark et al. 2009; Lee et al. 2009). We parameterize the stellar mass completeness of high-redshift surveys with a two-parameter fit, governing the amplitude (A) and redshift onset (z_c) of incompleteness:

$$\begin{equation} c_i(z) = 1 - \frac{A}{\exp (z_c-z)+1}. \end{equation} \tag{ 12 }$$

Because there is little reason to believe that surveys are missing a significant fraction of galaxies at z < 1, we require z_c > 0.8, and we set the completeness fraction to 1 for z < 1. This amounts to adopting the following formula for the completeness:

$$\begin{equation} c(z) = \left\lbrace {\begin{array}{@{}ll}1 & \hbox{if z < 1}\\ c_i(z) + (1-c_i(1)) & \hbox{if z > 1.}\end{array}} \right. \end{equation} \tag{ 13 }$$

We parameterize the effect of this incompleteness on the observed SFRs by a single parameter b, which sets the fraction of "bursty" versus "dusty" star formation. Fully bursty star formation (b = 1) would have no effect on the total observed CSFR, but the observed SSFRs would be boosted by the incompleteness factor compared to the SSFR averaged over longer periods of time. Fully dusty star formation (b = 0) would lower the observed CSFR by the same factor as the cosmic SM density. Presumably, a moderate amount of unobscured star formation would be enough to render a high-redshift galaxy observable—hence, dustiness will impact galaxy completeness only if star formation is largely obscured. Given that models for SFRs from observed galaxies already attempt to correct for dust obscuration, the main effect of dustiness on SSFRs is that the SSFRs for very dusty galaxies will not be observed. Therefore, in our current model, we assume that this has no effect on the average SSFR.

2.3. Star Formation Histories

Under our assumption of a Chabrier (2003) IMF and Bruzual & Charlot (2003) stellar evolution tracks, we use the FSPS package (Conroy et al. 2010; Conroy & Gunn 2010) to calculate the rate of stellar mass loss. We find the fraction of mass lost from a single stellar population as a function of the time since its formation to be well-fit by the following formula:⁶

$$\begin{equation} f_\mathrm{loss}(t) = 0.05 \ln \left(1 + \frac{t}{1.4\; \mathrm{Myr}}\right). \end{equation} \tag{ 14 }$$

2.4. Extragalactic Light and Merging Galaxies

From Equation (1), galaxies can build up stars either through mergers or internal star formation. Schematically, this implies that the unknown amount of stars deposited in mergers can be derived given separate constraints on the stellar mass growth of galaxies (i.e., from SMFs) and their SFRs:

$$\begin{equation} \Delta {\rm SM}_\mathrm{deposited} = \Delta {\rm SM}_\mathrm{central} - {\rm SFR}_\mathrm{central} \Delta t, \end{equation} \tag{ 15 }$$

where ΔSM_deposited is the amount of stellar mass deposited from merging satellites, ΔSM_central is the stellar mass growth of the central galaxy, and SFR_central is the SFR of the central galaxy over the chosen period Δt. (Note that the real equation is slightly more complicated due to stellar mass loss from passive evolution; see Appendix B).

When a galaxy merger occurs, the stars associated with the satellite galaxy may either be deposited onto the central galaxy or be ejected into the intrahalo light (IHL), also called the ICL for galaxy clusters. Some care is necessary when using these terms, because there exist two separate definitions of the ICL/IHL from simulations (i.e., stars not bound to the main galaxy) and from observations (stars not counted in the light profile of the main galaxy). In this work, we use "intracluster light" (and "ICL") to mean any remnant stars from past galaxy mergers which are not counted as part of the light profile from the main galaxy, regardless of their boundedness. We do not make a distinction between IHL and ICL; we use one term ("ICL") regardless of the size of the host halo, as this usage is by far more common in the literature.

Because we can constrain the total amount of stars in merging satellites through the halo merger rate (Appendix I) and the SMHM relation, we can also derive the amount of stellar mass deposited into the ICL (ΔICL):

$$\begin{equation} \Delta {\rm ICL} = \Delta {\rm SM}_\mathrm{incoming} - \Delta {\rm SM}_\mathrm{deposited}, \end{equation} \tag{ 16 }$$

where ΔSM_incoming is the total amount of stars in satellites that merge within Δt. Note that this quantity includes all stars ever ejected from merging galaxies. Some fraction of these ejected stars may have speeds higher than the escape velocity of the surrounding dark matter halo and will be scattered to very large distances; however, this fraction is expected to be on the order of a few percent or less (Behroozi et al. 2012).

This simple picture is made more complicated by inconsistencies in how galaxies and the ICL are separated. The total luminosity recovered and the fraction of it attributed to the galaxy are both dependent on surface-brightness limits, sky subtraction methods, and galaxy fitting methods, which are different for the different data sets we use and for different redshift ranges in the data sets themselves. Over most mass ranges and redshifts, this is not a problem: the rate of star formation in many galaxies results in an overall change in stellar mass between redshifts that is much larger than that attributable to a change in the luminosity modeling method. However, for massive galaxies (M_* > 10¹¹ M_☉) at low redshifts (z < 1), this is a more significant issue. For example, including 5% more light between z = 0.3 and z = 0.1 in a massive galaxy (0.05 mag) would result in an equal buildup in stellar mass as would a SSFR of about 2.5 × 10⁻¹¹ yr⁻¹ (the expected SSFR for such galaxies; Salim et al. 2007). As the true SSFR drops to 10⁻¹² yr⁻¹, the galaxy models must be self-consistent to about 0.2% (0.002 mag) to avoid introducing a comparable error, which is well beyond current calibration methods (Bernardi et al. 2010).

This said, changes in the mass of the galaxy because of definitional issues are completely degenerate with incoming mass from mergers in Equations (15) and (16): ΔSM_deposited is replaced by the sum of the definitional and deposited stellar mass changes (ΔSM_definitional + ΔSM_deposited) everywhere it appears. Notably, Equation (16) still gives us a way to robustly determine the buildup of stars in the ICL.

If we used Equations (16) and (15) directly, the SFR would have to be parameterized separately from the stellar mass growth. Instead, we rewrite Equation (15) in terms of the fraction of stellar mass growth that comes from the SFR:

$$\begin{equation} \Delta {\rm SM}_\mathrm{deposited} + \Delta {\rm SM}_\mathrm{definitional} = (1-f_{{\rm SFR}}) \Delta {\rm SM}_\mathrm{central} \end{equation} \tag{ 17 }$$

$$\begin{equation} f_{{\rm SFR}} = \frac{{\rm SFR}_\mathrm{central} \Delta t}{\Delta {\rm SM}_\mathrm{central}}. \end{equation} \tag{ 18 }$$

This is useful because there are strong physical priors on what f_SFR can be. Low-mass galaxies do not have enough incoming stellar mass in satellites to account for a significant part of their stellar mass growth, so f_SFR asymptotes to 100% for such galaxies. On the other hand, the highest-mass galaxies experience almost no internal star formation (Salim et al. 2007), meaning that f_SFR must approach 0 for such galaxies. Therefore, we can approximate the halo mass dependence of f_SFR using a double power law as

$$\begin{equation} f_{{\rm SFR}}(M_h) = \left[\left(\frac{M_h}{M_{h,{\rm ICL}}}\right)^\beta + 1\right]^{-1}. \end{equation} \tag{ 19 }$$

We allow M_{h, ICL} to be redshift-dependent using a two-parameter fit:

$$\begin{equation} M_{h,{\rm ICL}} = M_{h,{\rm ICL},0} + (a-1) M_{h,{\rm ICL},a}. \end{equation} \tag{ 20 }$$

Rather than add extra parameters for β, we make the assumption that f_SFR(10¹⁶ M_☉) is 1%. This is equivalent to the constraint

$$\begin{equation} \beta = \frac{\log _{10}(99)}{16 - \log _{10}(M_{h,{\rm ICL}})}. \end{equation} \tag{ 21 }$$

2.5. Methodology Summary

We use three types of observational constraints on star formation in galaxies: the SMF, the SSFR of galaxies (star formation rate per unit stellar mass), and the CSFR. SMFs are discussed in Section 3.1, CSFRs are discussed in Section 3.3, and SSFRs are discussed in Section 3.2.

To ensure a fair comparison between the different data sets, we convert each measurement to the same set of systematic assumptions. For the IMF, we convert all results to that of Chabrier (2003). We do not consider evolving IMFs, nor do we consider the uncertainties for a non-universal IMF, as both uncertainties are beyond the scope of this paper. For the stellar population synthesis and dust models, we convert all results to the models of Bruzual & Charlot (2003) and Blanton & Roweis (2007). In all cases, we convert authors' supplied units to physical units under the assumption that h = 0.7.

3.1. Stellar Mass Functions

We combine many overlapping data sets (Baldry et al. 2008; Moustakas et al. 2013; Pérez-González et al. 2008; Mortlock et al. 2011; Marchesini et al. 2009; Stark et al. 2009; Bouwens et al. 2011) to constrain the evolution of the SMF from z = 0 to z = 8. These are shown in Table 3.

Table 1. Table of Parameters

Symbol	Description	Equation	Section
M1	Characteristic halo mass	3	2.2.1
	Characteristic stellar mass to halo mass ratio	3	2.2.1
α	Faint-end slope of SMHM relation	3	2.2.1
δ	Strength of subpower law at massive end of SMHM relation	3	2.2.1
γ	Index of subpower law at massive end of SMHM relation	3	2.2.1
ν	Exponential cutoff of evolution of M*(Mh) with scale factor	4	2.2.1
ξ	Scatter in dex of true stellar mass at fixed halo mass	5	2.2.1
Mh, ICL	Characteristic halo mass at which half of stellar mass growth is due to mergers	18	2.4
μ	Systematic offset in stellar masses and SFRs	6	2.2.2
κ	Systematic offset in stellar masses for active galaxies	7	2.2.2
σ	Scatter in measured stellar mass at fixed true stellar mass	11	2.2.2
c	Galaxy detection completeness	13	2.2.2
b	Fraction of incompleteness due to burstiness (as opposed to dustiness)	C15	2.2.2
ρ	Correlation of SFR to stellar mass at fixed halo mass	C12	C.2

Download table as:  ASCIITypeset image

Table 2. Table of Priors

Symbol	Description	Equation	Prior
ξ0	Value of ξ at z = 0, in dex	5	G(0.16, 0.04)
ξa	Redshift scaling of ξ, in dex	5	G(0, 0.16)
μ0	Value of μ at z = 0, in dex	6	G(0, 0.14)
μa	Redshift scaling of μ, in dex	9	G(0, 0.22)
κ0	Value of κ at z = 0, in dex	7	G(0, 0.24)
κa	Redshift scaling of κ, in dex	10	G(0, 0.3)
σz	Redshift scaling of σ, in dex	11	G(0.05, 0.015)
A	Amplitude of galaxy detection completeness	12	U(0, 1)
zc	Onset redshift of galaxy detection completeness	12	zc > 0.8
b	Fraction of incompleteness due to burstiness (as opposed to dustiness)	⋅⋅⋅	U(0, 1)
ρ0.5	Correlation of SFR to stellar mass at fixed halo mass at a = 0.5	C12	U(0.23, 1.0)

Notes. G(x, y) denotes a Gaussian distribution with center x and width y. U(x₁, x₂) denotes a uniform distribution from x₁ to x₂. The remaining parameters have no explicit priors; M₁ and are explored in logarithmic space, whereas α, δ, and γ are explored in linear space.

Download table as:  ASCIITypeset image

Table 3. Observational Constraints on the Stellar Mass Function

Publication	Redshifts	Colors	Area	Notes
Baldry et al. (2008)	0.003–0.05	ugriz	4783 deg2	VBD
Moustakas et al. (2013)	0.05–1	UV–MIR	9 deg2	D
Pérez-González et al. (2008)	0.2–1.6	UV–MIR	0.184 deg2	ISD
Mortlock et al. (2011)	1–3	BVizH160	0.0125 deg2	VID
Marchesini et al. (2009)	1.3–4	B–MIR	0.142 deg2	ID
Marchesini et al. (2010)	3–4	UV–MIR	0.43 deg2	ID
Lee et al. (2012)	4–5	B–MIR	0.089 deg2	SD
Stark et al. (2009)	6	B–MIR	0.089 deg2	SD
Bouwens et al. (2011)	7–8	B–H160	0.0148 deg2	UD
Bradley et al. (2012)	8	MIR	0.076 deg2	UD

Notes. Letters correspond to the following corrections made to the published results: I, initial mass function; S, stellar population synthesis model; D, dust model; V, sample variance; B, surface brightness incompleteness; C, cosmology (specifically h); U, UV to stellar mass conversion according to González et al. (2011). The local results (z < 0.2) in the Moustakas et al. (2013) mass functions are taken from the Sloan Digital Sky Survey and cover an area of 2505 deg².

Download table as:  ASCIITypeset image

Different observational challenges present themselves in each redshift range. At low redshifts, the statistics are sufficient that systematic errors in calculating stellar masses are the largest source of uncertainty (Behroozi et al. 2010). At higher redshifts, this is no longer the case; statistical and sample variance errors can be of equal magnitude to estimated systematic errors. Nonetheless, under the assumption of a smoothly varying SMHM relation, our approach (which simultaneously constrains the SMHM relation against many SMFs at different redshifts) can constrain the statistical errors better than the individual error bars on SMFs might suggest. We refer the interested reader to Behroozi et al. (2010) for a discussion of the most common uncertainties affecting observation constraints on the SMF and to Section 2.2.1 for our parameterization of their effects. In this section, therefore, we limit ourselves to discussing special concerns relevant to using SMFs over this broad range of redshifts.

Questions about the accuracy of stellar mass estimates over this redshift range come from several previously published results that suggest that the evolution in the CSFR density is inconsistent with the estimated cosmic stellar mass density at z > 1 (Nagamine et al. 2006; Hopkins & Beacom 2006; Pérez-González et al. 2008; Wilkins et al. 2008). One possible explanation is an evolving or a non-universal IMF; a number of different lines of evidence exist in support of this theory (e.g., van Dokkum & Conroy 2010; Lucatello et al. 2005; Tumlinson 2007a, 2007b; van Dokkum 2008). Another possible reason is that star formation at high redshift could be bursty or dust-obscured, yielding an incomplete census of galaxies in Lyman-break surveys (Lee et al. 2009; Stark et al. 2009). On the other hand, Reddy & Steidel (2009) offer a simpler explanation. They appeal to luminosity-dependent reddening corrections in the ultraviolet luminosity functions at high redshift as well as steeper faint-end slopes for SMFs, and demonstrate that the purported discrepancy then largely vanishes.

In our analysis of recent CSFR literature (Section 3.2), we find that newer estimates of the CSFR are systematically lower at z > 3 than the data in Hopkins & Beacom (2006; see Figure 2). A similar conclusion is reached in Bernardi et al. (2010) and Moster et al. (2013). In addition, deep probes of the SMF at lower redshifts (z < 1.5) have resulted in findings of steep faint-end slopes (Baldry et al. 2008; Drory et al. 2009; Mortlock et al. 2011). We find this to be strong evidence in support of the explanation in Reddy & Steidel (2009). Nonetheless, we also consider simple models for the effects of bursty or dusty modes of star formation; an analysis of models for an evolving or non-universal IMF is beyond the scope of this paper.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The highest redshifts present unique challenges. Herein, for redshifts 7 < z < 8.5, we convert the UV luminosity function given in Bouwens et al. (2011) to an SMF according to the recipe in González et al. (2011). It has been suggested that this conversion is reasonable given the blue UV continuum slopes of such high-redshift galaxies, and follow-up observations from IRAC that imply that these early galaxies are nearly dust-free (Labbé et al. 2010). At the same time, concerns arise because the UV luminosity function presented in Bouwens et al. (2011) is steep (α ≈ −1.6 to −2.0, instead of the shallower α ≈ −1.5 to −1.7 at lower redshifts). One interpretation of this result is that assumptions for galaxy morphology (and, therefore, the ability to estimate the full UV flux or half-light radii) can bias the calculated UV luminosity function to have a steeper faint-end slope (Grazian et al. 2011). Indeed, any systematic effect that interferes more with brighter galaxies than with dimmer ones (e.g., dust) will tend to cause an overestimation of the slope of the UV luminosity function. González et al. (2011) assume that luminosity-dependent corrections for the UV-to-stellar-mass conversion do not evolve substantially from z = 4 to z = 7. This may in fact overestimate the dust corrections, an error that would lower the faint-end slope, partially correcting for a failure to account for galaxy morphology. In addition, recent work has shown that many estimates of stellar masses at z > 5 may have been overestimated due to nebular line emission (Stark et al. 2013); however, the conversions in González et al. (2011) do not share this problem (D. Stark 2013, private communication). As the use of UV-only stellar masses therefore remains controversial, we present results both with and without the converted Bouwens et al. (2011) SMFs for z = 7–8 in Appendix J.

3.2. Cosmic Star Formation Rates

3.3. Specific Star Formation Rates

We rely mainly on the high-resolution Bolshoi simulation, described in Klypin et al. (2011). Bolshoi follows 2048³ (≈8.6 billion) particles in a comoving, periodic box with side length 250 h⁻¹ Mpc from z = 80 to the present day. Its mass resolution (1.9 × 10⁸ M_☉) and force resolution (1 h⁻¹ kpc) make it ideal for studying the evolution of halos from 10¹⁰ M_☉ (e.g., satellites of the Milky Way) to the largest clusters in the universe (10¹⁵ M_☉). Bolshoi was run as a collisionless dark matter simulation with the Adaptive Refinement Tree Code (ART; Kravtsov et al. 1997; Kravtsov & Klypin 1999) assuming a flat, ΛCDM cosmology (Ω_M = 0.27, $\Omega _\Lambda = 0.73$, h = 0.7, σ₈ = 0.82, and n_s = 0.95). These cosmological parameters are consistent with results from both WMAP5 (Komatsu et al. 2009) and the latest WMAP7+BAO+H₀ results (Komatsu et al. 2011).

Halos in Bolshoi were identified using Rockstar, a seven-dimensional halo finder that uses phase space plus temporal information (Behroozi et al. 2013b). Merger trees were generated using a new algorithm presented in Behroozi et al. (2013c), which enforces physically consistent evolution of halo properties across timesteps to provide increased robustness against anomalies in the halo finder. Halo masses are calculated using spherical overdensities, according to the virial overdensity criterion of Bryan & Norman (1998).

To supplement results from Bolshoi where there exist concerns about statistics (i.e., for high-mass halos) or for initial conditions (as Bolshoi uses Zel'dovich/1LPT initial conditions), we also make use of two larger simulations. MultiDark follows the same number of particles in a much larger volume (1000 h⁻¹ Mpc) with 7 times worse force resolution and 64 times worse mass resolution, but using the identical simulation code and cosmology (Riebe et al. 2011). Consuelo follows 1400³ particles in a somewhat larger volume (420 h⁻¹ Mpc) than Bolshoi; it has eight times worse force resolution and fourteen times worse mass resolution (J. McBride et al., in preparation; see also Behroozi et al. 2013c). ⁷ It has a slightly different cosmology (Ω_m = 0.25, $\Omega _\Lambda = 0.75$, h = 0.7, n_s = 1.0, and σ = 0.8) and was run using the GADGET-2 code (Springel 2005). Notably, however, it uses the more accurate second-order Lagrangian Perturbation Theory (2LPT) initial conditions instead of Zel'dovich initial conditions, so it serves as a way to correct for inaccuracies in the high-redshift mass function for Bolshoi.

Results from these three simulations are used to calibrate the halo mass function (Appendix G), halo mass accretion histories (Appendix H), and subhalo disruption rates (Appendix I) used in this work.

The method presented above results in a posterior distribution for the set of parameters describing models that match observed SMFs, SSFRs, and CSFRs from z = 0 to z = 8. All data results in this paper are available for download online.⁸ Our best-fitting parameters with 1σ limits are as follows:

Intrinsic parameters:

$$\begin{eqnarray*} \nu &=& \exp (-4a^2)\nonumber \\ \log _{10}(\epsilon) &=& -1.777^{+0.133}_{-0.146} + \left({-}0.006^{+0.113}_{-0.361}(a-1)\right. \nonumber \\ &&\left. + \left({-}0.000^{+0.003}_{-0.104} \right)z \right)\nu -0.119^{+0.061}_{-0.012}(a-1)\nonumber \\ \log _{10}(M_1) &=& 11.514^{+0.053}_{-0.009} + \left({-}1.793^{+0.315}_{-0.330}(a-1)\right.\nonumber\\ &&\left. + \left({-}0.251^{+0.012}_{-0.125} \right)z \right)\nu \nonumber \\ \alpha &=& -1.412^{+0.020}_{-0.105} + \left(0.731^{+0.344}_{-0.296}(a-1) \right)\nu \nonumber \\ \delta &=& 3.508^{+0.087}_{-0.369} + \left(2.608^{+2.446}_{-1.261}(a-1)\right. \nonumber\\ &&\left. + \left({-}0.043^{+0.958}_{0.071} \right)z \right)\nu \nonumber \\ \gamma &=& 0.316^{+0.076}_{-0.012} + \left(1.319^{+0.584}_{-0.505}(a-1)\right. \nonumber\\ &&\left. + 0.279^{+0.256}_{-0.081}z \right)\nu \nonumber \\ \log _{10}(M_{h,{\rm ICL}}) &=& 12.515^{+0.050}_{-0.429} + \left({-}2.503^{-0.202}_{-2.078} \right)(a-1)\nonumber \\ \rho _{0.5} &=& 0.799^{+0.028}_{-0.355}.\nonumber \end{eqnarray*}$$

Systematic parameters:

$$\begin{eqnarray*} \mu &=& -0.020^{+0.168}_{-0.096} + 0.081^{+0.078}_{-0.036}(a-1)\nonumber \\ \kappa &=& 0.045^{+0.110}_{-0.051} + \left({-}0.155^{+0.133}_{-0.133} \right)(a-1)\nonumber \\ \xi &=& 0.218^{+0.011}_{-0.033} + \left({-}0.023^{+0.052}_{-0.068} \right)(a-1)\nonumber \\ \sigma &=& 0.070 + 0.061^{+0.017}_{-0.008}(z-0.1)\nonumber \\ c_i(z) &=& 0.273^{+0.103}_{-0.222} \left(1+\exp \left(1.077^{+3.502}_{-0.099}-z \right)\right)^{-1}\nonumber \\ b &=& 0.823^{+0.043}_{-0.629}.\nonumber \end{eqnarray*}$$

Our total χ² error for the best-fit model from all sources (observational and theoretical) is 245. For the number of observational data points we use (628), the nominal reduced χ² is 0.4. While the true number of degrees of freedom is not easily computed (note for example that covariance matrices are not available for most SMFs in the literature), this suggests that our best-fit is a reasonable match to the data. This is shown visually in Figure 3, which shows the evolution of the SMF and the posterior distribution for the observed CSFR and in Figure 4, which shows the posterior distribution of SSFRs.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

As discussed in Section 2, each model in the distribution contains complete information on the average stellar content and formation history of halos as a function of mass and redshift. Herein, we focus on a few interesting implications, leaving more extensive coverage to a follow-up paper. In Section 5.1, we present derived SFRs and SFHs as a function of halo mass and redshift. Next, in Section 5.2, we discuss constraints on the stellar mass to halo mass relation, both as a function of historical halo mass and as a function of halo mass at z = 0. In Section 5.3, we compare the trajectories of the stellar mass and ICL buildup in galaxies. Then, in Section 5.4, we discuss instantaneous baryon conversion efficiencies and how they relate to integrated baryon conversion efficiencies as well as stellar ages and formation times in Section 5.5. In Section 5.6, we compare our main results to those obtained by previous studies. We discuss the main effects of the uncertainties we have modeled in Section 5.7; finally, in Section 5.8, we present new fitting formulae for individual galaxy SFHs relevant to observers.

5.1. Star Formation Rates and Histories

We show derived SFRs as a function of halo mass in the left panel of Figure 5, and the corresponding SFHs for galaxies at z = 0 in the right panel. As discussed in Section 5.3, the contribution from merging galaxies to the central galaxy is small, so SFHs for galaxies at z = 0 trace the SFR as a function of their progenitor's halo mass. Also, for this reason, SFHs for galaxies at a given redshift z_g are nearly the same as for the z = 0 stellar populations, except truncated at z = z_g.⁹

Zoom In Zoom Out Reset image size

We show similar plots with 1σ uncertainties in Figure 6. The left-hand panel demonstrates that the SFR at fixed halo mass has been monotonically decreasing since very early redshifts. This rate of decrease is different for different halo masses. At moderate to high redshifts (z > 2), larger halo masses generically have larger average SFRs. However, at lower redshifts, the highest mass halos (M_h ≳ 10¹⁴ M_☉) become so inefficient that they have lower SFRs than group-scale (10¹³ M_☉) halos or Milky-Way sized (10¹² M_☉) halos.

Zoom In Zoom Out Reset image size

From the perspective of individual galaxies, it is more illuminating to look at the SFH in the right panel of Figure 6. Because halos continually gain mass over time, and do so more rapidly at early redshifts, the SFH for galaxies is not monotonically decreasing. Instead, it increases rapidly with time, approximately as a power law in time. Depending on present-day halo mass, the galaxy's SFR reaches a peak at a redshift between z = 0.5 to z = 2.5 (higher redshifts for higher halo masses) and then decreases until the present day. The rate of decrease depends again on the halo mass, with high halo masses shutting off more rapidly than lower halo masses. Cluster-scale (M_h ≳ 10¹⁴ M_☉) halos form most of their stars rapidly, at early times, whereas galaxies in Magellanic Cloud-scale halos (10¹¹ M_☉) form stars over an extended period of time (see also Section 5.4).

5.2. The Stellar Mass–Halo Mass Relation

We show constraints on SMHM relation from z = 0 to z = 8 in the left panel of Figure 7 and on the SMHM ratio in the right panel. As seen in our previous work (Behroozi et al. 2010), there is a strong peak in the stellar mass to halo mass ratio at around 10¹² M_☉ to at least z ∼ 4 and a weaker peak still visible to z ∼ 8. While the location of the peak appears to move to higher masses with increasing redshift (consistent with Leauthaud et al. 2012), the abundance of massive halos is also falling off with increasing redshift.

Zoom In Zoom Out Reset image size

In terms of dwarf galaxies (M_h ∼ 10¹⁰ M_☉), we only have constraints from observations at z = 0. These have been the subject of recent interest due to the finding of higher-than-expected stellar mass to halo mass ratios in dwarf galaxies around the Milky Way (Boylan-Kolchin et al. 2012). However, these expectations have been set largely by the assumption that the stellar mass to halo mass ratio remains a scale-free power law below 10¹¹ M_☉. As seen in Figure 7, the low mass power-law behavior is broken below 10¹¹ M_☉, corresponding with an upturn in the SMF below 10^8.5 M_☉ (Baldry et al. 2008; this result has also been seen by A. Kravtsov, in preparation). This underscores the danger of assuming that faint dwarfs obey the same physical scaling relations as Magellanic-Cloud-scale galaxies; moreover, it also is a strong argument against fitting the SMHM relation with a double power law (see discussion in Appendix D).

Concerning the range of allowed SMHM relations, the observational systematics are large enough that our results are marginally consistent with an unchanging SMHM relation from z = 6 to z = 0. Nonetheless, the feature with strongest significance is a gradual decrease in stellar mass in the median 10¹¹ M_☉ halo from z = 0 to z = 2, followed by an increase again for redshifts z > 6; also potentially indicated is an increase in stellar mass in the median M_h > 10¹³ M_☉ halo from z = 0 to z = 2. The best-fit SMHM relations at z = 7 and z = 8 are significantly different than at lower redshifts, with more stellar mass per unit halo mass. However, concerns about the reliability of the SMFs at those redshifts (see Section 3.1) urge caution in interpreting the physical meaning of this result.

A useful perspective on these results can be obtained by considering the historical stellar mass to halo mass ratio of halos, as shown in Figure 8. Despite the large systematic uncertainties, it is clear that halos go through markedly different phases of star formation. This evolution is most apparent for massive halos, as observations have been able to probe the properties of the progenitor galaxies all the way to z = 8. Specifically, high-redshift progenitors of today's brightest cluster galaxies (M_h ∼ 10¹⁴ M_☉ were relatively efficient in converting baryons to stars—comparable to the most efficient galaxies today. However, between redshifts 2–3, their efficiencies peaked, and thereafter they began to form stars less rapidly than their host halos were accreting dark matter. At the present day, such galaxies have an integrated star formation efficiency that is two orders of magnitude less than at their peak. The picture is less clear for progenitors of lower-mass galaxies because current observations cannot probe their progenitors as far back. Nonetheless, their apparent behavior in Figure 7 of rising to a peak efficiency and later falling is consistent with all available data.

Zoom In Zoom Out Reset image size

5.3. Stellar Mass and Intracluster Light Growth Histories

Our model constrains the buildup of stars in the ICL (see definition in Section 2.4) purely from observational galaxy data and measurements of the halo–halo merger rate in simulations (see also Watson et al. 2012 for an alternate method). In our best-fitting model, only 5% of stellar mass in mergers for 10¹⁴ M_☉ halos is allowed to be deposited onto the central galaxy since z = 1, and only 10% for 10¹³ M_☉ halos. For Milky-Way-sized and smaller halos, this number rises rapidly to 70%–80%. Yet, due to the sharply decreasing stellar mass to halo mass ratio for lower-mass halos, most of the incoming stellar mass will be in (rare) major mergers. Central galaxies are therefore relatively uncontaminated by stars from smaller recently merged satellites. At higher redshifts, however, larger fractions of the stellar mass in merging galaxies are allowed to be deposited onto the central galaxy.

In Figure 9, we show the amount of galaxies' present-day stellar mass and ICL/halo stars that was in place at a given redshift. The left-hand panel (stellar mass) shows that almost all stars in the central galaxy in present-day cluster-scale halos were in place at z = 2. However, since that time, a large number of their satellites have been disrupted into the ICL. Thus, for cluster scale halos, the stellar mass in the ICL exceeds the stellar mass in the central galaxy by a factor of four to five, consistent with observations (e.g., Gonzalez et al. 2005). We note that our model predicts what may seem to be a large ICL fraction for Milky Way-sized galaxies (10¹² M_☉). However, the Milky Way is a special case. It has not had a major merger for ∼10–11 Gyr (Hammer et al. 2007); however, as noted above, only major mergers can contribute substantially to the ICL. A major merger 11 Gyr ago would have, however, contributed less than 3% of the present-day stellar mass of the Milky Way into the ICL; allowing for passive stellar evolution, this would result in less than 2% of the luminosity of the Milky Way coming from the IHL, in excellent agreement with observations (Purcell et al. 2007).

Zoom In Zoom Out Reset image size

Finally, in Figure 10, we show the inferred fraction of stellar mass growth in galaxies coming from in situ star formation (as opposed to galaxy–galaxy mergers) as a function of halo mass and redshift. At all redshifts greater than 1, the vast majority of stellar mass growth is from star formation, and hence SMFs alone may be used to infer galaxy SFRs. However, the stellar mass growth rate of high-mass (M_h > 10¹³ M_☉), low-redshift halos is such that it cannot be explained entirely by in situ star formation. For these halos, SMFs alone cannot be used to infer their SFRs (see also the discussion in Section 6 and in Appendix J). Interestingly, in the local universe the transition between merger-dominated growth and star formation growth in the local universe is roughly the mass of the Milky Way—low-mass galaxies are dominated by star formation, and galaxies in halos more massive than 10¹² M_☉ are dominated by merging.

Zoom In Zoom Out Reset image size

5.4. Baryon Conversion Efficiencies

Assuming that the baryon accretion rate for a halo is equal to the universal baryon fraction (f_b, which is 0.17 in our cosmology) times the overall halo mass accretion rate, we can calculate the instantaneous baryon conversion efficiency implied by our models. In the left-hand panel of Figure 11, we show the instantaneous conversion efficiency at fixed halo mass. Remarkably, this efficiency is always highest for Milky-Way-sized halos (10¹² M_☉), corresponding well with the peak in the stellar mass to halo mass ratio seen in Figure 7. Moreover, this efficiency (20%–40%) is constant to within a factor of two over a remarkably large redshift range, suggesting that the efficiency is only a weak function of accretion rate in halos of this mass. For more massive halos, there is a decrease in the baryon conversion efficiency following redshift 2–3. This suggests that the ability of accreted material to cool onto star-forming regions becomes impaired at lower redshifts for such halos. However, this decrease is gradual, taking place over many Gyr (see also Behroozi et al. 2013a). Rather than an abrupt change in the character of infalling gas, this may suggest that the decreasing density and accretion rate of gas make it gradually more difficult for cold clumps to form, especially in the presence of an active galactic nucleus.

Zoom In Zoom Out Reset image size

On the right-hand side of Figure 11, we show the historical baryon conversion efficiency for progenitors of z = 0 halos. For massive halos at z = 0, their conversion efficiency climbs steeply toward higher redshifts as their halo mass falls to roughly 10¹² M_☉; once their halo mass drops below that value, they then become less and less efficient at very high redshifts. For less-massive halos, such as those that are 10¹² M_☉ at z = 0, their conversion efficiency has been increasing from early times to the present day.

On the left-hand side of Figure 12, we show extrapolated historical stellar mass to halo mass ratios—proportional to integrated baryon conversion efficiencies—similar to the bottom panel of Figure 7. Notably, stellar mass to halo mass ratios also peak when halos reach 10¹² M_☉, indicative of steeply falling star formation efficiencies at higher and lower masses. It would appear that the maximum integrated stellar mass efficiency is around 20%–40% at all redshifts.

Zoom In Zoom Out Reset image size

5.5. Stellar Ages

In the right-hand panel of Figure 12, we show the historical stellar mass in progenitors of halos at z = 0 relative to the present-day stellar mass. This should not be confused with Figure 9, which shows what amount of the currently remaining stellar mass was in place at a given redshift; these differ mainly because of passive stellar evolution (massive stars that formed sufficiently long ago will burn out by the present day).

These data allow us to derive stellar mass-weighted ages as a function of stellar mass, shown in the left-hand panel of Figure 13. These ages are consistent with the stellar populations in massive galaxies (and halos) forming at very early times and having little star formation continuing to the present day. Less-massive galaxies have younger average stellar ages, consistent with the ongoing star formation seen in the SFHs for such galaxies in Figure 6. While we have little information on the progenitors of galaxies less massive than 10⁹ M_☉, there is evidence that the average ages of the stellar populations may increase for such galaxies. Specifically, as shown in the right-hand panel of Figure 13, the time required to form 50% and 90% of the galaxy's stars increases toward lower stellar masses, indicating that the SFHs become more and more flat (constant), on average. A perfectly flat SFH would result in an average age of about half the age of the universe, or roughly 6.7 Gyr; this is somewhat higher than that seen for 10⁹ M_☉ galaxies. On the other hand, for individual dwarf galaxies, the SFH is likely to be stochastic, which may result in a large scatter around this average value.

Zoom In Zoom Out Reset image size

For the most massive galaxies today (∼10¹² M_☉), the time required to form 50% of their stars was extremely short—only 1–2 Gyr. This in turn translates to a high SFR, on the order of 200–1000 M_☉ yr⁻¹, which is consistent with observations of ultraluminous infrared galaxies (ULIRGs; Magnelli et al. 2012; Michałowski et al. 2010a, 2010b; Daddi et al. 2005; Chapman et al. 2004).

5.6. Comparison to Other Results

We show a comparison of our best-fit results for the stellar mass to halo mass ratio at z = 0.1 to previously published results in Figure 14, and we show comparisons at z = 1.0 and z = 3.0 in Figure 15. Where possible, conversions to our assumed cosmology and halo mass definition have been applied. The results in this work are almost identical to our previous results in Behroozi et al. (2010), with the exception of a deviation at low halo masses due to the updated SMFs used. We also compare to the new constraint of Reddick et al. (2012), which uses additional input from the correlation function and conditional SMF as measured by SDSS. There is a slight discrepancy between the two results for small masses due to the different satellite fractions obtained using M_peak and $v^\mathrm{peak}_\mathrm{max}$ (still within our systematic errors), but the Reddick et al. (2012) result is only well constrained by additional data in the range 10¹²–10¹⁴ M_☉.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Comparisons to other data sets at z ⩽ 1.0 (Moster et al. 2010; Guo et al. 2010; Wang & Jing 2010; Zheng et al. 2007a; Yang et al. 2009a; Hansen et al. 2009; Lin & Mohr 2004) are discussed extensively in Behroozi et al. (2010), so we do not repeat that discussion here. The z = 3 comparison is notable because it illustrates the large discrepancies that can occur when different SMFs are used at z ≳ 3. For example, Yang et al. (2012) perform modeling for two separate high-z SMFs (Pérez-González et al. 2008 and Drory et al. 2005). Part of the discrepancy between the two sets of results may be due to a somewhat restrictive redshift fit, and part of it may be due to systematic biases in stellar masses at high redshifts.

In Figure 16, we show a comparison of galaxy stellar mass histories between our results and those of Leitner (2012). A direct comparison is difficult because Leitner (2012) only considers star-forming galaxies. However, according to Equation (8), the active fraction for 10^9.5 M_☉ galaxies is 90%, and the active fraction for 10¹⁰ M_☉ galaxies is 65%—although most of the passive 10¹⁰ M_☉ galaxies have only recently become passive. Comparison of the stellar mass histories for those galaxies is thus somewhat more feasible. Our stellar mass histories agree remarkably well with those in Leitner (2012) except at recent times for 10¹⁰ M_☉ galaxies (where our model includes more passive galaxies) and at early times for 10^9.5 M_☉ galaxies, where observations do not constrain properties of their progenitors.

Zoom In Zoom Out Reset image size

5.7. Systematic Uncertainties

Several aspects of the allowed parameter range are notable. Most importantly, there is no necessity for a large number of galaxies to be missed (due to dustiness or burstiness) at high redshift; indeed, the allowed incompleteness is in general less than 30%, as shown in the top panel of Figure 17. This would imply that most high-redshift surveys are not missing a large fraction (>50%) of galaxies above their nominal detection thresholds.

Zoom In Zoom Out Reset image size

The allowed parameter range is consistent with having both μ = 0 and κ = 0 at z = 0; i.e., no systematic offsets in stellar mass are necessary at z = 0. This is partially by design: the parameterization of the SMHM was chosen so that it could fit the z = 0 constraints without need for μ or κ; moreover, the cosmic and specific SFRs are derivative constraints, meaning that they only constrain how the SMHM evolves with time. While the main fits favor some redshift evolution in the systematic parameters, it is still possible to reasonably fit all observations while disallowing all systematic corrections (i.e., assuming that the stellar masses and SFRs used in this paper are on average correct and that all the surveys were 100% complete). Details of this fit are presented in Appendix J. We also present discussions therein of how the allowed parameter space varies when certain data sets are not used to constrain the fit; specifically, excluding z > 6 data; excluding cosmic star formation constraints; and excluding all star formation constraints.

5.8. Fits to Star Formation Histories

The most commonly used SFH fitting formula is a declining exponential with time (SFH(t)∝e^−t/τ), used by almost all the observational papers cited in this work (see references in Tables 3 and 5). Various modifications to this form have been suggested (e.g., te^−t/τ or e^+t/τ in Maraston et al. 2010), and there has been increasing recent evidence that increasing SFHs are more appropriate for massive galaxies at high redshifts.

Our best constraints in this study are on average SFHs for halos at a given mass. Thus, we can place the strongest priors on the expected form of stacked galaxy SFHs. However, we can also place some weaker constraints on what individual galaxy SFHs should look like.

Individual galaxies' histories can differ from the average SFHs for three reasons. The first is stochasticity in the SFR on short timescales relative to the dynamical time of the halo, due to feedback effects as well as shot noise in the SFR. This variation cannot be modeled with our current approach, so we caution that the constraints on individual galaxy histories that we derive should be considered as smoothed over a period of at least 500 Myr. If enough color bands are available to avoid degeneracies with other fitting parameters, observers may consider adding the possibility of recent starbursts to their models to account for the fact that they can significantly alter observed galaxy colors.

The next reason for deviations is sustained deficits in the SFR for satellite galaxies compared to field galaxies (e.g., Wetzel et al. 2012). As shown in Figure 24 (see also, e.g., Kravtsov et al. 2004), the satellite fraction of halos falls both with increasing mass and increasing redshift. This is not a practical concern for z > 3, where the satellite fraction is less than 20% for all currently observable galaxies (M_* > 10^9.5 M_☉). At lower redshifts, the observation that satellites cease to form stars relatively quickly after accretion (Wetzel et al. 2012) means that declining SFHs may continue to be reasonable fits for satellite galaxies, even though most other galaxies of the same mass may be better fit by rising SFHs.

Finally, we note that differences can arise in individual galaxy SFHs because of the scatter in stellar masses at a given halo mass, as well as the scatter in mass accretion histories for halos (see Figure 18). We can attempt to model these effects by sampling random mass accretion histories from the Bolshoi simulation and then sampling random stellar mass growth histories as allowed by the scatter in stellar mass at fixed halo mass (i.e., choosing random stellar mass offsets at z = 8, z = 1.0, and z = 0 and using spline interpolation at intermediate times). Our results, expressed as a function of time since the beginning of the universe, are shown in Figure 19.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

These results have dramatically larger error bars than those in Figure 6 on account of the spread in mass accretion histories for halos. Nonetheless, some basic conclusions can be drawn. SFHs for z > 3 galaxies increase with time. Although there are different ways to parameterize these histories, a straightforward one is a direct power law form:

$$\begin{equation} {\rm SFH}(t) = A t^B\;\hbox{(for z>3)}. \end{equation} \tag{ 22 }$$

This is also equally capable of fitting the average SFRs in Figure 6 at z > 3. Our result is similar to that of Papovich et al. (2011); however, we find steeper slopes (B ∼ 3–4) than theirs (B ∼ 1.7) because they ignore all effects of mergers that occur from z = 8 to z = 3.

At lower redshifts, there is a mass-dependent turnover after which the SFH begins to decline. This happens at z ∼ 2–3 for 10¹¹ M_☉ galaxies, making a declining exponential at z = 0 a reasonable fit. However, it happens as late as z = 0.5 for 10^9.5 M_☉ galaxies—meaning that a declining exponential is never a good fit for these galaxies unless they are satellites. The best-fit in general for all the constraints on individual histories presented here is a double power law:

$$\begin{equation} {\rm SFH}(t) = A \left[\left(\frac{t}{\tau }\right)^B + \left(\frac{t}{\tau }\right)^{-C}\right]^{-1}; \end{equation} \tag{ 23 }$$

typical values of B and C range from 1 to 5. It may often be the case that the available data is insufficient to constrain all three shape parameters. In this case, a hybrid of the exponential decline and power law rise still provides a reasonable fit:

$$\begin{equation} {\rm SFH}(t) = A t^B \exp (-t/\tau). \end{equation} \tag{ 24 }$$

The existence of a "cold mode" of gas accretion that allows efficient star formation at high redshifts and shuts off for massive galaxies past z ∼ 2 has been predicted in hydrodynamical simulations for the past decade (Birnboim & Dekel 2003). In this work, we observe a transition after z ∼ 3, where galaxies in massive halos become progressively less efficient at forming stars, even after correcting for slowing accretion rates since that time (Figure 11, right panel). At fixed halo mass (Figure 11, left panel), the transition is less severe, but there has been a clear reduction in the conversion efficiency since z ∼ 3 to the present day. This suggests a picture in which dense infalling gas can cool to the galaxy disk efficiently at high redshifts; however, at lower redshifts, the infalling gas becomes sparse enough that it is susceptible either to shock heating or feedback from an active galactic nucleus.

Another important feature of our results is that merging galaxies' stars almost never reach the central galaxy in massive halos (<10% since z = 1 for halos larger than 10¹³ M_☉). This claim is not new (Purcell et al. 2007; Conroy et al. 2007), but it is reassuring that we duplicate this conclusion with high confidence. This is in stark contrast to what is asserted in (Moster et al. 2013; i.e., 80% of merging stars reach the central galaxy); however, their assertion is based on simulations of Milky-Way-sized halos, for which we find that indeed ∼70% of merging galaxies' stars reach the central galaxy. Thus, the fraction of merging stars that reach the central galaxy is likely a strong function of halo mass.

It is worth mentioning that the SFRs derived in Moster et al. (2013) appear reasonable despite this issue. We note that Moster et al. (2013) assumes that satellite stellar mass is fixed at the epoch of accretion, which results in less stellar mass in satellites than in this work. Additionally, Moster et al. (2013) assumes a higher value of the scatter in observed stellar mass versus true stellar mass (σ) than we do, which comes from comparing technique-to-technique scatter in recovering observed stellar masses. This results in a much larger allowed evolution for massive galaxies between z = 1 and z = 0, which is very sensitive to the value of σ (Behroozi et al. 2010). The combination of these two assumptions (less stellar mass in satellites and greater allowed evolution in the stellar mass of massive galaxies) roughly cancels the effect on the SFR of assuming that more stellar mass in satellites is allowed to merge into the central galaxy. We note, however, that our satellite stellar masses are already on the low side of what is allowed by clustering constraints (Reddick et al. 2012), leaving little room to reduce them further as is done in the Moster et al. (2013) analysis. In addition, the scatter in the difference between two stellar mass estimation techniques tends to overestimate the scatter between any one technique and the true value by a factor of ${\sim} \sqrt{2}$. This factor is consistent with other independent estimates of σ, for example, by evaluating the posterior distribution of allowed recovered stellar masses (Kajisawa et al. 2009; Conroy et al. 2009). Thus, both the actual evolution in the stellar mass of massive galaxies and the fraction of stars in mergers which reach the central galaxy may not be as large as the values implied by the Moster et al. (2013) analysis.

We note further that the amount of stellar mass that was contained in disrupted galaxies can easily reach five times the amount of stellar mass in the central galaxy for high-mass halos (Figure 9), consistent with observations in Gonzalez et al. (2005). Regardless of uncertainties in determining stellar masses, there is no possible way to funnel any significant fraction of this mass into the central galaxy without badly violating observed limits on the number density of massive galaxies. Rather, stars from merging galaxies must end up in a more extended distribution to avoid confusion with the central galaxy, a topic that will be explored in depth in a future paper (P. S. Behroozi et al., in preparation). From the opposite perspective, galaxy stellar populations contain mostly stars from the most-massive progenitor, with some contamination from major mergers. Surveys of stars within a galaxy's optical radius (such as for the Milky Way) thus capture a potentially limited range of different star-forming environments, at least for stars formed since z = 1. A wider range of environments can be sampled by considering the ICL and halo stars in general.

The most efficient halo mass for converting baryons into stars appears to be close to 10¹² M_☉, regardless of redshift. The progenitors of today's massive halos (≳ 10¹⁴ M_☉) passed through this mass quickly, and as a result formed most of their stars in a short time: a mere 1–2 Gyr (Figure 13). Given that baryon to star conversion becomes so inefficient above and below this halo mass, much of the shape of the SFHs seen in this paper (Figure 6) could be explained by considering the duration of time that a given halo's mass accretion history spends close to this most efficient halo mass. That is to say, the peakiness of the SFH depends on how quickly the halo mass accretion history passes through halos of mass ∼10¹² M_☉. As lower-mass halos have slower relative accretion rates, they have flatter SFHs; e.g., Milky-Way-sized halos have lingered in this mass range from z ∼ 1, and so have the current highest stellar mass to halo mass ratio, as well as SFHs that track their halo mass accretion rates well at late times.

Low-mass halos (10¹¹ M_☉) show an unusual behavior in terms of their baryon conversion efficiency. In the left panel of Figure 11, it would appear that the instantaneous baryon conversion efficiency of 10¹¹ M_☉ halos is high at high redshifts (z > 5), low at intermediate redshifts (2 < z < 5), and then high again at low redshifts. While this behavior matches what is seen in the integrated baryon conversion efficiency for such halos (Figure 12, left panel), the corresponding galaxies are often at the faint edge of the completeness limit for stellar mass surveys. For that reason, observational biases may be the most likely source of this behavior.

Our results have importance for future observational studies. The current levels of systematic errors in the inference of stellar masses and SFRs are the dominant systematic uncertainty in this work. As discussed in Section 3, SMFs at z > 3 can disagree by up to 0.5 dex, well beyond any errors expected from sample variance. Star formation rates are equally concerning, both at high redshifts and at low redshifts, especially in terms of SSFRs. As a result, systematic errors are the single most important aspect preventing better understanding of the formation histories of massive galaxies. As discussed in Section 5.8, we have provided new fitting formulae for galaxy SFHs, which should help reduce one of the largest sources of systematic error especially for low-mass and high-redshift galaxies. New studies of SMFs and SFRs from deep multiwavelength surveys like CANDELS (Koekemoer et al. 2011; Grogin et al. 2011), analyzed with self-consistent SFHs and assumptions about systematic errors, will provide a large step forward in constraining the model presented here.

We have presented a comprehensive analysis of average SFRs and SFHs in galaxies from z = 0 to z = 8, and their connection to the underlying growth and merging of dark matter halos, along with a treatment of the inherent uncertainties. Our approach provides a self-consistent picture of the growth of galaxies over this epoch in the cosmological context of ΛCDM structure formation. The model is able to match the evolution of CSFRs, SSFRs, and the SMF of galaxies over the last 13.2 billion years of cosmic time. Our main findings are as follows.

• 1.  
  Halos of mass ∼10¹² M_☉ appear to be the most efficient at forming stars at every epoch, with baryon conversion efficiencies of 20%–40% over nearly the entire redshift range of this study. Halos at higher and lower masses are less efficient by orders of magnitude, especially at low redshifts.
• 2.  
  The baryon conversion efficiency of more massive halos is still reasonably high (10%–20%) until z = 2–3, after which time it takes a steep downturn. This is consistent with expectations of reduced cold mode accretion onto the corresponding galaxies.
• 3.  
  We have characterized the fraction of stars in galaxies due to merging versus star formation. At z > 1, galaxy buildup in all but the most massive galaxies was strongly dominated by star formation. At the present day, the transition between merger-dominated buildup (at high mass) and star formation-dominated buildup (at low mass) occurs at roughly the mass of the Milky Way.
• 4.  
  We confirm previous results in support of most merging galaxies' stars being disrupted in the ICL for massive halos. We find, however, that the fraction of merging stars that reach the central galaxy may increase strongly with decreasing halo mass, to about 70% for Milky-Way-sized halos.
• 5.  
  For massive galaxies (M_* > 10¹¹ M_☉), we predict that at least as much stellar mass is present in the ICL as in the galaxy itself out to intermediate redshifts, possibly as high as z ∼ 1 for the largest clusters (although this depends on the exact cut between ICL and galaxy stellar mass).
• 6.  
  SFHs peak at z ∼ 3 for massive galaxies and are increasing nearly to the present day for small galaxies. We have presented new fitting formulae for galaxy SFHs: a featureless power law for high redshifts (Equation (22)), and more general formulae valid for all redshifts (Equations (23) and (24)).
• 7.  
  Upturns in the SMF at low stellar masses mean that dwarf galaxies (M_* < 10^8.5 M_☉) have a higher stellar mass to halo mass ratio than that predicted by extrapolating results from more massive galaxies. The upturn is most robustly observed at low redshifts, but there are also indications at higher redshifts of a steepening slope to the SMF at low stellar masses. The stellar mass to halo mass relation thus cannot be fit with a double power law alone.
• 8.  
  Systematic errors remain at high levels (∼0.2–0.3 dex in stellar masses, 0.15–0.3 dex in SFRs), and disagreement between galaxy abundances for high-redshift surveys (z > 3) is substantial.
• 9.  
  Galaxies with stellar masses below 10¹⁰ M_☉ have the least well-constrained properties in this study. This could be improved by future surveys and by incorporating information on the SFHs of local galaxies.

Support for this work was provided by an HST Theory grant. Program number HST-AR-12159.01-A was provided by NASA through a grant from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Incorporated, under NASA contract NAS5-26555. This research was also supported in part by the National Science Foundation under grant No. NSF PHY11-25915, through a grant to KITP during the program "First Galaxies and Faint Dwarfs," and through a ARCS Foundation Graduate Fellowship to P.S.B.

We thank Anatoly Klypin and Joel Primack for access to the Bolshoi and MultiDark simulations and Michael Busha and the LasDamas collaboration for access to the Consuelo simulation. The latter was run using computational resources at SLAC. We thank John Moustakas for access to the PRIMUS stellar mass functions prior to publication. This work has benefited from many interesting conversations with Tom Abel, Steve Allen, Andrew Benson, Andreas Berlind, Kevin Bundy, Joanne Cohn, Sandy Faber, Dusan Keres, David Koo, Andrey Kravtsov, Alexie Leauthaud, Avi Loeb, Yu Lu, Chung-Pei Ma, John Moustakas, Joel Primack, Rachel Reddick, Beth Reid, Daniel Stark, Chuck Steidel, Tommaso Treu, David Wake, Andrew Wetzel, Frank van den Bosch, Alexey Vikhlinin, and Naoki Yoshida, We thank Yu Lu, Ari Maller, Simone Wienmann, and Anatoly Klypin for helpful comments on an earlier draft, and Surhud More and Michal Michałowski for catching typographical mistakes. We additionally thank the anonymous referee for many suggestions which improved this paper.

As mentioned in the introduction to Section 2, there is a potential error in our method because we propagate average SFHs for halos at a given mass and redshift, as opposed to individual SFHs. We discuss effects from both a distribution in SFRs for galaxies at a given halo mass and from a distribution of halo mass accretion histories.

Bimodality has been observed in the SFRs of galaxies at fixed stellar mass, partially but not exclusively driven by quenched satellite galaxies (see, e.g., Wetzel et al. 2012). This may appear to be a problem for our method because the rate of specific mass loss (see Section 2.3) will be different depending on whether a given galaxy is active or passive. Remarkably, however, our method is robust to any distribution of SFRs in galaxies of the same stellar mass. Because our method only depends on the total stellar mass buildup in halos of a given mass M_h at any given time, the fact that some galaxies get a larger share of that buildup than others is not a concern. If certain galaxies sustain a higher growth rate for a long time, then the scatter in stellar mass at fixed halo mass will certainly increase; however, the possibility of scatter changing with time is an effect that we include explicitly in our model. Bimodality in SFRs is only a problem if it results in a distribution of stellar masses at fixed halo mass that cannot be reasonably modeled by a log-normal distribution (i.e., the distribution we assume in this work). So far, no evidence has emerged for any other distribution to be preferred (see comments in Section 2.2.1 where ξ is defined).

There is more possibility for error because we do not account for correlations between the SFH and the mass accretion history of a halo. For example, a rapidly accreting halo might have a more recent SFH than other halos at the same mass. In this case, when it transitions from one mass bin to another, we may use the incorrect SFH in Equation (B2). This effect is balanced by halos that accrete more slowly than other halos of the same mass, which cancels the effect to first order.

It is worth working through a more quantitative example to verify this assertion. Consider a population of galaxies all with the same stellar mass M_* = 10¹⁰ M_☉ today. Assume that each galaxy has had a constant SSFR for its entire lifetime (i.e., an exponential growth history), and moreover that the distribution of the SSFRs is log-normal with width 0.3 dex (similar to that seen in Salim et al. 2007) and median 10⁻⁹ yr⁻¹. The average (as opposed to median) SFR for this population is 12.7 M_☉ yr⁻¹; given the SFHs just specified, this translates into an average rate of change of stellar mass of 9.3 M_☉ yr⁻¹ over the next 100 Myr. If, on the other hand, we had simply assumed that all of the galaxies had SFHs corresponding to a constant SSFR of 10⁻⁹ yr⁻¹, we would interpret a buildup in stellar mass of 9.3 M_☉ yr⁻¹ as implying an instantaneous SFR of 13.4 M_☉ yr⁻¹, an error of 6%. Thus, we conclude that the errors on the average quantities considered here introduced by lumping together SFHs for halos with different mass accretion histories are strongly subdominant to existing observational uncertainties.

Because of passive stellar evolution, calculating the amount of stellar mass remaining in a galaxy (required for determining the SFR in Equation (1)), one must know the entire SFH. For a large simulation, keeping track of the individual stellar mass histories (covering hundreds of timesteps) for millions of galaxies repeated over several million realizations of the SMHM relation is inefficient. Instead, since we are interested in only the average SFR as a function of mass, we take an elegant (although approximate) shortcut.

Instead of keeping track of every halo/galaxy individually, we instead keep track of average SFHs at each timestep in bins of halo mass (which is a rank-3 tensor, ${\rm SFH}_{t,m}^{t_{{\rm form}}}$; t is the time now, t_form is the timestep at which the stellar mass was formed, and m is the halo mass). For central halos, we use the virial mass (as discussed in Section 4); for satellite halos, we use the mass of the largest progenitor (the peak mass) in the satellite's mass accretion history. At every timestep, we calculate the number of halos that transition between one mass bin and another ($T_{t,m_n}^{m_p}$), i.e., the number of halos that had a mass m_p at timestep t and a mass m_n at timestep t + 1 (see Appendices G and H). Similarly, we calculate the merger rates for subhalos at each timestep as a function of the subhalo peak mass and the target host halo mass ($M_{t,m_n}^{m_p}$), as discussed in Appendix I. Keeping track of SFHs then becomes an exercise in tensor multiplication, which is fast regardless of the number of halos in the simulation.

To convert Equation (1) into a full tensor equation, we need a small amount of additional notation. The stellar mass at a particular timestep t and halo mass m becomes SM_{t, m}, the SFR at a particular timestep and halo mass becomes SFR_{t, m}, and the number counts at a given timestep for halos in a given mass bin becomes N_{t, m}. In addition, we denote the difference in stellar mass from the previous timestep by ΔSM_{t, m}. Thus, Equation (1) formally becomes

$$\begin{eqnarray} \Delta {\rm SM}_{t,m} &=& \left({\rm SM}_{t,m} - \sum _{t_{{\rm form}}} {\rm SFH}_{t,m}^{t_{{\rm form}}} (1-C_\mathrm{loss}(t, t_{{\rm form}} - \Delta t, t_{{\rm form}}))\right)\nonumber \\ {\rm SFR}_{t,m} & = & \frac{f_{{\rm SFR}}(m,t)\Delta {\rm SM}_{t,m}}{\Delta t (1-C_\mathrm{loss}(t, t-\Delta t,t))}\nonumber \\ C_\mathrm{loss}(t_{{\rm now}}, t_a,t_b) & = & |t_b-t_a|^{-1} \int _{t_a}^{t_b} f_\mathrm{loss}(t_{{\rm now}}-x) dx, \end{eqnarray} \tag{ B1 }$$

where C_loss(t_now, t_a, t_b) is the average stellar mass loss fraction (a.k.a., "recycling" fraction) at the present time for stars that formed at a constant rate from t_a to t_b. (Note that the indices t and m are not meant to be summed over).

To calculate SFHs, we need to know both internal star formation and the amount of stars accreted in mergers. As discussed in Section 2.4, purely definitional changes in what constitutes the boundary between a galaxy and the ICL can effectively shift stellar mass from the ICL into the galaxy. Not knowing the relative fraction of each, we assume that merging stellar mass from satellites takes priority over stellar mass from the ICL. For convenience, we define ${\rm SM}^{\prime }_{t,m}$ and ${\rm ICL}^{\prime }_{t,m}$ as the total stellar mass in merging satellites and the total stellar mass in the main progenitors' ICL, respectively; we also define ${\rm ICLH}_{t,m}^{t_{{\rm form}}}$ to be the SFH of the ICL (analogously to ${\rm SFH}_{t,m}^{t_{{\rm form}}}$). The new SFHs will then be the sum of the old SFHs from all progenitors, with the addition of new stars formed at the latest timestep and material from mergers and the ICL:

$$\begin{equation} {\rm SFH}_{t+\Delta t,m}^{t_{{\rm form}}} = \sum _{m_p} \left[{\rm SFH}_{t,m_p}^{t_{{\rm form}}} \left(T_{t,m}^{m_p} + C^{{\rm merge}}_{t,m} M_{t,m}^{m_p}\right) + {\rm ICLH}_{t,m_p}^{t_{{\rm form}}} T_{t,m}^{m_p} C^{{\rm def}}_{t,m} \right] (N_{t+\Delta t,m})^{-1} + {\rm SFR}_{t,m} \delta _{t}^{t_{{\rm form}}} \Delta t \end{equation} \tag{ B2 }$$

$$\begin{equation} C^{{\rm merge}}_{t,m} = \min \left(1,\frac{(1-f_{{\rm SFR}}(m,t))\Delta {\rm SM}_{t+\Delta t,m}}{{\rm SM}^{\prime }_{t,m}}\right) \end{equation} \tag{ B3 }$$

$$\begin{equation} C^{{\rm def}}_{t,m} = \max \left(0, \frac{(1-f_{{\rm SFR}}(m,t))\Delta {\rm SM}_{t+\Delta t,m}-{\rm SM}^{\prime }_{t,m}}{{\rm ICL}^{\prime }_{t,m}}\right) \end{equation} \tag{ B4 }$$

$$\begin{equation} {\rm SM}^{\prime }_{t,m} = \sum _{m_p} {\rm SFH}_{t,m_p}^{t_{{\rm form}}} M_{t,m}^{m_p} (1-C_\mathrm{loss}(t+\Delta t, t_{{\rm form}} - \Delta t, t_{{\rm form}})) (N_{t+\Delta t,m})^{-1} \end{equation} \tag{ B5 }$$

$$\begin{equation} {\rm ICL}^{\prime }_{t,m} = \sum _{m_p} {\rm ICLH}_{t,m_p}^{t_{{\rm form}}} T_{t,m}^{m_p} (1-C_\mathrm{loss}(t+\Delta t, t_{{\rm form}} - \Delta t, t_{{\rm form}})) (N_{t+\Delta t,m})^{-1}. \end{equation} \tag{ B6 }$$

(Again, the indices t and m are not meant to be summed over; $\delta _{t}^{t_{{\rm form}}}$ is the standard identity tensor: $\delta _{t}^{t_{{\rm form}}} = 1$ if t = t_form and 0 otherwise). Note that, in this present analysis, we do not separate the SFHs of satellites and centrals. As discussed in Section 5.8, satellites at high redshifts (z > 1) tend to have short infall timescales, so they have roughly similar SFHs as centrals. At later times, this approximation is less valid; however, we can nonetheless correctly recover the average SFH across all halos (centrals and satellites) at a given halo mass, as discussed in Appendix A.

Calculating the ICL formation history as well as the total amount of stars in the ICL is similar:

$$\begin{equation} {\rm ICLH}_{t+\Delta t,m}^{t_{{\rm form}}} = \sum _{m_p} \left({\rm SFH}_{t,m_p}^{t_{{\rm form}}} \left(1-C^{{\rm merge}}_{t,m} \right) M_{t,m}^{m_p} + {\rm ICLH}_{t,m_p}^{t_{{\rm form}}} T_{t,m}^{m_p} \left(1-C^{{\rm def}}_{t,m} \right) \right) (N_{t+\Delta t,m})^{-1} \quad\quad \end{equation} \tag{ B7 }$$

$$\begin{equation} {\rm ICL}_{t,m} = \sum _{t_{{\rm form}}} {\rm ICLH}_{t,m_p}^{t_{{\rm form}}} (1-C_\mathrm{loss}(t, t_{{\rm form}} - \Delta t, t_{{\rm form}})). \end{equation} \tag{ B8 }$$

Note that the merging satellite ICL is not deposited into the central galaxy's ICL, which is the same approach taken by Conroy et al. (2007) and Purcell et al. (2007). Because the ICL of a merging satellite halo will get stripped long before the stars in its associated galaxy, the stripped satellite's ICL will end up in a much more extended, faint distribution than the host's ICL.

This formal specification is equivalent to a discrete differential equation for the stellar mass history. It is thus important to properly specify the boundary conditions; that is, the stellar mass history at t = 0. In the absence of constraints on SMFs above z ≈ 9, we extrapolate the fit to the SMHM relation to z = 15 and assume a uniform SFR prior to that redshift. While incorrect in detail at z = 15, this choice has little effect on the SFRs or SFHs of galaxies at later (and more visible) times.

Given the large number of equations in our framework, it is appropriate to describe exactly how a point in the parameter space of the cosmic stellar mass history (Equations (3)–(5)) and uncertainties (Equations (9)–(18)) gets converted into quantities comparable to published observations. The three observables in our case are SMFs, SSFRs, and the CSFR density. We deal with each in turn.

C.1. The Observed Stellar Mass Function

Equation (3) gives the median stellar mass as a function of halo mass for the entire redshift range we consider. If no scatter between stellar mass and halo mass were present, the conversion from the halo mass function to the fiducial SMF (ϕ_fid) would be a straightforward application of the chain rule:

$$\begin{equation} \phi _\mathrm{fid}(M_*) = \frac{dN}{d \log _{10} M_\ast } = \frac{dN}{d\log _{10} M_h(M_\ast)} \left(\frac{d\log _{10} M_\ast }{d\log _{10} M_h}\right)^{-1}. \end{equation} \tag{ C1 }$$

In other words, the SMF is equal to the halo mass function divided by the logarithmic derivative of Equation (3). As discussed in Behroozi et al. (2010), this may be convolved with the scatter in observed stellar mass at fixed halo mass ($P(\Delta _{M_\ast })$) to give the true SMF:

$$\begin{equation} \phi _\mathrm{true} (M_\ast) = \int _{-M_\ast }^{\infty } \phi _\mathrm{fid}(M_\ast + \Delta _{M_\ast }) P(\Delta _{M_\ast })d\Delta _{M_\ast }. \end{equation} \tag{ C2 }$$

In our model, both the scatter in true stellar mass at fixed halo mass (Equation (5)) and the scatter in observed stellar mass at fixed true stellar mass (Equation (11)) are log-normal distributions; hence $P(\Delta _{M_\ast })$ is a log-normal distribution with width equal to the sum in quadrature of the widths of the two sources of scatter ($\sqrt{\sigma ^2 + \xi ^2}$).

To convert from the true SMF to the observed SMF, we need to correct not only for observational completeness (Equation (13)), but also for systematic biases in the measured stellar masses (Equations (9) and (10)). This latter correction is made somewhat trickier by the fact that we only know the passive fraction as a function of the measured stellar mass. The true passive fraction, f_{passive, true}(M_*) is related to the observed passive fraction by

$$\begin{equation} \scriptsize{f_\mathrm{passive,obs}(M_{\ast,\mathrm{meas}}) = \frac{\phi _\mathrm{true}(M_{\ast,\mathrm{meas}}10^{-\mu }) f_\mathrm{passive,true}(M_{\ast,\mathrm{meas}}10^{-\mu })}{\phi _\mathrm{true}(M_{\ast,\mathrm{meas}}10^{-\mu }) f_\mathrm{passive,true}(M_{\ast,\mathrm{meas}}10^{-\mu }) + \phi _\mathrm{true}(M_{\ast,\mathrm{meas}}10^{-\mu -\kappa }) (1-f_\mathrm{passive,true}(M_{\ast,\mathrm{meas}}10^{-\mu -\kappa }))}.} \end{equation} \tag{ C3 }$$

In the case that κ = 0, then this equation reduces to f_{passive, true}(M_{*, meas}10^−μ) = f_{passive, obs}(M_{*, meas}). Otherwise, this equation becomes a recurrence relation:

$$\begin{equation} \frac{f_\mathrm{passive,true}(M_{\ast,\mathrm{meas}}10^{-\mu })}{1-f_\mathrm{passive,true}(M_{\ast,\mathrm{meas}}10^{-\mu -\kappa })} = \frac{\phi _\mathrm{true}(M_{\ast,\mathrm{meas}}10^{-\mu -\kappa })}{\phi _\mathrm{true}(M_{\ast,\mathrm{meas}}10^{-\mu })}\left[(1-f_\mathrm{passive,obs}(M_{\ast,\mathrm{meas}}))^{-1}-1\right]. \end{equation} \tag{ C4 }$$

Solving this equation is possible, but using the result has disadvantages. Specifically, calculating f_{passive, true} couples the values of ill-constrained regions of the SMF to well-constrained regions through this recurrence relation. The values for α—the faint-end slope of the SMHM relation—and κ will then be constrained by the MCMC algorithm attempting to improve fits to the well-constrained regions of the SMF, and will lose their physical significance.

We therefore compromise by decoupling the true passive fraction from the value of κ. We take f_{passive, true}(M_{*, meas}10^−μ) = f_{passive, obs}(M_{*, meas}), the exact solution for κ = 0. This results in a change to the definition of κ so that it no longer corresponds to an exact offset for active galaxy stellar masses; however, the correlation between κ and the average offset implied for active galaxies should still be high. With this approach, it is straightforward to calculate the measured SMF:

$$\begin{equation} \scriptsize{\phi _\mathrm{meas}(M_{\ast,\mathrm{meas}}) = f_\mathrm{passive,true}(M_{\ast,\mathrm{meas}}10^{-\mu })\phi _\mathrm{true}(M_{\ast,\mathrm{meas}}10^{-\mu }) + \phi _\mathrm{true}(M_{\ast,\mathrm{meas}}10^{-\mu -\kappa })(1-f_\mathrm{passive,true}(M_{\ast,\mathrm{meas}}10^{-\mu -\kappa })).} \end{equation} \tag{ C5 }$$

Or, in terms of the observed passive fraction,

$$\begin{equation} \scriptsize{\phi _\mathrm{meas}(M_{\ast,\mathrm{meas}}) = f_\mathrm{passive,obs}(M_{\ast,\mathrm{meas}})\phi _\mathrm{true}(M_{\ast,\mathrm{meas}}10^{-\mu }) + \phi _\mathrm{true}(M_{\ast,\mathrm{meas}}10^{-\mu -\kappa })(1-f_\mathrm{passive,obs}(M_{\ast,\mathrm{meas}}10^{-\kappa })).} \end{equation} \tag{ C6 }$$

Then, to correct for observational completeness, we multiply by the appropriate factor:

$$\begin{equation} \phi _\mathrm{obs}(M_\mathrm{\ast,\mathrm{meas}}) = c(z) \phi _\mathrm{meas}(M_\mathrm{\ast,\mathrm{meas}}). \end{equation} \tag{ C7 }$$

Finally, as surveys usually return the SMF over a range of redshifts (e.g., z₁ < z < z₂), we must perform one last integration to determine the survey mass function:

$$\begin{equation} \phi _\mathrm{survey}(M_{\ast,\mathrm{meas}}) = \frac{\int _{z_1}^{z_2} \phi _\mathrm{obs}(M_{\ast,\mathrm{meas}})dV_c(z)}{V_c(z_2) - V_c(z_1)}, \end{equation} \tag{ C8 }$$

where V_c(z) is the comoving volume out to redshift z.

C.2. Calculating the Observed Specific Star Formation Rates and Cosmic Star Formation Rate Density

Equation (B1) gives the average SFR for halos as a function of halo mass and redshift. Equation (3) gives the median stellar mass for halos as a function of halo mass and redshift. Ignoring scatter, the fiducial SSFR is simply

$$\begin{equation} {\rm SSFR}_\mathrm{fid}(M_\ast) = \frac{{\rm SFR}(M_h(M_\ast))}{M_\ast }. \end{equation} \tag{ C9 }$$

Scatter in the observed stellar mass at fixed halo mass results in contributions from lower- and higher-mass halos, weighted by number density. If there is no correlation between stellar mass and SFR beyond the correlation between halo mass and SFR, then the SSFR would be as follows:

$$\begin{equation} {\rm SSFR}(M_\ast)_\mathrm{nc} = M_\ast ^{-1} \phi _\mathrm{true}^{-1}(M_\ast)\int _{-M_\ast }^{\infty } \phi _\mathrm{fid}(M_\ast + \Delta _{M_\ast }) {\rm SFR}(M_h(M_{\ast }+\Delta _{M_\ast })) P(\Delta _{M_\ast })d\Delta _{M_\ast }. \end{equation} \tag{ C10 }$$

However, with a nonzero covariance (ρ) between SFR and stellar mass at fixed halo mass this equation is amended to become

$$\begin{eqnarray} {\rm SSFR}(M_\ast)_\mathrm{true} & =& M_\ast ^{-1} \phi _\mathrm{true}^{-1}(M_\ast)\int _{-M_\ast }^{\infty } \phi _\mathrm{fid}(M_\ast + \Delta _{M_\ast }) ({\rm SFR}(M_h(M_{\ast }+\Delta _{M_\ast }))10^{\rho \Delta _{M_\ast }^{\prime }}) P(\Delta _{M_\ast })d\Delta _{M_\ast }\nonumber \\ \Delta _{M_\ast }^{\prime } &=& \log _{10}\left(\frac{M_\ast + \Delta _{M_\ast }}{M_\ast \exp (0.5(\xi \ln 10)^2)}\right). \end{eqnarray} \tag{ C11 }$$

The formula for $\Delta _{M_\ast }^{\prime }$ is necessary because $M_h(M_\ast + \Delta _{M_\ast })$ gives the halo mass at which the median stellar mass is $M_\ast + \Delta _{M_\ast }$; however, the covariance between SFR and SM is calculated relative to the average stellar mass at fixed halo mass. The exponential in the denominator (exp (0.5(ξln 10)²)) is exactly the ratio of the average value to the median value of a log-normal distribution with standard deviation ξ.

There are restrictions on the plausible range of covariance coefficient ρ. At low redshifts, a multiple regression analysis on group catalog data from the SDSS (Wetzel et al. 2012) suggests that, at fixed halo mass, the average SFR rises by 0.23 dex per 1 dex increase in stellar mass. At high redshifts (e.g., z > 6), the covariance will be almost perfect (1 dex rise in SFR per 1 dex rise in stellar mass) because of the short times that galaxies have had to form stars. Between these two times, it is reasonable to assume that the covariance increases with increasing redshift, so we adopt the following formula to interpolate between z ∼ 0 and z ≳ 6:

$$\begin{equation} \rho (a) = 1 + (4\rho _{0.5} -3.23)a + (2.46-4\rho _{0.5})a^2. \end{equation} \tag{ C12 }$$

This formula for ρ(a) gives ρ(0) = 1, ρ(1) = 0.23, and ρ(0.5) = ρ_0.5, where ρ_0.5 is a free parameter in our analysis.

By comparison, the total (cosmic) SFR density is straightforward to calculate; it is simply the sum of all star formation in all halos:

$$\begin{equation} {\rm CSFR}_\mathrm{true} = \int _{-\infty }^{\infty } {\rm SFR}(M_h) \frac{dN}{d\log M_h} d\log M_h. \end{equation} \tag{ C13 }$$

To calculate the measured quantities (as opposed to true quantities), we adopt the same conversion between measured stellar mass and true stellar mass as for the SMF. However, we do not impose additional corrections for active galaxies; by definition, because recent star formation is more evident in such galaxies, it is somewhat easier to recover the SFR, as opposed to the overall stellar mass. Therefore, we only apply a uniform systematic correction and do so in the same way as for stellar masses:

$$\begin{equation} {\rm SFR}_\mathrm{meas}(M_h) = {\rm SFR}_\mathrm{true}10^{\mu }. \end{equation} \tag{ C14 }$$

To calculate SSFR_meas and CSFR_meas, we would therefore replace SFR_true with SFR_meas in Equations (C11) and (C13). There is one more complicating factor; namely the burstiness/dustiness of high-redshift galaxies. To account for this, we add additional factors to calculate the observed SSFR and CSFR:

$$\begin{eqnarray} {\rm SSFR}_\mathrm{obs}(M_\ast) & = & \frac{1}{1-b(1-c(z))} {\rm SSFR}_\mathrm{meas}(M_\ast)\nonumber \\ {\rm CSFR}_\mathrm{obs} & = & [1 - (1-b)(1-c(z))] {\rm CSFR}_\mathrm{meas}. \end{eqnarray} \tag{ C15 }$$

A double power law has been used to fit the stellar mass (or luminosity) to halo mass relation at least since Yang et al. (2003). While using it is simple and convenient, the quantity and quality of galaxy observations have improved dramatically since its first use. To test its continued usefulness, we have used abundance matching to extract the implied SMHM relation from our low-redshift SMFs (Moustakas et al. 2013; Baldry et al. 2008). To account for the presence of scatter in observed stellar mass at fixed halo mass, we use a method based on Richardson–Lucy deconvolution to extract the median stellar mass as a function of halo mass. We assume a combined scatter of 0.2 dex.

The deconvolved SMHM relation is shown in Figure 20. The figure demonstrates that a double power-law can approximate the relation, but can be discrepant by up to 0.1 dex. Adding one more parameter improves the fit substantially, to a mean offset of less than 0.025 dex in our case. The discrepancies in the SMHM relation translate directly to discrepancies in the assumed SMF. This may be seen in the right-hand panel of Figure 20; the fits for the SMHM relations have been used to populate the halo mass function and reconvolved with 0.2 dex scatter to determine the resulting SMFs. While the five-parameter fit that we use is largely within the error bars of the observational data, the double power-law fit gives results that are again discrepant by up to 0.1 dex or 4σ from the observed SMF.

Zoom In Zoom Out Reset image size

It is tempting to ignore this issue: after all, the size of the systematic uncertainties in stellar mass is significantly larger (∼0.25 dex). However, if it is indeed the case that the deviations from the power law are driven by systematic uncertainties, then part of the value of the physical interpretation of the double power law is lost: the high-mass and low-mass slopes are now contaminated by using the incorrect parameterization. In our five-parameter fit, the high-mass dependence and low-mass slope are somewhat more insulated from unparameterized systematic errors in the transition region, leading to a more fair estimation of the high- and low-mass behavior of the SMHM relation.

Naturally, no fit will perfectly reproduce the original function. Special caution is necessary when using a MCMC method where a parameterization to tightly constrained data is accompanied by additional systematic or nuisance parameters. This is because the MCMC method may try to recruit the additional parameters to improve the fit, rather than for their intended purpose (i.e., to cover the range of allowed uncertainties). As shown in Figure 20, our chosen parameterization can reproduce the local SMF to within 0.02 dex in the vast majority of cases. Hence, to avoid overfitting bias, we reduce the magnitude of errors by 0.02 dex at each data point before calculating the χ² value used with our MCMC method. At z < 0.2, this has the impact of decoupling the systematic uncertainty parameters from the detailed shape of the SMF; at higher redshifts, there is no effect, as the typical size scale of the errors for observational data is already 0.1–0.2 dex.

As a check of the accuracy of our analysis pipeline, it is instructive to compare with a direct abundance matching approach. In Figure 21, we show results for both direct abundance matching (left panel) and abundance matching including an assumed intrinsic (ξ) and measurement (σ) scatter that match our assumed fiducial defaults. Notably, several features present in our main analysis are evident in this simpler approach:

• 1.  
  The SMHM relation for low-mass halos (M_h < 10¹¹ M_☉) at z = 0 deviates from a power law.
• 2.  
  The halo mass at which the integrated star formation is most efficient rises by 0.5 dex from 10¹² M_☉ at z = 0 to 10^12.5 M_☉ at z = 2.
• 3.  
  Above z = 2 it is difficult to locate the halo mass at which the peak integrated efficiency occurs because surveys at high redshifts do not have sufficient volume to probe the exponential tail of the SMF (if one exists). If the UV luminosity function is any indication, however, the z = 8 results suggest that the peak may occur at lower masses above z = 2.
• 4.  
  There appears to be a turnaround in the SMHM relation near z ∼ 2. The integrated efficiency of M_h ∼ 10^11.5 M_☉ halos falls from z = 0 to z = 2, but then increases at higher redshifts.

Zoom In Zoom Out Reset image size

At high redshifts, it is difficult to tell what, if any, evolution there is in the shape of the SMHM relation. However, it is clear that either the efficiency or the characteristic halo mass scale is evolving. If we constrain the shape parameters in our fitting formula for the SMHM relation (α, δ, and γ) to the z = 0 values, we find that high redshifts strongly prefer a halo mass scale that scales as a power law with redshift, but have comparatively little implied evolution in the efficiency at high redshifts. For that reason, when we extrapolate to high redshifts unconstrained by data in our main analysis (z > 8.5), we assume that , α, δ, and γ do not evolve beyond z = 8.5, but that the characteristic halo mass M₁ continues to evolve according to the redshift scaling we have adopted.

For CSFRs, we selected papers from the astro-ph arXiv that

• 1.  
  Were posted within the past 6 yr (2006–2012).
• 2.  
  Contained the word "cosmic" and either the word "SFR" or the words "star" and "formation" in the title or abstract.
• 3.  
  Contained an estimate of the CSFR from observations.
• 4.  
  Were not superseded by later publications.
• 5.  
  Contained a reference to the IMF used.

We have in addition excluded constraints from long gamma-ray burst SFR estimates for this work (e.g., Kistler et al. 2009); we feel that not enough data points are available at this time to independently calibrate the evolution of known systematic biases (e.g., higher GRB rates in low-luminosity galaxies) with redshift. The selected papers are summarized in Table 4. From these, we calculate the interpublication variance in five redshift bins with 16 data points each. To take into account variation of the CSFR in bins, we subtract a double power-law fit to the data points in Figure 2:

$$\begin{equation} {\rm CSFR}(z) = \frac{C}{10^{A(z-z_0)} + 10^{B(z-z_0)}}, \end{equation} \tag{ F1 }$$

where A, B, C, and z₀ are constants given in Table 6; adding more parameters did not make the quality of the fit substantially better. Our results are given in Table 7; we find the average systematic error to range from 0.13 dex at z = 0 to 0.27 dex for z > 3.

For SSFRs we conducted a very similar literature search. In this case, the selection criteria on astro-ph were papers that were posted within the past 6 yr (2006-2012); contained the word "specific" and either the word "SFR" or the words "star" and "formation" in the title or abstract; contained an estimate of the SSFR as a function of stellar mass, both as calculated from observations; were not superseded by later publications; contained a reference to the IMF used; and did not limit themselves to a specific class of galaxies (submillimeter, ULIRGs, satellites, star-forming only, etc.). The selected papers are summarized in Table 5.

To estimate uncertainties, we interpolate data for each publication's tabulated redshift ranges to obtain (if possible) estimates of the SSFR at four stellar masses (10⁹, 10^9.5, 10¹⁰ and 10^10.5 M_☉). We divide the SSFR data into the same five redshift bins as for the CSFR data and estimate the systematic errors from the interpublication variance. To exclude the effects of evolution in the SSFR over a redshift bin, we again subtract a double power-law fit (Equation (F1)) with parameters determined separately for each stellar mass in Table 8; the quality of the fit is shown in Figure 22. It is worth remarking that Weinmann et al. (2011) conduct a similar analysis of SSFRs at M_* = 10^9.5 M_☉ and find a "plateau" at high redshifts; while the data are consistent with a plateau, they are also consistent with not having a plateau: indeed, any range of shallow slopes at high redshifts would be consistent with the data collected here.

Zoom In Zoom Out Reset image size

Our results for the interpublication variance as a function of stellar mass and redshift are shown in Table 9. Because only UV surveys can probe low stellar masses at high redshifts (z > 3), the published SSFRs for low stellar masses at high redshift are more in agreement (0.22 dex scatter) than for low stellar masses at low redshift (0.32 dex scatter)! This strongly suggests that the interpublication variance in these bins does not represent the true systematic uncertainties. Excluding high-redshift bins (z > 3), there are no obvious systematic trends in the interpublication variance, either with redshift or stellar mass. Indeed, a one-way analysis of variance test delivers no strong evidence that the interpublication variances differ across bins (p = 0.6). This conclusion is not affected by merging bins by stellar mass or by redshift, and the size of the effect is too small to be measurable without adding more data. We adopt the fiducial model of constant errors across stellar masses and redshifts, where the error is determined by merging all redshift and stellar mass bins except for the high-redshift bins (z > 3) which are known underestimates, we would estimate the average systematic error in SSFRs to be 0.28 dex.

Finally, we mention a commonly overlooked systematic error in optical SFRs, which comes from the bias introduced by separately fitting the dust model for each galaxy. As the dust model for each galaxy can never be derived exactly, individual galaxy SFRs will have some scatter relative to the true galaxy SFRs. One way to estimate the shape of this scatter is to compare two different SFR estimates with different assumptions; for example, Salim et al. (2007) finds log-normal scatter consistent with intrinsic scatter of 0.3 dex. In calibrating one estimator to another (or to observations), authors typically make sure that the median offsets are zero; however, they often fail to notice that the expected mean offset is nonzero because the log-normal distribution is skewed in linear space. Specifically, a log-normal distribution with median value 1 and scatter σ (in dex) has an expected mean value of

$$\begin{equation} \langle x \rangle = \exp (0.5(\sigma \ln 10)^2). \end{equation} \tag{ F2 }$$

For a log-normal distribution with 0.3 dex scatter, this corresponds to an overall bias of +0.1 dex. This correction does not apply to UV studies that apply a median dust correction to all galaxies or to those that already correct for this bias.

The halo mass function for distinct halos presented by Tinker et al. (2008) was calibrated for a wide array of simulations from z = 0 to z = 2.5. The present work requires a mass function from z = 0 to z > 8 including satellite halos, so some additional corrections must be made. In addition, as new studies (Reddick et al. 2012) suggest that abundance matching to halo properties at their peak value before stripping offers a better match to z = 0 clustering data, calibration for the quantities of interest using M_peak must also be applied.

To test the high-redshift calibration of the mass function in Tinker et al. (2008), we integrate the mass function in two simulations (Bolshoi and Consuelo) to 10^11.5 M_☉ and compare to the Tinker mass function for the appropriate cosmology. A mass-dependent correction is necessary to normalize the z = 0 mass function for Consuelo to the Tinker mass function to account for slight incompleteness in the Consuelo mass function at this level:

$$\begin{equation} \log _{10}\left(\frac{\phi _\mathrm{Consuelo}(M)}{\phi _\mathrm{Tinker}(M)}\right) = -0.0008 -0.042\left(\frac{M}{10^{11.51}\,M_{\odot }}\right)^{-0.69}. \end{equation} \tag{ G1 }$$

The net effect of this correction on Φ(10^11.5 M_☉) is an increase of 0.03 dex (7%) compared to the raw simulation result for Consuelo. No correction is necessary for the Bolshoi simulation. As shown in Figure 23, the redshift fit given in Tinker et al. (2008) matches Bolshoi to within 0.02 dex up to a redshift of 5–6. However, Consuelo diverges from the Tinker mass function by 0.06 dex (15%) already at z = 3. This discrepancy is due to the different initial conditions calculations used in the two simulations: Consuelo uses 2LPT, whereas Bolshoi and the Tinker MF use the Zel'dovich approximation (1LPT). The Zel'dovich approximation has been previously found to underestimate the early nonlinear collapse of overdense regions (Crocce et al. 2006; Jenkins 2010), which is reconfirmed by the comparison in Figure 23. For that reason, we prefer to use the correction from Consuelo, parameterized as follows:

$$\begin{equation} \log _{10}\left(\frac{\Phi _\mathrm{true}(10^{11.5}\,M_{\odot })}{\Phi _\mathrm{Tinker}(10^{11.5}\,M_{\odot })}\right) = \frac{0.144}{1+\exp \left[14.79(a-0.213)\right]}. \end{equation} \tag{ G2 }$$

As shown in Figure 23, this fit matches the correction from the Consuelo simulation to within 1%. However, as also shown in Figure 23, there is a slight mass dependence to the correction; the magnitude of the correction is larger for more massive halos, as their nonlinear collapse begins earlier. We parameterize the mass dependence as follows:

$$\begin{equation} \log _{10}\left(\frac{\Phi _\mathrm{true}(M)}{\Phi _\mathrm{Tinker}(M)}\right) = \log _{10}\left(\frac{\Phi _\mathrm{true}(10^{11.5}\,M_{\odot })}{\Phi _\mathrm{Tinker}(10^{11.5}\,M_{\odot })}\right) \left(\frac{M}{10^{11.5}\,M_{\odot }}\right)^\frac{0.5}{1+\exp (6.5a)}. \end{equation} \tag{ G3 }$$

This function is well behaved, even when extrapolated down to low masses, as shown in Figure 23.

Zoom In Zoom Out Reset image size

Parameterizing the corrections as we have done introduces several uncertainties. The use of only a single simulation to calibrate the mass function introduces sample variance uncertainties, but these are in fact small for a box as large as Consuelo (420 Mpc h⁻¹); at low redshifts, the median offset from the Tinker mass function is 0.2% (see Equation (G1)), and even at redshift 6, the Poisson and sample variance errors combined are on the order of 1% (Trenti & Stiavelli 2008) for Φ(10^11.5 M_☉). Much larger uncertainties come from the use of a mass-dependent incompleteness correction to the Consuelo mass function, because uncertainties in this correction are degenerate with the mass dependence of the normalization correction in Equation (G3). If the incompleteness correction needed to be larger at high redshift, the evidence for the mass-dependence in the normalization correction would weaken. We can take this as an estimate of the additional systematic uncertainties on the high-redshift mass function. In terms of the impact on the halos in which high-redshift galaxies (z = 7–8) are expected to reside, the number density of such objects would be underestimated by at most 0.05 dex (12%), according to Figure 23, if there were no mass dependence in the normalization correction. This is tiny in comparison to the uncertainties in the stellar masses and SFRs at high redshift. Indeed, even if the real universe corresponded exactly to the Tinker mass function, the total error introduced by using the corrections in this section (0.1 dex, or 26%) would again be dwarfed by uncertainties in measuring galaxy stellar content.

In terms of the satellite fraction, we determine fits using a simple, physically motivated model. It is an oft-observed feature of dark matter substructure that the number of satellite halos at a given mass scales approximately with the mass of the host—the number density of satellites in hosts of mass M_c is proportional to M_cϕ_c(M_c), where ϕ_c(M) is the mass function for central halos. Another feature of dark matter substructure is that the number of satellites at a given satellite-to-host mass ratio is largely self-similar across host halo masses; this requires that the proportionality must scale inversely with the satellite mass. Therefore, we can write down a simple formula for the satellite mass function, ϕ_s:

$$\begin{equation} \phi _s(M) \sim \frac{C(a)}{M}\! \int _M^\infty\! M_\mathrm{c} \phi _\mathrm{c}(M_\mathrm{c})d\log _{10}M_\mathrm{c}, \end{equation} \tag{ G4 }$$

where the proportionality factor, C(a), is a function of scale factor only. The low-mass limit of this function takes on an especially appealing form; as $\phi _c(M_\mathrm{c}) \sim \phi _0 M_c^{-1}$ for M_c is less than the exponential cutoff scale, M_cutoff, we find

$$\begin{equation} \frac{\phi _s(M)}{\phi _c(M)} \sim \frac{\frac{C(a)}{M} \int _M^{M_\mathrm{cutoff}} \phi _0 d\log _{10} M_c}{\phi _0 M^{-1}}, \end{equation} \tag{ G5 }$$

which can be simplified to

$$\begin{equation} \frac{\phi _s(M)}{\phi _c(M)} \sim C(a) \log _{10}\left(\frac{M_\mathrm{cutoff}(a)}{M}\right). \end{equation} \tag{ G6 }$$

Thus, the satellite fraction grows only logarithmically in the low-mass limit.

However, the picture in Equation (G4) is not entirely complete; because they have later assembly times, massive clusters are more likely to be undergoing major mergers than smaller halos. However, these major mergers also become absorbed more quickly, as the effects of dynamical friction are more pronounced in massive halos. These two (opposite) effects combine to give a satellite fraction slightly larger at the high-mass end than Equation (G4) would predict; moreover, this excess evolves nontrivially with redshift. For that reason, we use the slightly more flexible form of Equation (G6) and allow great freedom in the parameterization for both C(a) and M_cutoff(a).

Fitting the satellite fraction to a simulation is made especially tricky by the fact that simulations are incomplete for satellites at higher masses than for central halos, even before selecting satellites on M_peak or M_acc. For that reason, Consuelo is unsuitable for determining the fraction of all but the most massive satellites (M ∼ 10¹³ M_☉), which provides an especially poor handle on constraining the redshift evolution of the fit. Bolshoi fares better, being complete for satellites down to 10¹¹ M_☉; yet, given its use of ZA initial conditions, it may not seem especially suitable for high-redshift calibration, either!

However, there are two arguments in its favor. First, the satellite fraction is low at high redshift, so even large relative errors in the satellite fraction will translate into small errors on the overall mass fraction. Second, the influence of the initial conditions on the satellite fraction are largely through the effects of assembly bias; however, recent work has shown that satellite richness (when tallied above a fixed cut in accretion or peak mass) is almost uncorrelated with formation time (Wu et al. 2013). The impact of the initial conditions on the satellite fraction is thus minor compared to the effect of the mass function normalization.

Therefore, using the Bolshoi simulation, we fit the satellite fraction from z = 0 to z = 10; when satellite masses are determined by M_acc, we find

$$\begin{equation} \log _{10}(C(a)) = -2.69+11.68 a-28.88 a^2+29.33 a^3-10.56 a^4 \end{equation} \tag{ G7 }$$

$$\begin{equation} \log _{10}(M_\mathrm{cutoff}(a)) = 11.34+8.34 a-0.36 a^2-5.08 a^3+0.75 a^4 \end{equation} \tag{ G8 }$$

and find the following parameterization for the satellite fraction when satellite masses are determined by M_peak:

$$\begin{equation} \log _{10}(C(a)) = -1.91+6.23 a-15.07 a^2+15.02 a^3-5.29 a^4 \end{equation} \tag{ G9 }$$

$$\begin{equation} \log _{10}(M_\mathrm{cutoff}(a)) = 10.66+15.93 a-21.39 a^2+18.20 a^3-8.21 a^4. \end{equation} \tag{ G10 }$$

Both of these fits in comparison to satellite fractions in Bolshoi are shown in Figure 24; they are accurate to within the estimated uncertainties in the mass function (5%, Tinker et al. 2008).

Zoom In Zoom Out Reset image size

The most commonly used fitting formula for halo mass accretion rates is that proposed by Wechsler et al. (2002), which suggests that halo mass growth is an exponential function in redshift (M(z)∝exp (− cz)). As simulations improved, deviations from this simple formula were found, which led authors to propose more complicated formulae; e.g., M(z)∝(1 + z)^bexp (− cz) in McBride et al. (2009). However, little attention has been given to the high-redshift behavior of these formulae, with the predictable result that they fail to accurately capture median mass accretion histories for halos as recently as z ∼ 3 (see Figure 25, left panel).¹⁰

Zoom In Zoom Out Reset image size

This limitation a problem for this study. To motivate our choice of improved fitting formulae, we note an issue with all proposed fitting formulae to date, which is their exponential dependence at high redshifts. That is to say, the mass ratio of two halos M₁ and M₂ will asymptote to

$$\begin{equation} \frac{M_1(z)}{M_2(z)} \rightarrow \frac{M_1(0)}{M_2(0)} \exp ((c_2-c_1)z). \end{equation} \tag{ H1 }$$

This is problematic because more massive halos usually have steeper exponential dependencies (i.e., M₁ < M₂ implies c₁ < c₂). If M₁(z = 0) < M₂(z = 0), these formulae then predict that eventually the progenitor histories will cross and M₁(z) > M₂(z). Using the McBride et al. (2009) formula, for example, with the MultiDark simulation, we find that it predicts the progenitors of 10¹⁵ M_☉ halos to be less massive than the progenitors of 10^14.75 M_☉ halos as early as z = 5. However, this is not the case: in simulations, the hierarchy of average or median progenitor masses is maintained, even though the distance between them may decrease (see, for example, the right panel of Figure 25).

We therefore parameterize the mass accretion history of a single halo mass, and then parameterize the ratios of progenitor masses for other halo masses, which entirely avoids the problem of crossing progenitor mass histories. The mass bin we choose (10^13.15–10^13.40 M_☉) is well-resolved in all three dark matter simulations at our disposal, and sufficient statistics are available even in the smallest simulation (Bolshoi) to be able to reduce scatter to an acceptable level. Our resulting fit is

$$\begin{equation} M(M_0,z) = M_{13}(z) 10^{f(M_0,z)} \end{equation} \tag{ H2 }$$

$$\begin{equation} M_{13}(z) = 10^{13.276} (1+z)^{3.00} \left(1+\frac{z}{2}\right)^{-6.11}\exp (-0.503z)\,M_{\odot } \end{equation} \tag{ H3 }$$

$$\begin{equation} f(M_0,z) = \log _{10}\left(\frac{M_0}{M_{13}(0)}\right)\frac{g(M_0,1)}{g(M_0,\frac{1}{1+z})} \end{equation} \tag{ H4 }$$

$$\begin{equation} g(M_0,a) = 1 + \exp (-4.651(a-a_0(M_0)) \end{equation} \tag{ H5 }$$

$$\begin{equation} a_0(M_0) = 0.205 - \log _{10}\left[\left(\frac{10^{9.649}\,M_{\odot }}{M_0}\right)^{0.18} + 1\right], \end{equation} \tag{ H6 }$$

where M(M₀, z) gives the median virial mass (Bryan & Norman 1998) for progenitors of halos with mass M₀ at z = 0. Despite its complexity, this seven-parameter fit (three parameters for M₁₃(z), one parameter for g(M₀, a), and three parameters for a₀(M₀)) provides an excellent fit to mass accretion histories for a wide range of halo masses and redshifts, as shown in Figure 25.

Halo mergers are important for our study because they impact the rate at which stellar mass is transferred from low-mass galaxies to high-mass galaxies as well as the rate at which stars are ejected into the ICL. While there have been many studies of the rate at which halos become subhalos (see, e.g., Fakhouri et al. 2010) in dark matter simulations, these studies are less relevant to our work because it is expected that satellite galaxies should last at least as long as the corresponding subhalo that hosts them (e.g., Yang et al. 2012; Reddick et al. 2012). Studying the subhalo destruction rate is more difficult because of issues with halo finding completeness when subhalos pass close to the centers of their hosts (Behroozi et al. 2013b, 2013c), and relatively few authors have attempted to correct for this issue in their work (c.f., Wetzel & White 2010).

However, with the advantage of a phase-space halo finder and a merger tree algorithm that corrects for such inconsistencies, we are in a unique position to provide a calibration of the merger rate (or disruption rate) of subhalos into host halos all the way back to the earliest redshifts in simulations. We find that the functional form proposed in Fakhouri et al. (2010) is an excellent fit to the subhalo disruption rate; however, the parameter values and redshift evolution are different due to the different merger definition, cosmology, and mass definition used.

Letting M_h be the mass of the host halo and θ be the ratio of subhalo peak mass to host mass, we find that the disruption rate (expressed in terms of the decrement in number of subhalos per unit host halo per unit redshift per log interval in subhalo mass ratio) is

$$\begin{equation} {-}\frac{d^2N}{dz d\log _{10}\theta }(M,\theta,z) = A(M,(1+z)^{-1})\theta ^{-1}\exp (b(M,\theta)) \end{equation} \tag{ I1 }$$

$$\begin{equation} A(M,a) = 0.0316 \left(1 + \frac{0.20\ln (a)}{1+a}\right) \left(\frac{M}{10^{12}\,M_{\odot }}\right)^{0.03+0.05a} \end{equation} \tag{ I2 }$$

$$\begin{equation} b(M,\theta) = \left[10^{1.929} \left(\frac{M}{10^{12}\,M_{\odot }}\right)^{0.11}\theta \right]^{0.2586}. \end{equation} \tag{ I3 }$$

If θ is instead the ratio of subhalo accretion mass to host halo mass, we find instead a slightly simpler dependence:

$$\begin{equation} A_\mathrm{acc}(M,a) = 0.0139 \left(1 + \frac{0.25\ln (a)}{1+a}\right) \left(\frac{M}{10^{12}\,M_{\odot }}\right)^{0.052+0.072a} \end{equation} \tag{ I4 }$$

$$\begin{equation} b_\mathrm{acc}(M,\theta) = [10^{3.091} \theta ]^{0.1965}. \end{equation} \tag{ I5 }$$

We show the quality of the fits for the peak subhalo mass in Figure 26. For ease of comparing simulation data to the fit, we have scaled out host mass and redshift dependencies according to Equations (I1)–(I3), so that all plotted lines are directly comparable to the subhalo disruption rate for M = 10¹² M_☉ at z = 0.

Zoom In Zoom Out Reset image size

Figure 26 shows that scaling in both mass and redshift of the subhalo disruption rate has been successfully captured in Equations (I1)–(I3). The remaining deviations are due to Poisson noise in the disruption rates for halo masses where Bolshoi has few halos (especially in the highest mass bin) and stochasticity in halo mass estimations for low-mass halos. Specifically, as the peak mass captures the highest mass ever recorded in a halo's mass accretion history, greater variance in the halo mass estimator translates to artificially increased peak masses for low-mass halos. This manifests as a slight upturn in the subhalo disruption rate close to the halo mass resolution limit of the simulation.

We consider four models with alternate sets of assumptions to the main analysis in this paper.

• 1.  
  A model ("No Systematics") that excludes all nuisance parameters related to errors in converting luminosities into stellar masses (i.e., μ = κ = 0, η_cal = ∞) and to incompleteness (A = 0).
• 2.  
  A model ("No z > 6 Constraints") that excludes all data constraints at high redshift (z > 6), including data in Bouwens et al. (2011), Bradley et al. (2012), Schaerer & de Barros (2010), McLure et al. (2011), and the data points above z = 6 in Bouwens et al. (2012).
• 3.  
  A model ("No CSFR Constraints") that excludes all constraints from observed CSFRs.
• 4.  
  A model ("No SFR Constraints") that excludes all constraints from observed SFRs (both CSFRs and SSFRs) as well as excluding all systematic nuisance parameters as in the first alternate model ("No Systematics").

Comparisons with observed data are presented in Figure 27. For ease of comparison across mass ranges, we show comparisons between derived baryon conversion efficiencies and stellar mass histories for the main model in this paper and the alternate models in Figures 28 and 29.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

First, it may be verified that the model including no systematics succeeds remarkably well in matching all three kinds of observations; while it is potentially somewhat high at low redshifts for the CSFR, this is not unexpected on account of issues with the galaxy/ICL definition (Section 2.4). The No Systematics model has much tighter error bars on the stellar mass to halo mass ratios (not shown) than the full model at all redshifts, which suggests that systematics remain the single largest source of errors in determining stellar masses. Figure 28 is comparable to Figure 5 for the No Systematics model, except that at high redshifts and low halo masses, there is an upturn in the SFR that is not present in the full model.

The model excluding z > 6 data gives almost identical results as the full model for z ⩽ 6, which is encouraging. The "predictions" for high-redshift SMFs are somewhat lower than in the best-fitting model, which results in a higher SFR (i.e., rate of stellar mass growth) over the period from z ∼ 7 to z ∼ 6.

The model excluding all cosmic star formation constraints is remarkably similar to the full model in our analysis. This is as expected: the specific SFR in combination with the SMF gives a constraint on the total amount of star formation. Yet, the CSFR data has tighter error bars than the SSFR in combination with the SMF, at least for z < 1 (see Sections 3.2 and 3.3), so it provides slightly better constraints at those redshifts.

Excluding all SFR data results in a good fit to the CSFR at high redshifts, where growth in stellar mass is rapid—i.e., where SMFs at successive redshifts are significantly different even despite large errors. At lower redshifts, the predicted CSFR for the best-fitting model is discrepant from observations. At these redshifts, the growth in stellar mass is poorly constrained because the intrinsic SMFs are changing comparatively little, especially from z = 0.5 to z = 0. A small difference in the choice of SMF evolution at low redshifts therefore has a much larger impact on the relative error in SFRs. Thus, although the stellar mass to halo mass ratios are similar to the No Systematics model (and the abundance matching model in Appendix E), the SFRs derived using this model are unreliable.