How Are Galaxies Assigned to Halos? Searching for Assembly Bias in the SDSS Galaxy Clustering

For the simulations used in this work, we make use of the Rockstar (Behroozi et al. 2013b) halo catalogs in the z = 0 snapshot of the Small MultiDark of Planck cosmology (hereafter SMDP; Klypin et al. 2016). This high-resolution N-body simulation⁵ was carried out using the GADGET-2 code (see Klypin et al. 2016 and references therein), following the Planck ΛCDM cosmological parameters Ω_m = 0.307, Ω_b = 0.048, Ω_Λ = 0.693, σ₈ = 0.823, n_s = 0.96, h = 0.678. It has a box size of (0.4 h⁻¹ Gpc)³ with 3840³ simulation particles with m_p = 9.6 × 10⁷ h⁻¹ M_⊙ and gravitational softening length of  = 1.5 h⁻¹ kpc.

The Rockstar algorithm, in SMDP, can reliably resolve halos with ≥100 particles, which corresponds to ${M}_{\mathrm{vir}}\geqslant 9.6\,\times {10}^{9}\,{h}^{-1}\,{M}_{\odot }$ (Rodríguez-Puebla et al. 2016). The SMDP simulation satisfies both the size and resolution requirements of studying the galaxy–halo connection in a wide range of luminosity thresholds. For fainter galaxy samples, the faintest galaxies reside in lower mass halos, which requires high resolution. Meanwhile for luminous galaxy samples, their lower number densities requires a large comoving volume.

Furthermore, since we are studying the higher order halo-occupation statistics, the concentration dependence in particular, it is important to use a simulation that can resolve the internal structure of halos. We only keep the Rockstar halos with more than 1000 particles, for which the halo concentration can be determined reliably. This corresponds to a minimum viral mass of ${M}_{\mathrm{vir}}\geqslant 9.6\times {10}^{10}\,{h}^{-1}\,{M}_{\odot }$. Based on this halo mass cut, we exclude the faintest DR7 galaxies samples from our analysis.

In the context of Subhalo-Abundance-Matching models, which require subhalo completeness in the low mass limit, the SMDP simulation has been used to model the faintest galaxy samples in the SDSS data (see Guo et al. 2016). The added advantage of using the SMDP simulation over some of the other industrial simulation boxes commonly used in the literature, such as Bolshoi (Klypin et al. 2011, 2016) simulation, is its larger comoving volume. Larger volume makes this simulation more suitable for performing inference with L_⋆ (corresponding to M_r ∼ −20.44, see Blanton et al. 2003) and more luminous than L_⋆ galaxy samples that occupy larger comoving volumes.

2.2.1. Standard Model without Assembly Bias

For our standard HOD model, we assume the HOD parameterization from Zheng et al. (2007). According to this model, a dark matter halo can host a central galaxy and some number of satellite galaxies. The occupation of the central galaxies follows a nearest-integer distribution, and the occupation of the satellite galaxies follows a Poisson distribution. The expected number of centrals and satellites as a function of the host halo mass of M_h are given by the following equations

$$\begin{eqnarray}&&\langle {N}_{{\rm{c}}}| {M}_{h}\rangle =\displaystyle \frac{1}{2}\left[1+\left(\displaystyle \frac{\mathrm{log}{M}_{h}-\mathrm{log}{M}_{\min }}{{\sigma }_{\mathrm{logM}}}\right)\right],\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&\langle {N}_{{\rm{s}}}| {M}_{h}\rangle ={\left(\displaystyle \frac{{M}_{h}-{M}_{0}}{{M}_{1}}\right)}^{\alpha }.\end{eqnarray} \tag{ 2 }$$

For populating the halos with galaxies, we follow the procedure described in Hahn et al. (2017), and Hearin et al. (2016). The central galaxies are assumed to be at the center of the host dark matter halos. We assume that the central galaxies are at rest with respect to the bulk motion of the halos and their velocities are given by the velocity of the center of mass of their host halo. Note that this assumption is shown to be violated in brighter than L_⋆ galaxy samples (see Guo et al. 2015b). But since we are not considering the redshift-space 2PCF multipoles in our study, we do not expect this velocity bias to impact our inference. We place the satellite galaxies within the virial radius of the halo following a Navarro–Frenk–White profile (hereafter NFW; Navarro et al. 2004). This approach is different from other simulation-based halo-occupation modeling techniques (see Guo et al. 2016; Zheng & Guo 2016) in that the positions of the satellites are not assigned to the dark matter particles in the N-body simulation. The concentration of the NFW profile comes directly from the halo catalog. The velocities of the satellite galaxies are given by two components. The first component is the velocity of the host halo. The second component is the velocity of the satellite galaxy with respect to the host halo, which is computed following the solution to the NFW profile Jeans equations (More et al. 2009). We refer the readers to Hearin et al. (2016) for a more comprehensive and detailed discussion of the forward modeling of the galaxy mock catalogs.

