CONSTRAINTS ON THE RELATIONSHIP BETWEEN STELLAR MASS AND HALO MASS AT LOW AND HIGH REDSHIFT

In order to link the stellar mass of a galaxy m to the mass of its dark matter halo M, we need to specify the SHM ratio. A direct comparison of the halo mass function n(M) and the galaxy mass function ϕ(m) helps to constrain the SHM function. If we assume that every host (sub) halo contains exactly one central (satellite) galaxy and that each system has exactly the same SHM ratio m/M, the galaxy SMF can be derived trivially from the halo mass function and has the same features. The galaxy mass function derived for m/M = 0.05 is compared to the observed Sloan Digital Sky Survey (SDSS) galaxy mass function in Figure 1. The observed galaxy mass function is steeper for high masses and shallower for low masses than the one derived from the halo mass function. Thus, for a constant SHM ratio there will inevitably be too many galaxies at the low and high mass end.

Zoom In Zoom Out Reset image size

This implies that the actual SHM ratio m/M is not constant, but increases with increasing mass, reaches a maximum around m* and then decreases again. Hence, we adopt the following parameterization, similar to the one used in Yang et al. (2003):

$$\begin{equation} \frac{m(M)}{M} = 2 \left(\frac{m}{M}\right)_0 \left[\left(\frac{M}{M_1}\right)^{-\beta } + \left(\frac{M}{M_1}\right)^{\gamma }\right]^{-1}. \end{equation} \tag{ 2 }$$

It has four free parameters: the normalization of the SHM ratio (m/M)₀, a characteristic mass M₁, where the SHM ratio is equal to (m/M)₀, and two slopes β and γ which indicate the behavior of m/M at the low and high mass ends, respectively. We use the same parameters for the central and satellite populations, since—unlike luminosity—the stellar mass of satellites changes only slightly after they are accreted by the host halo.

Note that though both β and γ are expected to be positive, they are not restricted to be so. The SHM relation is therefore not necessarily monotonic.

Having set up the model, we now need to constrain the four free parameters M₁, (m/M)₀, β, and γ. To do this, we populate the halos in the simulation with galaxies. The stellar masses of the galaxies depend on the mass of the halo and are derived according to our prescription (Equation (2)). The positions of the galaxies are given by the halo positions in the N-body simulation.

Once the simulation box is filled with galaxies, it is straightforward to compute the SMF Φ_mod(m). As we want to fit this model mass function to the observed mass function Φ_obs(m) by Panter et al. (2007), we choose the same stellar mass range (10^8.5 M_☉–10^11.85 M_☉) and the same bin size. The observed SMF was derived using spectra from the Sloan Digital Sky Survey Data Release 3 (SDSS DR3); see Panter et al. (2004) for a description of the method.

Furthermore, it is possible to determine the stellar mass dependent clustering of galaxies. For this, we compute projected galaxy CFs w_p,mod(r_p, m_i) in several stellar mass bins which we choose to be the same as in the observed projected galaxy CFs of Li et al. (2006). These were derived using a sample of galaxies from the SDSS DR2 with stellar masses estimated from spectra by Kauffmann et al. (2003).

We first calculate the real space CF ξ(r). In a simulation, this can be done by simply counting pairs in distance bins:

$$\begin{equation} \xi (r_i) = \frac{dd(r_i)}{N_p(r_i)}-1, \end{equation} \tag{ 3 }$$

where dd(r_i) is the number of pairs counted in a distance bin and N_p(r_i) = 2πN²r²_iΔr_i/L³_box, where N is the total number of galaxies in the box. The projected CF w_p(r_p) can be derived by integrating the real space correlation function ξ(r) along the line of sight:

$$\begin{equation} w_p(r_p) = 2 \int _{0}^{\infty } dr_{||} \xi \left(\sqrt{r_{||}^2+r_p^2}\right) = 2 \int _{r_p}^{\infty } dr \frac{r \xi (r)}{\sqrt{r^2-r_p^2}} \;, \end{equation} \tag{ 4 }$$

where the comoving distance (r) has been decomposed into components parallel (r_||) and perpendicular (r_p) to the line of sight. The integration is truncated at 45 Mpc. Due to the finite size of the simulation box (L_box = 100 Mpc), the model correlation function is not reliable beyond scales of r ∼ 0.1L_box ∼ 10 Mpc.

In order to fit the model to the observations, we use Powell's directions set method in multidimensions (e.g., Press et al. 1992) to find the values of M₁, (m/M)₀, β, and γ that minimize either

$$\begin{eqnarray*} &&\chi _r^2 = \chi _r^2 (\Phi)= \frac{\chi ^2(\Phi)}{N_{\Phi }} \end{eqnarray*}$$

(mass function fit) or

$$\begin{eqnarray*} &&\chi _r^2 = \chi _r^2 (\Phi)+\chi _r^2 (w_p) = \frac{\chi ^2(\Phi)}{N_{\Phi }} + \frac{\chi ^2(w_p)}{N_r\;N_m} \end{eqnarray*}$$

(mass function and projected CF fit) with N_Φ and N_r the number of data points for the SMF and projected CFs, respectively, and N_m the number of mass bins for the projected CFs.

In this context, χ²(Φ) and χ²(w_p) are defined as follows:

$$\begin{eqnarray*} \chi ^2(\Phi) &=& \sum _{i=1}^{N_{\Phi }} \left[ \frac{\Phi _{\rm mod}(m_i)-\Phi _{\rm obs}(m_i)}{\sigma _{\Phi _{\rm obs}(m_i)}} \right]^2,\\ \chi ^2(w_p) &=& \sum _{i=1}^{N_m} \sum _{j=1}^{N_r} \left[ \frac{w_{p,\rm mod}(r_{p,j},m_i)-w_{p,\rm obs}(r_{p,j},m_i)}{\sigma _{w_{p,\rm obs}(r_{p,j},m_i)}} \right]^2\;, \end{eqnarray*}$$

with $\sigma _{\Phi _{\rm obs}}$ and $\sigma _{w_{p,{\rm obs}}}$ the errors for the SMF and projected CFs, respectively. Note that for the simultaneous fit, by adding the reduced χ²_r, we give the same weight to both data sets.

In order to obtain estimates of the errors on the parameters, we need their probability distribution prob(A|I), where A is the parameter under consideration and I is the given background information. The most likely value of A is then given by: A_best = max(prob(A|I)).

