A COMPREHENSIVE ANALYSIS OF UNCERTAINTIES AFFECTING THE STELLAR MASS–HALO MASS RELATION FOR 0 < z < 4

In this paper, we analyze seven main sources of uncertainties in stellar mass functions applicable to photometric surveys, listed below in rough order of importance.

• 1.  
  Choice of stellar IMF. The luminosity of stars scales as their mass to a large power (dlnL/dlnM ∼ 3.5), while the IMF, defined as the number of stars formed per unit mass, scales approximately as ξ ∝ M^−2.3 (Salpeter 1955). Thus, the light from a galaxy is dominated by the most massive stars, while the total stellar mass is dominated by the lowest mass stars. This fact, which has been known for decades (e.g., Tinsley & Gunn 1976), implies that the assumed form of the IMF at M ≲ 0.5 M_☉ will have no effect on the integrated light from galaxies, but will have a large effect on the total stellar mass. Nonetheless, constraints on the total dynamical mass of spheroidal systems provide a valuable independent check on the form of the IMF at low masses. Cappellari et al. (2006) find, for example, that the IMF proposed by Chabrier (2003) for the solar neighborhood is consistent with dynamical constraints on masses of nearby ellipticals.We do not marginalize over IMFs in our error calculations for the reason that it is relatively simple to convert from our choice of IMF (Chabrier 2003) to other IMFs. For our purposes, the IMF becomes a serious source of systematic uncertainty only if the IMF varies with galaxy properties or if it evolves. It should be noted, however, that the IMF can also introduce more complicated systematic effects associated with the inferred star formation rate (SFR), which may in turn impact stellar masses in non-trivial ways (e.g., Conroy et al. 2010).
• 2.  
  Choice of SPS model. SPS modeling efforts have grown substantially in sophistication in the past decade (e.g., Leitherer et al. 1999; Bruzual & Charlot 2003; Le Borgne et al. 2004; Maraston 2005); yet, significant uncertainties remain (e.g., Charlot et al. 1996; Charlot 1996; Yi 2003; Lee et al. 2007; Conroy et al. 2009). Treatments of convection vary, leading to different main sequence turn off times for intermediate mass stars. Advanced stages of stellar evolution—including blue stragglers, thermally pulsating AGB stars (TP-AGB), and horizontal branch stars—are poorly understood, both observationally and theoretically. The theoretical spectral libraries contain known deficiencies, especially for M giants, where the effective temperatures are low and where hydrodynamic effects become important. Empirical stellar libraries to some extent circumvent these issues, although they are plagued by incomplete coverage in the Hertzsprung–Russell diagram and difficulties associated with deriving stellar parameters. See Conroy et al. (2009) and Percival & Salaris (2009) for recent reviews.Different SPS models treat these issues differently, which can result in large systematic differences in the derived stellar mass. For instance, the model of Maraston (2005) compared to Bruzual & Charlot (2003) has systematic differences of 0.1 dex in stellar mass (Salimbeni et al. 2009; Pérez-González et al. 2008). However, Salimbeni et al. (2009) reports that the model in Bruzual (2007; with a revised treatment of TP-AGB stars) yields systematic differences in stellar mass relative to Bruzual & Charlot (2003), which ranges from 0.05 dex for 10¹¹ M_☉ galaxies to 0.3 dex for 10^9.5 M_☉ galaxies. Conroy et al. (2009) show that uncertainty in the luminosity of the TP-AGB phase can shift stellar masses by as much as ±0.2 dex.
• 3.  
  Parameterization of SFHs. In order to estimate stellar masses, model libraries are constructed with a large range in SFHs, dust attenuation (see below), and, often but not always, metallicity. The adopted functional form of the SFH is another source of systematic uncertainty, as typically very simple functional forms are assumed. Several authors have investigated various aspects of this problem. When attempting to model observed photometry, Pérez-González et al. (2008) found that a single stellar population model (in particular, star formation proportional to e^−t/τ, which is a commonly used parameterization) systematically underpredicts stellar mass by 0.18 dex compared to a double stellar population model (exponentially decaying star formation followed by a later starburst). The particular parameterization of SFHs may also lead to systematic differences as a function of stellar mass. Lee et al. (2009) analyzed a sample of mock Lyman-break galaxies at z ∼ 4–5 and found that simple SFHs produced best-fit stellar masses that were under- or overestimated by ∼±50% depending on the rest-frame galaxy color. This bias was attributed to the chaotic SFHs of the mock galaxies.
• 4.  
  Choice of dust attenuation model. Because dust reddens starlight, it is difficult to separate the effects of dust from stellar population effects, especially when fitting optical photometry. The effects of dust are also known to change depending on galaxy inclination (e.g., Driver et al. 2007). Hence, the choice of dust attenuation law has a non-trivial effect on the inferred stellar population ages and, consequently, SFHs derived from photometric and even spectroscopic surveys (Panter et al. 2007). In terms of the effect on stellar masses, Pérez-González et al. (2008) compared the dust models of Calzetti et al. (2000) and Charlot & Fall (2000), finding a systematic difference of 0.10 dex. Panter et al. (2007) found a similar difference between Calzetti et al. and models based on extinction curves from the Small and Large Magellanic Clouds. The effects of varying dust attenuation models have also been explored recently by Marchesini et al. (2009) and Muzzin et al. (2009), with similar results. We have used the stellar population fitting procedure described in Conroy et al. (2009) to compare the Calzetti et al. dust attenuation law to the dust model used in the kcorrect package (Blanton & Roweis 2007). We find a median offset of 0.02 dex but also a systematic trend such that two galaxies whose stellar masses are estimated with Calzetti dust attenuation and are separated by 1.0 dex will have kcorrect masses separated by (on average) only 0.92 dex.
• 5.  
  Statistical errors in individual stellar mass estimates. Stellar mass estimates for each galaxy are subject to statistical errors due to uncertainty in photometry as well as uncertainty in the SPS parameters for a given set of model assumptions. We treat this here as a random statistical error. While it may seem that random scatter in individual stellar masses should on average have no systematic effect, it in fact introduces a systematic error analogous to the Eddington bias (Eddington 1940) observed in luminosity functions. As the stellar mass function drops off steeply beyond a certain characteristic stellar mass, there are many more low stellar mass galaxies that can be upscattered than there are high stellar mass galaxies that can be downscattered by errors in stellar mass estimates. This asymmetric scatter implies that the drop-off in number density at high masses becomes shallower in the presence of scatter. We discuss this effect in detail in Section 3.1.2 (also see the Appendix in Cattaneo et al. 2008).
• 6.  
  Sample variance. Surveys of finite regions of the universe are susceptible to large-scale fluctuations in the number density of galaxies. This is no longer a dominant source of uncertainty for the volumes probed at low redshift by the SDSS, but it is an important consideration for higher-redshift surveys, which cover much smaller comoving volumes. Most authors who consider sample variance attempt to minimize it by averaging over several fields (e.g., Pérez-González et al. 2008; Marchesini et al. 2009). Very few surveys at z>0 attempt to estimate the magnitude of the error except by computing the field-to-field variance, which is often an underestimate when insufficient volume is probed (Crocce et al. 2010). We detail a more accurate method based on simulations to model the error arising from sample variance in Section 3.2.4.