2.2.2. Model with Assembly Bias

Now let us provide a brief overview of HOD modeling with Heaviside assembly bias (referred to as the decorated HOD) introduced in Hearin et al. (2016). At a fixed halo mass M_h, halos are split into two populations: a population of halos with the 0.5 percentile of highest concentration and a population of halos with the 0.5 percentile of lowest concentration. For simplicity, we call the first population "type-1" halos, and the second population "type-2" halos. In the decorated HOD model, the expected number of central and satellite galaxies at a fixed halo mass M_h in the two populations are given by

$$\begin{eqnarray}&&\langle {N}_{c,i}| {M}_{h},c\rangle =\langle {N}_{c}| {M}_{h}\rangle +{\rm{\Delta }}{N}_{c,i},\,i=1,2\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&\langle {N}_{s,i}| {M}_{h},c\rangle =\langle {N}_{s}| {M}_{h}\rangle +{\rm{\Delta }}{N}_{s,i},\,i=1,2\end{eqnarray} \tag{ 4 }$$

where $\langle {N}_{c}| {M}_{h}\rangle$ and $\langle {N}_{s}| {M}_{h}\rangle$ are given by Equations (1) and (2), respectively, and we have ΔN_s,1 + ΔN_s,2 = 0 and ΔN_c,1 + ΔN_c,2 = 0. These two conditions ensure the conservation of HOD. At a given host halo mass M_h, the central occupation of the the two populations follows a nearest-integer distribution with the first moment given by Equation (3) and the satellite occupation of the the two populations follows a Poisson distribution with the first moment given by Equation (4).

In this occupation model, the allowable ranges that quantities ΔN_c,1 and ΔN_s,1 can take are given by

$$\begin{eqnarray}&&{\rm{\Delta }}{N}_{c,i,\min }\leqslant {\rm{\Delta }}{N}_{c,i}\leqslant {\rm{\Delta }}{N}_{c,i,\max },\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{\rm{\Delta }}{N}_{c,,\min }=\max \{-\langle {N}_{c}| {M}_{h}\rangle ,\langle {N}_{c}| {M}_{h}\rangle -1\},\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&{\rm{\Delta }}{N}_{c,i,\max }=\min \{\langle {N}_{c}| {M}_{h}\rangle ,1-\langle {N}_{c}| {M}_{h}\rangle \},\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&-\langle {N}_{s}| {M}_{h}\rangle \leqslant {\rm{\Delta }}{N}_{s,i}\leqslant \langle {N}_{s}| {M}_{h}\rangle .\end{eqnarray} \tag{ 8 }$$

Afterward, the assembly bias parameter ${ \mathcal A }$ is defined in the following way:

$$\begin{eqnarray}&&{\rm{\Delta }}{N}_{\alpha ,1}({M}_{h})=| {{ \mathcal A }}_{\alpha }| {\rm{\Delta }}{N}_{\alpha ,1}^{\max }({M}_{h})\,\mathrm{if}\,{{ \mathcal A }}_{\alpha }\gt 0,\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&{\rm{\Delta }}{N}_{\alpha ,1}({M}_{h})=| {{ \mathcal A }}_{\alpha }| {\rm{\Delta }}{N}_{\alpha ,1}^{\min }({M}_{h})\,\mathrm{if}\,{{ \mathcal A }}_{\alpha }\lt 0,\end{eqnarray} \tag{ 10 }$$

where the subscript α = c, s stands for the centrals and satellites respectively, and ${\rm{\Delta }}{N}_{\alpha ,1}^{\max }({M}_{h})$, ${\rm{\Delta }}{N}_{\alpha ,1}^{\min }({M}_{h})$ are given by Equations (5) and (8).

For a given ${{ \mathcal A }}_{\alpha }$, once ΔN_α,1 is computed using Equation (9)—if ${{ \mathcal A }}_{\alpha }\gt 0$—or Equation (10)—if ${{ \mathcal A }}_{\alpha }\lt 0$—, ΔN_α,2 = 1 − ΔN_α,1 is computed. At a fixed halo mass M_h, once the first moments of occupation statistics for the type-1 and type-2 halos are determined, we perform the same procedure described in 2.2.1 to populate the halos with mock galaxies.

