THE CONCENTRATION DEPENDENCE OF THE GALAXY–HALO CONNECTION: MODELING ASSEMBLY BIAS WITH ABUNDANCE MATCHING

This study uses several cosmological boxes, as listed in Table 1. The c250-2048 box comes from the "Chinchilla" suite, run with the L-Gadget code, a variant of Gadget (Springel 2005). The "Chinchilla" suite spans a wide range of box sizes and resolutions with the same cosmology (M. R. Becker et al. 2016, in preparation; some boxes in this series were described in Mao et al. 2015 and were also used in Desmond & Wechsler 2015). Bolshoi and BolshoiP have the same box size and resolution as c250-2048, but have different cosmologies and were run with the ART N-body code (Klypin et al. 2011). MDPL is part of the "MultiDark" suite (Klypin et al. 2016), and was also run with Gadget. The three DarkSky boxes of different sizes are smaller boxes that accompany the $8\,\mathrm{Gpc}\,{h}^{-1}$ box from the "Dark Sky" Simulations (Skillman et al. 2014), run with the 2HOT code (Warren 2013). Here we refer to these boxes as DarkSky-250 (ds14_j_2560), DarkSky-400 (ds14_i_4096), and DarkSky-Gpc (ds14_b). The particles used to build the halo catalogs and merger trees for the DarkSky-Gpc box were down-sampled (1/32 particles) from a high-resolution box run with 10240³ particles.

Table 1.  Cosmological and Simulation Parameters for Boxes Used in This Study

Box Name	Side Length	Particle	h	${{\rm{\Omega }}}_{M}$	ns	${\sigma }_{8}$	Particle Mass	Code
	$(\mathrm{Mpc}\,{h}^{-1})$	number					$({M}_{\odot }\,{h}^{-1})$
c250-2048	250	20483	0.7	0.286	0.96	0.82	1.44 × 108	L-Gadget
Bolshoi	250	20483	0.7	0.27	0.95	0.82	1.35 × 108	Art
BolshoiP	250	20483	0.678	0.295	0.968	0.823	1.49 × 108	Art
MDPL	1000	38403	0.678	0.307	0.96	0.823	1.51 × 109	L-Gadget
DarkSky-250	250	25603	0.688	0.295	0.968	0.834	7.63 × 107	2Hot
DarkSky-400	400	40963	0.688	0.295	0.968	0.834	7.63 × 107	2Hot
DarkSky-Gpc	1000	102403a	0.688	0.295	0.968	0.834	4.88 × 109b	2Hot

Notes.

^aHalo catalogs and merger trees are constructed with 1/32 of the total particle number. ^bEffective mass of the down-sampled particles.

Download table as:  ASCIITypeset image

For each of these boxes, we use the halo catalog generated by the Rockstar halo finder (Behroozi et al. 2013a) and the Consistent Trees merger tree builder (Behroozi et al. 2013b). We use the virial overdensity (${{\rm{\Delta }}}_{\mathrm{vir}}$) as our halo mass definition (Bryan & Norman 1998).

The mock galaxy catalogs used in this work are generated with the abundance matching technique. We follow the procedure of Behroozi et al. (2010) and Reddick et al. (2013) in order to implement abundance matching with scatter in luminosity at fixed halo proxy. First, we deconvolve the scatter from the luminosity function. We then abundance match luminosity with the halo proxy, producing a catalog of galaxy luminosities. Finally, we replace the scatter by adding a log-normal scatter to the catalog.

We make measurements of the projected two-point correlation function, ${w}_{p}({r}_{p})$, from the mock catalogs as follows. We use the plane-parallel approximation in redshift-space and integrate along one of the axes (i.e., the line of sight), with an integration depth of $2{z}_{\max }=80\,\mathrm{Mpc}\,{h}^{-1}$. Redshift-space distortions are applied along the integration axis before integration. We account for the periodic boundary conditions of the cosmological boxes when computing the projected correlation function, and hence we can evaluate the expected number of uniformly random pairs explicitly without using an estimator.

