The Spitzer South Pole Telescope Deep Field Survey: Linking galaxies and haloes at z=1.5

# The Spitzer South Pole Telescope Deep Field Survey: Linking galaxies and halos at z=1.5

## Abstract

We present an analysis of the clustering of high-redshift galaxies in the recently completed 94 deg Spitzer-SPT Deep Field survey. Applying flux and color cuts to the mid-infrared photometry efficiently selects galaxies at in the stellar mass range , making this sample the largest used so far to study such a distant population. We measure the angular correlation function in different flux-limited samples at scales (corresponding to physical distances Mpc) and thereby map the one- and two-halo contributions to the clustering. We fit halo occupation distributions and determine how the central galaxy’s stellar mass and satellite occupation depend on the halo mass. We measure a prominent peak in the stellar-to-halo mass ratio at a halo mass of , 4.5 times higher than the value. This supports the idea of an evolving mass threshold above which star formation is quenched. We estimate the large-scale bias in the range and the satellite fraction to be , showing a clear evolution compared to . We also find that, above a given stellar mass limit, the fraction of galaxies that are in similar mass pairs is higher at than at . In addition, we measure that this fraction mildly increases with the stellar mass limit at , which is the opposite of the behavior seen at low-redshift.

###### keywords:
cosmology: observations — galaxies: evolution — galaxies: high-redshift — galaxies: halos — large scale structure of universe
12

## 1 Introduction

Many observational studies have measured dark matter halo masses in order to find correlations with the properties of the galaxies they host. Various works have utilized gravitational lensing of background objects (Mandelbaum et al., 2006; Gavazzi et al., 2007; Bolton et al., 2008; Auger et al., 2010; Cacciato et al., 2009, 2013; Velander et al., 2011), virial temperatures derived from X-rays (Lin et al., 2003; Lin & Mohr, 2004; Peterson & Fabian, 2006; Hansen et al., 2009) and dynamics of satellites (More et al., 2009, 2011). These methods have achieved high accuracy, but are also observationally expensive to carry out on large samples and for small haloes, which limits the statistical strength and range of application. A less direct but more comprehensive method of linking galaxies to haloes is abundance matching (Conroy et al., 2006; Vale & Ostriker, 2006; Moster et al., 2010; Guo et al., 2010; Behroozi et al., 2010; Moster et al., 2010, 2013), which uses the merger trees from -body dark matter simulations as input and assumes that the halo mass is the main determinant of galaxy luminosity and stellar mass. The basic idea is to cumulatively match observed galaxy luminosity functions and halo mass functions by placing progressively less luminous galaxies in less massive haloes. By design, this method reproduces the luminosity (or stellar mass) function, and is able to predict the clustering of galaxies in many cases (Conroy et al., 2006; Conroy & Wechsler, 2009; Moster et al., 2010).
Direct measurements of galaxy clustering are another powerful way to connect galaxies with the underlying dark matter distribution. As a function of physical separation , clustering is commonly measured in the form of the two-point spatial correlation function (SCF; Peebles 1980). The relation between the distributions of galaxies and dark matter can be parametrized through the galaxy bias (Kaiser, 1984; Coles, 1993; Fry & Gaztanaga, 1993; Mo & White, 1996; Kauffmann et al., 1997; Sheth & Tormen, 1999; Tinker et al., 2005, 2010), which is given by the scaling between the SCFs of these two fields:

 ξg(r,z)=ξm(r,z)b2g(r,z). (1)