2.2.3. Redshift-space Distortion

As described in Hearin et al. (2016b) and Hahn et al. (2017), this approach makes no appeal to the fitting functions used in the analytical calculation of the 2PCF. The accuracy of these fitting functions is limited (Tinker et al. 2008a, 2010; Watson et al. 2013). Our approach also does not face the known issues of the treatment of halo exclusion and scale-dependent bias that can lead to potential inaccuracies in halo-occupation modeling (see van den Bosch et al. 2013).

The projected 2PCF w_p(r_p) can be computed by integrating the 3D redshift-space 2PCF ξ(r_p, π) along the line of sight (where r_p and π denote the projected and line-of-sight separation of galaxy pairs, respectively):

$$\begin{eqnarray}&&{w}_{p}({r}_{p})=2{\int }_{0}^{{\pi }_{\max }}\xi ({r}_{p},\pi )\,d\pi \end{eqnarray} \tag{ 11 }$$

For our 2PCF calculations, we use the w_p measurement functionality of the fast and publicly available pair-counter code CorrFunc⁶ (Sinha 2016). To be consistent with the SDSS measurements described in Section 3, w_p(r_p) is obtained by the line-of-sight integration to ${\pi }_{\max }=40\,{h}^{-1}\,\mathrm{Mpc}$. Note that w_p(r_p) is measured in units of h⁻¹ Mpc. To be consistent with Guo et al. (2015b), we use the same binning (as specified in Section 3) to measure w_p. In addition to the projected 2PCF, we use the number density given by the number of mock galaxies divided by the comoving volume of the SMDP simulation.

Note that a full forward model of the data requires running the simulation at different redshifts, generation of light cones, accounting for the complex survey geometry and systematic errors such as fiber collisions. Using the z = 0 output of the SMDP simulation in our forward model of the spatial distribution of galaxies is only an approximation. This approximation can be justified by the small redshift range of the SDSS DR7 main galaxy sample. As described in Zehavi et al. (2011), using random catalogs with angular window function of the data in measurements of galaxy clustering accounts for the geometry of the data. As described in Section 3, the fiber collision correction method of Guo et al. (2012) is applied to the SDSS clustering measurements used in this study. Therefore we do not account for that effect in our forward model.

Given the SDSS measurements described in Section 3, we aim to constrain the HOD model without assembly bias (described in 2.2.1), and the HOD model with assembly bias (described in 2.2.2) for each luminosity-threshold sample, by sampling from the posterior probability distribution $p(\theta | d)\propto p(d| \theta )\pi (\theta )$ where θ denotes the parameter vector and d denotes the data vector. In the standard HOD modeling θ is given by

$$\begin{eqnarray}&&\theta =\{\mathrm{log}{M}_{\min },\,{\sigma }_{\mathrm{logM}},\,\mathrm{log}{M}_{0},\,\alpha ,\,\mathrm{log}{M}_{1}\},\end{eqnarray} \tag{ 15 }$$

and in the HOD modeling with assembly bias we have

$$\begin{eqnarray}&&\theta =\{\mathrm{log}{M}_{\min },\,{\sigma }_{\mathrm{logM}},\,\mathrm{log}{M}_{0},\,\alpha ,\,\mathrm{log}{M}_{1},\,{{ \mathcal A }}_{\mathrm{cen}},\,{{ \mathcal A }}_{\mathrm{sat}}\}.\end{eqnarray} \tag{ 16 }$$

Furthermore, data (denoted by d) is the combination of [n_g, w_p(r_p)]. The negative log-likelihood (assuming negligible covariance between n_g and w_p(r_p)) is given by

$$\begin{eqnarray}\begin{array}{rcl}-2\mathrm{ln}p(d| \theta ) & = & \displaystyle \frac{{\left[{n}_{g}^{\mathrm{data}}-{n}_{g}^{\mathrm{model}}\right]}^{2}}{{\sigma }_{n}^{2}}+{\rm{\Delta }}{w}_{{\rm{p}}}^{{\rm{T}}}\widehat{{C}^{-1}}{\rm{\Delta }}{w}_{{\rm{p}}}\\ & & +\ \mathrm{const}.,\end{array}\end{eqnarray} \tag{ 17 }$$