In this study, we use the luminosity function (for abundance matching) and the two-point clustering measurements (for comparison) extracted by Reddick et al. (2013). These measurements were made on the volume-limited samples from the New York University Value Added Galaxy Catalog (NYU-VAGC; Blanton et al. 2005), based on Data Release 7 from the Sloan Digital Sky Survey (Padmanabhan et al. 2008; Abazajian et al. 2009). We note that these measurements are quite consistent with the measurements of Zehavi et al. (2011), but here a consistent sample was used to determine both the luminosity function and clustering measurements. We refer the readers to Section 2 and Appendix C of Reddick et al. (2013) for details on these measurements. In this work, we focus primarily on constraining our models with galaxies of luminosity $\sim {L}_{* }$ and brighter in order to be conservative about the resolution requirements for the complete halo and subhalo samples needed for abundance matching, but we present results from dimmer samples in Section 4.4.

In Section 4, when we compare the SDSS data to the predicted ${w}_{p}({r}_{p})$ obtained from the mock catalogs with a ${\chi }^{2}$ statistic, we need to take the covariance matrix into account. The covariance matrix has three contributions: (1) the sample variance (jackknife covariance) of the SDSS data, (2) the sample variance (jackknife covariance) of the mock catalog, and (3) the covariance due to the stochastic scatter in the mock catalog. We sum up these three contributions quadratically to obtain the final covariance matrix.

First, the sample variance in the SDSS measurements is estimated by jackknifing the SDSS data set, as detailed in Reddick et al. (2013). Second, we estimate the sample variance in the predictions of abundance matching due to the finite volume of the N-body simulations used in this work. We employ a jackknife procedure in order to estimate the contribution to the covariance matrix from this effect. Each N-body box is divided into smaller square "sky areas." Each patch has side lengths of $25\,\mathrm{Mpc}\,{h}^{-1}$, $25\,\mathrm{Mpc}\,{h}^{-1}$, and the original box side length. We then omit one patch at a time in the jackknifing process (i.e., omitting everything along the line of sight in the square patches), and compute the jackknife covariance. Third, the final contribution to the covariance comes from the scatter in abundance matching. Since we apply log-normal random scatter in luminosity directly to the catalogs, multiple catalogs generated with the same abundance matching parameters produce slightly different predictions for ${w}_{p}({r}_{p})$. Thus, from each set of abundance matching parameters, we generate 40 catalogs, compute ${w}_{p}({r}_{p})$ for each, and calculate the covariance on ${w}_{p}({r}_{p})$ due to this random variation.

We note that the estimate of the covariance of the mock ${w}_{p}({r}_{p})$ has a direct impact on the goodness of fit, and hence on the derived constraints on the abundance matching parameters. Nevertheless, Norberg et al. (2009) find that the jackknife method does not typically underestimate the covariance.

A key simplifying assumption in the abundance matching framework is that some galaxy property (for both centrals and satellites) is tightly correlated with some dark matter (sub)halo property. Physically, this halo property should be connected to, for example, the mass, velocity, or potential of the dark matter halo. However, at present, we do not understand the details of galaxy formation well enough to specify the exact connection. Hence, in practice, the halo proxy used in abundance matching could be a function of multiple halo properties instead of a single halo property.

In particular, to test whether, or how strongly, galaxy properties within halos at a fixed mass should be dependent on halo concentration or assembly history, here we present an interpolation scheme that generalizes the conventional abundance matching model to allow continuously adjustable concentration dependence. To build such a scheme, we adopt the parameterization used in Mao et al. (2015), defining a new generalized proxy to be used in abundance matching:

$$\begin{eqnarray}&&{v}_{\alpha }:= {v}_{\mathrm{vir}}{\left(\displaystyle \frac{{v}_{\max }}{{v}_{\mathrm{vir}}}\right)}^{\alpha },\end{eqnarray} \tag{ 1 }$$

where ${v}_{\max }$ is the maximal circular velocity and

$$\begin{eqnarray}&&{v}_{\mathrm{vir}}:= {\left(\displaystyle \frac{{{GM}}_{\mathrm{vir}}}{{R}_{\mathrm{vir}}}\right)}^{1/2}={\left(\displaystyle \frac{4\pi }{3}{{\rm{\Delta }}}_{\mathrm{vir}}{\rho }_{\mathrm{crit}}{G}^{3}\right)}^{1/6}{M}_{\mathrm{vir}}^{1/3},\end{eqnarray} \tag{ 2 }$$