The SCF of dark matter depends on the cosmology, and can be prescribed analytically given those parameters (Eisenstein & Hu, 1999; Smith et al., 2003). Thus, the bias of a galaxy sample is directly determined by its SCF. In general, the bias depends on the spatial scale and redshift (Fry, 1996; Moscardini et al., 1998; Tinker et al., 2005; Moster et al., 2010), since galaxies and dark matter do not evolve in the exact same manner in time or space. The measurement of can therefore reveal a precise description of the connection between galaxies and dark matter.
Many studies up to intermediate redshifts () have investigated galaxy clustering with samples selected in different ways (Phleps et al., 2006; Zheng et al., 2007; Coil et al., 2008; Blake et al., 2008; Brown et al., 2008; McCracken et al., 2008; Meneux et al., 2008, 2009; Simon et al., 2009; Ross et al., 2010; Foucaud et al., 2010; Abbas et al., 2010; Zehavi et al., 2011; Matsuoka et al., 2011; Wake et al., 2011; Jullo et al., 2012; Leauthaud et al., 2012; Hartley et al., 2013; Mostek et al., 2013; de la Torre et al., 2013; Donoso et al., 2014). The most common conclusion is that clustering strength is correlated with luminosity, red color and morphology (towards early-type). Galaxies on the extreme of these properties are highly biased and therefore they live in massive haloes.
These conclusions can be obtained just by analyzing the overall amplitude of the bias. However, the precise form of this observable as a function of spatial separation contains more information about the inner structure of the haloes. The halo occupation distribution (HOD) is a simple parametric framework to accurately model the bias (Ma & Fry, 2000; Seljak, 2000; Peacock & Smith, 2000; Scoccimarro et al., 2001; Cooray & Sheth, 2002; Berlind & Weinberg, 2002; Berlind et al., 2003; Kravtsov et al., 2004; Zheng et al., 2005). It considers galaxies to be either centrals or satellites, and the number of these that a halo can host is fully determined by the halo mass.
One of the advantages of the HOD framework is that its parameters have a clear physical meaning, and thus when fitting the clustering one can gain a deeper insight into the connection between the galaxies and their host haloes. For example, the HOD framework can directly relate the average stellar mass of the central galaxies to a particular halo mass. As shown in many studies at , the ratio of these masses is highest around a halo mass of (Zheng et al., 2007; Yang et al., 2012; Zehavi et al., 2012; Leauthaud et al., 2012; Reddick et al., 2013; Behroozi et al., 2013; Moster et al., 2013; Wang et al., 2013). This implies that there is a characteristic halo mass where galaxy formation has been more efficient. The qualitative explanation for this is that at low halo masses the gravitational potential is not deep enough to halt the expulsion of gas due to stellar winds (Benson et al., 2003), while high-mass haloes have heated up the intra-halo medium by gravitational heating and AGN feedback (Croton et al., 2006; Bower et al., 2006; van de Voort et al., 2011a) so that infalling gas gets heavily shocked and cannot easily cool and condense (Birnboim & Dekel, 2003; Dekel & Birnboim, 2006; Kereš et al., 2005, 2009). These two trends can be reduced to a comparison between dynamical and gas cooling times in haloes, such that for low masses and for high-masses. A possible consequence is that the peak halo mass is related to a characteristic quenching mass that sets (Neistein et al., 2006) and marks a transition between star forming and quenched haloes. Indeed, massive red galaxies with little star formation have been shown to live in massive haloes (Coil et al., 2008; Zehavi et al., 2011), supporting the idea of the red sequence of galaxies arising when they become quenched (Croton et al., 2006; Bower et al., 2006). This blue/red dichotomy is present in the nearby Universe (Kauffmann et al., 2003, 2004; Baldry et al., 2004), and starts its build-up around (Bell et al., 2004; Cooper et al., 2006; Muzzin et al., 2013; Wang et al., 2013). Thus, when haloes become large enough, they quench their star formation. A consequence of this is that the most massive galaxies today have no significant ongoing star formation. This effect has been called archeological downsizing (Cowie et al., 1996; Juneau et al., 2005; Conroy & Wechsler, 2009), and is also inferred from the lack of evolution in the massive end of the stellar mass function (Pérez-González et al., 2008; Marchesini et al., 2009; Muzzin et al., 2013). HOD models have shown that the stellar-to-halo mass ratio (SHMR) evolves in the sense that the peak moves to lower halo masses with increasing time, at least since (Coupon et al., 2012; Leauthaud et al., 2012). This trend has been predicted to persist up to by extensions of HOD that use conditional stellar mass functions (Yang et al., 2003, 2012; Wang et al., 2013) and abundance matching studies (Moster et al., 2013; Behroozi et al., 2013). A possible mechanism for this would involve evolution in , which is supported by the idea that the universal gas fraction drops with time and therefore star formation becomes more difficult with time at fixed halo mass (van de Voort et al., 2011a, b). However, this is still a matter of debate (Conroy & Ostriker, 2008; Tinker & Wetzel, 2010). For instance, Leauthaud et al. (2012) present evidence in favor of this evolution being set by quenching below a critical galaxy-halo mass ratio instead of a critical halo mass. Such a mechanism would also shift the SHMR toward lower masses with time.
We have described the basic processes that can determine , based on the comparison of and as a function of halo mass. This basic model can be extended to include modes of galactic outflows, which are then directly constrained by the observed slope of the SHMR. The stellar mass growth of a galaxy is heavily regulated by the expulsion of gas, which could be mainly sourced by supernovae feedback (Murray et al., 2005). The stellar mass loss rate, , can be broken down in two contributions: pressure-supported energy injection (energy-driven winds) and coherent momentum transfer (momentum-driven winds). The energy and momentum deposition rates, and , can be related from first principles to the mass via a proxy of the kinetic velocity field, : and . This suggests that galaxies with low velocity fields, and therefore low masses, may have their outflows dominated by energy-driven winds (Dutton & van den Bosch, 2009). Thus, a larger contribution from this type of wind would result in a steeper low mass slope of the SHMR.
At high masses (and high ), these arguments would point to a dominance of momentum-driven winds. However, the winds in this regime are also sourced by radiative AGN feedback, which is expected to have a strong contribution (Vogelsberger et al., 2013). In addition, a large merger rate between central galaxies will result in a flattening of the SHMR (Leauthaud et al., 2012). With all these processes at play, the high-mass slope is less straightforward to interpret than the low-mass one, but it can still offer important constraints on this combination of mechanisms.
Galaxy clustering combined with HOD modeling provide particularly solid measurements of the SHMR whenever the selection of galaxies spans the relevant range of stellar masses. At , such measurements have proven to be very difficult given the lack of large volume-limited samples. Wake et al. (2011) use the 0.25 deg NEWFIRM survey (van Dokkum et al., 2009), but the low number statistics made it difficult to map the turnover of the SHMR. In this study, we use a 94 deg mid-infrared survey to select galaxies with stellar masses ranging from and fit an HOD model to the angular correlation function. We present the most robust measurement to date of the peak of the SHMR at .
In addition, the HOD yields particularly strong constraints on the satellite population of a given galaxy sample. We determine what fraction of the galaxies are satellites, and how the abundance of these depends on the halo mass. Moreover, we measure a proxy for the occurrence of galaxy pairs of similar mass, and find that it mildly decreases toward high luminosities. Although we do not achieve a robust detection, this represents the opposite trend to what is seen at low-redshift. The processes that produce this relationship are strongly tied to the accretion and merger events between galaxies and haloes, as well as the quenching of star formation in satellites.
The paper is organized as follows. In Section 2, we describe all datasets that are used. In Section 3 we describe how we adapt redshift and stellar mass distributions from a reference optical + mid-IR survey. In Section 4 we define the two-point clustering statistic and the method used to compute it. In Section 5, we describe the model that links galaxies to haloes. In Section 6, we explain the fitting procedure of the HOD to the observed clustering. In Sections 7, 8 and 9 we discuss the results obtained regarding the SHMR, the satellite galaxies and the large-scale bias, respectively. We end with a short summary in Section 10. For the reader that is only interested in the results, we recommend reading Sections 7 and beyond.
Additionally, we include several appendices where many of the details are covered. Appendix A presents a calibration of systematic effects in the photometry. Appendix B compares the results obtained from using different reference catalogs to draw redshift and stellar mass distributions. Appendix C calculates the systematic offset in the clustering amplitude due to the geometry of the survey. Appendix D presents the formalism of the halo model. Appendix E investigates the removal of low-redshift sources from the sample using optical data. Appendix F explores different choices of free parameters used in the HOD fits to the clustering.
Throughout this paper we use the following cosmology: , and Mpc. All magnitudes are in the Vega system and masses are in units of .

## 2 datasets

Our main dataset is the Spitzer South Pole Telescope Deep-Field Survey (SSDF; Ashby et al. 2013b), a 93.8 deg photometric survey using the IRAC 3.6 and 4.5m bands (hereafter [3.6] and [4.5]). The mosaics have a nominal integration time of 120 seconds. We used Source Extractor (SExtractor; Bertin & Arnouts 1996) in dual image mode, detecting galaxies in [4.5] and extracting the flux from fixed 4 apertures in both IRAC channels. These aperture fluxes were then corrected to total fluxes using growth curves from isolated point sources found in the mosaics. A detailed description of the survey and a public photometric catalog are presented in Ashby et al. (2013b). However, here we use a deeper private catalog and account for faint-end photometric bias and detection completeness (see Appendix A). We determine the 5 limit in [4.5] to be 18.19 mag, in agreement with Ashby et al. (2013b).
We use the near-infrared 2MASS Point Source Catalog (Skrutskie et al., 2006) to identify and remove sources brighter than mag, most of which are likely to be stars. In addition, we visually inspected some of these sources in the IRAC mosaics and determined an empirical relation between their -band magnitude and the maximum radius where their 4.5m flux caused a clear suppresion in the detection of nearby sources. This relation was then applied to the rest of the -selected sample and the resulting radii were used to mask all SSDF sources enclosed within from the main catalog. For reference, the radii corresponding to and sources were 41 and 8.4  respectively. We also masked out low coverage gaps in the survey, yielding a final effective area of 88.8 deg.
Finally, in order to better understand the redshift distribution of our IRAC-selected sample in the SSDF, we use public catalogs in two other regions of the sky: the COSMOS-UltraVista field (hereafter COSMOS, Muzzin et al. 2013) and the Extended Groth Strip (hereafter EGS, Barro et al. 2011a,b). These two surveys have publically accessible IRAC photometry, photometric redshifts, and stellar masses. In the following Section we describe how we used these catalogs to infer the redshift and mass distributions of SSDF samples.

## 3 Control Sample

This study requires knowledge of the redshift and stellar mass distribution of the SSDF galaxy sample. However, our main dataset is too limited to obtain reliable values for these observables. Therefore, the strategy is to import this information from a reference survey that contains optical data and IRAC photometry with a higher accuracy. We consider the catalogs from COSMOS and EGS, which include photometric redshifts and stellar masses. We will adopt COSMOS as the fiducial dataset because it is larger and has better statistics, and in Appendix B we show how our results do not change significantly when using EGS instead. The reference catalog is degraded to become a control sample whose photometric errors match those of the SSDF. Then, applying the same IRAC selection in both SSDF and the control sample allows us to match the derived distributions of redshift and mass. A brief description of the COSMOS photometry can be found in Appendix B.