where Δ${w}_{{\rm{p}}}$ is a two-dimensional vector, ${\rm{\Delta }}{w}_{{\rm{p}}}({r}_{{\rm{p}}})={w}_{{\rm{p}}}^{\mathrm{data}}({r}_{{\rm{p}}})\,-{w}_{{\rm{p}}}^{\mathrm{model}}({r}_{{\rm{p}}})$, and $\widehat{{C}^{-1}}$ is the estimate of the inverse covariance matrix that is related to the inverse of the total covariance matrix (Equation (14)) ${\widehat{C}}^{-1}$, following Hartlap et al. (2007):

$$\begin{eqnarray}&&\widehat{{C}^{-1}}=\displaystyle \frac{N-d-2}{N-1}\,{\widehat{C}}^{-1},\end{eqnarray} \tag{ 18 }$$

where N = 400 is the number of the jackknife samples, and d = 12 is the length of the data vector w_p. Another important ingredient of our analysis is specification of the prior probabilities π(θ) over the parameters of the halo-occupation models considered in this study. For both models, we use uniform flat priors for all the parameters. The prior ranges are specified in the Table 1. Note that a uniform prior between −1 and 1 is chosen for assembly bias parameters since these parameters are, by definition, bounded between −1 and 1.

For sampling from the posterior probability, given the likelihood function (see Equation (17)) and the prior probability distributions (see Table 1), we use the affine-invariant ensemble MCMC sampler (Goodman & Weare 2010) and its implementation emcee (Foreman-Mackey et al. 2013). In particular, we run the emcee code with 20 walkers and we run the chains for at least 12,000 iterations. We discard the first one-third part of the chains as burn-in samples and use the reminder of the chains as production MCMC chains. Furthermore, we perform an auto-correlation convergence test (Goodman & Weare 2010) to ensure that the MCMC chains have reached convergence.

Once the posterior probability distributions are obtained, we estimate the corresponding minimum χ² (or, equivalently, the maximum likelihood) of each HOD model for each luminosity-threshold sample. Finding the minimum χ² is done by using the scipy implementation of the BFGS algorithm (Byrd et al. 1994). The BFGS parameter bounds are given by the 68% confidence intervals of the posterior probability distributions.

In this section, we present the constraints derived for the two assembly bias parameters: the satellite assembly bias parameter (${{ \mathcal A }}_{\mathrm{sat}}$) and the central assembly bias parameter ${{ \mathcal A }}_{\mathrm{cen}}$. As shown in Figure 1, for all the five luminosity-threshold samples in the SDSS DR7 data, our constraints on the parameter ${{ \mathcal A }}_{\mathrm{sat}}$ are consistent with zero. On the other hand, our constraints on the parameter ${{ \mathcal A }}_{\mathrm{cen}}$, albeit not tightly constrained, show a trend can be summarized as the following. In the most luminous galaxy samples, i.e., M_r < −21.5 and M_r < −21, ${{ \mathcal A }}_{\mathrm{cen}}$ is poorly constrained and the constraints are equivalent to zero. As we investigate less luminous samples, M_r < −20.5, −20, −19.5, our constraints on ${{ \mathcal A }}_{\mathrm{cen}}$ shift toward positive values, with the M_r < −20.5 sample favoring the highest values for ${{ \mathcal A }}_{\mathrm{cen}}$.

Zoom In Zoom Out Reset image size

According to Hearin et al. (2016), the underlying theoretical considerations for explaining assembly bias (as formulated in the decorated HOD model) of the central and satellite galaxies are different. The large-scale clustering, or the two-halo term in the galaxy clustering, is mainly governed by the clustering of the central galaxies. The central galaxy clustering can be thought of as the weighted average over the halo clustering. The large-scale bias of the dark matter halo clustering depends not only on mass, but also on the other properties of halos beyond mass, such as concentration (Wechsler et al. 2006; Gao & White 2007; Miyatake et al. 2016), spin (Gao & White 2007), formation time (Gao & White 2007; Li et al. 2008), and maximum circular velocity of the halo V_max (Sunayama et al. 2016).

In particular, findings of Wechsler et al. (2006) and Sunayama et al. (2016) have demonstrated that for halos with mass below the collapse mass (M ≤ M_col ≃ 10^12.5 M_⊙) the large-scale bias of high-V_max (or equivalently high-c halos at a fixed halo mass) is larger than that of the low-V_max (low-c halos). This signal reverses and weakens for the high mass halos (M ≥ M_col). Note that the halo concentration traces the maximum circular velocity V_max such that halos with higher V_max have higher concentration and vice versa (see Prada et al. 2012). For halos described by NFW profile, at a fixed halo mass, halos with higher values of concentration have higher values of V_max.