with ${{\rm{\Delta }}}_{\mathrm{vir}}$ being the virial overdensity and ${\rho }_{\mathrm{crit}}$ the critical density.

This generalized proxy captures the continuously varying dependence on concentration through the parameter α because the ratio ${v}_{\max }/{v}_{\mathrm{vir}}$ can be viewed as a proxy for halo concentration. In principle, this ratio can be replaced by f(c) with a general function f. Nevertheless, using this ratio facilitates comparisons with other proxies that have been used in the literature. In particular, when $\alpha =0$, the dependence on concentration is turned off, as ${v}_{\alpha =0}={v}_{\mathrm{vir}}\propto {M}_{\mathrm{vir}}^{1/3}$, and when $\alpha =1$, this proxy reduces to the maximal circular velocity ${v}_{\alpha =1}={v}_{\max }$.

Note that Equation (1) can be written as

$$\begin{eqnarray}&&\mathrm{log}{v}_{\alpha }=\alpha \mathrm{log}{v}_{\max }+(1-\alpha )\mathrm{log}{v}_{\mathrm{vir}}.\end{eqnarray} \tag{ 3 }$$

Hence, this new proxy can be considered as a linear combination of $\mathrm{log}{v}_{\max }$ and $\mathrm{log}{v}_{\mathrm{vir}}$ on a logarithmic scale. As illustrated in Figure 1, on a log–log plot of ${v}_{\max }$ and ${v}_{\mathrm{vir}}$, the value of α affects the direction of the halo ranking. In this figure, the arrow represents the direction of descending rank when the halos are ranked by v_α, and the slope of the arrow is $\alpha /(1-\alpha )$. Hence, different values of α effectively rank the halos with different slopes. As a result, at a given number density, different values of α select out different halos. In particular, a larger value of α selects out more low-mass, high-concentration halos, and also more subhalos.

Zoom In Zoom Out Reset image size

We note that the specific choice of the parameterization of the concentration dependence should not impact our results significantly, as the essence of our model is to vary how much we weight the concentration of halos when ranking halos by their masses in the abundance matching procedure. However, one could instead weight other halo properties, such as halo formation time, in order to study the dependence on other properties in abundance matching. In this work, we only study the dependence on concentration. Nevertheless, we expect that qualitatively similar results would also apply to other proxies that are highly correlated with concentration.

So far, we have only discussed how to model the concentration dependence in our new proxy. In abundance matching, the choice of epoch at which the ranking proxy is evaluated also significantly impacts the results (Reddick et al. 2013; Chaves-Montero et al. 2016). For example, if the proxy is evaluated at the present day, subhalos are usually ranked lower due to stripping, and the resulting mock catalog is less clustered. Conroy et al. (2006) argues that the time at which a subhalo enters the virial radius of its parent halo is a natural time at which to set proxies. Reddick et al. (2013) further uses the peak values of those proxies (e.g., ${M}_{\mathrm{peak}}$ and ${v}_{\mathrm{peak}}$) throughout history. In this work, we limit our discussion to a single choice of epoch. We evaluate the value of v_α for each halo at the epoch when ${M}_{\mathrm{vir}}$ reaches its peak value, and let ${\hat{v}}_{\alpha }$ denote this quantity. In follow-up work, we will explore this choice of proxy epoch in detail.

Since ${\hat{v}}_{\alpha }$ is evaluated at the time of peak mass for each halo, ranking with ${\hat{v}}_{\alpha =0}$ and ${\hat{v}}_{\alpha =1}$ is equivalent to ranking with ${M}_{\mathrm{peak}}$ and ${v}_{\max }$ at ${M}_{\mathrm{peak}}$ respectively. (However, the former is only approximately true in our case because ${{\rm{\Delta }}}_{\mathrm{vir}}$ in Equation (2) has a weak dependence on the scale factor, and for different halos, ${M}_{\mathrm{peak}}$ occurs at different scale factors.) Our choice of evaluating the abundance matching proxy at the scale when ${M}_{\mathrm{peak}}$ rather than when ${v}_{\mathrm{peak}}$ occurs was motivated by the finding that halos at the largest circular velocities may be out of dynamical equilibrium (Ludlow et al. 2012); e.g., Behroozi et al. (2014) showed that ${v}_{\mathrm{peak}}$ is commonly set by major mergers, and hence may not represent the physical time at which the subhalo started to be stripped. Evaluating the proxy at the scale factor of ${M}_{\mathrm{peak}}$ then avoids this unphysical epoch probed by ${v}_{\mathrm{peak}}$, and is similar to using the relaxation criterion proposed by Chaves-Montero et al. (2016). Nevertheless, for the purpose of abundance matching, the difference between matching to ${v}_{\mathrm{peak}}$ and to ${v}_{\max }({a}_{\mathrm{mpeak}})$ is minimal, as the rank orders are very similar when halos are ranked by these two proxies. As a result, the clustering signals with these two proxies are also similar.