For every source in the reference catalog, we have magnitudes, colors, photometric redshifts and stellar masses. The goal is to infer the SSDF distributions of these parameters by degrading the reference photometry, which is done using the SSDF photometric errors. We calculate the scatter in SSDF magnitudes and colors as a function of these same variables, using the results from the photometric simulations described in Appendix A. These scatter profiles are shown in Figure 1. At fixed [4.5] magnitude, the scatter in color increases for larger colors since these imply fainter [3.6] magnitudes. In the case of the reference sample, since it is 2 magnitudes deeper than SSDF (see Appendix B), we can safely consider its photometric scatter as negligible in comparison.
The degradation of the reference catalog into a control sample consists of transforming the specific values (e.g., magnitude) of each source in the catalog into Gaussian probability density functions (PDFs). These PDFs are defined in the parameter space of apparent magnitude [4.5] (), color () and photometric redshift (): . The centroids are given by , which correspond to the parameter vectors of the sources in the reference catalog. The standard deviations are . The first two components in are the functions and , which are shown by the curves in Figure 1. The redshift component does not have a counterpart in the SSDF catalog, but we apply a variable redshift smoothing kernel equivalent to 100 comoving Mpc, in order to filter the effect of large-scale structure. This amounts to within our redshift range. However, we find that this redshift filtering has a minimal effect in the results, varying the clustering amplitude and galaxy number density at the level.

### 3.1 Main redshift distribution

If we consider galaxies with apparent magnitudes within some bracket , we can compute the distribution in color and redshift space of the control sample:

 K(C,z)=1NrefNref∑j=1∫ΔMdmP(m,C,z;→μj,→σj). (2)

Here, we have marginalized each individual PDF over and summed them in the resulting space of , using the reference catalog (subscript “ref”). A similar procedure to derive full redshift distributions based on the Bayesian combination of individual redshift likelihood functions was performed by Brodwin et al. (2006a, b). The normalization of is the total number of sources in the reference catalog, . Figure 2 shows the application of Equation (2) for , which are the limits for our full SSDF sample (see Section 3.2). The top panel correspond to the color versus redshift distribution from the raw reference catalog. The lower panel shows the control sample, which is how SSDF sources are expected to be distributed. For comparison, we have also plotted a galaxy evolutionary track for a single stellar population with solar metallicity and formation redshift of , computed using the Bruzual & Charlot (2003) models with Chabrier (2003) IMF.
There is a clear correlation between color and redshift at . This occurs because going from to , the IRAC bands map the galaxy spectrum across the stellar bump at rest-frame -band. This results in a monotonic change in observed color within . An insightful description of this phenomenology can be found in Muzzin et al. (2013). We can take advantage of this effect to select galaxies in redshift using a color cut. A lower color threshold needs to be high enough to reject galaxies (see Figure 2), while also keeping a number of higher redshift sources that is large enough to measure a robust clustering signal. An upper threshold is also necessary, since very red colors are characteristic of active galactic nuclei (Stern et al., 2005). The best compromise is a color cut of , as shown in Figure 2.
With a given , we can derive the redshift distribution of sources:

 ϕ(z)=1NrefNref∑j=1∫ΔM∫ΔCdmdcP(m,c,z;→μj,→σj). (3)

Note that is equal to 1 only when and represent the full ranges spanned by the reference sources. We denote such distribution as , while the one corresponding to the color selection is denoted as .
We are assuming that this color cut selects a representative sample of the galaxy population. However, it is important to check whether such a selection is biased toward young or old galaxies. It has been shown that older galaxies exhibit a higher clustering amplitude than their younger counterparts (Skibba et al., 2014, and references therein). In our case, we are tracing a part of the spectral energy distribution that is much less sensitive to star formation history. To illustrate this point, we use the EzGal package (Mancone & Gonzalez, 2012) to compare the colors for passive and star-forming galaxies using assorted stellar population models (Bruzual & Charlot, 2003; Maraston, 2005; Conroy et al., 2009). For the passive galaxies we assume a single burst model with formation redshift . For the star-forming galaxies we run models with exponentially declining star formation, using  Gyr and an initial formation redshift . The difference in is . For the galaxy sample used in our analysis, this difference is comparable to or smaller than the photometric errors (see Fig. 1). Thus, the associated systematic bias due to the color selection will be small compared to the photometric errors and intrinsic scatter in galaxy colors, and can be neglected for the current analysis.

### 3.2 Definition of subsamples

Our main science sample of SSDF galaxies is determined by the apparent magnitude and color cuts of and . The first of these cuts imposes an upper magnitude limit at the 80% completeness level (see Appendix A.3), and the second is tuned to select galaxies at high-redshift while avoiding AGN. A density map of this selection can be seen in Figure 3, representing a slice of the Universe at .
We further split the main sample into 13 subsamples, with faint limits over the range in steps of 0.2 mag. The bright limit is in all of them. We do this instead of a selection within differential magnitude bins because the halo occupation framework presented in Section 5 requires cumulative samples in order to link halo masses and galaxy masses. We note that this approach carries the drawback of producing a correlation between the different samples. This correlation is strong between neighboring samples, but not dominant otherwise. Due to the steep variation of the stellar mass function, any given sample is mostly comprised by galaxies close to its low-mass threshold, making their clustering less sensitive to the most massive population (see Matsuoka et al. 2011).
The photometric scatter increases for fainter samples. Thus, we calculate the redshift distribution (see Equation 3) for each sample, obtaining sets of , , where is the sample index (going from brightest to faintest). We can also define the number density completeness as , which determines the fraction of galaxies as a function of redshift that the color cut retains. Figure 4 shows a comparison of for (the largest sample). At the peak of the color-cut distribution we have that 0.3, and we will use this factor to scale up and correct the number density (see below). Figure 5 shows for the smallest and largest samples (i.e., brighter and fainter thresholds, =1,13), where each curve is shown normalized to 1. We also derive cosmic variance errors using the prescriptions from Moster et al. (2011), which are based on analytical predictions of dark matter structure given a particular survey geometry (see also Brodwin et al. 2006a). The peak in these redshift distributions is consistently around in all samples. In general, the samples consist of a population that has approximately the same absolute luminosity and stellar mass (see Section 3.3), plus a contribution of “contaminant” galaxies that are intrinsically much less luminous. These contaminants represent 12% (37%) of all galaxies in our full (brightest) sample. When setting a brighter flux threshold, the high-redshift population becomes less dominant since these galaxies are closer to the turnover of the luminosity function. The consequence of this is the clear trend where brighter samples have a stronger low-redshift bump. The contribution of the latter to the clustering is modeled in the following Sections. Alternatively, we show in Appendix E that our results remain unchanged if instead we employ shallow optical data to remove most of the low-redshift sources.
With these redshift distributions, we can calculate the spatial number density of observed galaxies at the pivot redshift . Here, we use the SSDF survey area and the effective number of observed sources in the subsamples, . This number is derived by summing the inverse of the completeness value for all galaxies, using the relation from Figure 17. Then, the true number of galaxies within can be written as

 Nktrue=NkobsfkN(zp)ϕkcut(zp)∫ϕkcut(z′)dz′δz. (4)

 V=dV(zp)dzδz=cΩχ2(zp)H(zp)δz, (5)

where is the comoving radial distance, is the Hubble function, is the speed of light and steradians is the solid angle subtended by the survey. Hence, the number density at results in

 nkg=NktrueV. (6)

Note that this quantity is the result of combining the SSDF observed number counts (via ) and the color fractions of the control sample. Table 1 shows the values of these number densities for all samples.

### 3.3 Stellar masses

The stellar masses in the reference catalog are also retrieved to construct our control sample. We use those based on Bruzual & Charlot (2003, hereafter BC03) stellar grids, Chabrier (2003) IMF and Calzetti et al. (2000) dust extinction. Unless otherwise noted, all stellar masses are given under these prescriptions.
For the purposes of this paper, we need to calculate stellar masses for two different selections of galaxies. One is the median mass of all galaxies within each sample, derived at every redshift bin, . This mass will be used to derive a redshift scaling of the galaxy bias in Section 5.2. The other, , is the median mass of the galaxies at the pivot redshift () and at the magnitude limit of each sample. This is the stellar mass that will be linked in Section 7 to a particular halo mass.