Furthermore, investigation of the scale dependence of halo assembly bias has shown that the ratio of the bias of high-V_max halos to low-V_max halos has a bump-like feature in the quasi-linear scales ∼0.5 Mpc h⁻¹–5 Mpc h⁻¹. From the theoretical standpoint, this phenomenon has been attributed to a population of distinct halos with $M\sim {10}^{11.7}\,{h}^{-1}\,{M}_{\odot }$ at the present time that are close to the most massive groups and clusters (Sunayama et al. 2016, and see More et al. 2016 for the observational investigation of this signal by means of galaxy–galaxy lensing). The clustering of these population of halos is therefore dictated by that of the massive halos. Note that this scale-dependent bias feature vanishes in high mass halos.

Consequently, at a fixed halo mass less than M_col, assignment of more central galaxies to the high-c halos (higher expected number of central galaxies in the high-c halos) gives rise to a boost in the galaxy clustering in the linear scales as well as in regimes corresponding to the one-halo to two-halo transition. For the more massive halos (M ≥ M_col), we expect the large-scale clustering boost to reverse sign, and the quasi-linear bump feature to vanish.

Figure 2 demonstrates the 68% and 95% posterior predictions for the projected 2PCF w_p from the occupation model without assembly bias (shown in red) and the occupation model with assembly bias (shown in blue) for all five luminosity-threshold samples. For the brightest galaxies, M_r < −21.5 and M_r < −21.0, the posterior prediction of w_p from the two models are consistent with one another. Note that these galaxies reside in the most massive halos (M > M_col) for which the scale dependence of the halo assembly bias and the difference between the large-scale bias of the high-c and low-c halos is negligible.

Zoom In Zoom Out Reset image size

Figure 3 shows the fractional difference between the 68% and 95% posterior predictions of w_p and the SDSS data. It is evident from Figure 3 that some improvement on modeling the clustering of the samples of L_⋆ and slightly less brighter than L_⋆ galaxies can be achieved by employing the more complex halo-occupation model with assembly bias. As a result of apportioning more central galaxies to the high-c halos relative to the low-c halos, in the samples with luminosity thresholds of M_r < −20.5, −20, −19.5, the posterior predictions for w_p are slightly improved in the intermediate scales (1 ∼ 2 Mpc h⁻¹) and large scales. This can be also noted in significantly lower χ² values, at the cost more model flexibility and higher degrees of freedom, achieved by the assembly bias model in these luminosity-threshold samples (see Table 2).

Zoom In Zoom Out Reset image size

The luminosity-dependent trend in the constraints on the central assembly bias for the six dimmest samples can be attributed to the fact that the halo concentration is highly correlated with the maximum circular velocity V_max, which is a tracer of the potential well of dark matter halos (Prada et al. 2012). In a dark matter halo described by an NFW profile, the depth of the gravitational potential well of dark matter halos can be directly measured by the maximum circular velocity V_max (van den Bosch et al. 2014). In particular, the magnitude of the potential well at the center of an NFW halo, where the central galaxy is assumed to reside, scales as ${V}_{\max }^{2}$. More specifically, we have:

$$\begin{eqnarray}&&{\rm{\Phi }}(r=0)=-{\left(\displaystyle \frac{{V}_{\max }}{0.465}\right)}^{2},\end{eqnarray} \tag{ 19 }$$

where Φ(r = 0) is the central potential of an NFW profile. Note that V_max is also the quantity often used in the abundance-matching technique in which the luminosity of galaxies is monotonically matched to V_max (see for example Reddick et al. 2013; Lehmann et al. 2017; Guo et al. 2016; Rodríguez-Puebla et al. 2016). The trend between the constraint on ${{ \mathcal A }}_{\mathrm{cen}}$ and the luminosity threshold of the samples may suggest that, at a fixed halo mass, the central galaxies in the samples (M_max < −19.5, −20, −20.5), have a tendency to reside in dark matter halos with deeper gravitational potential well.