We first demonstrate the impact of α, as defined in Section 3.1, on clustering. Figure 2 shows the wide range of clustering predictions that can be produced by varying α. We find that changing α can significantly change the clustering, and that a higher value of α produces a more clustered catalog.

Zoom In Zoom Out Reset image size

There are two effects that contribute to this result. First, at a given halo mass, on average, subhalos have higher concentrations than host halos. Hence, when a higher value of α is used, subhalos are more likely to make it through the threshold cut, resulting in a more clustered sample. Effectively, increasing α increases the difference between the luminosity–halo mass relation of host halos and that of subhalos. This effect impacts both the one- and two-halo terms, and also boosts the satellite fraction.

Second, when we use a higher value of α, high-concentration halos are ranked higher in the catalog and are more likely to make it through the threshold cut. Since, in this mass regime, high-concentration halos are more clustered due to halo assembly bias (Wechsler et al. 2006), the resulting catalog is also more clustered. This effect impacts mostly the two-halo term, and is less significant in brighter samples. In the sections below, we discuss these two effects in detail.

It is known that increasing the scatter in abundance matching would decrease clustering strength because it brings in lower-mass halos (e.g., Reddick et al. 2013). Thus, there exists a degeneracy between α and the scatter. This degeneracy is demonstrated in Figure 3, which shows the correlation function for several values of α and scatter. The clustering strength decreases with decreasing α and also with increasing scatter. Nevertheless, the scatter has a stronger effect on the brighter samples, while α has a stronger effect on the fainter samples. This implies that samples of different thresholds are likely to give different constraints on α and scatter, and might be able to break the degeneracy between α and scatter.

Zoom In Zoom Out Reset image size

Here, we define the satellite fraction to be the fraction of satellites in bins of luminosity. In this study, we did not apply a group finder on the mock galaxies, so galaxies labeled as satellites are exactly the same as those labeled as subhalos in the initial halo catalog. That is, the satellite fraction we measured here is actually the fraction of galaxies assigned to subhalos in each luminosity bin. A subhalo is defined as any halo whose center falls within another larger halo. We refer to a halo that is not a subhalo as a host halo.

Figure 4 shows the satellite fraction as a function of luminosity for several values of α. As expected, increasing α increases the satellite fraction, since subhalos are, on average, more concentrated than host halos of the same mass; hence, subhalos are ranked higher when α is larger. This is especially true at the faint end because the ratio ${v}_{\max }/{v}_{\mathrm{vir}}$ differs more between subhalos and host halos for low-mass halos.

Zoom In Zoom Out Reset image size

Applying scatter to abundance matching increases the satellite fraction at the bright end because more satellites in the fainter luminosity bins are scattered up to the brighter luminosity bins. Applying scatter does not significantly change the satellite fraction for samples fainter than ${M}_{r}=-21$.

Another way to demonstrate this change in the satellite fraction is to look at the difference between the luminosity–halo mass relation of host halos (central galaxies) and that of subhalos (satellite galaxies). Previous studies have explored the case in which the stellar mass–halo mass relations of central and satellite galaxies differ from each other (e.g., Neistein et al. 2011; Rodríguez-Puebla et al. 2012, 2013). Here, using the α parameter, we can evaluate this difference quantitatively. Figure 5 shows the luminosity–halo mass (L − M_h) relations for all halos, only host halos, and only subhalos, for different values of α. We see that changing α changes the overall L–M_h relations very little, but changes the difference between the halo and subhalo L–M_h relation. Specifically, increasing α effectively more strongly differentiates the L–M_h relations for halos and subhalos, while maintaining the overall L–M_h relation.