Ideally, in order to calculate we would compute the median mass of those galaxies within the given selection range in the parameter space of redshift, magnitude and color. However, classifying galaxies on whether they fall in that range is not straightforward, since each galaxy is represented by an extended probability distribution in the parameter space. A more insightful approach is to compute the probability that a galaxy’s true parameter vector falls within the specified range, which reads as

 ϱj(z)=∫ΔM∫ΔCdmdcP(m,c,z;→μj,→σj), (7)

where and is an index that identifies every galaxy in the reference catalog. Thus, we can calculate as the weighted median stellar mass over all galaxies in the reference catalog, where the set of act as weight coefficients to the individual stellar masses . By definition, the weighted median mass represents the mass value where the weighted integral of the mass distribution above and below that value is the same, and we write it in the following condensed form

 Missing \left or extra \right (8)

where we have omitted the implicit dependence on . Conveniently, the weighted distribution of masses follows closely a log-normal distribution. Thus, Equation 8 returns almost the same value as the weighted mean of , which allows us to adopt the standard deviation from the latter distribution:

 ΣM=∑(logMj⋆−log¯Mfull⋆)2ϱj∑ϱj, (9)

which is typically dex. Even though this scatter is rather large, these log-mass distributions are single-peaked and approximately symmetric, so their mean value is well-defined and physically meaningful. We use the scatter to estimate the error in the mean as

 ¯ΣM=ΣM√Nind (10)

with

 Nind=(∑ϱj)2∑ϱj2. (11)

Here, represents the effective number of independent elements in the ensemble. This number is proportional to the sum of contributing weights and inversely proportional to their scatter. It equals the total number of elements in the reference catalog in the limit of .
Figure 6 shows for all samples. The errors are from Equation (10) and the solid curve is a 5th order polynomial fit to the points of the faintest sample. We use that curve plus an offset to fit the data from the rest of the samples, since it becomes noisier at brighter limits. Here, we take advantage of the fact that stellar mass scales with flux approximately in a linear manner. It is clear from the figure that the mass is tightly correlated with the redshift of observation within . Beyond that, the relation flattens out significantly. The reason for this is that at , the [4.5] band samples the rising spectral slope of the stellar bump (Muzzin et al., 2013). This offsets the k-correction in a way that galaxies of a certain intrinsic near-infrared luminosity have a similar apparent [4.5] magnitude across a range of redshift. A consequence of this is that any [4.5] limited sample becomes roughly stellar mass limited at (see also Figure 14 in Barro et al. 2011b). Nonetheless, we do not attempt to take advantage of this effect by averaging stellar masses at high-redshift. The modeling in this work is based on well-defined median masses as a function of redshift, independent of the form of that redshift dependence. However, the flattening of this curve does benefit our study to some extent. Since there is an inherent uncertainty in how well represented the SSDF data is with the control sample, it is convenient that the stellar masses are naturally more constrained than a case where they had a strong redshift dependence.
We can calculate at in an analogous way, considering a selection within . The corresponding weights are

 κj(M)=∫ΔCdcP(M,c;zp,→μj,→σj), (12)

and replacing with in Equations (7-10) gives us the median mass , along with its error. An example of the stellar mass histograms that are obtained using these weights can be seen in Figure 7. They are shown for the brightest and faintest magnitude limits of our samples. These distributions are symmetric and have a well-defined mean. The scatter is generally large, with values around 0.2 dex. However, the error in the logarithmic mean (see Equation (10)) is typically quite small 0.03 dex. The values of these limiting masses are displayed in Table 1. Thus, we have calculated the median mass of all galaxies in our samples as a function of redshift, and the median mass of galaxies around the pivot redshift at each sample magnitude limit.

## 4 two-point clustering

Given a population of galaxies in a three-dimensional space, one can define the joint probability of finding two such objects in volume elements , separated by a distance (Peebles, 1980; Phillipps et al., 1978):

 δP(r)=¯N2(1+ξg(r))δV1δV2. (13)

Here, is the density of galaxies and is the SCF, which quantifies the clustering strength of the field as a function of . The SCF can also be interpreted as the differential probability of finding two objects separated by a given distance, with respect to the case of a random distribution.
The SCF for a galaxy population can be directly computed if the individual distances (redshifts) to those galaxies are known. However, in our case we are limited to individual sky positions and the ensemble redshift distribution. Therefore, we are interested in the angular correlation function (ACF), which is the projection of the SCF onto the 2D sphere. Analogously to the SCF, the ACF represents the differential probability with respect to a random distribution of finding two galaxies separated by a particular angle. The ACF is related to the SCF through the Limber projection (Limber, 1953), which integrates the SCF along the line of sight using the normalized redshift distribution as a weight kernel (Phillipps et al., 1978; Coupon et al., 2012):

 ω(θ)=2c∫∞0dzH(z)ϕ2(z)∫∞0dyξg(r=√y2+D2c(z)θ2), (14)

where is the radial comoving distance, is the Hubble function, is the speed of light and is the angular separation given in radians.
In order to measure we use the estimator presented in Hamilton (1993), which counts the number of galaxy pairs with respect to those of a random sample distributed in the same geometry:

 ^ω(θ)=RR(θ)GG(θ)GR2(θ)−1, (15)

where GG, GR and RR are total number of galaxy-galaxy, galaxy-random and random-random pairs separated by an angle . We have also tested the estimator from Landy & Szalay (1993), which returns results that are practically indistinguishable from those using Equation 15.
In order to account for the completeness shown in Figure 17, we make a small generalization of Equation 15. Instead of counting all pairs with values of 1, we use a weighted scheme where each pair of sources is counted as a product of weights . Random sources have , and galaxies have weights equivalent to the inverse of the completeness value at its apparent magnitude. We counts pairs by brute force in discrete angular bins using the graphics processing unit (GPU) on a desktop computer. We have developed our own code, which yields computation times of the order of 1000 times faster than using a CPU-based run with 16 cores. Our code is written in PyCUDA3, which is a Python wrapper of the CUDA, the programming language that interfaces with the device.
The estimator in Equation 15 implicitly assumes that the average galaxy density of the survey is the same as the all-sky value. However, since the survey is a small fraction of the sky, its density is higher (structures cluster more toward smaller scales) and this results in a systematic suppression of . We correct for this effect, even though it is not significant for our results. Details can be found in Appendix C. The values of the corrected ACF for all samples are displayed in Table 4.

### 4.1 Error estimation

We estimate errors with the jackknife technique, which uses the observed data and is very effective in recovering the covariance of between different scales. First, the entire sample is divided into spatial regions of equal size. Then, the correlation is run times, each one excluding one of those regions from the sample. The value of the estimator is the average of those iterations and the covariance between angular bins is given by (Scranton et al., 2002)

 Cjk=N−1NN∑i=0[^ωi(θj)−¯ω(θj)][^ωi(θk)−¯ω(θk)]. (16)

We also compare the jackknife errors with those obtained from mock simulations, which are described in Appendix A. We find that both sets of errors have a good agreement, with differences around 20%. Although our mock simulations only cover large-scales, the systematic differences between mock and jackknife errors are not expected to vary significantly across different scales for a projected statistic like (Norberg et al., 2009).

## 5 placing galaxies in haloes

The galaxy bias (see Equation 1) encodes all the information that can be extracted from the two-point galaxy distribution, given a particular cosmology. Thus, our aim is to construct a precise model of and adjust the resulting correlation function to match the observed clustering of galaxies. The main idea behind this model is to assume a halo distribution and place galaxies in haloes according to a set of simple rules, as explained below.