As we have to assume that all our parameters are coupled, we can only compute the probability for a given set of parameters. This probability is given by:

$$\begin{eqnarray*} &&{\rm prob(}M_1,(m/M)_0,\beta,\gamma \vert I) \propto \exp (-\chi ^2). \end{eqnarray*}$$

In a system with four free parameters A, B, C, and D one can calculate the probability distribution of one parameter (e.g., A) if the probability distribution for the set of parameters is known, using marginalization:

$$\begin{eqnarray*} {\rm prob(}A\vert I) &=& \int _{-\infty }^{\infty }{\rm prob(}A,B\vert I) dB\\ &=& \int _{-\infty }^{\infty }{\rm prob(}A,B,C,D\vert I) d\!B d\hspace{-1pt}C d\!D. \end{eqnarray*}$$

Once the probability distribution for a parameter is determined, one can assign errors based on the confidence intervals. This is the shortest interval that encloses a certain percentage X of the area under the posterior probability distribution. For the 1σ error X = 68%, while for the 2σ error X = 95%. Assuming that the probability distribution has been normalized to have unit area we seek A₁ and A₂ such that

$$\begin{eqnarray*} &&\int _{0}^{A_1} {\rm prob(}A\vert I) d\!A = \int _{A_2}^{\infty } {\rm prob(}A\vert I) dA = \frac{1-X}{2}. \end{eqnarray*}$$

Finally, the parameter A is given as $A = A_{\rm best}\mbox{ }_{-\sigma _{-}}^{+\sigma _{+}}$ with σ₊ = A₂ − A_best and σ₋ = A_best − A₁. The errors derived in this way only include sources that have been considered when computing χ². The calculation of the errors applies for uncorrelated data points. Since in our case the data points are correlated, the values of the errors are slightly modified. Also errors caused by cosmic variance are not included.

First, we fit to the SDSS SMF and use the derived best-fit parameters to calculate the model projected correlation functions. Note that for now, we do not take into account any possible scatter in the m(M) relation. We will consider scatter in Section 4.5.

We see in Figure 1 that our fit produces excellent agreement with the observed SMF. Using the approach described above we also compute the errors on the parameters. The results are summarized in Table 1.

Having derived the best-fit parameters, we can predict the projected CFs. We present the results both including and not including orphan galaxies, where we have fitted to the SMF for each case.

Figure 2 shows a comparison between our model and the SDSS projected correlation functions in five stellar mass bins ranging from log  m/ M_☉ = 9.0 to log  m/ M_☉ = 11.5 with a bin size of 0.5 dex. The correlation function that has been derived without orphans is too low at small scales and can be regarded as a lower limit. Neglecting these galaxies results in an underprediction of satellite galaxy clustering. As on small scales, the projected CF depends mainly on the one-halo term this results in the underprediction of w_p(r_p). This effect weakens for the clustering of more massive galaxies as they are more likely to be central galaxies and thus not effected by tidal stripping at all.

Zoom In Zoom Out Reset image size

The agreement with the observationally derived w_p(r_p) for the catalog including orphaned galaxies is very good, which is also reflected in the low value of χ²_r(w_p) = 3.83. Note that this value has been calculated with the parameters from the mass function fit given above and does not correspond to a fit to the projected CFs.

Note that we plot the projected CFs only up to 20 Mpc. Because of the finite box size, the clustering of host halos and thus central galaxies is underpredicted at large scales independent of mass. Additionally, due to the lack of long-wavelength modes, massive halos, and galaxies can be underproduced leading to an underprediction of w_p for the massive objects, independent of scale. However, the latter effect is very small, since the abundance of the massive halos in our simulation agrees very well with the predicted average (Sheth & Tormen 1999).

As a test we also used the present mass instead of the maximum mass for subhalos. We then found that the projected CF was particularly underpredicted on small scales. This effect is due to tidal stripping of subhalos and is thus strongest at small scales where the subhalo contribution dominates.

We now investigate whether we can improve the agreement between the model and the observed projected CFs by performing a combined fit as described above. We obtain the same parameters as those we derived from the fit to the SMF alone. This seems surprising, but on further inspection we find that this is due to χ²(m) being a lot more sensitive to changes of the parameters than χ²(w_p). This means that if one changes the parameters a little in order to improve the fit to the projected correlation functions, one can get a slightly better agreement between the model and the observed projected CFs only at the cost of a large disagreement between the model and the observed SMFs. In other words: χ²(w_p) is much flatter around its minimum than χ²(w_p), as shown in Figure 3.

Zoom In Zoom Out Reset image size

This means that, assuming that both central and satellite galaxies follow the same SHM relation, the model that matches the SMF can reproduce the correct clustering. However, if subhalos have a different SHM ratio there is an infinite number of solutions that match the SMF but produce very different correlation functions. The only way to constrain the SHM relations then is to take the clustering data into account. By adopting different SHM relations for central and satellite populations it is even possible to produce a slightly better fit to the correlation functions (Wang et al. 2006).

On the other hand, if one wants to predict clustering as a function of stellar mass (e.g., at higher redshift) then one has to make an assumption about how the SHM ratios of central and satellite galaxies are related. We made the very simple assumption, that the relation between the stellar mass of central galaxies and the virial mass of their host halo and the relation between the stellar mass of satellite galaxies and the mass of the subhalo at the time of accretion is the same, and have shown that this leads to very good predictions for the mass dependent clustering. We conclude that under this simple assumption we can use our model to predict clustering as a function of stellar mass.

The upper panel of Figure 4 shows the derived stellar mass as a function of halo mass. The light shaded area gives the 68% confidence interval, while the dark and light shaded areas together give the 95% confidence interval. These have been derived using a set of different models computed on a mesh, as described in Section 3.3.

Zoom In Zoom Out Reset image size

For the SHM ratio, we apply the same procedure. The result is shown in the lower panel of Figure 4. We see that the SHM ratio has the form we expected: it increases with increasing halo mass, reaches its maximum value around M₁ and then decreases again.

We now explore the effects of changing each parameter in order to understand how they affect the SMF. If we keep M₁, β, and γ fixed and only vary (m/M)₀, this corresponds to changing the stellar mass of the galaxy that lives inside each halo by a constant factor. This has no impact on the form of the SMF. Its shape stays the same, while only the position on the stellar mass axis changes. Due to the monotonic form of the SMF, this directly determines the value of the normalization ϕ*. For a larger value of (m/M)₀, we get a larger value of ϕ*.