Zoom In Zoom Out Reset image size

In Figure 4, we also compare our results with the observed satellite fraction. The observed satellite fraction measurements are taken from Reddick et al. (2013), who used a group finder (Tinker et al. 2011) applied to the same sample used to make the clustering measurements. Since we did not apply the same group finding procedure on our mock catalogs, this comparison is subject to the systematic errors introduced by the group finder. The gray band shown in Figure 4 is the sum of the error from SDSS data and the estimated systematic error introduced by the group finder; the latter was estimated by taking the one-sided difference between the satellite fractions before and after the catalog was processed with a group finder, shown in the left panel of Figure 21 in Reddick et al. (2013). We see that these systematic errors increase significantly at the bright end, due to the fact that the group finder does not always select the right galaxy as the central. However, both scatter and group finding have much smaller impacts at luminosities dimmer than ${M}_{r}=-21$, which is also where α has a larger impact. Up to the systematic errors, the observed satellite fraction agrees well with the model prediction when $\alpha \sim 0.5$. We show in Section 4.2 that this is also consistent with the inference from galaxy clustering.

In our model, α also controls how much halo assembly bias can manifest in the mock catalogs as galaxy assembly bias. To quantify this, we need to separate the effects of the satellite fraction and halo assembly bias. To that end, we shuffle our mock catalogs to remove halo assembly bias, but leave the satellite fraction intact. Here we adopt the same shuffling procedure as described by Zentner et al. (2014). We divide the catalogs into bins of halo masses, with a bin width of 0.1 dex. For each bin, we first shuffle the central galaxies, and then independently shuffle the satellites (while retaining their relative positions to the central galaxies). This procedure preserves the halo occupation and the satellite fraction in the catalogs by construction.

Figure 6 shows the relative difference in ${w}_{p}({r}_{p})$ between the shuffled and unshuffled catalogs. Since the shuffling procedure preserves the satellite fraction, the difference seen in this figure comes from halo assembly bias alone. We see that the difference is larger for fainter samples and for larger values of α. This behavior is expected: halo assembly bias impacts the fainter samples more strongly, and catalogs with a larger value of α have stronger halo assembly bias and are more clustered. We note here that although models with more concentration dependence have stronger assembly bias, there is still some assembly bias in the models with α = 0, because the relationship between ${M}_{\mathrm{vir}}$ and ${M}_{\mathrm{peak}}$ has some dependence on formation time and/or halo concentration.

Zoom In Zoom Out Reset image size

Figure 6 also shows the scale dependence of assembly bias for each value of α. We find that halo assembly bias impacts both the two-halo term and the transition regime around $1\mbox{--}2\,\mathrm{Mpc}\,{h}^{-1}$, in agreement with the findings of Sunayama et al. (2016). At the smallest scales, the original catalog and the shuffled catalog exhibit similar clustering, since the clustering at small scales is dominated by the change in satellite fraction. This implies that our ${\hat{v}}_{\alpha }$ model is distinct from merely introducing halo assembly bias to a non-biased catalog (e.g., modeling the HOD). In particular, varying α simultaneously changes the amount of halo assembly bias and the satellite fraction.

In the previous section, we show how this α parameter, which controls the concentration dependence in abundance matching, affects the clustering signals in the mock catalog. Given this finding, here we investigate whether the current galaxy clustering measurement can already provide constraints on the this α parameter. Since the effect of the α parameter on the clustering signals and that of the scatter in abundance matching are degenerate, here we present the joint constraints on α and scatter using the SDSS galaxy catalog.

We compute the ${\chi }^{2}$ statistic to evaluate the goodness-of-fit for a set of values in the $(\alpha ,\mathrm{scatter})$ parameter space for each threshold. We also compute the ${\chi }^{2}$ statistic for several different cosmological boxes to determine whether the constraint on $(\alpha ,\mathrm{scatter})$ varies significantly between boxes using different cosmologies and codes.

The ${\chi }^{2}$ statistic is computed as

$$\begin{eqnarray}&&{\chi }^{2}=\displaystyle \sum _{i}\displaystyle \sum _{j}d({r}_{p}^{i})d({r}_{p}^{j}){C}^{-1}({r}_{p}^{i},{r}_{p}^{j}),\end{eqnarray} \tag{ 4 }$$