• 7.  
  Redshift errors. Photometric redshift errors blur the distinction between GSMFs at different redshifts. While a galaxy may be scattered either up or down in redshift space, volume-limited survey light cones will contain larger numbers of galaxies at higher redshifts, meaning that the GSMF as reported at lower redshifts will be artificially inflated. Moreover, as galaxies at earlier times have lower stellar masses, surveys will tend to report artificially larger faint-end slopes in the GSMF. However, as these errors are well known, it is easy to correct for their effects on the stellar mass function, as has been done for the data in Pérez-González et al. (2008; see the Appendix of Pérez-González et al. 2005 for details on this process).

For completeness, we remark that galaxy–galaxy lensing will also result in systematic errors in the GSMF at high redshifts because galaxy magnification will result in higher observed luminosities. However, from ray-tracing studies of the Millennium simulation (Hilbert et al. 2007), the expected scatter in galaxy stellar masses from lensing is minimal (e.g., 0.04 dex at z = 1) compared to the other sources of scatter above (e.g., 0.25 dex from different model choices). For that reason, we do not model galaxy–galaxy lensing effects in this paper.

Recently, it has become clear that current estimates of the evolution in the cosmic SFR density are not consistent with estimates of the evolution of the stellar mass density at z>1 (Nagamine et al. 2006b; Hopkins & Beacom 2006b; Pérez-González et al. 2008; Wilkins et al. 2008a). The origin of this discrepancy is currently a matter of debate. One solution involves allowing for an evolving IMF with redshift (Davé 2008; Wilkins et al. 2008a). While such a solution is controversial, a number of independent lines of evidence suggest that the IMF was different at high redshift (Lucatello et al. 2005; Tumlinson 2007a, 2007b; van Dokkum 2008). Reddy & Steidel (2009) offer a more mundane explanation for the discrepancy. They appeal to luminosity-dependent reddening corrections in the ultraviolet luminosity functions at high redshift and demonstrate that the purported discrepancy then largely vanishes.

In contrast to results at z>1, there does seem to be an accord that for z < 1 both the integrated SFR and the total stellar mass are in good agreement if one assumes (as we have) a Chabrier (2003) IMF (see Wilkins et al. 2008b; Pérez-González et al. 2008; Hopkins & Beacom 2006a; Nagamine et al. 2006a; Conroy & Wechsler 2009).

Because of the discrepancy between reported SFRs and stellar masses in the literature, it is clear that estimates of uncertainties in GSMFs and SFRs at z>1 tend to underestimate the true uncertainties; for this reason, we separately analyze results for z < 1 in Section 4 and z>1 in Section 5 of this paper.

Most studies on the GSMF report Schechter function fits as well as individual data points; many also provide statistical errors. However, even when systematic errors are reported (either in Schechter parameters or at individual data points), the systematic error estimates are of limited value unless one is also able to model shifts in the GSMF caused by such errors.

Fortunately, based on the discussion in Section 2.1.1, there seem to be two main classes of systematic errors causing shifts in the GSMF.

• 1.  
  Over/underestimation of all stellar masses by a constant factor μ. This appears to cover the majority of errors, including most differences in SPS modeling, dust attenuation assumptions, and stellar population age models.
• 2.  
  Over/underestimation of stellar masses by a factor which depends linearly on the logarithm of the stellar mass (i.e., depends on a power of the stellar mass). This covers the majority of the remaining discrepancies between different SPS models and different stellar age models.

Both forms of error are modeled with the equation

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

Without loss of generality, we may take M₀ = 10^11.3 M_☉ (the fixed point of the variation between the Bruzual 2007 and Bruzual & Charlot 2003 models found by Salimbeni et al. 2009), allowing the prior on M₀ to be absorbed into the prior on μ.

For the prior on μ, we consider four contributing sources of uncertainty. We adopt estimates of the uncertainty from the SPS model (≈0.1 dex), the dust model (≈0.1 dex), and assumptions about the SFH (≈0.2 dex) from Pérez-González et al. (2008) as detailed in Section 2.1.1. Additionally, we have the variation in κ log₁₀(M₀) (at most 0.1 dex, as |κ| ≲ 0.15—see below). Assuming that these are statistically independent, they combine to give a total uncertainty of 0.25 dex, which is consistent with the accepted range for systematic uncertainties in stellar mass (Pérez-González et al. 2008; Kannappan & Gawiser 2007; van der Wel et al. 2006; Marchesini et al. 2009). For lack of adequate information (i.e., different models) to infer a more complicated distribution, we assume that μ has a Gaussian prior. As more studies of the overall systematic shift, μ, become available, our assumptions for the prior on μ and the probability distribution will likely need corrections. We remark, however, that our results can easily be converted to a different assumption for μ, as μ simply imparts a uniform shift in the intrinsic stellar masses relative to the observed stellar masses.

For the prior on κ, the result of Salimbeni et al. (2009) would suggest |κ| ≲ 0.15. As mentioned in Section 2.1.1, we found that |κ| ≈ 0.08 between the Blanton & Roweis (2007) and Calzetti et al. (2000) models for dust attenuation. Li & White (2009) finds |κ| ≲ 0.10 between Blanton & Roweis (2007) and Bell et al. (2003) stellar masses. Without a large number of other comparisons, it is difficult to robustly determine the prior distribution for κ; however, motivated by the results just mentioned, we assume that the prior on κ is a Gaussian of width 0.10 centered at 0.0.

We remark that some authors have considered much more complicated parameterizations of the systematic error. For example, Li & White (2009) consider a four-parameter hyperbolic tangent fit to differences in the GSMF caused by different SPS models, as well as a five-parameter quartic fit. However, we do not consider higher-order models for systematic errors for several reasons. First, given that second- and higher-order corrections will result only in very small corrections to the stellar masses in comparison to the zeroth-order correction (μ ≈ 0.25 dex), the corrections will not substantially affect the systematic error bars. Second, we do not know of any studies which would allow us to construct priors on the higher-order corrections. Finally, with higher-order models, there is the serious danger of over-fitting—that is, with very loose priors on systematic errors, the best-fit parameters for the systematic errors will be influenced by bumps and wiggles in the stellar mass function due to statistical and sample variance errors. Hence, the interpretive value of the systematic errors becomes increasingly dubious with each additional parameter.

In addition to the systematic effects discussed in the previous section, measurement of stellar masses is subject to statistical errors. Even for a fixed set of assumptions about the dust model, SPS model, and parameterization of SFHs, stellar masses will carry uncertainties because the mapping between observables and stellar masses is not one-to-one. This additional source of uncertainty has unique effects on the GSMF. Observers will see a GSMF (ϕ_meas) which is the true or "intrinsic" GSMF (ϕ_true) convolved with the probability distribution function of the measurement scatter. For instance, if the scatter is uniform across stellar masses and has the shape of a certain probability distribution P, we have