### 5.1 The Halo Occupation Distribution

The distribution of dark matter haloes under the CDM paradigm has been well-studied both phenomenologically and through simulations (Ma & Fry, 2000; Cooray & Sheth, 2002; Berlind & Weinberg, 2002), leading to a halo model where the halo mass function, the bias and the halo density profile are determined by the halo mass. The Halo Occupation Distribution (HOD) is a statistical framework that has been developed to link the halo model with the distribution of galaxies (Berlind & Weinberg, 2002; Cooray & Sheth, 2002; Kravtsov et al., 2004). The HOD is mainly described with the probability that a halo of a given virial mass hosts galaxies; one central and satellites distributed according to a NFW profile. All galaxies are linked to some halo, and the occupation is independent of their formation history and environment (Zentner et al., 2005). This assumption is generally valid, since the induced changes in the galaxy bias due to environment are expected to be only at the 5% level (Croton et al., 2007; Zu et al., 2008), while the overall uncertainties in galaxy clustering studies are typically larger. For our work in particular, the main source of error arises from the uncertainty in the shape of the redshift distribution, which is explored in Appendix B by comparing results from the use of COSMOS and EGS as reference catalogs. The variations in galaxy bias are around 10%  and they do not alter qualitatively any of the final conclusions. Therefore, given that the environmental effects in the galaxy bias are expected to be smaller, we consider them negligible for the current purposes.
The average distribution of central galaxies as a function of halo mass can written as (Zheng et al., 2005, 2007):

 Nc(M)=12[1+erf(logM−logMminσlogM)]. (17)

This implies that . Thus, sets a step-like transition where half of the haloes above this mass will host a central galaxy, and this transition is smoothed by the scatter . The number of satellites galaxies is drawn from a Poisson distribution with mean

 Ns(M)=Nc(M)(M−M0M′1)α, (18)

and are assumed to follow a NFW (Navarro et al., 1997) density profile from the halo center. The factor accounts for the constraint that only haloes with a central galaxy may host satellites. Equation 18 represents a power law, where sets the steepness, defines the typical mass scale for this distribution being close to unity and represents the mass below which the power-law is cut off. In addition, one can derive the characteristic mass where a halo hosts exactly one satellite on average, , by imposing and noting that generally . In the case where it reduces simply to , and when then . The occupation distribution of the total number of galaxies in a halo can be expressed as the sum of the central and satellite terms:

 N(M)=Nc(M)+Ns(M). (19)

Other HOD derived quantities are the effective galaxy bias

 beffg=1ng∫dMdn(M)dMN(M)bh(M), (20)

and the fraction of satellite galaxies

 fsat=1ng∫dMdn(M)dMNs(M). (21)

Further details about the halo model used here can be found in the Appendix D.

### 5.2 Redshift scaling

We aim to fit the halo model at the pivot redshift . However, the galaxies in our samples have redshift distributions that are too broad to be neglected or averaged over (see Fig. 5). Thus, our approach is to produce , and scale it using a simple prescription to generate at all other redshifts. We then use Equation 14 to make a redshift projection of onto .
The scaling we apply is based on how the large-scale clustering (represented by the two-halo term , see Equation 49) varies with redshift. This change in amplitude is driven by the growth factor of the dark matter and the galaxy bias . Hence, we can write

 ξg(z)=G2(z)G2(zp)b2g(z)beffg2ξpg, (22)

where and are set at by construction. Here, we have made the approximation that the entire correlation function can be scaled with a single factor. However, the relative amplitude of the one and two-halo terms is known to evolve (Conroy et al., 2006; Watson et al., 2011), in the sense that typically the one-halo term is more prominent at higher redshift. We have tested how would change if we allow for some differential redshift scaling between the one and two-halo terms of . This was done applying a linearly redshift-dependent factor to the one-halo term, in addition to the general scaling from Equation 22. In this way, below and above the one-halo becomes reduced and boosted, respectively. We find that this has a very little effect on the resulting ACF. This is because the redshift distributions of our galaxies are more or less symmetric, so that the relative scaling of the one-halo above and below is almost canceled when these contributions are summed together. In reality, this relative scaling might have a more complex dependence on redshift, but we believe that the linear representation we considered here is adequate given the symmetric and peaked forms of our redshift distributions. Thus, we find that our model is not sensitive to the particular evolution of the one-halo term and do not incorporate it in the determination of our results.
Our general approach is to calculate the bias as a function of the evolving median stellar mass, . For this purpose, we make use of the galaxy bias as a function of stellar mass and redshift presented in Moster et al. (2010, hereafter M10), . However, we do not use their bias values directly since we need to enforce that , i.e., the bias function has to match the HOD bias at the redshift of the fit. Our bias function is normalized to hold that constraint, but the scaling at other redshifts is adopted from M10 (for a given stellar mass). To accomplish this, first we define the stellar mass where . Ideally, would be equal to the median mass of the sample from Section 3.3, , but they differ. This is not surprising, since the modeling in M10 is based on abundance matching, which is different from our clustering approach and can potentially yield differing values of the bias. In addition, some variations are expected given the differences in models and codes used to derive stellar masses in M10 and our reference sample. However, the M10 masses by themselves are not relevant to us, and they simply represent a quantity or label that links brighter populations of galaxies with a larger bias. Therefore, it is sufficient to assume a priori that all these masses hold a monotonic relationship with sample luminosity, which has been proven correct a posteriori. In other words, does scale monotonically with across all samples. So, for a given sample, what we calculate is the offset at . In this way, we are able to “convert” our stellar masses into M10 masses. Our bias function then becomes

 bg(z)=bM10g(log¯Mfull⋆(z)+ΔlogM⋆,z). (23)

In Section 3.3 we calculated for all samples. For the brightest ones, the data becomes a bit noisy due to the low number density (see Figure 6). This mass distribution is consistent with having a constant shape and varying it by some normalization that scales with the magnitude limit of the sample. Thus, we adopt the functional form of the largest sample , which is given by the polynomial fit shown in Figure 6. The normalization of this function does not need to be taken into account, since it will be implicitly incorporated in .
The stellar mass dependence of the bias has an important effect on . Because the stellar mass of our samples is larger at high-redshift (about 10 times larger at than at ; see Figure 6), the bias will place a stronger weight there than at low-redshift. This helps to minimize the contribution of the undesired low-redshift bump at . Additionally, there are other functions weighing in the Limber projection, which is proportional to , as inferred from Equations 14 and 22. Figure 8 shows the comparison between the redshift distribution and the full Limber kernel. It is shown for the brightest sample because it is the one with the highest fraction of low-redshift contaminants. In the end, the low-redshift contribution to the clustering is minimized due to the decrease in , which suppresses . This effect is convenient for our analysis, since it makes the clustering properties of our samples highly representative of the Universe. Moreover, in Appendix E we investigate how the final results are impacted by the use of available optical data in the SSDF field to remove low-redshift sources. We find that the changes in the results are negligible compared to keeping these sources and modeling their weak contribution to the clustering, as done in this Section.
A possible concern at this point is that the results from M10 would be “built-in” to ours through the coupling with Equation 23. However, the normalization of the bias is set by our own data, and it is the redshift modulation that we incorporate from these authors. In addition, we have explored variations of and determined that our model is not very sensitive to such changes. As seen in Figure 8, the redshift modulation plays a role in weighing galaxies at vs . Basically, any function that down-weights the low-redshift bump will do so in a manner that it becomes quickly subdominant. We just need a function that reflects approximately the variation of the bias with stellar mass and redshift, which is precisely what is provided by scalings from M10. Our model does not strongly depend on the detailed form of this function, and we have verified that our results are not pre-set in a significant way by those in M10.