Varying only M₁, we find that the shape of the SMF changes drastically. For a higher M₁ than our best-fit value, we get too many massive galaxies and too few low mass galaxies, while for a lower value of M₁ we get too few massive galaxies and too many low mass galaxies. This is because M₁ is the characteristic mass corresponding to the highest SHM ratio. In the SMF, this corresponds to the knee and we get a SMF which has its knee at the stellar mass corresponding to M₁. For a larger M₁, the knee is shifted to a higher stellar mass. Together, M₁ and the maximum SHM ratio (m/M)₀ determine the normalization of the SMF ϕ and the characteristic mass m*.

Changing β affects mainly the low-mass slope of the SMF. For larger values of β, the slope becomes shallower. As β influences mainly the slope of the low mass end of the SMF, it is strongly related to the parameter α of the Schechter function. A small value of β corresponds to a high value of α.

If we change γ, this mainly impacts the slope of the massive end of the SMF. For larger values of γ than for its best-fit value, the slope of the massive end becomes steeper. As γ affects mainly the slope of the massive end of the SMF, it is not coupled to a parameter of the Schechter function though it is related to the high-mass cutoff, assumed to be exponential in a Schechter function.

Figure 5 shows the contours of the two-dimensional probability distributions for the parameters pairs. We see a correlation between the parameters [M₁, γ] and [(m/M)₀, γ] and an anti-correlation between [β, γ], [β, M₁], and [(m/M)₀, M₁]. There does not seem to be a correlation between [β, (m/M)₀].

Zoom In Zoom Out Reset image size

Up until now we have assumed that there is a one-to-one, deterministic relationship between halo mass and stellar mass. However, in nature, we expect that two halos of the same mass M may harbor galaxies with different stellar masses, since they can have different halo concentrations, spin parameters, and merger histories.

For each halo of mass M, we now assign a stellar mass m drawn from a lognormal distribution with a mean value given by our previous expression for m(M) (Equation (2)), with a variance of σ²_m. We assume that the variance is a constant for all halo masses, which means that the percent deviation from m is the same for every galaxy. This is consistent with other halo occupation models, SAMs and satellite kinematics (Cooray 2006; van den Bosch et al. 2007; More et al. 2009b).

Assuming a value of σ_m = 0.15 dex and fitting the SMF only, we find the values given in Table 2. These values lie within the (2σ) error bars of the best-fit values that we obtained with no scatter. The largest change is on the value of γ, which controls the slope of the SHM relation at large halo masses. The SMF and the projected CFs for the model including scatter are shown in Figures 1 and 2, respectively, and show very good agreement with the observed data.

In Figure 6, we compare our model without scatter with the model including scatter. We have also included the relation between halo mass and the average stellar mass. Especially at the massive end scatter can influence the slope of the SMF, since there are few massive galaxies. This has an impact on γ and as all parameters are correlated scatter also affects the other parameters. We thus see a difference between the model without scatter and the most likely stellar mass in the model with scatter in Figure 6.

Zoom In Zoom Out Reset image size

In order to specify the CMF Φ(m|M), we divide the galaxy population into a central and a satellite part, as in the updated CLF formalism (Zheng et al. 2005; Zehavi et al. 2005; Cooray 2006; Yang et al. 2008; Cacciato et al. 2009). The central part is Φ_c(m|M) and the satellite part is Φ_s(m|M). Then the total CMF is the sum of both parts:

$$\begin{equation} \Phi (m \vert M) = \Phi _c(m \vert M) + \Phi _s(m \vert M). \end{equation} \tag{ 7 }$$

Note that both Φ_c(m|M) and Φ_s(m|M) are statistical functions and should not be regarded as the mass functions of galaxies living in a given individual halo.

For the central population, we expect the CMF to have a peak around the stellar mass m_c that corresponds to the host halo's virial mass M in the SHM relation (Equation 2). Due to the halo mass bin size this distribution gets smeared out, because halos in the interval [M₁, M₂] contain central galaxies of stellar masses m ∈ [m₁(M₁), m₂(M₂)]. Thus, $\tilde{\Phi }(m)/\Delta n$ will be finite inside the interval [m₁(M₁), m₂(M₂)] and zero elsewhere with a normalization such that the number of central galaxies per halo equals one. This can be regarded as scatter σ_bin due to the binning. If we add intrinsic scatter σ_m to Equation (2), we expect Φ_c(m|M) to be a lognormal with a maximum around m_c(M) and a variance of σ²_m. To this scatter, the binning scatter σ_bin adds in quadrature (assuming that σ_bin and σ_m are uncorrelated), resulting in a total scatter of σ²_c = σ²_m + σ²_bin. For both cases (σ_m = 0 and σ_m ≠ 0), we use a lognormal distribution:

$$\begin{equation} \Phi _c(m \vert M) = \frac{1}{\sqrt{2\pi }\ln 10\;m\;\sigma _c} \exp \left[-\frac{\log \, ^2(m/m_c)}{2\sigma _c^2}\right]\;, \end{equation} \tag{ 8 }$$

where the mean m_c(M) and width σ²_c(M) are parameterized functions of the halo mass M.

For the satellite population, we adopt a Schechter function with a steeper slope for the massive end. This is done by squaring the argument of the exponential function in the Schechter function:

$$\begin{equation} \Phi _s(m \vert M) = \frac{\Phi _s^*}{m_s} \left(\frac{m}{m_s}\right)^{\alpha _s} \exp \left[ -\left(\frac{m}{m_s}\right)^2 \right]\;. \end{equation} \tag{ 9 }$$

Also here the parameters Φ*_s(M), m_s(M), and α_s(M) are functions of the host halo mass M. They are the normalization, the characteristic mass, and the low-mass slope of the satellite population of host halos of mass M.

We populate the halos and subhalos in our simulation with central and satellite galaxies according to the prescription in Section 3. Then we choose halo mass bins between log  M/ M_☉ = 10.2 and log  M/ M_☉ = 15.0 with a bin size of ΔM = 0.4 dex. In every halo mass bin, we seek all galaxies which live in a host halo with a mass in that bin, which we divide between central and satellite galaxies. For these populations, we then compute two separate SMFs which we normalize such that the number of central galaxies per host halo equals one. This procedure then yields for every halo mass bin a central and a satellite distribution $(d\tilde{n}_g / d\log \, M) \Delta n_h$.

Using Equation (6), we can now relate the SMF in a halo mass bin to the CMF:

$$\begin{eqnarray} \frac{d\tilde{n}_g(m)}{d\log \, M}\frac{1}{\Delta n_h} &=& \frac{\ln {10}}{\Delta n_h} \; M \; \frac{d\tilde{n}_g(m)}{dM} \nonumber \\ &=& \ln {10} \; M \; \frac{\tilde{\Phi }(m)}{\Delta n_h} \nonumber \\ &\approx & \ln {10} \; M \; \Phi (m\vert M). \end{eqnarray} \tag{ 10 }$$

Now we can fit the five parameters m_c(M), σ_c(M), m_s(M), Φ*_s(M), and α_s(M) to the SMFs in each halo bin. We compute and fit the central and the satellite parts separately.

The left panels of Figure 7 show the CMF in a subsample of halo mass bins running from log  M/ M_☉ = 10.2 ± 0.2 to 15.0 ± 0.2, where we have not included intrinsic scatter in the SHM relation. For the satellite part, only galaxies with a mass above the completeness limits for each halo mass bin (as indicated in Figure 7) have been used in the fit.

Zoom In Zoom Out Reset image size

In low-mass halos (log  M/ M_☉ < 11.0), the contribution from satellite galaxies is very small and the central contribution dominates until log  M/ M_☉ = 12.0. For massive halos (log  M/ M_☉ > 13.0), the satellite contribution dominates by number. The mean of the lognormal fit to the central contribution also increases with halo mass as stipulated by the model derived in Section 3. The characteristic mass scale of the satellite contribution also increases with halo mass meaning that the most massive satellite galaxies have a mass which is comparable to the mass of the central galaxy.

The scatter of the central contribution σ_c(M) decreases with halo mass. As we did not include any scatter in the model, this scatter reflects the width (0.4 dex) of the halo mass bins (σ_bin). The halo mass dependence of σ_c(M) arises because a fixed halo mass bin is mapped to a smaller galaxy mass bin for larger halo mass due to the shape of the SHM relation. Another feature of the CMF is the slope for low mass satellite galaxies α_s(M) which becomes shallower with increasing halo mass.

In this section, we investigate the halo mass dependence of the five parameters of the CMF: m_c(M), σ_c(M), m_s(M), Φ*_s(M), and α_s(M). They have been fixed by fitting to the SMFs in each halo mass bin. We introduce a parameterization in order to describe the dependence on halo mass and constrain these by a fit to each parameter. The results are presented in Table 3. This provides a complete description of the CMF.

As we have already determined the mean relation between the stellar mass of a galaxy and the mass of its halo, the form of m_c(M) has to be the same and can thus be described by Equation (2):

$$\begin{equation} m_c(M) = 2 \; M \; \left(\frac{m_c}{M}\right)_0 \left[\left(\frac{M}{M_{1c}}\right)^{-\beta _c} + \left(\frac{M}{M_{1c}}\right)^{\gamma _c}\right]^{-1}. \end{equation} \tag{ 11 }$$

This yields four parameters (m_c/M)₀, M_1c, β_c, and γ_c.

In the upper left panel of Figure 8, m_c(M) is plotted as a function of halo mass. Note that by construction, it has the same form as the SHM relation.

Zoom In Zoom Out Reset image size

The scatter of the central galaxy contribution is high for low halo masses and decreases for more massive halos. The middle left panel of Figure 8 shows σ_c(M) as a function of halo mass. As one can see, σ_c(M) goes to a constant value both for low and high halo masses while it decreases with halo mass. We therefore choose the following parameterization:

$$\begin{equation} \sigma _c(M) = \sigma _{\infty } + \sigma _1 \left[ 1-\frac{2}{\pi } \arctan \left(\xi \log \, \frac{M}{M_2} \right) \right]. \end{equation} \tag{ 12 }$$

This yields four more parameters σ_∞, σ₁, ξ, and M₂. Here, σ_∞ sets the high mass limit of σ_c(M) while σ₁ sets the difference between the low and high mass limits of σ_c(M). The parameter M₂ determines the mass scale at which the transition occurs and ξ sets the strength. For a large (small) value of ξ, the transition occurs in a small (large) interval around M₂.

The specific shape of σ_c(M) can be explained by the form of the SHM relation (Equation (2)). As we have not included any scatter in this relation (σ_m = 0), the width of the lognormal function of the central galaxy distribution arises from the width of the halo mass bin (σ_c = σ_bin). A halo mass interval [M₁, M₂] contains only central galaxies with stellar masses of m ∈ [m₁(M₁), m₂(M₂)]. The lower left panel of Figure 8 illustrates this by showing how halo mass bins affect the bin size of the stellar mass. If we choose the same bin size for low and high mass halos, we get different bin sizes for low and high mass galaxies, due to the changing slope of m(M). Therefore, the transition occurs where the slope of m(M) changes which is around M₁, so the value of M₂ is very close to that value.

As Figure 7 shows that the satellite contribution falls off around the mean mass of the central galaxy, we expect the characteristic mass of the modified Schechter function m_s(M) to follow m_c(M). We therefore describe m_s(M) with the same function we used for the parameterization of m_c(M):

$$\begin{equation} m_s(M) = 2 \; M \; \left(\frac{m_s}{M}\right)_0 \left[\left(\frac{M}{M_{1s}}\right)^{-\beta _s} + \left(\frac{M}{M_{1s}}\right)^{\gamma _s}\right]^{-1}. \end{equation} \tag{ 13 }$$

This function yields four parameters (m_s/M)₀, M_1s, β_s, and γ_s.

The upper right panel of Figure 8 plots m_s(M) as a function of halo mass. We see that the shape is similar to that of m_c(M). Note that m_s(M) is always lower than m_c(M), while the deviation increases with increasing halo mass. This implies that for high halo masses the satellite contribution to the CMF falls off before the mean mass of the central galaxy.

The normalization of the modified Schechter function is small for low halo masses and increases with the mass of the host halo. The middle right panel of Figure 8 shows Φ*_s(M) as a function of halo mass. We see that Φ*_s(M) can be described by a power law and choose the following parameterization:

$$\begin{equation} \Phi _s^*(M) = \Phi _0 \left(\frac{M}{\,{M}_\odot }\right)^{\lambda }. \end{equation} \tag{ 14 }$$

We get two more parameters, Φ₀ and λ. The normalization of Φ*_s(M) is given by Φ₀ and the slope by λ. The shape of Φ*_s(M) implies that the probability for a host halo to harbor satellite galaxies (in a given stellar mass range) increases with increasing halo mass.