The satellite assembly bias can only significantly alter the galaxy clustering at small to- intermediate scales. Assigning more satellite galaxies to lower (or higher) concentration halos affects the one-halo term through increasing the satellite–satellite pair counts $\langle {N}_{{\rm{s}}}{N}_{{\rm{s}}}\rangle$ (Hearin et al. 2016). This results in boosting the small-scale clustering. But, as pointed out by Hearin et al. (2016), the amount by which small-scale clustering increases also depends on the sign of the central assembly bias parameter ${{ \mathcal A }}_{\mathrm{cen}}$. Formation history of the halos can lead to the dependence of the abundance of subhalos on halo concentration (Zentner et al. 2005; Mao et al. 2015) at fixed halo mass, and since the occupation of the satellite galaxies is related to the abundance of subhalos, the satellite occupation may depend on halo concentration.

However, our results suggest that. for all the luminosity-threshold samples considered in this study, the satellite assembly bias parameter is largely unconstrained and consistent with zero. We do not expect the galaxy clustering data to be a sufficient statistics for obtaining constraints on the satellite assembly bias. Group statistics probes the high mass end of the galaxy–halo connection and is sensitive to the parameters governing the satellite population (see Hearin et al. 2013; Hahn et al. 2017). Therefore these measurements may shed some light on potential presence of assembly bias in satellite population.

It is important to note that the halo mass range in which the central assembly bias ${{ \mathcal A }}_{\mathrm{cen}}$ affects the central galaxy population is the mass range in which the condition $0\lt \langle {N}_{{\rm{c}}}| M\rangle \lt 1$ is met. Consequently, a larger scatter parameter ${\sigma }_{\mathrm{logM}}$ increases the dynamical mass range in which assembly bias affects the galaxy clustering. Note that, in the luminosity regimes for which we obtain tighter constraints on ${{ \mathcal A }}_{\mathrm{cen}}$, the best-estimate values of the scatter parameter ${\sigma }_{\mathrm{logM}}$ appear to be higher in the model with assembly bias. This is evident in Figure 4. The model with assembly bias tends to push ${\sigma }_{\mathrm{logM}}$ to higher values. This can be attributed to the tendency of this model to increase the effective dynamical mass range of central assembly bias (Hearin et al. 2016).

Zoom In Zoom Out Reset image size

As shown in Table 2, in the HOD model with assembly bias, the constraints found on the scatter parameter are not tight. This is in keeping with the results of Guo et al. (2015b), which uses the same SDSS measurements and finds that the scatter parameter remains largely unconstrained when only n_g and w_p are used as observables. Note that scatter is better constrained for the most luminous galaxy samples (this is attributed to the steep dependence of the halo bias and halo mass function on halo mass in the high mass end). But since these samples live in the most massive halos, we do not expect the tighter constraints on scatter to help us constrain the central assembly bias parameter. Guo et al. (2015b) shows that, by employing additional measurements such as the monopole (ξ₀), quadruple (ξ₂), and hexadecapole (ξ₄), one can obtain tighter constraints on the scatter parameter. Tightening the constraints on the scatter parameter can lead to more precise inference of the central assembly bias parameter.

As shown in Table 2, our constraints on the underlying standard HOD model obtained from the model with assembly bias and the model without assembly bias are in good agreement. The only cases in which there are mild tensions between the constraints found from the two models on the underlying HOD parameters are the M_r < −20 and the M_r < −20.5 samples. However, these tensions are still within 1σ level. For instance, Figure 4 shows that in the M_r < −20.5 sample, the constraint on the parameter α found from the model without assembly bias favors slightly higher values than the constraint found from the model with assembly bias. Also the scatter parameter ${\sigma }_{\mathrm{logM}}$ is more tightly constrained in the standard HOD model. However, it is important to emphasize that these constraints are still in agreement with each other within a 1σ level.

Zentner et al. (2014) shows that, in the mock catalogs that exhibit significant levels of assembly bias, using a simple mass-only occupation model can lead to considerable biases in inference of the galaxy–halo connection parameters. Although we cannot rule out moderate levels of assembly bias in our findings, we do not find any considerable discrepancy between the two models in terms of estimating the underlying HOD parameters.

A few galaxy–halo connection methods have been proposed in the literature that give rise to assembly bias in the galaxy population. Zentner et al. (2014) demonstrates that the abundance-matching techniques based on V_max (Conroy et al. 2006; Hearin & Watson 2013; Reddick et al. 2013) exhibit some levels of assembly bias. We aim to provide a comparison between the impact of assembly bias on clustering in these mock catalogs and the mock catalogs predicted from our constraints on the decorated HOD model for L_⋆-type galaxies. In particular, we consider the abundance-matching catalogs produced by Hearin & Watson (2013). These catalogs have been extensively studied for examining potential systematic effects of galaxy assembly bias on cosmological (McEwen & Weinberg 2016) and halo-occupation (Zentner et al. 2014) parameter inferences.