## 6 HOD Model fits

The fitting procedure is based on maximizing the likelihood of the model given the observable modobs, with

 χ2=N∑i=0N∑j=0[ωm(θj)−¯ω(θj)]C−1ij[ωm(θk)−¯ω(θk)]. (24)

Here, and are the model and observed ACFs, is the covariance matrix from Equation 16 and is the number of angular bins. The halo occupation model we consider has a total of 5 parameters: , , , and . Even though the signal-to-noise of our ACFs is very good ( with respect to the null hypothesis), the fact that it is the result of projecting the SCF across a wide redshift distribution reduces our constraining power on the HOD model. Thus, to avoid over-ftting the data, we choose to fix a number of parameters. We have run sets of Monte Carlo Markov chains to explore the sensitivity of the model to different choices of constraints. To evaluate this sensitivity, we use the Akaike criterion (Akaike, 1974), which states that an extra free parameter is justified only when the new best-fit is reduced by an amount larger than 2. For either and , this criterion is not fullfilled. Thus, we follow Conroy et al. (2006) and set . We also fix , following a number of studies that support typical values (More et al., 2009; Behroozi et al., 2010; More et al., 2011; Wake et al., 2011; Moster et al., 2013; Reddick et al., 2013; Behroozi et al., 2013). In the case of , we have that , which would mildly favor setting it free. However, this parameter has an intrinsic degeneracy with and when left free to float, the best fit values show a significant stochastic component in their behavior with respect to sample luminosity. It cannot be constrained as well as , and thus we decide to fix it to a common choice in the literature that is also supported by simulations, (Kravtsov et al., 2004; Zentner et al., 2005; Tinker et al., 2005; Zheng et al., 2005; Zehavi et al., 2011; Wake et al., 2011; Leauthaud et al., 2012). None of the final conclusions in this work change whether or not we allow to vary freely. Additionally, Equation 40 fixes through the observed galaxy number density, leaving as the only parameter left in the fit. In Appendix F we comment on how the resulting HOD model changes if we leave nearly all parameters free in the fit.
Obtaining the best-fit value of is straightforward. The error in the fit can be estimated from the width of the likelihood distribution, but it does not account for departures arising from cosmic variance. To account for that, we perform a set of 100 random realizations of the redshift distribution and number density at , which we call and . We find the best HOD fit each time, along with the corresponding derived parameters. Each redshift -bin in is drawn from a normal distribution with mean and standard deviation as in Section 3.1 (see Figure 5). The value for is produced in a similar manner; using a normal distribution with mean and standard deviation equal to the value and error of in Table 1. The scatter in all parameters from the random realizations is clearly dominant over that arising from the width of the likelihood in the fiducial fit, especially due to the variations in . We can therefore approximate the final errors as those from the random realizations.
Figure 9 shows the observed ACF and the model fits for a few samples. The values and errors for all relevant parameters are given in Table 2.

## 7 The stellar-to-halo mass ratio

Haloes of masses equal to host on average 0.5 central galaxies with luminosities greater than the sample threshold (Equation 17). Zheng et al. (2007) showed analytically that central galaxies living in these particular haloes have a median luminosity corresponding to the limit of the sample. This links halo masses with luminosities, albeit with some scatter dex (Zehavi et al., 2011; Coupon et al., 2012). However, we are interested in the connection with stellar masses, which also have a well-defined mean and scatter at fixed luminosity (see Section 3.3). Hence, we can link to , with a scatter () that ought to be close to the quadratic sum of the scatters from the luminosity- and luminosity- relations. We measure the latter to be around 0.2 dex, and the former is expected to be similar. Thus, the fact that we fix might seem an underestimation. However, as we will discuss in Appendix F, an unconstrained HOD fit does not prefer larger values for this parameter. Also, the final results do not change significantly by increasing it to larger values as 0.4 dex. We thus retain our choice and proceed.
The values of and for our samples can be found in Table 1. Their ratio yields the SHMR, which is plotted as a function of halo mass in Figure 10 for our different sets of stellar masses. The vertical error bars are a combination of the halo mass uncertainty and the error in the median stellar masses (equation 10), i.e., it does not represent the scatter in stellar mass at fixed halo mass. It is interesting that the error bars do not get notably bigger for brighter samples, even though the ACFs of those are much noisier and the stellar mass errors are indeed larger. The reason is that becomes progressively less sensitive to the HOD fit at higher luminosities. The fit is based on , which falls close to the steep drop of the halo mass function (equation 34) in the bright samples and makes the overall HOD model be weakly affected by the satellite occupation (e.g., Equation 40). Thus, the error contribution from is minimized, and that from and increases, keeping the total error roughly constant across the different samples.

### 7.1 Comparison to other results at z=1.5

There are several studies that have tried to constrain the SHMR at , based on abundance matching (Moster et al., 2013; Behroozi et al., 2013), HOD modeling (Zheng et al., 2007; Wake et al., 2011; Coupon et al., 2012) and extensions using conditional luminosity functions (Yang et al., 2012; Wang et al., 2013). Some of these works also provide their own parametric form for the SHMR as a function of halo mass, and we will use three of them to fit our points. These are the forms from Yang et al. (2012), Moster et al. (2013) and Behroozi et al. (2013) (hereafter Y12, M13, B13), which read:

 SY(m)=MY0(mMYp)αY+βY(1+mMYp)−βY, (25)
 SM(m)=2NM⎡⎢⎣(mMMp)−βM+(mMMp)γM⎤⎥⎦−1 (26)

and

 logSB(m)=log(ϵBMBp)+f(logmMBp)−f(0)−log(m) (27)

with

 Extra open brace or missing close brace (28)