The slope of the modified Schechter function for the satellite contribution becomes shallower for more massive halos. The lower right panel of Figure 8 shows α_s(M) as a function of halo mass and shows that α_s(M) goes to a constant value for both low and high halo masses. Similar to σ_c(M), we choose the parameterization:

$$\begin{equation} \alpha _s(M) = \alpha _{\infty } + \alpha _1 \left[ 1-\frac{2}{\pi } \arctan \left(\zeta \log \, \frac{M}{M_3} \right) \right]. \end{equation} \tag{ 15 }$$

This yields four more parameters α_∞, α₁, ζ, and M₃. Here, α_∞ sets the high mass limit of α_c(M) while α₁ sets the difference between the low and high mass limits of α_c(M). The mass scale at which this transition occurs is determined by M₃ and ζ sets its strength. The transition occurs in a small (large) interval around M₃ for a large (small) value of ζ.

Until now, we have used the SHM Equation (2) without any intrinsic scatter. In this section, we investigate how the CMF and the parameters change if we include a scatter σ_m as described in Section 4.5. This scatter is again assumed to be constant with host halo mass.

The right panels of Figure 7 show the resulting CMF in a subsample of halo mass bins for an intrinsic scatter of σ_m = 0.15. The central part is now no longer near-constant in the interval [m(M − ΔM/2), m(M + ΔM/2)] as in the left panels of Figure 7 (where σ_m = 0.0) but has the form of a lognormal with a broader distribution for bigger σ_m. As the scatter has been taken from a lognormal distribution, the central galaxy contribution to the CMF is distributed in the same way. Hence, σ_c(M) changes with respect to the model that does not include artificial scatter. We notice that at the massive end the binning scatter σ²_bin and the intrinsic scatter σ²_m add to the total scatter σ²_tot. At the low mass end, however, the total scatter is less than what has been obtained by using no intrinsic scatter. This shows that the two forms of scatter do not add in quadrature and indicates that they are correlated.

We compare m_c(M), σ_c(M), m_s(M), Φ*_s(M), and α_s(M) for σ_m = 0 and σ_m = 0.15 and show the resulting parameters in Table 3 (Columns 4 and 5) and in Figure 8. The mean mass of the central galaxy m_c(M) does not change much if artificial scatter is introduced. The most likely stellar mass of a central galaxy is still given by the SHM relation, so the mean of the Gaussian in logarithmic space stays the same. Also the parameters of the satellite population (m_s(M),Φ*_s(M), and α_s(M)) do not change significantly.

In order to compare our results to other HOD models, it is useful to compute the average number of galaxies per halo 〈N〉, as this is the main prediction of the HOD approach. To compute 〈N〉(M) from the CMF, we simply integrate Φ(m|M) over the desired stellar mass range:

$$\begin{equation} \langle N\rangle (M) = \int _{m_1}^{m_2} \Phi (m\vert M) dm. \end{equation} \tag{ 16 }$$

As we have divided Φ(m|M) into a central galaxy contribution Φ_c(m|M) and a satellite galaxy contribution Φ_s(m|M), we can compute separate occupation numbers for central and satellite galaxies:

$$\begin{eqnarray*} \langle N\rangle (M) &=& \int _{m_1}^{m_2} \Phi _c(m\vert M) dm + \int _{m_1}^{m_2} \Phi _s(m\vert M) dm\\ &=& \langle N_c\rangle (M) + \langle N_s\rangle (M). \end{eqnarray*}$$

The average number of central galaxies per halo 〈N_c〉(M) is given by

$$\begin{equation} \langle N_c\rangle (M) = \case{1}{2}\left[{\rm erf}(\eta _2)-{\rm erf}(\eta _1)\right] \;, \end{equation} \tag{ 17 }$$

with the error function erf(x) and the integration boundaries

$$\begin{equation*} \eta _1 = \frac{\log \, (m_1/m_c)}{\sqrt{2}\sigma _c} \quad {\rm and} \quad \eta _2 = \frac{\log \, (m_2/m_c)}{\sqrt{2}\sigma _c} \;. \end{equation*}$$

The average number of satellite galaxies per halo 〈N_s〉(M) is

$$\begin{equation} \langle N_s\rangle (M) = \frac{\Phi _s}{2}\left[\Gamma \left(\frac{\alpha _s}{2}+\frac{1}{2},\kappa _1\right)-\Gamma \left(\frac{\alpha _s}{2}+\frac{1}{2},\kappa _2\right)\right] \;, \end{equation} \tag{ 18 }$$

with the upper incomplete gamma function Γ(a, x) and the integration boundaries

$$\begin{equation*} \kappa _1 = (m_1/m_s)^2 \quad {\rm and} \quad \kappa _2 = (m_2/m_s)^2 \;. \end{equation*}$$

Figure 9 shows the resulting occupation numbers for the values of the CMF parameters that were derived in Section 5.3 (using a scatter of σ_m = 0.15). The five lines in each panel correspond to different stellar mass bins.

Zoom In Zoom Out Reset image size

The left panel shows the average number of central galaxies per halo 〈N_c〉(M) as a function of halo mass. In the middle panel, the average number of satellite galaxies per halo 〈N_s〉(M) as a function of halo mass is shown. The right panel plots the average number of all galaxies per halo 〈N_tot〉(M) as a function of halo mass. A galaxy of a low stellar mass can thus either be a central galaxy of a low-mass halo, or a satellite galaxy of a massive halo. It is not likely to live in a halo of intermediate mass.

As it is common in the literature to plot occupation numbers not for stellar mass intervals, but for galaxy samples with a mass above a given threshold, we need to adjust Equations (17) and (18). The stellar mass threshold is then given by m₁ while m₂ → ∞. This yields for the average number of central galaxies

$$\begin{equation} \langle N_c\rangle (M, m_1) = \frac{1}{2}\left[1 - {\rm erf}\left(\frac{\log \, (m_1/m_c)}{\sqrt{2}\sigma _c}\right)\right]\;, \end{equation} \tag{ 19 }$$

since erf(x → ∞) → 1, and for the average number of satellite galaxies

$$\begin{equation} \langle N_s\rangle (M, m_1) = \frac{\Phi _s}{2} \; \Gamma \left[\frac{\alpha _s}{2}+\frac{1}{2},\left(\frac{m_1}{m_s}\right)^2\right],\; \end{equation} \tag{ 20 }$$