This abundance-matching catalog was built based on the Bolshoi N-body simulation (Klypin et al. 2011) using the adaptive refinement tree code (ART Kravtsov et al. 1997). The box size for this simulation is 250 h⁻¹ Mpc, the number of simulation particles is 2048³, the mass per simulation particle m_p is $1.35\times {10}^{8}\,{h}^{-1}\,{M}_{\odot }$, and the gravitational softening length is 1 h⁻¹ kpc. The halos and subhalos in this simulation are identified using the Rockstar algorithm (Behroozi et al. 2013b). The Hearin & Watson (2013) catalogs make use of V_peak (maximum V_max throughout the assembly history of halo) as the subhalo property to be matched to galaxy luminosity.

As noted by Zentner et al. (2014) and McEwen & Weinberg (2016), these galaxy mock catalogs show significant levels of assembly bias in the central galaxy population. This has been demonstrated by investigating the difference in w_p between the randomized mock catalogs and the original mock catalogs. Randomization is performed in a procedure described in Zentner et al. (2014), which we briefly summarize here. First, halos are divided into different bins of halo mass with width of 0.1 dex. Then all central galaxies are shuffled among all halos within each bin. Once the centrals have been shuffled, within each bin, the satellite systems are shuffled among all halos in that mass bin, preserving their relative distance to the center of halo. This procedure preserves the HOD, but erases any dependence of the galaxy population on the assembly history of halos. Therefore, assembly bias is erased in the randomized galaxy catalog.

For L_⋆ galaxies, the difference in w_p between the randomized and the original catalogs of Hearin & Watson (2013) is shown in Figure 5 with the red curves. As demonstrated in Figure 5 (and as previously noted by Zentner et al. 2014; McEwen & Weinberg 2016), the relative difference in w_p is only significant in relatively large scales (r_p > 1 Mpc h⁻¹). This implies that, in these catalogs, only the central occupation is affected by assembly bias. This is in agreement with our findings.

Zoom In Zoom Out Reset image size

Furthermore, we investigate whether the impact of assembly bias on galaxy clustering predicted by our findings are in agreement with the abundance-matching catalogs of Hearin & Watson (2013). First, we make random draws from the 68% confidence intervals of the posterior probability distribution function over the parameters of the model with assembly bias. Then, we create mock catalogs with these random draws and then we compute the difference in w_p between the randomized catalogs and the original catalogs. The relative difference in w_p predicted from our constraints are shown with blue curves in Figure 5.

We note that our findings follow the same trend. That is, we see negligible difference in the small-scale clustering and more considerable differences in w_p on larger scales (${r}_{p}\gt 1\,\mathrm{Mpc}\,{h}^{-1}$). Similar to findings of Zentner et al. (2014), Lehmann et al. (2017), and McEwen & Weinberg (2016), we see that the impact of assembly bias on galaxy clustering becomes less significant in brighter galaxy samples. Furthermore, we notice that our mock catalogs favor more moderate changes in galaxy clustering as a result of assembly bias.

Using the constraints on the model with assembly bias and without, we want to address the question of whether assembly bias is supported by galaxy clustering observations. Within the standard HOD framework, the distribution of the galaxies is modeled using a simple description based on the mass-only ansatz: $P(N| M)$. The decorated HOD model, however, provides a more complex description of the data by adding a secondary halo property (halo concentration in this study) and a more flexible occupation model: $P(N| M,c)$.

In order to investigate whether observations demand a higher level of model complexity, we compare the models with and without assembly bias. In particular, we make use of two simple methods for model comparison: the Akaike Information Criterion (AIC, Akaike 1974, see Gelman et al. 2014 for detailed discussion on AIC), and the Bayesian Information Criterion (BIC, Schwarz 1978). BIC and AIC are more computationally tractable than alternative approaches such as computing the fully marginalized likelihood. The underlying assumption of these information criteria is that models that yield higher likelihoods are more preferable but, at the same time, models with more flexibility are penalized.

Suppose ${{ \mathcal L }}^{\star }$ is the maximum likelihood achieved by the model, N_par is the number of free parameters in the model, and N_data is the number of data points in the data set. Then we have