The superscript labels refer to the author names. The position of the peak is mostly modulated by the pivot mass . In Y12 and M13, the low mass logarithmic slopes are set by , and the high-mass slopes by , , respectively. In the case of B13, the link between the slopes and the parameters is less straightforward, but the low and high-mass regimes are mostly modulated by and . tunes the high-mass behavior of going from logarithmic at to power-law at (see Behroozi et al. 2013). We set all of these parameters free when perform orthogonal regression fits of to our measurements of the SHMR. However, we do enforce and (see Yang et al. 2012), which are limits by definition.
These authors mainly use measures of the stellar mass functions at different redshifts to build a redshift evolution model of the SHMR. They provide explicit redshift dependence for all parameters in Equations 26-27. Thus, we use them to compare our measurements to the predicted SHMR of these authors at the redshift of our survey. This is shown in Figure 10, where we plot their predictions at low and high-redshift. In addition, the specific parameter values of the curves, for both the predictions and the fits to our data, are displayed in Table 2. The normalization values , and are also fitted, although we disregard any interpretation of them because there are important systematic uncertainties in the stellar masses between different authors. A thorough examination of these to allow a meaningful comparison is beyond the scope of this paper. For the current purposes, we simply assume that the differences in stellar masses are due to a simple logarithmic offset. This assumption holds well when comparing different sets of masses in the COSMOS and EGS catalogs. In addition, M13 and B13 use stellar masses based on BC03 and Chabrier IMF, which matches our fiducial choice of masses. Y12 use masses produced with the Fioc & Rocca-Volmerange (1997) models and Kroupa IMF, but we still do not expect a significant deviation from a constant offset when compared to our masses (Barro et al., 2011b). We have checked this based on the masses from this particular model that are also available in the EGS control catalog. We also note that we use the SHMR in Y12 that is based on fits “CSMF/SMF1”, where only stellar mass functions are utilized.
We limit the comparison between all measurements to the centroid position and slopes of the SHMR. There are some discrepancies when comparing our results to the predictions from the other authors, but these are not dramatic (see below). The centroid of the SHMR is computed as the actual peak position in the parametric relations, and the fits of these models to our data show (see Table 2). We did not compute confidence intervals for the predictions because it requires knowledge of the explicit covariance between their parameters fits. Our value is larger than these predictions; based only on our errors, it lies 0.8 above Y12, 1.2 above M13 and 2.7 above B13. Because their errors are not being taken into account, these offsets should not be treated as absolute levels of inconsistency with respect to our study.
For the slopes, there are also some slight discrepancies. To better visualize this comparison, we have plotted in Figure 11 the prediction and the fits to our data for each parametric model. All curves in each panel are scaled in the x-axis to match the peak of the prediction, and scaled in the y-axis to set all peak heights to zero. The idea is to fix the peak position (in both axes) of all curves to better compare the slopes on either side. In this case, the slopes are the approximate power-law index at either side of the peak, and is not necessarily linked to a parameter in a unique manner (except for M13, where the slopes are independently controlled by , ). In comparison to M13, our low mass slopes are steeper (higher ) and the high-mass slopes are shallower (lower ) than their predictions. With respect to B13, the low mass slopes are in agreement but our high-mass slopes are shallower. In the case of Y12, our slopes are steeper at both low and high-mass.
As explained in Section 1, the low-mass slope can be directly related to the importance of energy- versus momentum-driven winds. In general, we find a steeper low-mass slope than the predictions, which favors energy-driven winds. At high-masses, the interpretation of the slope is less clear, since AGN feedback and galaxy mergers should also have an important contribution.
Another important study to compare our measurements with is Wake et al. (2011). It is based on HOD modeling of stellar mass limited galaxies at , which makes it similar to our work. These authors had the advantage of using data with accurate photometric redshifts and stellar masses, but also the drawback of sampling a small region of the sky (NEWFIRM survey, 0.25 deg). They had very few galaxies around and therefore it was not possible to map the full peak of the SHMR. Their data is shown in Figure 10, where we have scaled the stellar masses by 50% to roughly transform them from the M05 model to BC03. These authors performed a parametric fit and found a peak at (an estimated uncertainty was not provided), which lies 2.5 above our result of .

### 7.2 Evolution with redshift

At this point, we can compare our result for the peak in the SHMR with other studies at different epochs and trace its evolution with redshift (Figure 12). We include HOD results from Zehavi et al. (2011), Zheng et al. (2007), Leauthaud et al. (2012), Coupon et al. (2012) and Wake et al. (2011), as well as predictions from M13, B13 and Y12. As mentioned earlier, our peak lies above the predictions and below the value inferred by Wake et al. (2011). Looking at the trend with values at other epochs, the peak mass seems to have evolved in a monotonic and quasi-linear way with redshift. Our data supports a change of log through . This means that the halo mass scale that is most efficient at forming and accreting stars to the central galaxy has decreased by a factor of 4.5 during this redshift range. Thus, the downsizing trend of galaxies has continued steadily during the last 10 Gyrs. Low mass galaxies have grown faster than their haloes, while the opposite trend happened for high-mass galaxies.

## 8 Satellite galaxies

### 8.1 Satellite fraction

The satellite fractions for all of our samples are displayed on Table 1 and plotted in the lower panel of Figure 13. Note that the satellites making up this fraction are above the sample flux limit, i.e., does not refer to the total fraction of satellites that a central galaxy at the flux limit has. The satellite fraction clearly decreases towards the brighter end, which is a manifestation of the drop in the halo mass function. Based upon the model we use, is the scale that sets the occupation number of satellites in a halo of a given mass, and the number of such haloes is given by the mass function. If approaches the cutoff scale of the mass function, then the satellite contribution to the total density will be reduced compared to that of central galaxies. This effect is seen in most studies (Zheng et al., 2007; Zehavi et al., 2011; Wake et al., 2011; Coupon et al., 2012; Tinker et al., 2013).
Our faint-limit value is (see Figure 13). Compared to the results of obtained at by Zehavi et al. (2011), our value is suggestive of a mild increase in the satellite fraction with cosmic time. A similar conclusion was also reached by Coupon et al. (2012) based on their comprehensive study of samples at , and such evolution is predicted by some simulations (e.g., Wetzel et al. 2009, 2013). We caution, however, that the sample in (Zehavi et al., 2011) extends to fainter absolute magnitudes that our data set ( versus ), and so the evidence from this comparison is suggestive rather than conclusive.

### 8.2 The M1/Mmin relation

A deeper insight into the relationship between haloes and their satellites is given by the ratio. As mentioned in previous Sections, haloes typically become occupied by a central galaxy at and gain an additional satellite at . Thus, at fixed , lowering would directly increase the overall satellite fraction. However, holds further clues in relation to the galaxies that occupy these haloes. Because of the decline in the halo mass function towards the massive end, most of the haloes are small and have masses around . These will typically host galaxies that are also small, with stellar mass close to the sample limit . The satellites considered have masses that are also near this limit and living in haloes near , where the central galaxy can have a mass much larger than . However, if approaches , then its central galaxy will have a mass closer to . In the case of , the satellite will have a stellar mass around and the central will be slightly more massive than that. Thus, when this ratio is smaller, there is an increased fraction of centrals that have a satellite of similar stellar mass.
In the local Universe, (Zehavi et al., 2011; Beutler et al., 2013). On the other hand, we measure at . A decline of this ratio with redshift had been predicted by simulations (Kravtsov et al., 2004; Zentner et al., 2005) and measured by abundance matching (Conroy et al., 2006) and other HOD studies. Ratios at different redshifts and fixed number density are shown in the top panel of Figure 12, where it can be seen that there is a general increase towards later times. The reason why both and are higher at low-redshift is due to the evolution of the mass function. At low-redshift, there are very large haloes that can host multiple satellites, which helps increase the average satellite fraction. However, the fraction of galaxies that are satellites with masses close to the limit of the sample is still larger at high-redshift. This happens because the halo infall timescale is lower, which enhances their accretion rate onto other structures and reduces the gap between and (see below and Conroy et al. 2006).
However, the most interesting result we derive is the trend of with sample luminosity. The middle panel of Figure 13 shows indications of a slight rise at high luminosities, which is not obviously expected. At , this ratio has been observed to be constant or decrease with increasing luminosity at fixed redshift (Zheng et al., 2007; Blake et al., 2008; Abbas et al., 2010; Zehavi et al., 2011; Matsuoka et al., 2011; Leauthaud et al., 2012; Beutler et al., 2013). Simulations also predict that the accretion rate is larger for more massive haloes at all times (Zentner et al., 2005; Wetzel et al., 2009; McBride et al., 2009; Fakhouri & Ma, 2008; Fakhouri et al., 2010), which would lower at the bright end. We measure the opposite trend at . Interestingly, Wake et al. (2011), the only other study at this redshift that measured HOD for stellar mass selected samples, also obtained a slight increase of this ratio with sample mass. However, those authors did not explore this effect in depth. The results from Coupon et al. (2012) also hint a similar trend at in haloes of mass , although they are consistent with a constant ratio.
In order to better compare the results from a few different authors, we plot as a function of cumulative number density in Figure 14. For visual clarity, we show in the left (right) panel those results that follow an increasing (decreasing) trend with density, along with our data. A caveat in this comparison is that, in reality, the number density of a given population does not remain constant through redshift. However, the data from Zheng et al. (2007) and Coupon et al. (2012) do not follow the same trends, even though they sample similar redshifts. Thus, from the observational side, the data do not offer a consensus regarding the trend with luminosity. At , our results and those from Wake et al. (2011) do support a minimal rise in with luminosity.
As shown in Figure 14, two families of curves can be defined. Our results, Wake et al. (2011) and Coupon et al. (2012) show a similar shape, offset in the y-axis according to redshift. Meanwhile, Zehavi et al. (2011) and Zheng et al. (2007) are similar to one another. One possible effect leading to the disparity between the two sets of results may be the particular selection of galaxy samples. Those in Zehavi et al. (2011), Zheng et al. (2007) and Coupon et al. (2012) are limited by absolute magnitude in the optical. Wake et al. (2011) selects directly in stellar mass, and we make a luminosity selection that is later matched to a stellar mass limited sample. Thus, we find no clear explanation for the existence of these two families of curves regarding sample selection. In addition, all these authors (including us) use a very similar form of the HOD, and variations in the chosen cosmology do not have such a strong impact. Regarding possible systematics effects in our modeling, we test different possibilities in Appendices F and E and find nothing that would alter our conclusions.