since Γ(a, x → ∞) → 0.

Figure 10 shows occupation numbers for different stellar mass thresholds. The left panel shows the average number of central galaxies per halo 〈N_c〉(M) as a function of halo mass. The middle panel plots the average number of satellite galaxies per halo 〈N_s〉(M) as a function of halo mass. It is similar to the middle panel of Figure 9, while it is larger at a given halo mass. In the right panel, the average number of all galaxies per halo 〈N_tot〉(M) as a function of halo mass is shown.

Zoom In Zoom Out Reset image size

Numerous variations on halo occupation models have been presented in the literature. In this section, we describe some of the most popular ones and compare them to our model. As many authors use different IMFs and definitions of halo masses, we convert all results to the conventions that we have used in this work (Kroupa IMF and virial overdensity).

In the non-parametric model (Vale & Ostriker 2006; Conroy et al. 2006; Shankar et al. 2006), galaxy properties, such as luminosity and stellar mass, are monotonically related to the mass of dark matter halos. Using the observed galaxy SMF, the most massive halo is matched to the most massive galaxy:

$$\begin{equation} n_g(>m_i) = n_h(>M_i). \end{equation} \tag{ 21 }$$

In this way, the observed SMF is automatically reproduced. Applying this procedure and fitting the parameters of the SHM relation to the result, we have derived the values given in Table 4. These are in good agreement with the parameters of our model, except for β. We find that this is due to the shape of the SHM ratio for low masses. For the non-parametric model, m(M < M₁) cannot be perfectly described by a single power law, as is assumed in our model.

Adding an additional parameter and assuming a fitting function with five free parameters, we are able to fit the SHM relation predicted by the non-parametric model quite precisely. The fifth parameter accounts for the deviation from the power law at high and low masses. Using the parameterization

$$\begin{equation} m(M)=m_0 \; \frac{(M/M_1)^{\gamma _1}}{\left[1+(M/M_1)^{\beta }\right]^{(\gamma _1-\gamma _2)/\beta }}, \end{equation} \tag{ 22 }$$

we determine the values given in Table 5. Figure 11 shows the results of four- and five-parameter fits to the SHM relation derived via the non-parametric method, compared with our usual model. In the range where we applied the mass function fit, the non-parametric model lies within our error bars.

Zoom In Zoom Out Reset image size

In Wang et al. (2006), a model similar to ours is used to constrain the SHM ratio. The halo catalog is taken from the Millennium simulation (Springel et al. 2005); halos are identified using a FoF group finder while substructure is found using the SUBFIND algorithm of Springel et al. (2001). As observational constraints, the authors use a SMF which they compute from the SDSS DR2 data using the mass estimates of Kauffmann et al. (2003) and the projected CFs of Li et al. (2006).

The parameterization they use is similar to ours, with four free parameters that can easily be converted to M₁, (m/M)₀, β, γ, and an unconstrained scatter. These are fixed by generating a grid of models and the best-fit model is defined as the one for which χ² = χ²(Φ) + χ²(w_p) is minimal. They find that their fit improves if they take a different set of parameters for central and satellite galaxies. In Table 4, we compare our best-fit parameters with their central galaxy best-fit parameters which have been updated in Wang et al. (2007). We show these results in Figure 11.

The values of M₁ and (m/M)₀ are in very good agreement with our values, but the slopes are both higher, resulting in fewer massive and fewer low mass galaxies. The reason for the difference in the low mass end is the different simulation used. As the resolution of the simulation in our model is higher, the low mass end can be constrained more tightly. For the massive end, the difference in γ can be explained by the additional unconstrained scatter that is used in Wang et al. (2006). As the mass function is steep at high masses and shallow for low masses, a change in the scatter will influence the number of massive galaxies strongly, while it will have only a small effect on the low mass end. As the other three parameters, M₁, (m/M)₀, and β, are coupled to the Schechter function parameters, there are two parameters to constrain the slope of the massive end of the SMF. This degeneracy can cause the difference in γ between the two models. The fact that in the Millennium simulation the cosmology is different to that of our simulation also affects the value of the parameters.

The relation between stellar mass and halo mass can be observationally constrained using galaxy–galaxy lensing. Gravitational lensing induces shear distortions of background objects around foreground galaxies, allowing the mass of the dark matter halo to be estimated. Mandelbaum et al. (2005, 2006) have used SDSS data to calibrate the predicted signal from a halo model which has been derived from a dissipationless simulation. They have extracted the mean halo mass as a function of stellar mass. The lensing data for combined early- and late-type galaxies (R. Mandelbaum 2008, private communication) are shown in Figure 11 and are in excellent agreement with our model.

As we discussed in the introduction, SAMs of galaxy formation attempt to predict the relationship between dark halo mass and stellar mass by a priori modeling of physical processes, such as the growth of structure, cooling, star formation, and stellar and active galactic nucleus (AGN) feedback. We compare our results with predictions from the latest version of the SAMs of Somerville & Primack (1999); see Somerville et al. (2008). For this, we compute the mean stellar mass of central galaxies as a function of the mass of the host halo in halo mass bins. The results are shown in Figure 11 and are in good agreement with our model. This is not surprising, as the physical parameters in the model of Somerville et al. (2008) have been tuned to match the observed SMF at z = 0.

In Wang et al. (2006), the authors use the SAM of Croton et al. (2006) and link galaxy properties, such as the stellar mass, to the mass of the halo in which the galaxy was last a central object M_infall. They fit the same four-parameter function that they used for their empirical model (described above) to obtain the parameter estimates from the SAM. We summarize these results in Table 4, and show them in Figure 11.

The two slopes are in very good agreement with our results. However, the normalization in the Croton et al. (2006) SAM is ∼25% higher and the characteristic mass is ∼25% lower than what we found and what Wang et al. (2006) found for their model. This is because the SAM of Croton et al. (2006) does not produce a perfect fit to the observed SMF.

Another direct way of studying galaxy properties as a function of halo mass is using the SDSS group catalog presented in Yang et al. (2007). In this approach, galaxies are first linked together into "groups" using a FoF algorithm. Each group is then assigned a total halo mass by matching to the theoretical dark matter halo mass function. Yang et al. (2008) present the relation between the mean stellar mass of the central galaxy and the host halo mass. We fit the parameters of Equation (2) to their relation and present the results in Table 4.