$$\begin{eqnarray}&&\mathrm{BIC}=-2\,\mathrm{ln}{{ \mathcal L }}^{\star }+{N}_{\mathrm{par}}\,\mathrm{ln}{N}_{\mathrm{data}},\end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray}&&\mathrm{AIC}=-2\,\mathrm{ln}{{ \mathcal L }}^{\star }+2{N}_{\mathrm{par}}.\end{eqnarray} \tag{ 21 }$$

Given a data set, models with lower values of AIC and BIC are more desired. That is, higher model complexity (larger N_par) is favored only if ${{ \mathcal L }}^{\star }$ is sufficiently higher.

Figure 6 shows the comparison between the BIC and AIC scores for the model with assembly bias and the model without assembly bias. The more complex decorated HOD model is favored by the information criteria if ΔBIC or ΔAIC < 0, and strongly favored if ΔBIC or ΔAIC < −5. With the exception of the M_r < −20.5 sample, the model without assembly bias is still preferred by both information criteria. Although including assembly bias improves the fit to the observed clustering, overall these improvements are not enough to justify adding the extra assembly bias parameters to the model.

Zoom In Zoom Out Reset image size

In the M_r < −20.5 sample where both AIC and BIC scores improve for the assembly bias model, we have the greatest improvement in terms of the minimum χ². This is also the sample with the tightest constraints on the parameters ${{ \mathcal A }}_{\mathrm{cen}}$ and ${{ \mathcal A }}_{\mathrm{sat}}$. While the information criteria favor the model with assembly bias in this sample, both ΔAIC and ΔBIC are only marginally negative. Therefore, even in this sample, the improvement in the achieved ${{ \mathcal L }}^{\star }$ by additional parameters is not enough to make the decorated model strongly favored by AIC and BIC.

Given the SDSS clustering measurements described in Section 3, We repeat the inference of the assembly bias parameters ${{ \mathcal A }}_{\mathrm{sat}}$ and ${{ \mathcal A }}_{\mathrm{cen}}$ with the BolshoiP simulation (Klypin et al. 2016). This N-body simulation is carried out with a setting similar to the Bolshoi simulation with the exception that, in the BolshoiP simulation, Planck cosmology is adapted and the mass per simulation particle is $1.49\times {10}^{8}\,{h}^{-1}\,{M}_{\odot }$.

The summary of constraints are shown in Figure 7. In Figure 7, the constraints from the SMDP and the BolshoiP simulations are shown with circles and crosses, respectively. Additionally, the upper and lower bounds on the inferred parameters reported by Zentner et al. (2016) are shown in shaded blue regions. In the case of central assembly bias, our constraints based on the SMDP and the BolshoiP simulations are consistent. For the luminosity-threshold samples M_r < −20.5, −20, −19.5, where the central assembly bias parameters are strongly positive, the constraints obtained from the SMDP simulation are tighter. For the three faintest luminosity-threshold samples, our constraints are fully consistent with those of Zentner et al. (2016). For the brightest samples (M_r < −21.5, −21), however, the Zentner et al. (2016) constraints on ${{ \mathcal A }}_{\mathrm{cen}}$ are slightly higher than ours.

Zoom In Zoom Out Reset image size

In the case of satellite assembly bias, our constraints based on the SMDP and the BolshoiP simulations are fully consistent with each other and consistent with ${{ \mathcal A }}_{\mathrm{sat}}$ = 0. For the luminosity-threshold samples M_r < −20.5, −20, where the estimated central assembly bias parameters are strongly positive, the ${{ \mathcal A }}_{\mathrm{sat}}$ constraints obtained from the SMDP simulation are tighter. For all the luminosity-threshold samples except for the brightest sample, the Zentner et al. (2016) constraints on ${{ \mathcal A }}_{\mathrm{sat}}$ are slightly higher than ours. But given the large uncertainties on ${{ \mathcal A }}_{\mathrm{sat}}$, our constraints on this parameter seem to be consistent with those of Zentner et al. (2016).

In terms of the effect of assembly bias on galaxy clustering, note that mock catalogs created using the inferred parameters with the SMDP simulation show the same behavior as what we observe in the abundance-matching catalogs presented in Zentner et al. (2014) and Lehmann et al. (2017). That is, the difference in w_p between the mock catalogs and the randomized catalogs is mostly on large scales where the clustering is governed by the central galaxies. That is, the impact of assembly bias on the satellite occupation is negligible and only the central occupation is affected.