$$\begin{equation} \phi _\mathrm{meas}(M) = \int _{-\infty }^\infty \phi _\mathrm{true}(10^y) \, P(y - \log _{10}(M)) dy, \end{equation} \tag{ 2 }$$

where y is the integration variable, in units of log₁₀ mass. As derived in Appendix B, the approximate effect of the convolution is

$$\begin{equation} \log _{10} \left(\frac{\phi _\mathrm{meas}(M)}{\phi _\mathrm{true}(M)}\right) \approx \frac{\sigma ^2}{2}\ln (10)\,\left(\frac{d\log \phi _\mathrm{true}(M)}{d\log M}\right)^2, \end{equation} \tag{ 3 }$$

where σ is the standard deviation of P. That is to say, the effect of the convolution strongly depends on the logarithmic slope of ϕ_true. Where the slope is small (i.e., for low-mass galaxies), there is almost no effect. Above 10¹¹ M_☉, where the GSMF becomes exponential, there can be a dramatic effect, with the result that ϕ_true is more than an order of magnitude less than ϕ_meas because it becomes far more likely that stellar mass calculation errors produce a galaxy of very high perceived stellar mass than it is for there to be such a galaxy in reality (see, for example, Cattaneo et al. 2008).

For the observed z ∼ 0 GSMF, we take the probability distribution P to be log-normal with 1σ width 0.07 dex from the analysis of the photometry of low-redshift luminous red galaxies (LRGs; Conroy et al. 2009). Kauffmann et al. (2003) found similar results regarding the width of P. This function only accounts for the statistical uncertainties mentioned above and does not include additional systematic uncertainties. In light of Equation (3), we use LRGs to estimate P because LRGs occupy the high stellar mass regime where measurement errors are most likely to affect the shape of the observed GSMF. However, the single most important attribute of the distribution P is its width; the main results do not change substantially if an alternate distribution with non-Gaussian tails beyond the 1σ limits of P is used.

For higher redshifts, we scale the width of the probability distribution to account for the fact that mass estimates become less certain at higher redshift (e.g., Conroy et al. 2009; Kajisawa et al. 2009):

$$\begin{equation} P(\Delta \log _{10} M_\ast, z) = \frac{\sigma _0}{\sigma (z)} P_0\left(\frac{\sigma _0}{\sigma (z)} \Delta \log _{10} M_\ast \right), \end{equation} \tag{ 4 }$$

where P₀ is the probability distribution at z = 0 (as discussed above), σ₀ is the standard deviation of P₀, and σ(z) gives the evolution of the standard deviation as a function of redshift.

Conroy et al. (2009) did not give a functional form for σ(z), but they calculate for a handful of massive galaxies that σ(z = 2) is ≈0.18 dex, as compared to σ(z = 0)≈ 0.07 dex. Kajisawa et al. (2009) performed a similar calculation (albeit with a different SPS model) of the distribution in several redshift bins; their results show gradual evolution for σ(z) out to z = 3.5 for high stellar mass galaxies consistent with a linear fit:

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

The results of Kajisawa et al. (2009) suggest that σ_z = 0.03–0.06 dex for LRGs. As this is consistent with the value of σ_z = 0.05 dex which would correspond to Conroy et al. (2009), we adopt the linear scaling of Equation (5) with a Gaussian prior of σ_z = 0.05 ± 0.015 dex.

Note that the effect of this statistical error on the stellar mass function is minimal below 10¹¹ M_☉, and therefore does not affect the SM–HM relation for halos below ∼10¹³ M_☉, as discussed in Section 4.2. While this scatter does have an effect on the shape of the stellar mass function for high-mass galaxies, the qualitative predictions we make from this analysis are generic to all types of random scatter.

For the principal simulation in this study (L80G), we used a pure dark matter N-body simulation based on Adaptive Refinement Tree (ART) code (Kravtsov et al. 1997; Kravtsov & Klypin 1999). The simulation assumed flat, concordance ΛCDM (Ω_M = 0.3, Ω_Λ = 0.7, h = 0.7, and σ₈ = 0.9) and included 512³ particles in a cubic box with periodic boundary conditions and comoving side length 80 h⁻¹ Mpc. These parameters correspond to a particle mass resolution of ≈3.2 × 10⁸ h⁻¹ M_☉. For this simulation, the ART code begins with a spatial grid size of 512³; it refines the grid up to 8 times in locally dense regions, leading to an adaptive distance resolution of ≈1.2 h⁻¹ kpc (comoving units) in the densest parts and ≈0.31 h⁻¹ Mpc in the sparsest parts of the simulation.

In this simulation, halos and subhalos were identified using a variant of the Bound Density Maxima algorithm (Klypin et al. 1999). Halo centers are located at peaks in the density field smoothed over a 24 particle SPH kernel (for a minimum resolvable halo mass of 7.7 × 10⁹ h⁻¹ M_☉). Nearby particles are classified as bound or unbound in an iterative process; once all the locally bound particles have been found, halo parameters such as the virial mass M_vir and maximum circular velocity V_max may be calculated. (See Kravtsov et al. 2004 for complete details on the algorithm.) The simulation is complete down to V_max ≈ 100 km s⁻¹, corresponding to a galaxy stellar mass of 10^8.75 M_☉ at z = 0.

The ability of L80G to track satellites with high mass and force resolution gives it several uses. Merger trees from L80G inform our prescription for converting analytical-central-only halo mass functions to mass functions which include satellite halos (see Section 3.2.2). Additionally, the merger trees allow for evaluation of different models of satellite stellar evolution with full consistency (see Section 3.3.2). Finally, the knowledge of which satellite halos are associated with which central halos allows for estimates of the total stellar mass (in the central and all satellite galaxies)–halo mass relation (see Section 4.3.6).