where $d({r}_{p}):= {w}_{p}^{\mathrm{mock}}({r}_{p})-{w}_{p}^{\mathrm{SDSS}}({r}_{p})$, $C({r}_{p}^{i},{r}_{p}^{j})$ is the covariance, and $i,j$ denote the indices of the bins of r_p. Note that here C already includes the covariance from the SDSS data, as well as the covariance from jackknifing the mock catalog and from multiple realizations of abundance matching.

Note that in our analysis, we compute the p-value using the ${\chi }^{2}$ test independently for each set of parameter values. Hence the p-values presented in the figures below are not comparisons between different parameter values, but always are with respect to the null hypothesis.

Figure 7 shows the constraints from the three DarkSky boxes and the MDPL box for four different luminosity thresholds separately. Note that the two $1\,\mathrm{Gpc}\,{h}^{-1}$ boxes do not have the resolution to generate a complete sample below roughly ${M}_{r}=-21$, and hence we omit the lowest luminosity panels for these boxes in Figure 7. We will discuss detailed resolution requirements for abundance matching in upcoming work.

Zoom In Zoom Out Reset image size

We see several interesting features here. First, the degeneracy between α and scatter is most visible in the samples of ${M}_{r}\lt -21.5$ and −21. In both cases, we see the degeneracy as expected: a larger α requires a larger scatter to balance the additional clustering since more highly concentrated halos are included.

Second, as expected, larger boxes provide stronger constraints, indicating that the constraining power of most previous studies, which have almost exclusively used boxes of $\sim 250\,\mathrm{Mpc}\,{h}^{-1}$ on a side, have been dominated by sample variance. This is especially true for the brighter samples because the numbers of galaxies in those samples are small. While the sample of ${M}_{r}\lt -22$ from DarkSky-250 and -400 provide little constraint on α and scatter, the samples from the $1\,\mathrm{Gpc}\,{h}^{-1}$ boxes give clear constraints on scatter, and exclude zero scatter in this range of α at $p\lt 0.001$.

Third, on the faint end, we obtain a much stronger constraint on α. The luminosity dependence of halo bias is significantly weaker in this regime, and thus these galaxies do not provide strong constraints on scatter. However, with DarkSky-400, this sample excludes both $\alpha =0$ and 1 at $p\lt 0.001$.

If we assume that α and scatter are constant with respect to luminosity, then the samples at different thresholds can be combined to produce an overall constraint on α and scatter. Here, we also assume that the constraints from samples at different thresholds are independent. This assumption is only approximately correct for two reasons. First, although the sample at each threshold is dominated by the fainter galaxies, it does include galaxies from higher thresholds. Second, for a given simulation, the clustering signals at different thresholds are also correlated. Here we assume the independence for simplicity, and because the effects are both small, we do not expect that they will significantly impact our conclusions.

The combined constraint from four thresholds (${M}_{r}\lt -20.5$, −21, −21.5, and −22) for DarkSky-400 is shown in Figure 8. This combined constraint excludes both $\alpha =0$ (resembling ${M}_{\mathrm{peak}}$) and 1 (resembling ${v}_{\mathrm{peak}}$) at $p\lt 0.001$, and also excludes zero scatter at $p\lt 0.001$. The best-fit values for DarkSky-400 are $\alpha ={0.57}_{-0.27}^{+0.20}$ and scatter $={0.17}_{-0.05}^{+0.03}$ dex. This value of α is consistent with the value that best matches the observed satellite fraction shown in Section 3.4.

Zoom In Zoom Out Reset image size

The ${w}_{p}({r}_{p})$ corresponding to this best-fit model is shown in Figure 9. We find that, with this new ${\hat{v}}_{\alpha }$ proxy, we can reproduce the ${w}_{p}({r}_{p})$ observed in the SDSS data at all four luminosity thresholds very closely, with a fixed value of α and scatter. In the same figure, the best-fit ${w}_{p}({r}_{p})$ for MDPL is also shown. The large size of the MDPL box results in much smaller errors on the mock ${w}_{p}({r}_{p})$, yet we still obtain excellent agreement with observations. We note that the agreement is good down to the small scales measured by Zehavi et al. (2011).