We note that the characteristic mass and the normalization derived from the group catalog are both higher than our model parameters. The high-mass slope of the SHM relation in the group catalog is shallower than that of our model. The low-mass slope is also shallower, however, the constraints on the low-mass slope in the group catalog are weak, since the lowest halo masses are log(M/ M_☉) ∼ 11.7. This can also be seen in Figure 11 where we show the SHM relation of the group catalog for comparison.

In order to constrain the SHM relation, we have to first select observational SMFs at the redshifts we want to investigate. Because of the trade-off between surveying large areas and obtaining deep samples, measurements of the SMF at high redshift tend to suffer from limited dynamic range. Therefore, it is important to think about how the constraints on our four SHM function parameters arise from the observations.

The characteristic mass M₁ and the maximum SHM ratio (m/M)₀ mostly depend on galaxies and halos of intermediate mass. The high-mass slope γ is fixed by the number of massive galaxies since these live in the massive halos. On the other hand, the low-mass slope β is set by the number of low mass galaxies since these live in the low-mass halos.

For a survey with a fixed area on the sky, the observed volume is smaller for low redshifts (z ≲ 1) than for high redshifts. In order to compute the SMF at high galaxy masses, the observed volume has to be relatively large, as massive galaxies are rare. Thus for low redshifts, one has to choose a wide survey (large area) to determine the SMF for massive galaxies and properly constrain γ. Constraining the SMF at the low mass end requires a high level of completeness for low mass galaxies, which are very faint objects. Hence, we have to choose a deep survey that can detect faint galaxies in order to constrain β.

Taking these considerations into account, we choose the SMFs presented in Drory et al. (2004) to constrain the parameters M₁, (m/M)₀, and γ at low redshifts. The authors derive the SMFs using MUNICS which is a wide area, medium-deep survey selected in the K band. The detection limit is K ≈ 19.5 and the subsample the authors use covers 0.28 deg². We apply our method using these mass functions and take the three parameters from that analysis.

However, the MUNICS survey is not deep enough to detect galaxies that are fainter than the characteristic mass of the SMF (the knee) and thus is not sufficient to constrain the parameter β. To constrain β, we choose the SMFs derived in Fontana et al. (2006). This work is based on the GOODS-MUSIC sample, a multicolor catalog extracted from the survey conducted over the Chandra Deep Field South. The catalog is selected in the z₈₅₀ and K bands, covers an area of 143.2 arcmin², and is complete to a typical magnitude of K ≈ 23.5. We apply our method using the SMFs computed with the z₈₅₀ band selected sample and take the parameter β from that analysis.

For high-redshift (z ≳ 1), we use the SMFs presented in Fontana et al. (2006) to constrain all four parameters. For high redshifts, the volume of a redshift bin becomes large enough to sample massive galaxies, and therefore the GOODS-MUSIC sample is sufficient to constrain γ.

We convert all SMFs which use a Salpeter IMF to the Kroupa/Chabrier IMF.

Having selected the observational SMFs for a set of different redshifts, we fit the four free parameters M₁, (m/M)₀, β, and γ to the observations. The errors on the parameters are derived in a similar way as explained in Section 3.3, but instead of using confidence intervals we have fitted a Gaussian to the probability distributions of M₁, (m/M)₀, and γ and a lognormal to the probability distribution of β.

Figure 12 shows the observed and the model SMFs for different redshifts (indicated at the top of each panel). The values of the resulting four parameters for the different redshifts are shown in Table 6 and the redshift evolution is plotted in Figure 13. The characteristic mass M₁ grows with increasing redshift, while the normalization of the SHM ratio (m/M)₀ becomes smaller with increasing redshift. This means that there is less stellar content in a halo of a given mass at a higher redshift.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The high-mass slope γ can be constrained only weakly. This is due to the limitation of the available galaxy surveys. As the area of the survey is small, the volume in which galaxies are detected is limited, and thus massive galaxies are very rare. This results in large error bars for the SMF for massive galaxies which propagate into the error bars of γ. The situation improves slightly for higher redshifts as the volume of higher redshift bins is larger and thus more massive galaxies can be observed. The value of γ decreases with increasing redshift. For higher redshifts (z > 1), the error bars on γ become very large because of the limited area covered by the available deep surveys (in this case, GOODS). We leave it up to the reader to assess the reliability of our results at z > 1 based on our quoted error bars.

The low-mass slope β seems to increase with redshift until z ≈ 2 and then drops to a low value. However, as the redshift increases it becomes more and more difficult to observe low mass galaxies which are very faint. Thus, the high-redshift values for β are not very well constrained and perhaps not to be fully trusted. We therefore assume that β grows with increasing redshift.

As we explained in Section 4.4, β is strongly related to the parameter α of the Schechter function. A small value of β corresponds to a large absolute value of α, while a large value of β results in a low absolute value of α. This would mean that for higher redshifts the SMF would become shallower, in contradiction with observations (e.g., Fontana et al. 2006 show that the absolute value of α increases with redshift). However, one has to remember that the halo mass function also changes with redshift and becomes steeper. Thus, the halo mass function steepens more than the SMF, so β has to increase in order to compensate.

With the derived parameter values, it becomes possible to interpolate and find the SHM ratio at any redshift. This is done by choosing a redshift parameterization for each of the parameters.

As M₁ and (m/M)₀ do not change much above a redshift of z > 1.5, we choose power laws for the redshift dependence:

$$\begin{equation} \log \, M_1(z) = \log \, M_1\vert _{z=0} \cdot (z+1)^{\mu } \; \end{equation} \tag{ 23 }$$

and

$$\begin{equation} \left(\frac{m}{M}\right)_0(z) = \left(\frac{m}{M}\right)_{z=0} \cdot (z+1)^{\nu } \; \end{equation} \tag{ 24 }$$

with the normalizations M₀ and (m/M)_{z = 0} and the slopes μ and ν.

To parameterize γ over redshift, a linear dependence would lead to a negative γ at a certain redshift. Though this is not forbidden, it leads to a SHM ratio which would be increasing monotonically with halo mass which is inconsistent with feedback processes at the massive end. Hence, we also choose a power-law parameterization for γ:

$$\begin{equation} \gamma (z) = \gamma _0 \cdot (z+1)^{\gamma _1} \; \end{equation} \tag{ 25 }$$

with the normalization γ₀ and the slope γ₁.

From Figure 13, we are not able to infer whether β converges to a constant value. Thus, we adopt a simple linear parameterization:

$$\begin{equation} \beta (z) = \beta _1 \cdot z + \beta _0 \;. \end{equation} \tag{ 26 }$$

Note that we have also tried other parameterizations (constant β, decreasing β) but could not reproduce the observed SMFs. Using the linear parameterization for β and the power laws for the other parameters, we were able to compute SMFs that are in good agreement with the observed ones.

A fit to the derived values presented in Table 6 yields the parameters given in Table 7. As we do not fully trust the derived values of β for z ≳ 2, we neglect these two values and fit a line to the remaining values of β.

Having developed a redshift-dependent model of the SHM relation, we now test this model by computing interpolated SMFs for different redshifts. For this, we use the method described in Section 3. However, now we do not use the parameters that have been derived at each redshift by fitting the model to the observations, but we use the eight parameters of the redshift-dependent SHM relation that have been derived in the previous section.

The resulting interpolated SMFs are compared to the observations (and the fitted mass functions) in Figure 12. For z ≲ 2, we see excellent overall agreement, the interpolated mass functions mostly overlap with the error bars of the observations.

The SMFs for the high-redshifts z ≳ 2 are too low. The deviations are largest at the low mass end. However, if we look at Figure 12, we see that β is higher than the derived value for the two highest redshifts which results in a low-mass slope that is too shallow.

To compare the relation at different redshifts, we use the redshift-dependent SHM relation with the eight parameters that have been derived in the previous section. Figure 14 plots stellar mass versus halo mass for different redshifts. The plot shows that at a fixed low halo mass (e.g., M = 10¹¹ M_☉), galaxies that live in such halos are more massive at low redshift (m ∼ 10⁹ M_☉ for z = 0) than galaxies that live in a halo of the same mass at a higher redshift (m ∼ 10⁸ M_☉ for z = 2). In contrast, massive halos contain more massive galaxies at high redshift, while at low redshifts the galaxies in massive halos have less mass. However, as halos also become more massive over time, one cannot identify a halo of a certain mass at high redshifts with a halo of the same mass at low redshifts. Thus, the fact that at a given (high) halo mass the mass of the central galaxy is lower at present than at an earlier epoch does not imply that individual galaxies lose mass during their evolution. This only means that large halos accrete dark matter faster than large galaxies grow in stellar mass, while the growth of low-mass halos is slower than that of the central galaxies they harbor (see also Conroy & Wechsler 2009). Because of its statistical nature, our model is not suitable for following the evolution of an individual galaxy through cosmic time. We also note that the SHM relation at the massive end (M ≳ 10¹³ M_☉) undergoes very little evolution, which has also been found by Brown et al. (2008).

Zoom In Zoom Out Reset image size

Having determined the SHM relation as a function of redshift, we are now able to populate halos with galaxies at any redshift. We choose a set of redshifts and populate the halos with galaxies, deriving the stellar masses from the redshift-dependent SHM relation. We divide these galaxies into six samples of different stellar mass between log  m/ M_☉ = 8.5 and 11.5. For each of these samples, we compute the real space CF ξ(r) by counting pairs in distance bins (Equation (3)). This leads to six CFs for every selected redshift.

Figure 15 shows the CFs for six different redshifts as a function of stellar mass. We also plot the correlation function of dark matter at the respective redshifts for comparison. For all redshifts, we see that massive galaxies are clustered more strongly than low mass galaxies. The higher the redshift, the more the CFs for different stellar masses differ. For high redshift, there are very few massive galaxies in our limited volume simulation box, and so the error bars become larger.

Zoom In Zoom Out Reset image size

At low redshift (z ≲ 1), observational measurements of stellar mass dependent galaxy clustering have recently been published using the VIMOS-VLT Deep Survey (VVDS) and the zCOSMOS Survey (Meneux et al. 2008, 2009). In order to compare our model predictions to these data, we compute correlation functions for the same stellar mass bins and thresholds and convert these to projected correlation functions as described in Section 3.2. Figure 16 plots the observed projected correlation functions (symbols) and the model predictions (lines) for different stellar mass bins or thresholds in three redshift bins for the zCOSMOS Survey and one redshift bin for the VVDS. There is good general agreement between the model and observations. The zCOSMOS clustering amplitude agrees very well with the model for r_p < 1 Mpc, but for z < 0.8 deviates at larger distances and becomes higher than the prediction. As suggested by Meneux et al. (2009), this may be because the COSMOS field represents an overdense volume at these redshifts. In contrast, the VVDS clustering amplitudes are lower than those predicted by our model, leading to the speculation that perhaps the VVDS represents an underdense region.

Zoom In Zoom Out Reset image size

The bias of any object may be defined as the square root of the ratio between the CF of the object ξ_o(r) and the CF of dark matter particles ξ_dm(r):

$$\begin{equation} b(r)=\sqrt{\frac{ \xi _o(r) }{ \xi _{dm}(r) }}. \end{equation} \tag{ 27 }$$

Here, we focus on the galaxy two-point CF ξ_gg(r, m, z), which in addition to the distance between the galaxies also depends on the redshift and the stellar mass of the galaxies:

$$\begin{equation} b(r,m,z)=\sqrt{\frac{ \xi _{gg}(r,m,z) }{ \xi _{dm}(r,z) }}. \end{equation} \tag{ 28 }$$

From our predicted galaxy CFs, we compute the bias for every redshift and stellar mass by averaging between r = 2 Mpc and 10 Mpc, where b(r) is roughly a constant (as one can see from Figure 15, the scale dependence of the bias is quite weak). Figure 17 shows the redshift dependence of the bias. The symbols represent the averaged value of the bias, while the solid lines correspond to a fit to the symbols. For this, we have used a power-law form:

$$\begin{equation} b(z)=b_0 (z+1)^{b_1}+b_2 \quad, \end{equation} \tag{ 29 }$$

where the parameters b₀, b₁, and b₂ are functions of stellar mass. The fit parameters are given in Table 8.

Zoom In Zoom Out Reset image size

This shows that the bias at a fixed stellar mass increases with increasing redshift. Massive galaxies are biased more strongly than galaxies of lower mass at any redshift. We find that the bias of massive galaxies evolves more rapidly than that of low mass ones (cf. White et al. 2007; Brown et al. 2008). Since the bias of massive halos evolves more rapidly than that of low mass galaxies, this seems to be a feature of any model in which the SHM relation is monotonically increasing (i.e., the most massive galaxies reside in the most massive halos).