### 8.3 Physical mechanisms for a mass-dependent evolution

The rise of with luminosity is not clearly detected. However, Zehavi et al. (2011) and Zheng et al. (2007) very clearly measure the opposite behavior at redshifts and , respectively, so that even if our data follows a flat trend at , it would imply that evolution has taken place. Interestingly, there is no obvious mechanism that could be responsible for this change, and we speculate with some possibilities in what follows. The dynamical processes at play can be reduced to a competition between accretion and destruction of satellites. Regarding the former, big structures have recently assembled a larger fraction of their mass than smaller counterparts, at all times (Wechsler et al., 2002; Zentner et al., 2005; Fakhouri et al., 2010). In other words, the specific growth rate of haloes is an increasing function of mass. Regarding the latter, the dynamical time in bigger haloes is larger, contributing to a slower destruction of accreted satellites. These effects yield a larger number of recently accreted and undisrupted satellites in larger haloes, which would produce a decrease in towards higher masses.
So, what additional mechanism can reverse this trend at high-redshift? This mechanism could involve the ratio of destruction to accretion being larger at high-masses, which is possible if the dynamical timescale decreases considerably with mass. However, a caveat in these scenarios is that we are implicitly considering that galaxies are accreted or disrupted in the same way as haloes, which does not have to hold. What we are really tracking are galaxies, since is inversely proportional to the occurrence of galaxy pairs with masses close to the sample limit. Thus, there could be a star formation dependent process that drives the trend we see with stellar mass. For example, Wetzel et al. (2013) show that the star formation in satellites fades at the same rate as the central galaxy for a few Gyrs after accretion, but then undergoes a rapid quenching period. They find that quenching timescale is shorter for more massive satellites. Thus, if the central galaxy outgrows the satellites in a way proportional to its own mass, this would produce a lower fraction of similar mass pairs and play in favor of our trend. In addition, such a mechanism would need to become milder toward low-redshift, so that the trend becomes inverted. This allows us to restate our previous question: what physical process would more efficiently quench satellites in similar mass pairs and is more important at high redshift? We do not have a plausible answer for this question.

## 9 Galaxy bias

At small scales, the complex baryonic processes of galaxy formation break the homology between the spatial distribution of galaxies and dark matter. However, at large-scales, the gravitational effects of dark matter dominate the dynamics and the overdensity of some selection of galaxies is expected to match that of dark matter multiplied by a scaling factor, the galaxy bias . In the HOD models, this quantity is described as a number-weighted average of the halo bias (see Equation 20) and ideally would match the square-root of the ratio between the large-scale SCF of galaxies and dark matter (equation 1).
Our measurements of the effective galaxy bias are shown in the top panel of Figure 13, where we plot against apparent magnitude threshold. Bright (massive) galaxies have a larger bias than faint (small) ones, a trend that has been determined in many other studies (Benoist et al., 1996; Norberg et al., 2001; Tegmark et al., 2004; Zehavi et al., 2005; Brodwin et al., 2008; Brown et al., 2008; Foucaud et al., 2010; Zehavi et al., 2011; Matsuoka et al., 2011; Coupon et al., 2012; Beutler et al., 2013; Jullo et al., 2012; Mostek et al., 2013) and is expected because luminous galaxies reside on average in more massive haloes, which are more biased with respect to dark matter (White et al., 1987; Kauffmann et al., 1997). In addition to , we also fit the large-scale bias directly to our measured ACF at deg ( comoving Mpc at ). This fit does not depend on the HOD, and is performed by scaling the dark matter SCF in a similar way to the procedure in Section 5.2, but leaving the bias as a free parameter. The inputs from the galaxy population are the redshift distributions and the evolution in the median mass to modulate the bias across redshifts, but not the galaxy number density. The shaded regions in the top panel of Figure 13 represent the confidence intervals for the direct bias fit, which is in good agreement with the HOD values. Thus, we find that the HOD modeling of our data makes a good description of the large-scale bias.
Nonetheless, we note that this description is not perfect. In Appendix F we comment on how a HOD fit with all parameters allowed to vary freely makes float down to unphysical values , trying to maintain a high bias that otherwise would yield a smaller value due small shifts in the fitted and . The overall HOD is not very sensitive to and therefore this is not a significant problem. However, it points towards the amplitude of our observed ACFs being slightly too large to be perfectly reproduced by the combination of halo bias and halo mass functions.

### 9.1 Comparison to Wake et al. (2011)

We find a slight bias excess in our data, a result that has been noted to a larger extent in other HOD studies of stellar mass limited samples at . Matsuoka et al. (2011) and Wake et al. (2011) find that their ACFs are too strong to be reproduced by a halo model with the observed density of galaxies. Those fits to the clustering plus density were compared to fits to the clustering only, where the number density was not fixed to the observed value. The latter fit was able to reproduce the ACFs, but with a bias about 50% higher than the clustering plus density fit in their most massive and distant samples. The samples of Wake et al. (2011) are directly comparable to our study, since they are defined by lower stellar mass limits. In the lower panel of Figure 15 we show our standard HOD bias measurements as a function of stellar mass and the bias results from those authors. To make the comparison more direct, we plot our results for the Maraston (2005; hereafter M05) evolutionary models with the Kroupa (2001) IMF. Here we comment on the two types of HOD fits in Wake et al. (2011), and how they compare to our results:

• Clustering only: Wake et al. (2011) find the effective bias from this fit to be the closest to a direct measurement of the large-scale bias. However, these values are high compared to our findings. The 0.2 dex offset relative to our work could be due to a difference in stellar mass estimates (Figure 15). Given that both studies employ very similar stellar mass models (M05 stellar grids, Kroupa IMF, and Calzetti 2000 extinction), this possibility seems unlikely. Another explanation would be sample variance due to the small size of the survey in Wake et al. (2011), which could lead to an excess in the clustering signal. Our survey is almost 400 times larger in area and therefore significantly less impacted by this effect.

• Clustering + density: Wake et al. (2011) find biases from this fit that are also not fully consistent with ours. Their bias is larger (smaller) than our values in the low (high) mass end. However, the observed number densities are discrepant in the opposite way, and the stellar mass at which the densities and biases match is roughly the same, . In any HOD model, the effective bias is anticorrelated with number density if the rest of the parameters are held fixed. Thus, if Wake et al. (2011) and our study had the same observed densities, the HOD bias from both surveys could perhaps be in full agreement.

Thus, we speculate that the high clustering amplitude in Wake et al. (2011) might be predominantly a consequence of cosmic variance (as also suggested by those authors), and their clustering + density HOD fits would be consistent with ours if the observed comoving number densities were the same.

### 9.2 Comparison to other studies

We also compare our bias results with other measurements at different redshifts based on stellar masses in the top panel of Figure 15. These studies include M10, Foucaud et al. (