Zoom In Zoom Out Reset image size

We repeat our study on the clustering with the other simulations listed in Table 1 to test the robustness of our results and to investigate their cosmology dependence. We use four $250\,\mathrm{Mpc}\,{h}^{-1}$ boxes with approximately the same mass resolution but with different cosmologies, three different N-body codes, and different initial conditions.

Figure 10 shows the $p\gt 0.05$ regions in $(\alpha ,\mathrm{scatter})$ from these four boxes, and also DarkSky-400 for reference. Despite the difference between these boxes, the p = 0.05 contours agree with one another very well, and the best-fit points are all in proximity in this parameter space. This result demonstrates the robustness of our analysis. It also suggests that, within the range of cosmologies tested here (all modern ΛCDM cosmologies but with a range of values of, e.g., ${{\rm{\Omega }}}_{M}$ and ${\sigma }_{8}$), the cosmology dependence is weak enough that it cannot be distinguished in these $250\,\mathrm{Mpc}\,{h}^{-1}$ boxes.

Zoom In Zoom Out Reset image size

Since abundance matching models are based on the halo catalogs of N-body simulations, they suffer from the same limitations due to finite resolution. Particularly, for dimmer samples, abundance matching models tend to underpredict small-scale clustering (Guo et al. 2016b). To avoid possible impact of the limited resolution and to obtain unbiased constraints on α and scatter, we only used galaxy samples brighter than ${M}_{r}=-20.5$ in the main results. We demonstrate that α can already be constrained even with only the bright samples.

Nevertheless, the model presented here can also provide good fits to dimmer galaxy samples given its flexibility. Even with the best-fit values of α and scatter used in Figure 9 (i.e., solely from the bright samples), we can obtain reasonably good matches to the clustering of dimmer samples, as shown in Figure 11. We note that for these four dimmer samples with ${M}_{r}\gt -20.5$, observed galaxies are more clustered at small scales and less so at large scales when compared to the model prediction with these particular parameter values. This hints at larger values of both α and scatter, and hence at the mass dependence of α and scatter. However, this hint could be a result of the bias due to resolution limit, and simulations of higher resolution are needed to obtain a definitive conclusion on this possible mass dependence of α and scatter.

Zoom In Zoom Out Reset image size

We note that Reddick et al. (2013) are able to get reasonable fits to the clustering measurements by abundance matching to ${v}_{\mathrm{peak}}$. However, the amount of scatter required in the ${v}_{\mathrm{peak}}$ case is large (0.22 dex) compared to other constraints in the literature (e.g., More et al. 2009). Additionally, this ${v}_{\mathrm{peak}}$ model requires the exclusion of subhalos whose current mass is less than some fraction of the peak mass (using the parameter ${\mu }_{\mathrm{cut}}:= {M}_{\mathrm{vir},\mathrm{now}}/{M}_{\mathrm{peak}}$), and we do not find this to be required with ${\hat{v}}_{\alpha }$. Furthermore, the ${v}_{\mathrm{peak}}$ model did not provide a good fit to the brightest samples when matching to luminosity (see the top left panel of Figure 26 of Reddick et al. 2013), nor did it fit the satellite fraction without excluding halos of low ${M}_{\mathrm{vir},\mathrm{now}}/{M}_{\mathrm{peak}}$ (i.e., with ${\mu }_{\mathrm{cut}}=0$ in Reddick et al. 2013, Figure 22).

The analysis in our present paper uses a larger box with about four times the volume, and thus provides more constraining power. In Figure 10, one can see that given the degeneracy between scatter and α, the smaller Bolshoi box does allow for a region with $\alpha =1$ (corresponding to ${v}_{\mathrm{peak}}$) with higher scatter (≳0.2). This region is consistent with the best-fit result of Reddick et al. (2013), but is ruled out here with the larger DarkSky-400 box.

The core idea of abundance matching is based on two key assumptions: (1) all galaxies live in dark matter density peaks, and (2) galaxy properties are well correlated with halo properties. However, abundance matching should not be viewed as a "parameter-free" model, but instead, can be viewed as a flexible description of the galaxy–halo connection whose parameters can be constrained by observations.