We also make use of a secondary simulation from the Large Suite of Dark Matter Simulations (LasDamas Project, http://lss.phy.vanderbilt.edu/lasdamas/) in our sample variance calculations. The L80G simulation is too small for use in calculating the sample variance between multiple independent mock surveys, but the larger size of the LasDamas simulation (420 h⁻¹ Mpc, 1400³ particles) makes it ideal for this purpose. However, the LasDamas simulation has poorer mass resolution (a minimum particle size of 1.9 × 10⁹ M_☉) and force resolution (8 h⁻¹ kpc), making it unable to resolve subhalos (particularly after accretion) as well as L80G. The LasDamas simulation assumes a flat, ΛCDM cosmology (Ω_M = 0.25, Ω_Λ = 0.75, h = 0.7, and σ₈ = 0.8) which is very close to the WMAP5 best-fit cosmology (Komatsu et al. 2009). Collisionless gravitational evolution was provided by the gadget-2 code (Springel 2005). Halos are identified using friends of friends with a linking length of 0.164. The subfind algorithm Springel (2005) is used to identify substructure.

As mentioned, the primary use of the LasDamas simulation is in sampling the halo mass functions in mock surveys to model the effects of sample variance on high-redshift pencil-beam galaxy surveys. The mock surveys are constructed so as to mimic the observations in Pérez-González et al. (2008). In each mock survey, three pencil-beam light cones (matching the angular sizes of the three fields in Pérez-González et al. 2008) with random orientations are sampled from a random origin in the simulation volume out to z = 1.3. Thus, by comparing the halo mass functions in individual mock surveys to the mass function of the ensemble, the effects of sample variance may be calculated with full consideration of the correlations between halo counts at different masses.

The analytic mass functions of Tinker et al. (2008) are used to calculate the abundance of halos in several cosmological models. We calculate mass functions defined by M_vir, using the overdensity specified by Bryan & Norman (1998).⁵ This results in an overdensity (compared to the mean background density) Δ_vir which ranges from 337 at z = 0 to 203 at z = 1 and smoothly approaches 180 at very high redshifts. Following Tinker et al. (2008), we use spline interpolation to calculate mass functions for overdensities between the discrete intervals presented in their paper.

The mass functions in Tinker et al. (2008) only include central halos. We model the small (≈20% at z = 0) correction to the mass function introduced by subhalos to first order only, as the overall uncertainty in the central halo mass function is already of order 5% (Tinker et al. 2008). In particular, we calculate satellite (mass at accretion) and central mass functions in our simulation (L80G) and fit Schechter functions to both, excluding halos below our completeness limit (10^10.3 M_☉). Then, we plot the difference between the Schechter parameters (the difference in characteristic mass, Δlog₁₀M_*; the difference in characteristic density, Δlog₁₀ϕ₀; and the difference in faint-end slopes, Δα) as a function of scale factor (a). This gives the satellite mass function (ϕ_s) as a function of the central mass function (ϕ_c), which allows us to use this (first-order) correction for central mass functions of different cosmologies:

$$\begin{equation} \phi _s(M) = 10^{\Delta \log _{10} \phi _0} \, \left(\frac{M}{M_0 \cdot 10^{\Delta \log _{10} M_0}}\right)^{-\Delta \alpha } \, \phi _c(M / 10^{\Delta \log _{10} M_0}). \end{equation} \tag{ 6 }$$

From our simulation, we find fits as shown in Figure 1:⁶

$$\begin{eqnarray} \Delta \log _{10} \phi _0(a) &=& {-}0.736 - 0.213 \, a, \nonumber \\ \Delta \log _{10} M_0(a) &=& 0.134 - 0.306 \, a,\\ \Delta \alpha (a) &=& {-}0.306 + 1.08 \, a - 0.570 \, a^2. \nonumber \end{eqnarray} \tag{ 7 }$$

The mass function used here may be easily replaced by an arbitrary mass function, as detailed in Appendix A.

Zoom In Zoom Out Reset image size

Our fiducial results are calculated assuming WMAP5 cosmological parameters. In order to model uncertainties in cosmological parameters, we have sampled an additional 100 sets of cosmological parameters from the WMAP5+BAO+SN MCMC chains (from the models in Komatsu et al. 2009) and generated mass functions for each one according to the method in the previous section. Hence, to determine the variance in the derived SM–HM relation caused by cosmology uncertainties, we recalculate the relation for each sampled mass function according to the method described in Appendix A.

Large-scale modes in the matter power spectrum imply that finite surveys will obtain a biased estimate of the number densities of galaxies and halos as compared to the full universe. That is to say, matching observed GSMFs measured from a finite survey to the halo mass function estimated from a much larger volume will introduce systematic errors into the resulting SM–HM relation. These errors cannot be corrected unless one has knowledge of the halo mass function for the specific survey in question, which is in general not possible.

However, we can still calculate the uncertainties introduced by the limited sample size. While we cannot determine the true halo mass function for the survey, we can calculate the probability distribution of halo mass functions for identically shaped surveys via sampling light cones from simulations. If we rematch galaxy abundances from the observed GSMF to the abundances of halos in each of the sampled light cones, then the uncertainty introduced by sample variance is exactly captured in the variance of the resulting SM–HM relations.

In detail, we create our distribution of halo mass functions by sampling one thousand mock surveys from the LasDamas simulation (see Section 3.2.1) corresponding to the exact survey parameters used in Pérez-González et al. (2008). We fit Schechter functions to the halo mass functions of each mock survey (over all redshifts), and we calculate the change in Schechter parameters (Δlog₁₀ϕ₀, Δlog₁₀M₀, and Δα) as compared to a Schechter fit to the ensemble average of the mass functions. Using the distribution of the changes in Schechter parameters, we may mimic to first order the expected distribution of halo mass functions for any cosmology. In particular, we use an equation exactly analogous to Equation (6) to convert the mass function for the full universe (ϕ_full) and the distribution of Δlog₁₀ϕ₀, Δlog₁₀M₀, and Δα into a distribution of possible survey mass functions (ϕ_obs):

$$\begin{eqnarray} \phi _{\rm obs}(M) &=& 10^{\Delta \log _{10} \phi _0} \, \left(\frac{M}{M_0 \cdot 10^{\Delta \log _{10} M_0}}\right)^{-\Delta \alpha } \nonumber\\ &&\times \phi _\mathrm{full}(M / 10^{\Delta \log _{10} M_0}). \end{eqnarray} \tag{ 8 }$$

Hence, to obtain the variance in the SM–HM relation caused by finite survey size, we recalculate the relation for each one of the survey mass functions thus computed according to the method described in Appendix A.

An important uncertainty in the abundance matching procedure is introduced by intrinsic scatter in stellar mass at a given halo mass. Suppose that M_*(M_h) is the average (true) galaxy stellar mass as a function of host halo mass. For a perfect monotonic correlation between stellar mass and halo mass, i.e., without scatter between stellar and halo mass, it is straightforward to relate the true or "intrinsic" stellar mass function (ϕ_true) to the halo mass function (ϕ_h) via

$$\begin{equation} \frac{dN}{d\log _{10} M_\ast } = \frac{dN}{d\log _{10} M_h} \, \frac{d\log M_h}{d\log M_\ast }, \end{equation} \tag{ 9 }$$

where N is the number density of galaxies, so that

$$\begin{equation} \phi _\mathrm{true}(M_\ast (M_h)) = \phi _h(M_h) \, \left(\frac{d\log M_\ast (M_h)}{d\log M_h}\right)^{-1}. \end{equation} \tag{ 10 }$$

Intuitively, as the halos of mass M_h get assigned stellar masses of M_*(M_h), the number density of galaxies with mass M_*(M_h) will be proportional to the number density of halos with mass M_h. The above equations are simply a mathematical representation of the traditional abundance matching technique.

Equation (10) remains useful in the presence of scatter. If we know the expected scatter about the mean stellar mass, say in the form of a probability density function P_s(Δlog₁₀M_*|M_h), then we may still relate ϕ_true to ϕ_h via an integral similar to a convolution:

$$\begin{eqnarray} \phi _\mathrm{true}(x) & =& \int _{0}^\infty \phi _h(M_h(M_\ast)) \, \frac{d \log M_h(M_\ast)}{d\log M_\ast } \nonumber \\ && {\times}\, P_s\left(\log _{10} \frac{x}{M_\ast } \bigg| M_h(M_\ast)\right) d\log _{10} M_\ast, \end{eqnarray} \tag{ 11 }$$

where M_h(M_*) is the inverse function of M_*(M_h).

This similarity to a convolution is no coincidence—mathematically, it is analogous to how we model random statistical errors in stellar mass measurements in Section 3.1.2. Namely, if we define ϕ_direct to equal the right-hand side of Equation (10),

$$\begin{equation} \phi _\mathrm{direct}(M_\ast) \equiv \phi _h(M_h(M_\ast)) \, \frac{d \log M_h}{d\log M_\ast }, \end{equation} \tag{ 12 }$$

and if we assume a probability density distribution independent of halo mass (i.e., scatter in stellar mass at fixed halo mass is independent of halo mass), then ϕ_true is exactly related to ϕ_direct by a convolution:

$$\begin{equation} \phi _\mathrm{true}(M_\ast) = \int _{-\infty }^\infty \phi _\mathrm{direct}(10^y) \, P_s(y - \log _{10} M_\ast) dy, \end{equation} \tag{ 13 }$$

which is mathematically identical to Equation (2) in Section 3.1.2.

Then, if one calculates ϕ_direct from ϕ_true, one may find M_h(M_*) via direct abundance matching. Namely, integrating Equation (12), we have

$$\begin{equation} \int _{M_h(M_\ast)}^\infty \phi _h(M)d\log _{10} M = \int _{M_\ast }^\infty \phi _\mathrm{direct}(M_\ast)d\log _{10} M_\ast. \end{equation} \tag{ 14 }$$

Equivalently, letting $\Phi _h(M_h) \equiv \int _{M_h}^\infty \phi _h(M) d\log _{10} M$ be the cumulative halo mass function, and letting $\Phi _\mathrm{direct}(M_\ast) \equiv \int _{M_\ast }^\infty \phi _\mathrm{direct}(M_\ast) d\log _{10} M_\ast$ be the cumulative "direct" stellar mass function, we have

$$\begin{equation} M_h(M_\ast) = \Phi _h^{-1}(\Phi _\mathrm{direct}(M_\ast)), \end{equation} \tag{ 15 }$$

and one may similarly find M_*(M_h) by inverting this relation.

Our approach in all equations except for Equation (13) allows a halo mass-dependent scatter in the stellar mass, but to date the data appears to be consistent with a constant scatter value. For example, using the kinematics of satellite galaxies, More et al. (2009) find that the scatter in galaxy luminosity at a given halo mass is 0.16 ± 0.04 dex, independent of halo mass. Using a catalog of galaxy groups, Yang et al. (2009b) find a value of 0.17 dex for the scatter in the stellar mass at a given halo mass, also independent of halo mass. Here, we thus assume a fixed value for the scatter in stellar mass at fixed halo mass, ξ, to specify the standard deviation of P_s(Δlog₁₀M_*). As the Yang et al. (2009b) value is consistent with the More et al. (2009) value, we set the prior using the More et al. (2009) value and error bounds on ξ. We assume a Gaussian prior on the probability distribution for ξ, and we assume that the scatter itself is log-normal.

When a galaxy is accreted into a larger system, it will likely be stripped of dark matter much more rapidly than stellar mass because the stars are much more tightly bound than the halo. It has been demonstrated that various galaxy clustering properties compare favorably to samples of halos where satellite halos—i.e., subhalos—are selected according to their halo mass at the epoch of accretion, M_acc, rather than their current mass (e.g., Nagai & Kravtsov 2005; Conroy et al. 2006; Vale & Ostriker 2006; Berrier et al. 2006). These results support the idea that satellite systems lose dark matter more rapidly than stellar mass.

As commonly implemented (e.g., Conroy et al. 2006), the abundance matching technique matches the stellar mass function at a particular epoch to the halo mass function at the same epoch, using M_acc rather than the present mass for subhalos. As M_acc remains fixed as long as the satellite is resolvable, the standard technique implies that the satellite galaxy's stellar mass will continue to evolve in the same way as for centrals of that halo mass. Therefore, a subtle implication of the standard technique is that satellites may continue to grow in stellar mass, even though M_acc remains the same. A different model for satellite stellar evolution (e.g., in which stellar mass which does not evolve after accretion) would therefore involve different choices in the satellite matching process.

The fiducial results presented here use the standard model where satellites are assigned stellar masses based on the current stellar mass function and their accretion-epoch masses. However, we also present results for comparison in which satellite masses are assigned utilizing the stellar mass function at the epoch of accretion, corresponding to a situation in which satellite stellar masses do not change after the epoch of accretion. In order to maintain self-consistency for the latter method, we use full merger trees (from L80G, the simulation described in Section 3.2.1) to keep track of satellites and to assure that, e.g., mergers between satellites before they reach the central halo preserve stellar mass.

Finally, we note that any specific halo-finding algorithm may introduce artifacts in the halo mass function in terms of when a satellite halo is considered absorbed/destroyed. This can have a small effect on satellite clustering as well as number density counts. Wetzel & White (2010) suggest an approach that avoids some of the problems associated with resolving satellites after accretion. Namely, they suggest a model where satellites remain in orbit for a duration that is a function of the satellite mass, the host mass, and the Hubble time, after which time they dissolve or merge with the central object. Although we have not modeled this explicitly, our satellite counts are consistent with their recommended cutoff—they suggest considering a satellite halo absorbed when its present mass is less than 0.03 times its infall mass; in our simulation, only 0.1% of all satellites fall below this threshold.

In abundance matching, knowledge of the halo mass function and the stellar mass function uniquely determines the SM–HM relation. Hence, parameterizing the stellar mass function yields an indirect parameterization for the SM–HM relation as well. Numerous papers (e.g., Cole et al. 2001; Bell et al. 2003; Panter et al. 2004; Pérez-González et al. 2008) have found that the GSMF is well approximated by a Schechter function:

$$\begin{equation} \phi (M_\ast,z) = \phi ^\star (z)\left(\frac{M_\ast }{M(z)}\right)^{-\alpha (z)}\exp \left(-\frac{M_\ast }{M(z)}\right), \end{equation} \tag{ 16 }$$

where the Schechter parameters ϕ^⋆(z), M(z), and α(z) evolve as functions of the redshift z. In many previous works on abundance matching (e.g., Conroy et al. 2009), it is the Schechter function for the stellar mass function that sets the form of the SM–HM relation.

More recently, however, several authors have noted that the GSMF cannot be matched by a single Schechter function for z < 0.2 to within statistical errors (e.g., Li & White 2009; Baldry et al. 2008), in part because of an upturn in the slope of the GSMF for galaxies below 10⁹ M_☉ in stellar mass. It is possible that a conspiracy of systematic errors causes the observed deviations, but there is no fundamental reason to expect the intrinsic GSMF to be fit exactly by a Schechter function (see discussion in Appendix C). In any case, our full parameterization–either the stellar mass function or the error parameterization—must be able to capture all the subtleties of the observed stellar mass function. Hence, we are inclined to adopt a more flexible model than the Schechter function of Equation (16). Other authors, wrestling with the same problem, have chosen to adopt multiple Schechter functions, including the eleven-parameter triple piecewise Schechter-function fit used by Li & White (2009). While accurate, these models often add complication without increasing intuition.

Rather than attempting to parameterize the stellar mass function, we could use abundance matching directly to derive the SM–HM relation for the maximal-likelihood stellar mass function, and then find a fit which can parameterize the uncertainties in the shape of the relation. This process is complicated by the various errors which we must take into account. Recall from Equations (2) and (13) that

$$\begin{equation} \phi _\mathrm{meas}(M_\ast) = \phi _\mathrm{direct}(M_\ast) \circ P_s(\Delta \log _{10} M_\ast) \circ P(\Delta \log _{10} M_\ast), \end{equation} \tag{ 17 }$$

where "○" denotes the convolution operation, P_s is the probability distribution for the scatter in stellar mass at fixed halo mass, and P is the probability distribution for errors in observed stellar mass at fixed true stellar mass. However, if we obtain ϕ_direct by deconvolution of the observed stellar mass function ϕ_meas, we may use direct abundance matching (Equation (15)) to determine the maximum likelihood form of M_h(M_*).

Figure 2 shows the result of calculating M_h(M_*) at z ∼ 0.1 via deconvolution and direct matching of the stellar mass function as described in the previous section. We choose the maximum-likelihood value for the distribution function P_s (namely, 0.16 dex log-normal scatter), and we use WMAP5 cosmology for the halo mass function ϕ_h in the derivation.

Zoom In Zoom Out Reset image size

While deconvolution plus direct abundance matching gives an unbiased calculation of the relation, there are several problems which prevent it from being used directly to calculate uncertainties.

• 1.  
  Deconvolution will tend to amplify statistical variations in the stellar mass function—that is, shallow bumps in the GSMF will be interpreted as convolutions of a sharper feature.
• 2.  
  Deconvolution will give different results depending on the boundary conditions imposed on the stellar mass function (i.e., how the GSMF is extrapolated beyond the reported data points), the effects of which may be seen at the edges of the deconvolution in Figure 2.
• 3.  
  Deconvolution becomes substantially more problematic when the convolution function varies over the redshift range, as it does for our higher-redshift data (z>0.2).
• 4.  
  Deconvolution cannot extract the relation at a single redshift; instead, it will only return the relation averaged over the redshift range of galaxies in the reported GSMF.

For these reasons, we choose to find a fitting formula instead. In the discussion that follows, we fit M_h(M_*) (the halo mass for which the average stellar mass is M_*) rather than the more intuitive M_*(M_h) (the average stellar mass at a halo mass M_h) primarily for reasons of computational efficiency. From Equations (12) and (17), the calculation of what observers would see (ϕ_meas) for a trial SM–HM relation requires many evaluations of M_h(M_*) and no evaluations of M_*(M_h). If we had instead parameterized M_*(M_h), and then inverted as necessary in the calculation of ϕ_meas, our calculations would have taken an order of magnitude more computer time.

It is well known from comparing the GSMF (or the luminosity function) to the halo mass function that high-mass (M_* ≳ 10^10.5 M_☉) galaxies have a significantly different SM–HM scaling than low-mass galaxies, which is usually attributed to different feedback mechanisms dominating in high-mass versus low-mass galaxies. The transition point between low-mass and high-mass galaxies—seen as a turnover in the plot of M_h(M_*) around M_* = 10^10.6 M_☉ in Figure 2—defines a characteristic stellar mass (M_*,0) and an associated characteristic halo mass (M₁). Hence, we consider functional forms which respect this general structure of a low stellar mass regime and a high stellar mass regime with a characteristic transition point:

$$\begin{eqnarray*} \log _{10}(M_h(M_\ast)) &= &\log _{10}(M_1) \; \hbox{[characteristic halo mass]}\\ &+&\!\! f_{\rm low}(M_\ast / M_{\ast,0})\; \hbox{[low-mass functional form]}\\ & +&\!\! f_{\rm high}(M_\ast / M_{\ast,0}) \; \hbox{[high-mass functional form],} \end{eqnarray*}$$

where f_low and f_high are dimensionless functions dominating below and above M_*,0, respectively.

For low-mass galaxies (M_* < 10^10.5 M_☉), we find the SM–HM relation to be consistent with a power law:

$$\begin{eqnarray} \frac{M_h(M_\ast)}{M_1} & \approx & \left(\frac{M_\ast }{M_{\ast,0}}\right)^\beta,\hbox{ or} \nonumber \\ \log (M_h(M_\ast)) & \approx & \log (M_1) + \beta \log \left(\frac{M_\ast }{M_{\ast,0}}\right). \end{eqnarray} \tag{ 18 }$$

For high-mass galaxies, we find the SM–HM relation to be inconsistent with a power law. In particular, the logarithmic slope of M_h(M_*) changes with M_*, with dlog M_h/dlog M_* always increasing as M_* increases. This may seem like a small detail; after all, by eye, it appears that a power law could be a reasonable fit for high-mass galaxies in Figure 2. In addition, because previous authors (e.g., Moster et al. 2010; Yang et al. 2009a) have used power laws, it may not seem necessary to use a different functional form.

In order to explore this issue, we tried a general double power-law functional form for M_h(M_*) which parameterized a superset of the fits used in Moster et al. (2010) and (Yang et al. 2009a; in particular, the same form as in Equation (C2) in Appendix C). We found that this approach had two major problems common to any such power-law form.

• 1.  
  As the logarithmic slope of M_h(M_*) increases with increasing M_*, the best-fit power law for high-mass galaxies will depend on the upper limit of M_* in the available data for the GSMF. Thus, the best-fit power law will depend on the number density limit of the observational survey used rather than on any fundamental physics. Moreover, for studies such as this one which consider redshift evolution, the different number densities probed at different redshifts result in a completely artificial "evolution" of the best-fit power law.
• 2.  
  The best-fit power-law will not depend on the highest-mass galaxies alone; instead, it will be something of an average over all the high-mass galaxies. Because the logarithmic slope is increasing with M_*, this means that the best-fit power law for M_h(M_*) will increasingly underestimate the true M_h(M_*) at high M_*. Namely, the fit will underestimate the halo mass corresponding to a given stellar mass, and therefore (as lower-mass halos have higher number densities) will result in stellar mass functions systematically biased above observational values. However, a systematic bias in our functional form will influence the best-fit values of the systematic error parameters. The systematic bias caused by assuming a power law form turns out to be most degenerate with the scatter in stellar mass at fixed halo mass (ξ). As a result, for the MCMC chains which assumed a double power law form for M_h(M_*), the posterior distribution of ξ was 0.09 ± 0.02 dex, which just barely lies within 2σ of the constraints from More et al. (2009).

These problems are not as significant if one only considers the stellar mass function at a single redshift, or if one does not allow for the systematic errors which change the overall shape of the stellar mass function (κ, ξ, and σ(z)). However, we find that the issues listed above exclude the use of a power law for our purposes. Instead, we find that M_h(M_*) asymptotes to a sub-exponential function for high M_*: namely, a function which climbs more rapidly than any power-law function, but less rapidly than any exponential function. We find that high–mass galaxies (M_*>10^10.5 M_☉) are well fit by the relation

$$\begin{eqnarray} M_h(M_\ast) & \stackrel{\sim }{\propto }& 10^{(\frac{M_\ast }{M_{\ast,0}})^\delta }, \quad\hbox{ or} \nonumber \\ \log _{10}(M_h(M_\ast)) & \rightarrow & \log _{10}(M_1) + \left(\frac{M_\ast }{M_{\ast,0}}\right)^\delta, \end{eqnarray} \tag{ 19 }$$

where δ sets how rapidly the function climbs; δ → 0 would correspond to a power law, and δ = 1 would correspond to a pure exponential. Typical values of δ at z = 0 range from 0.5to0.6. It is not obvious what physical meaning can be directly inferred from the choice of a sub-exponential function—after all, the stellar mass of a galaxy is a complicated integral over the merger and evolution history of the galaxy—but it could suggest that the physics driving the M_h(M_*) relation at high mass is not scale-free.

Although this form now matches the asymptotic behavior for the highest and lowest stellar mass galaxies, one additional parameter is necessary to match the functional form of the deconvolution. That is to say, galaxies in between the extremes in stellar mass will lie in a transition region, as they may have been substantially affected by multiple feedback mechanisms. The width of this transition region will depend on many things, e.g., how long galaxies take to gain stellar mass, how much of the stellar mass present came from quiescent star formation as opposed to mergers, and the degree of interaction between multiple feedback mechanisms. Hence, instead of having M_h(M_*) become suddenly sub-exponential for galaxies larger than M_*,0, we allow for a slow "turn-on" of the more rapid growth. The behavior of M_h(M_*) is best fit by modifying the previous equation to

$$\begin{equation} \log _{10}(M_h(M_\ast)) \rightarrow \log _{10}(M_1) + \frac{\big(\frac{M_\ast }{M_{\ast,0}}\big)^\delta }{1+ \big(\frac{M_{\ast }}{M_\ast,0}\big)^{-\gamma }}. \end{equation} \tag{ 20 }$$

The denominator, 1 + (M_*/M_*,0)^−γ, is large for M_* < M_*,0, and it falls to unity for M_*>M_*,0 at a rate controlled by γ. A larger value of γ implies a more rapid transition between the power-law and sub-exponential behavior (typical values for (γ) at z = 0 are 1.3–1.7). As the non-constant piece of M_h(M_*) in Equation (20) is $\frac{1}{2}$ for M_* = M_*,0, we add a final factor of $-\frac{1}{2}$ to compensate so that M_h(M_*,0) = M₁.

To summarize, our resulting best-fit functional form has five parameters:

$$\begin{eqnarray} \log _{10}(M_h(M_\ast)) &=& \log _{10}(M_1) + \beta \,\log _{10}\left(\frac{M_\ast }{M_{\ast,0}}\right)\nonumber\\ && + \frac{\big(\frac{M_\ast }{M_{\ast,0}}\big)^\delta }{1 + \big(\frac{M_{\ast }}{M_{\ast,0}}\big)^{-\gamma }} - \frac{1}{2}, \end{eqnarray} \tag{ 21 }$$

where M₁ is a characteristic halo mass, M_*,0 is a characteristic stellar mass, β is the faint-end slope, and δ and γ control the massive-end slope. The best fit using this functional form is shown in Figure 2, and it achieves excellent agreement over the entire range of stellar masses.

Deconvolving the GSMF at higher redshifts does not suggest that anything more than linear evolution in the parameters is necessary, at least out to z = 1. While the characteristic mass of the GSMF and the characteristic mass of the halo mass function certainly evolve, the change in the shapes of the two functions is relatively slight. As we wish for the functional form to have a natural extension to higher redshifts, we parameterize the evolution in terms of the scale factor (a):

$$\begin{eqnarray} \log _{10}(M_1(a)) & = & M_{1,0} + M_{1,a} \, (a-1), \nonumber \\ \log _{10}(M_{\ast,0}(a)) & = & M_{\ast,0,0} + M_{\ast,0,a} \, (a-1), \nonumber \\ \beta (a) & = & \beta _0 + \beta _a \, (a-1),\\ \delta (a) & = & \delta _0 + \delta _a \, (a-1), \nonumber \\ \gamma (a) & = & \gamma _0 + \gamma _a \, (a-1), \nonumber \end{eqnarray} \tag{ 22 }$$

where a = 1 is the scale factor today.

By far the largest contributor to the error budget of the SM–HM relation is the systematic error parameter μ. As the effect of μ is to multiply all stellar masses by a constant factor, and as the width of the error bars in Figure 5 corresponds almost exactly to the prior on μ, we may conclude that reducing the error on the systematic shifts in stellar mass calculations would represent the single largest improvement in our understanding of the shape of the SM–HM relation. Figure 6 shows the substantially smaller error bars that result if systematic errors (μ and κ) in the stellar mass calculations are neglected.

The effect of ignoring scatter in stellar mass at fixed halo mass (i.e., setting ξ = 0) is shown at two redshifts in Figure 7. We find that the change is insignificant below halo masses of 10¹² M_☉, and is within statistical error bars below 10¹³ M_☉ for z = 1. This is a result of the fact that the slope of the stellar mass function below 10^10.5 M_☉ in stellar mass (corresponding to 10¹² M_☉ in halo mass) is not steep enough for scatter to have significant impact (also see Tasitsiomi et al. 2004). Because ξ>0 results in high stellar mass galaxies being assigned to lower-mass halos than they would be otherwise (due to the higher number density of lower-mass halos), the effect is that higher-mass halos contain fewer stars on average than they would for ξ = 0. The effect of setting ξ = 0 exceeds systematic error bars only for the very highest mass halos, above 10^14.5 M_☉.

Zoom In Zoom Out Reset image size

We note that our posterior distribution constrains ξ to be less than 0.22 dex at the 98% confidence level. Higher values for ξ would result in GSMFs inconsistent with the steep falloff of the Li & White (2009) GSMF (also see discussion in Guo et al. 2010).

The significance of including or excluding random statistical errors in stellar mass calculations, σ(z), is also shown in Figure 7. The effect of this type of scatter on the SM–HM relation is mathematically identical to the effect of scatter in stellar mass at fixed halo mass. As σ(z = 0) (∼0.07 dex) is much smaller than the expected value of ξ (∼0.16 dex), the convolution of the two effects is only marginally different from including ξ alone at z = 0; this results in only a minor effect on the SM–HM relation. The effect becomes more pronounced at z = 1 for the reason that σ(z = 1) (∼0.12 dex) becomes more comparable to ξ; and so, including the effects of statistical errors in stellar mass becomes as important as modeling scatter in stellar mass at fixed halo mass.

In Figure 8, we show a comparison of best fits for the stellar mass fraction using abundance matching with three different halo mass functions: analytic prescriptions for WMAP5 and WMAP1 (see Section 3.2.2), as well as the mass function taken directly from the L80G simulation (see Section 3.2.1). The difference between the L80G simulation and the analytic WMAP1 mass function is slight, as the L80G simulation uses WMAP1 initial conditions (h = 0.7, Ω_m = 0.3, Ω_Λ = 0.7, σ₈ = 0.9, n_s = 1); the difference is consistent with sample variance for the relatively small (80 h⁻¹ Mpc) size of the simulation. The difference between SM–HM relations using WMAP1 and WMAP5 cosmologies is within the systematic errors at all masses. When systematic errors are neglected, the two cosmologies yield SM–HM relations that are noticeably different only at low halo masses (M < 10¹² M_☉).

Zoom In Zoom Out Reset image size

Figure 9 shows the results of including uncertainties in the WMAP5 cosmological parameters. As described in Section 3.1, this is done using halo mass functions calculated with parameters resampled from the cosmological parameter chains provided by the WMAP team. Only at z ∼ 0 are the changes in error bars significant enough to justify mention. Here, the uncertainty in cosmology begins to exceed other sources of statistical error for halos below 10¹² M_☉ due to the small errors on the GSMF at the stellar masses associated with such halos (Li & White 2009). However, the cosmology uncertainties are still well within the systematic error bars.

Zoom In Zoom Out Reset image size

Because of the large volume of the SDSS, sample variance contributes insignificantly to the error budget for the SM–HM relation below z = 0.2. Above that redshift, the comparatively limited survey volume of Pérez-González et al. (2008) results in sample variance becoming an important contributor to the statistical error for halos below 10¹² M_☉ (Poisson noise dominates for larger halos). If the effects of sample variance were ignored, the statistical error spreads for our derived SM–HM relations at z = 1 would shrink from 0.12 dex to 0.09 dex for 10¹¹ M_☉ halos, and from 0.05 dex to 0.04 dex for 10^12.25 M_☉ halos. As with other types of errors, these considerations are well below the limits of the systematic error bars.

We caution that our error bars including sample variance at z>0 have a very specific meaning. Namely, they include the standard deviation in our fitting form which might be expected if the survey in Pérez-González et al. (2008) had been conducted on alternate patches of the sky. Sample variance at redshifts z>0 impacts only the linear evolution of the SM–HM relations we derive, as the large volume probed by the SDSS constrains the SM–HM relation very well at z ∼ 0. Because our fit is matched to the ensemble of reported data between 0 < z < 1, it is less vulnerable to the effects of sample variance in individual redshift bins. Instead, it is affected most by overall shifts in the number densities reported for the entire high-redshift survey. While this means that our fit gives a more robust SM–HM relation at all redshifts, some caution must be used when comparing our relation to results derived from the GSMF in a single redshift bin (e.g., 0.2 < z < 0.4). These will have much larger uncertainties due to sample variance than SM–HM relations derived (like ours) from GSMFs along a light cone probing a large redshift range.

To demonstrate the credibility of our approach for calculating the appropriate error bars including sample variance, we repeated our analysis of the SM–HM relation using GSMFs from Drory et al. (2009; which appeared as we were completing this work) instead of Pérez-González et al. (2008) for 0.2 < z < 1 and retaining the GSMF from Li & White (2009) for z < 0.2. Although the COSMOS survey in Drory et al. (2009) covers a much larger area (∼9× the area in Pérez-González et al. 2008), the fact that it is a single field means that the expected sample variance is only slightly smaller than for the combined fields in Pérez-González et al. (2008). As might be expected, the SM–HM relation at z = 0.1 using the Drory et al. (2009) GSMF is identical to the result in Figure 6 because of the strong constraining power of the SDSS data sample. At z = 1, the SM–HM relation generated by using the Drory et al. (2009) GSMF is within our quoted statistical and sample variance errors, as shown in Figure 12. This may not be surprising unless one considers that clustering results suggest an overdensity at the 2σ–3σ level in the COSMOS redshift bin z = 1 (Meneux et al. 2009). However, because our method fits the Drory et al. (2009) GSMFs across the entire redshift range, the excess at z = 1 is partially offset by an underdensity at z = 0.5. This demonstrates the robustness of our fitting method to the effects of sample variance except on the scale of the entire survey.

Finally, we consider the changes in both the stellar mass fraction and the total stellar mass fraction (total stellar mass in central galaxy and all satellites/total halo mass) induced by different satellite evolution models (see Section 3.3.2 for details on the two models). As shown in Figure 10, fixing satellite stellar mass at the redshift of accretion (lines labeled as SMF_acc) has virtually no effect on either fraction as compared to allowing satellite stellar mass to evolve the same way as centrals with the same mass (labeled as SMF_now). Because comparison of different satellite evolution models requires tracking satellites through merger trees, Figure 10 shows results only for satellites in the L80G simulation.

Zoom In Zoom Out Reset image size

The treatment of satellites may have a somewhat larger impact on the total stellar mass fraction, including the stellar mass of both central and satellite galaxies within a halo. This is shown for both models in Figure 10. Because of the steep falloff in stellar mass for low-mass galaxies, the total stellar mass fraction has only minimal contribution from satellites for low-mass halos and deviates significantly from the stellar mass fraction for centrals only at halo masses M_h>10^12.5 M_☉. At cluster-scale masses (M_h ∼ 10¹⁴ M_☉), accreted satellites have on average a higher ratio of stars to dark matter than the central galaxy, and the total stellar mass fraction can be many times the central stellar mass fraction. However, the impact of the two models for satellite treatment on this ratio is small. Profiles of satellite galaxies in clusters should be able to distinguish better between such models.

Systematic stellar mass offsets resulting from modeling choices result in the single largest source of uncertainties (∼0.25 dex at all redshifts). The contribution from all other sources of error is much smaller, ranging from 0.02 to 0.12 dex at z = 0 and from 0.07 to 0.16 dex at z = 1. On the other hand, this statement is only true when all contributing sources of scatter in stellar masses are considered. Models that do not account for scatter in stellar mass at fixed halo mass will overpredict stellar masses in 10^14.25 M_☉ halos by 0.13–0.19 dex, depending on the redshift. Models that do not account for scatter in calculated stellar mass at fixed true stellar mass will overpredict stellar masses in 10^14.25 M_☉ halos by 0.12 dex at z = 1. Hence, it is important to take both these effects into account when considering the SM–HM connection either at high masses or at high redshifts.