By introducing this new interpolation scheme with the parameter α, we demonstrate that the abundance matching technique is more flexible than the version that was originally proposed. This interpolation scheme also provides a novel interpretation of the matching proxy. Traditionally, when we compare the performance of two abundance matching proxies, we tend to overemphasize the physical meaning of the proxy that performs better. With this α parameter, we demonstrate that, under the framework of abundance matching, there is indeed nothing special about the maximal circular velocity. It is only that observations of clustering statistics favor more concentration dependence than using simply halo mass as a proxy.

On a different note, the α parameter affects the galaxy clustering in the resulting catalog by changing the satellite fraction and the amount of assembly bias. However, we also note that, within the framework of abundance matching, these two effects (assembly bias and satellite fraction) are linked in a specific way when one changes the parameter α. This link is physically justified if all galaxies live in resolved halos and if galaxy and halo properties can be effectively rank matched with one of the proxies considered. On the contrary, the model in Hearin et al. (2016) does not assume this link, and the two effects can be adjusted separately. Nevertheless, with the clustering statistics we tested here, there is no evidence that this link, implicitly assumed when one uses abundance matching, needs to be broken.

This linked feature also enables us to constrain α with only the two-halo clustering. In fact, when we exclude small scales in our analysis, we obtain a consistent, though weaker, constraint on α. This is promising due to the more difficult nature of modeling the smallest scales, which can be impacted by fiber collisions in the data, and by resolution and baryonic effects in the simulations. At present, our best constraint on α still comes from scales in the one-halo regime, but stronger large-scale constraints will be possible as data samples become larger. This result suggests that many of the key details of the galaxy–halo connection may be constrained even without the smallest scales.

It is also important to note that, in addition to the concentration dependence discussed in this work, there is still a rich set of parameters that can potentially be included in abundance matching without breaking the core assumptions mentioned above, such as using non-constant or non-Gaussian scatter, evaluating the matching proxy at different epochs, and adopting different treatments for central and satellite galaxies. With future simulations that have larger volumes and higher resolutions, we can constrain these potential abundance matching parameters, and in return obtain insights on the physical processes of galaxy formation.

Several other probes can provide complementary constraining power on the α parameter. Although, in this paper, we have not completed a full analysis of satellite fractions, Figure 4 already demonstrates that the satellite fraction as a function of luminosity can provide independent constraints on α. Similarly, other group statistics, such as the conditional luminosity function, should also provide additional constraints on α and scatter.

As an example, R. M. Reddick et al. (2016, in preparation) have studied the conditional luminosity function of galaxies in the redMaPPer cluster sample (Rykoff et al. 2014). This sample consists of a very large number of photometrically identified clusters, and hence allows for very small statistical errors on the parameters. This work finds that for models with lower scatter, data require a stronger anti-correlation between satellite occupation and central luminosity. Since satellite occupation is also anti-correlated with host halo concentration (Zentner et al. 2005; Mao et al. 2015), the result of R. M. Reddick et al. (2016, in preparation) implies an anti-correlation between scatter and α (i.e., the concentration dependence of luminosity). This result would then be complementary to the clustering results presented here, since the latter finds a positive correlation between scatter and α, provided that the correlation between scatter and α behaves the same in both luminosity-selected and redMaPPer samples.

Although we do not investigate this directly here, measurements of galaxy voids may be able to put further constraints on the amount of assembly bias (Tinker et al. 2008; Tinker & Conroy 2009). Combining clustering results with measurements of galaxy–galaxy lensing may be able to put further limits on the scatter and on the concentration dependence (e.g., Tasitsiomi et al. 2004; Mandelbaum et al. 2006; Neistein & Khochfar 2012). In addition, data that have more information on redshift dependence, such as the pseudo-multipole correlation function, can provide more constraints on these parameters (Reid et al. 2014; Guo et al. 2016b; Saito et al. 2016).

Another way to put a physical prior on the parameters in empirical models is to compare the model predictions with hydrodynamic simulations. For example, Chaves-Montero et al. (2016) evaluated the galaxy–halo connection in the EAGLE simulation with various abundance matching models with different epochs at which the matching proxy is evaluated. In this work, we establish that the galaxy luminosity has at least some dependence on halo concentration. It will be interesting to fully understand whether and to what extent such a luminosity dependence arises in modern hydrodynamic simulations, and what physical parameters it depends on.