ON DETECTING HALO ASSEMBLY BIAS WITH GALAXY POPULATIONS

Similar to previous observational studies of assembly bias, we assume that if the formation history of galaxies in a halo could somehow be related to that of the host halo, it is the central galaxy that best traces the halo formation history. Therefore, the starting point of our analysis is a catalog of central galaxies, for which we use the galaxy group catalog of Yang et al. (2007, hereafter Y07). The version of the catalog used is based on data release seven (DR7; Abazajian et al. 2009) of the Sloan Digital Sky Survey (SDSS; York et al. 2000). The Y07 group-finding algorithm is essentially an adaptive matched filter applied to the spectroscopic sample of SDSS. Its unique features include (1) the identification of groups of all richness, ranging from massive clusters down to systems with just one member, (2) the designation of central and satellite galaxies among the group members, (3) the assignment of the total mass of the groups via an abundance matching method (see Section 1; in this study we adopt the mass based on the luminosity content ranking). The DR7 group catalog contains 374,052 systems¹² over 7748 deg² out to z = 0.2 (Wang et al. 2014).

After selecting the galaxies identified as centrals from the Y07 group catalog, we seek measurements of sSFR and SFH from public data sets. For the former, we adopt the measurements from the DR7 version of the MPA/JHU value-added galaxy catalog (Brinchmann et al. 2004), which are derived from emission lines as well as the strength of the 4000 Å break.

For the resolved SFH, we utilize the results of the VESPA algorithm (Tojeiro et al. 2009), which are available for query in the WFCAM science archive. In short, VESPA considers combinations of bursts of star formation of different metallicity, and chooses the one that best fits the observed SDSS spectrum. For each galaxy, it provides the stellar mass formed in 16 temporal (age) bins, where the first bin spans the period of 9–14 Gyr before the redshift of the galaxy (and other bins are closer to the present epoch). In addition to using two of the popular stellar population synthesis models (Bruzual & Charlot 2003; Maraston 2005) to model the spectrophotometric evolution of the galaxies, two different dust extinction models are also considered in VESPA. Therefore, for each galaxy, four SFHs are provided (Tojeiro et al. 2009). For simplicity, we adopt the SFHs based on the Bruzual & Charlot (2003) model, and demand the consistency between the SFHs derived from the two dust models (that is, a galaxy should be classified as early- or late-forming in both models).

In the following, unless otherwise stated, an early-forming galaxy is defined as having formed 50% of its final stellar mass in the first temporal bin, or roughly by z = 1.9 for a galaxy observed at z = 0.1, which is close to the mean redshift of the samples we use in Section 4. In contrast, a late-forming galaxy would have formed 50% of its final mass after the first temporal bin (or at $z\lt 1.9$ if it is observed at z = 0.1).

To measure the large-scale bias of the galaxy samples, we first calculate the redshift space, two-point auto-correlation function via the standard estimator (Landy & Szalay 1993)

$$\begin{eqnarray}&&{\xi }_{s}=\displaystyle \frac{\langle \mathrm{DD}\rangle -2\langle \mathrm{DR}\rangle +\langle \mathrm{RR}\rangle }{\langle \mathrm{RR}\rangle },\end{eqnarray} \tag{ 1 }$$

where $\langle \mathrm{DD}\rangle$, $\langle \mathrm{DR}\rangle$, and $\langle \mathrm{RR}\rangle$ are the normalized numbers of data–data, data–random, and random–random pairs in a given separation bin, respectively. In practice, to bypass complications due to redshift space distortions, we measure ${\xi }_{s}$ in two dimensions, both perpendicular to and along the line of sight (denoted as r_p- and π-directions, respectively), then compute the projected correlation function by integrating ${\xi }_{s}$ over the π direction

$$\begin{eqnarray}&&{w}_{p}({r}_{p})=2{\displaystyle \int }_{0}^{{\pi }_{{\rm{up}}}}d\pi {\xi }_{s}({r}_{p},\pi ).\end{eqnarray} \tag{ 2 }$$

Following common practice, the integration upper limit is chosen to be ${\pi }_{{\rm{up}}}=60\;{h}^{-1}$ Mpc, so that the results are not affected by redshift space distortion below ${r}_{p}\sim 30\;{h}^{-1}$ Mpc (e.g., Padmanabhan et al. 2007). As will be described below, our galaxy samples have sizes of several tens of thousands. We therefore use one million random points over the DR7 footprint, generated with the mask provided by the NYU value-added galaxy catalog (Blanton et al. 2005). The covariance matrix of the correlation function is calculated using the jackknife resampling method. We have divided the DR7 footprint into 100 equal-area subregions, and have constructed the jackknife samples by omitting each of the subregions in turn. The covariance matrix is then estimated as

$$\begin{eqnarray}&&{\rm{C}}({w}_{p,i},{w}_{p,j})=\displaystyle \frac{N-1}{N}\displaystyle \sum _{{\ell }=1}^{{\ell }=N}({w}_{p,i}^{{\ell }}-\bar{{w}_{p,i}})({w}_{p,j}^{{\ell }}-\bar{{w}_{p,j}}),\end{eqnarray} \tag{ 3 }$$

where N = 100, ℓ is the index of the jackknife samples, $i,j$ are indices of r_p bins, and $\bar{{w}_{p,i}}$ is the mean from all the jackknife samples. Note that N is much larger than the number of radial bins in our correlation function measurements (see, e.g., Appendix D of Hirata et al. 2004, which includes a method of simulating the impact of noise on the jackknife-based ${\chi }^{2}$ that we use here to estimate p-values), and the size of the jackknife subregions is much larger than the maximum size of the clustering used for the analysis ($\lesssim 35\;{h}^{-1}$ Mpc), so that the samples nearly satisfy the i.i.d. (independent and identically distributed) requirement. These conditions validate the use of jackknife resampling.

Galaxy–galaxy weak lensing, the coherent tangential shear of background galaxies due to foreground lens galaxies, provides a simple way to probe the connection between the lens galaxies and matter via their cross-correlation function ${\xi }_{{gm}}$. This cross-correlation can be related to the projected surface density¹³ via

$$\begin{eqnarray}&&{\rm{\Sigma }}({r}_{p})=\bar{\rho }\int [1+{\xi }_{{gm}}(\sqrt{{r}_{p}^{2}+{\pi }^{2}})]d\pi .\end{eqnarray} \tag{ 4 }$$

The surface density is then related to the observable quantity for lensing,

$$\begin{eqnarray}&&{\rm{\Delta }}{\rm{\Sigma }}({r}_{p})={\gamma }_{t}({r}_{p}){{\rm{\Sigma }}}_{c}=\bar{{\rm{\Sigma }}}(\lt {r}_{p})-{\rm{\Sigma }}({r}_{p}).\end{eqnarray} \tag{ 5 }$$

This observable quantity can be expressed as the product of two factors, a tangential shear ${\gamma }_{t}$ and a geometric factor

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{c}=\displaystyle \frac{{c}^{2}}{4\pi G}\displaystyle \frac{{D}_{{\rm{S}}}}{{D}_{{\rm{L}}}{D}_{\mathrm{LS}}{(1+{z}_{L})}^{2}},\end{eqnarray} \tag{ 6 }$$

where D_L and D_S are angular diameter distances to the lens and source, D_LS is the angular diameter distance between the lens and source, and the factor of ${(1+{z}_{L})}^{-2}$ arises due to our use of comoving coordinates.

As sources to measure the galaxy–galaxy lensing signal, we use a catalog (Reyes et al. 2012) of 1.2 background galaxies per arcminutes² with weak lensing shears estimated using the re-Gaussianization method (Hirata & Seljak 2003) and photometric redshifts from Zurich Extragalactic Bayesian Redshift Analyzer (ZEBRA; Feldmann et al. 2006). The catalog is characterized in detail in several papers (see Reyes et al. 2012; Mandelbaum et al. 2012, 2013; Nakajima et al. 2012).

The lensing measurement begins with the identification of background galaxies around each lens (with photometric redshift larger than the lens spectroscopic redshift). Inverse variance weights are assigned to each lens-source pair, including both shape noise and measurement error terms in the variance:

$$\begin{eqnarray}&&{w}_{\mathrm{ls}}=\displaystyle \frac{1}{{{\rm{\Sigma }}}_{{\rm{crit}}}^{2}({\sigma }_{e}^{2}+{\sigma }_{{\rm{SN}}}^{2})},\end{eqnarray} \tag{ 7 }$$

where ${\sigma }_{e}^{2}$ is the shape measurement error due to pixel noise, and ${\sigma }_{{\rm{SN}}}^{2}$ is the rms intrinsic ellipticity (both quantities are per component, rather than total; the latter is fixed to 0.365 following Reyes et al. 2012). Use of photometric redshifts, which have nonzero bias and significant scatter, gives rise to a bias in the signals that can be easily corrected using the method from Nakajima et al. (2012). This bias is a function of lens redshift, and is calculated including all weight factors for each lens sample taking into account its redshift distribution. For typical low-redshift lens samples in this work, the bias is of the order of 1%, far below the statistical errors.

ΔΣ in each radial bin can be computed via a summation over lens-source pairs "ls" and random lens-source pairs "rs":

$$\begin{eqnarray}&&{\rm{\Delta }}{\rm{\Sigma }}({r}_{p})=\displaystyle \frac{{\displaystyle \sum }_{\mathrm{ls}}{w}_{\mathrm{ls}}{e}_{t}^{(\mathrm{ls})}{{\rm{\Sigma }}}_{{\rm{crit}}}({z}_{l},{z}_{s})}{2{ \mathcal R }{\displaystyle \sum }_{\mathrm{rs}}{w}_{\mathrm{rs}}},\end{eqnarray} \tag{ 8 }$$

where r_p is the comoving projected radius from the lens, e_t is the tangential ellipticity component of source galaxy with respect to the lens position, ${ \mathcal R }\approx (1-{e}_{{\rm{rms}}}^{2})$ is the shear responsivity (Bernstein & Jarvis 2002) that converts from the ensemble average distortion to shear, and ${e}_{{\rm{rms}}}$ is the rms distortion per component. The division by $\sum {w}_{\mathrm{rs}}$ accounts for the fact that some of our "sources" are physically associated with the lens, and therefore not lensed by it (see, e.g., Sheldon et al. 2004). Finally, we subtract off a similar signal measured around random points with the same area coverage and redshift distribution as the lenses, to subtract off any coherent systematic shear contributions (Mandelbaum et al. 2005b); this signal is statistically consistent with zero for all scales used in this work.

To calculate the error bars, which are dominated by shape noise, we use the jackknife resampling method. The maximum scale used for the fits in the lensing analysis is $1\;{h}^{-1}$ Mpc, which for a typical lens redshift is far below the typical size of each jackknife-resampled region. Thus, the jackknife method is a reasonable approach to getting the covariance matrix for the projected mass profile.

All of the projected mass around lens galaxies contributes to the galaxy–galaxy lensing signal. This includes contributions from the host dark matter halo in which the lens galaxy resides ("one-halo term"), and from other dark matter halos ("two-halo term") that are part of the large-scale structure associated with the lens. For central galaxies, the one-halo term simply corresponds to the ${\rm{\Delta }}{\rm{\Sigma }}$ for the host dark matter halo. For satellite galaxies, there are two contributions to the one-halo term: on small scales, a contribution from the satellite subhalo, and on larger scales (0.3–$2\;{h}^{-1}$ Mpc) a contribution from the host halo itself. See, e.g., Mandelbaum et al. (2005a) for illustrations of these contributions to the galaxy–galaxy lensing signal. The distinctive shape of the satellite contribution in the lensing signal will be a diagnostic of possible contamination of our "central galaxy" sample by satellites in this work.

We fit the observed lensing profiles to the prediction for a pure Navarro et al. (1997, hereafter NFW) profile. The NFW profile in principle has two parameters: the concentration and the total mass within some fiducial radius, for which we use a spherical overdensity of $200{\rho }_{\mathrm{crit}}$. We fix the relationship between concentration c₂₀₀ and mass M_200c using a fitting formula from Diemer & Kravtsov (2015).

For a given mass (and therefore concentration), the ${\rm{\Delta }}{\rm{\Sigma }}$ is predicted via direct integration. Finally, we fit for the mass via ${\chi }^{2}$ minimization using the Levenberg–Marquardt algorithm. This fit is carried out twice for each sample. In the first case, we fit from 0.04 to $0.3\;{h}^{-1}$ Mpc, and compare the best-fit profile with the observed signal out to larger scales, $1\;{h}^{-1}$ Mpc. If the sample is contaminated by satellites, this will be evident in a substantial excess in the observed signal for ${r}_{p}\gt 0.3\;{h}^{-1}$ Mpc compared to the theoretical one fit to ${r}_{p}\lt 0.3\;{h}^{-1}$ Mpc and then extrapolated to larger scales. We use this test as a way to identify samples that are not truly central galaxies and that therefore cannot be used for our analysis (see Section 3). If we do not see any sign of satellite contamination, then we redo the fit using ${\rm{\Delta }}{\rm{\Sigma }}$ from 0.04 to $1\;{h}^{-1}$ Mpc (as in Section 4). As shown directly by Mandelbaum et al. (2005a), using fits to lensing signals from simulated galaxy samples based on N-body simulations, the best-fitting mass will lie between the true median and mean halo mass for the sample, and is most easily interpreted for samples with relatively narrow mass distributions.

To compare our measurements of assembly bias against theoretical expectations, we use three sets of cosmological simulations. The first one is a subset of N-body simulations presented in Diemer & Kravtsov (2014, specifically the L0250, L0500, and L1000 boxes). Each of these simulations follows 1024³ dark matter particles. Combined with the box sizes of 250, 500, and $1000\;{h}^{-1}$ Mpc on a side, the resulting particle mass resolutions are $1.1\times {10}^{9}$, $8.7\times {10}^{9}$, and $7\times {10}^{10}\;{h}^{-1}\;{M}_{\odot }$, respectively. The force resolution lengths are 5.8, 14, and 33 comoving ${h}^{-1}$ kpc, respectively. The cosmological parameters used in these simulations (${{\rm{\Omega }}}_{m}\;=1-{{\rm{\Omega }}}_{{\rm{\Lambda }}}=0.27$, h = 0.7, ${\sigma }_{8}=0.82$, ${M}_{{\rm{nl}}}=3.2\times {10}^{12}\;{h}^{-1}\;{M}_{\odot }$ at z = 0) are close to the ones adopted in the current study. Halos and merger trees were extracted using the Rockstar (Behroozi et al. 2013a) and consistent-trees (Behroozi et al. 2013b) codes. For more details of the simulations, we refer the reader to Diemer & Kravtsov (2014).

The main quantity we wish to derive from this set of simulation data is the redshift of formation ${z}_{{\rm{form}}}$ for the halos, for which we consider two commonly used definitions: (1) z₅₀, the epoch when a halo has acquired 50% of its final mass, and (2) ${z}_{{\rm{mah}}}$, which is obtained from the mass accretion history of halos, following the prescription of Wechsler et al. (2006). We fit the mass accretion history of halos by the form

$$\begin{eqnarray}&&M(z)\propto \mathrm{exp}(-\alpha z),\end{eqnarray} \tag{ 9 }$$

with the Levenberg–Marquardt algorithm for the minimization of the least squares, and define ${z}_{{\rm{mah}}}=2/\alpha -1$. Since ${z}_{{\rm{mah}}}$ is determined by the overall merger history, it is less sensitive to individual events in halo assembly than other definitions of ${z}_{{\rm{form}}}$ (Wechsler et al. 2006). In Section 4.1 we will compare the clustering properties of central galaxies selected by their SFH with that of halos, selected by either definition of ${z}_{{\rm{form}}}$.

The second set of simulations we use is in the form of a mock galaxy catalog, which is made publicly available by Watson et al. (2015) and is produced by applying the age-matching model to the Bolshoi simulation (Klypin et al. 2011). The Bolshoi simulation box size is the same as that of the L0250 simulation mentioned above, but with much better mass resolution ($1.35\times {10}^{8}\;{h}^{-1}\;{M}_{\odot }$). The cosmology adopted by the Bolshoi simulation is the same as that of Diemer & Kravtsov (2014).

To create the mock catalog, Watson et al. (2015) first use the standard abundance matching technique to assign stellar mass to the subhalos (by associating stellar mass with maximum circular velocity a subhalo has ever attained). Then, for subhalos within a given stellar mass range, a sSFR is further assigned to each subhalo (following the observed distribution of sSFR from galaxies of the same stellar mass range), under the simplifying assumption that older subhalos have lower sSFR. Assembly bias is therefore maximally built into the properties of the mock galaxies via a monotonic relationship between age and sSFR. In the resulting mock catalog that contains 207,950 galaxies, we have information such as the position within the simulation box, the designation as a central or satellite, the stellar mass and sSFR, and the host halo mass, for each galaxy. In Section 4.2 we will make use of this catalog and compare the clustering properties of central galaxies selected by sSFR with the predictions from the age-matching model.

As our third set of simulations, we make use of the catalog from the semi-analytic model of Guo et al. (2011), which is built upon the halo and subhalo merger history extracted from the Millennium Run simulation (Springel et al. 2005). The latter follows 2160³ dark matter particles in a ${(500{h}^{-1}\mathrm{Mpc})}^{3}$ box, with the mass and force resolutions of $8.6\times {10}^{8}\;{h}^{-1}\;{M}_{\odot }$ and $5\;{h}^{-1}$ kpc, respectively. The cosmology adopted differs from our default one mainly on ${\sigma }_{8}$: ${{\rm{\Omega }}}_{m}=1-{{\rm{\Omega }}}_{{\rm{\Lambda }}}=0.25$, h = 0.73, ${\sigma }_{8}=0.9$. The semi-analytic model provides, as a function of redshift, properties such as the total halo mass, stellar mass, broadband color, sSFR, and mean mass-weighted stellar age for the galaxies. We have extracted information for central galaxies whose present-day halo mass is ${M}_{200c}\geqslant {10}^{11}\;{h}^{-1}\;{M}_{\odot }$, and computed z₅₀ and ${z}_{{\rm{mah}}}$ for the stellar mass assembly history. In Section 4.1 we will compare the prediction of assembly bias from this model with our observations.

2.4.1. Magnitude of Assembly Bias

It is informative to examine the expected magnitude of assembly bias as a function of halo mass. We have used the suite of simulations of Diemer & Kravtsov (2014) to calculate the ratio of real space correlation functions between the early- and late-forming halos in several halo mass bins. In each mass bin, we designate halos that have ${z}_{{\rm{form}}}$ that is higher (lower) than the mode of the distribution as early-forming (late-forming). For ${z}_{{\rm{form}}}$, we consider both z₅₀ and ${z}_{{\rm{mah}}}$. Note that the way halos are split into early- and late-forming ones here is slightly different from what we will employ when comparing to actual galaxy samples (Section 4). The results are shown in Figure 1. From bottom to top, each row represents a mass bin in order of increasing halo mass. The left (right) column is obtained when ${z}_{{\rm{form}}}={z}_{{\rm{mah}}}$ (z₅₀). For our definition of halo formation time and the halo sample selection, we expect the magnitude of assembly bias to be small—only appreciable at $\mathrm{log}({M}_{200c}/{h}^{-1}\;{M}_{\odot })\lt 13$. This informs our decision to focus on halos around $\mathrm{log}({M}_{200c}/{h}^{-1}\;{M}_{\odot })\approx 12$ in Sections 3 and 4. It is interesting to note that while the amplitude of assembly bias is similar for both definitions of ${z}_{{\rm{form}}}$, the detailed halo mass dependence is somewhat different; we see clearly the sign change at $\mathrm{log}({M}_{200c}/{h}^{-1}\;{M}_{\odot })\approx 14$ when using ${z}_{{\rm{mah}}}$ (that is, younger halos are more strongly clustered; see also Wechsler et al. 2006; Dalal et al. 2008), but no evidence for it with z₅₀ for all the mass bins we have examined, which is consistent with the results of Li et al. (2008).

Zoom In Zoom Out Reset image size

Finally, from Figure 1 we see that the magnitude of assembly bias at large scales appears to depend on scales. This behavior is seen in both the Diemer & Kravtsov (2014) and Millennium simulations. Such scale dependence of the assembly bias potentially has very significant implications for cosmological constraints using large-scale clustering and will be investigated in a future publication.

2.4.2. Robustness of Halo Mass Inference from Lensing Measurements

The most important finding from Section 3 is that, after identifying samples with the same Y07 halo mass estimates and splitting by either the SFH or sSFR, the mean halo mass of the resulting early- and late-forming samples would actually be quite different. This makes such a procedure an invalid way to test for assembly bias. This then motivates us to develop a two-step approach: we can start with two central galaxy samples (denoted as samples 1 and 2) that are selected by certain halo mass proxy, according to which the mean masses ${M}_{1}\lt {M}_{2}$. Referring to the mean halo masses of the resulting early- and late-forming subsamples after splitting the samples by either the SFH or sSFR as M_1e, ${M}_{1l}$, M_2e, and ${M}_{2l}$ (so that in general ${M}_{1e}\gt {M}_{1l}$ and ${M}_{2e}\gt {M}_{2l}$), which are then measured by weak lensing, we can achieve ${M}_{1e}\approx {M}_{2l}$ (e.g., within $1\sigma$) by adjusting the chosen range of mass proxy of sample 2 (M₂) relative to that of sample 1 (M₁).

It is also found in Section 3 that there are non-negligible satellite contamination in the Y07 central galaxy sample. As these satellites are typically residing in halos more massive than the ones that host our central galaxies, their large-scale bias reflects that of their hosts, and therefore the inclusion of these satellites would dilute any assembly bias signal we are after. Below we describe details of the steps that lead to our final samples.

4.1.1. Sample Construction

Our starting point is a halo mass proxy that can help define samples 1 and 2. For the SFH-selected samples, we choose to use the central galaxy stellar mass–halo mass relationship for this purpose, and adopt the measurements by More et al. (2011, hereafter M11), which are done separately for red and blue centrals (see the solid curves in Figure 4), using satellite kinematics. Using these relations as an approximate guideline, we can select central galaxies in certain stellar mass ranges so that they might live in halos of similar masses. As suggested by Figure 4, we have to combine blue centrals of higher stellar mass with red ones of lower stellar mass.¹⁴ The VESPA-based SFH is then used to classify the galaxies into early- and late-forming ones, for which mean halo masses are measured by galaxy–galaxy lensing.

Zoom In Zoom Out Reset image size

In implementing the above procedure, we follow M11 and estimate the stellar mass via

$$\begin{eqnarray}\mathrm{log}\left(\displaystyle \frac{{M}_{{\rm{star}}}}{{h}^{-2}\;{M}_{\odot }}\right) & = & -0.406+1.097(g-r)\\ & & -0.4({M}_{r}-5\mathrm{log}h-4.64)\end{eqnarray} \tag{ 10 }$$

(see also Bell et al. 2003). Here M_r is the SDSS r-band absolute magnitude, $(g-r)$ is estimated in the rest-frame; for both of these we use the SDSS Petrosian magnitudes with k-corrections taken from the NYU value-added galaxy catalog (Blanton et al. 2005). After examination of our galaxies in the color-stellar mass space, we adopt a red/blue division line that is somewhat different from that used in M11:

$$\begin{eqnarray}&&{(g-r)}_{{\rm{div}}}=0.67+0.03\mathrm{log}[{M}_{{\rm{star}}}/({10}^{10}{h}^{-2}\;{M}_{\odot })].\end{eqnarray} \tag{ 11 }$$

The next step in our procedure is the removal of satellites. We have found that satellite contamination is more serious in red galaxy samples (i.e., more enhanced one-halo term in the galaxy-clustering signal). Thus, in the above process, after a red central sample is selected, we apply a FOF algorithm (see Jian et al. 2014 for a description of the code) to identify groups of galaxies, using comoving linking lengths along and perpendicular to the line of sight of 1 and 0.08 (in units of the mean galaxy separation)¹⁵ , respectively. We only keep the most massive member in the groups identified. During this step, we consider not only the red galaxy sample in question, but also other red centrals in the Y07 central catalog, as well as blue centrals in halos with $\mathrm{log}({M}_{200b}/{h}^{-1}\;{M}_{\odot })\gt 12.5$ (according to the M11 relations). Then, after combining the blue central sample with the satellite-trimmed red central sample, the FOF algorithm is applied again to further remove remaining satellites. Only after these steps do we separate the sample into early- and late-forming ones. About 7%–15% of galaxies are removed this way.

It should also be noted that not every galaxy can be unambiguously identified as early- or late-forming from the VESPA SFH, as we have opted for purity and demanded consistency between the two dust models (Section 2.1). The early- and late-forming sample sizes would therefore be further reduced (by 5%–17%) compared to the parent samples.

We have thus obtained a pair of early- and late-forming samples that we believe are dominated by central galaxies in their halos and have similar mean halo masses. The way they are constructed is depicted in Figure 4. In the figure, two horizontal bands (shaded regions) are shown around $\mathrm{log}({M}_{{\rm{200b}}}/{h}^{-1}\;{M}_{\odot })\approx 12.4$. Each band denotes the initial halo mass range (M₁ and M₂ in the notation used above) guessed from the central stellar mass–halo mass relations. The green one represents sample 2 from which a late-forming sample is derived. Similarly, the yellow band represents sample 1 from which an early-forming sample is obtained. More specifically, the log stellar mass ranges for which the red and blue galaxies are selected are 9.9–10.2 and 10.2–10.45 (10.35–10.6 and 10.52–11.1) for sample 1 (2). These ranges are delineated by the vertical dashed and dotted lines. The resulting early- and late-forming samples contain 18,200 and 26,071 galaxies, respectively. Using galaxy–galaxy lensing, it is found that these samples have mean masses (M_200c) that are statistically consistent: ${M}_{1e}=({9.5}_{-2.0}^{+2.5})\times {10}^{11}\;{h}^{-1}\;{M}_{\odot }$ and ${M}_{2l}=({8.4}_{-1.8}^{+2.2})\times {10}^{11}\;{h}^{-1}\;{M}_{\odot }$ (see the bottom panel of Figure 5). Below we refer to these as our SFH samples.

Zoom In Zoom Out Reset image size

In the top and middle panels of Figure 5, we show the comparison of the correlation functions for these samples. It is clear that the power at small scales (e.g., $\lesssim 1\;{h}^{-1}$ Mpc) is much reduced compared to Figures 2 and 3, due to the removal of satellites with our FOF procedure. (The lensing signal shown in the bottom panel matches the NFW fit for ${r}_{p}\lesssim 2\;{h}^{-1}$ Mpc, which further supports the conclusion that we do not have an appreciable satellite population in the final galaxy samples.) At large scales (${r}_{p}\gtrsim 5\;{h}^{-1}$ Mpc), we see in the middle panel the (square of) relative bias scatters around unity (indicated by the horizontal dotted line), suggesting that there is not a strong difference in the large-scale clustering of the two samples. Over the scales $5\mbox{--}35\;{h}^{-1}$ Mpc, the square root of the mean ratio is 1.00 ± 0.12.

4.1.2. Comparison with Theoretical Expectations

To properly interpret our findings, we first compare the observed relative bias with that obtained from the z = 0 output of the L0250 simulation (Section 2.4). Note that observationally we can only robustly measure the mean halo mass of our galaxy samples, but not the form of their halo mass distribution. To proceed, we thus make the assumption that the distribution of halo mass of the observed galaxy samples follows a log-normal form,¹⁶ with parameters ${M}_{{\rm{cen}}}$ and ${\sigma }_{\mathrm{log}M}$ representing the mean and standard deviation of the Gaussian in log space. To find out the values of ${M}_{{\rm{cen}}}$ and ${\sigma }_{\mathrm{log}M}$ appropriate for our galaxy samples, we consider a grid of combinations of these parameters, with $\mathrm{log}{M}_{{\rm{cen}}}$ ranging from 11.4 to 12.5 (with an interval of 0.1 dex), and ${\sigma }_{\mathrm{log}M}$ from 0.05 to 0.45 (with an interval of 0.05). Given an observed galaxy sample, for each combination of $({M}_{{\rm{cen}}},{\sigma }_{\mathrm{log}M}$), we predict the corresponding ${\rm{\Delta }}{\rm{\Sigma }}$ signal and obtain ${\chi }^{2}$ by comparing with the observed profile. We consider all models that satisfy ${\chi }^{2}\leqslant {\chi }_{{\rm{min}}}^{2}+2.3$ as plausible (corresponding to the 68% probability interval for two parameters¹⁷ ), where ${\chi }_{{\rm{min}}}^{2}$ is given by the model with minimum ${\chi }^{2}$ on the grid. To further split the halo samples into early- and late-forming ones, we proceed as follows. Given the median redshift ${z}_{{\rm{med}}}$ of an observed galaxy sample, we compute the redshift ${z}_{{\rm{div}}}$ corresponding in lookback time 9 Gyr prior to ${z}_{{\rm{med}}}$, and use it as a proxy for the boundary of the first temporal bin in VESPA. That is, an early-forming (late-forming) galaxy sample typically has formed 50% of its final stellar mass before (after) ${z}_{{\rm{div}}}$. ${z}_{{\rm{div}}}$ for our early- and late-forming galaxy samples are 1.84 and 2.11, respectively. For each of the models that can represent an observed early-forming (late-forming) galaxy sample, we select halos with ${z}_{{\rm{form}}}\gt {z}_{{\rm{div}}}$ (${z}_{{\rm{form}}}\lt {z}_{{\rm{div}}}$) from the simulation and compute the correlation function. Our theoretical expectation is then the average over the correlation functions from all these models. The uncertainties are estimated from jackknife resampling after splitting the simulation box into 125 sub-cubes. We note that the acceptable models show a degeneracy between ${M}_{{\rm{cen}}}$ and ${\sigma }_{\mathrm{log}M}$ (in the sense that the lower the mass, the higher ${\sigma }_{\mathrm{log}M}$), and have ${M}_{{\rm{cen}}}$ ranging from 11.6 to 12.0, and ${\sigma }_{\mathrm{log}M}$ from 0.05 to 0.45, and the spread in correlation functions from these models is much smaller (5%–6%) than the mean values. The fraction of early-forming halos among the halos used in the models is 0.37, fairly close to the observed sample (0.41).

In the top panel of Figure 5, we show as a magenta (cyan) short dashed curve the projected correlation function from the early-forming (late-forming) halo sample. Here we have adopted ${z}_{{\rm{mah}}}$ for ${z}_{{\rm{form}}}$, but using z₅₀ leads to similar conclusions. The ratio ${b}_{{\rm{th}}}^{2}({r}_{p})$ of the early-forming-to-late-forming halo correlation functions from the simulation is shown as the dashed curve in the middle panel of the Figure, with the (pink) shaded region enclosing the $1\sigma$ uncertainties from the model. Using the covariance matrices built from the ratio of the early- and late-forming projected correlation functions and their associated jackknife samples for both the observations and models (O and T, respectively), we infer the probability of the observed and theoretical relative bias squared to be consistent with each other from

$$\begin{eqnarray}{\chi }^{2} & = & \displaystyle \sum _{{ij}}({w}_{e}({r}_{i})/{w}_{l}({r}_{i})-{b}_{{\rm{th}}}^{2}({r}_{i})){D}_{{ij}}^{-1}\\ & & \times ({w}_{e}({r}_{j})/{w}_{l}({r}_{j})-{b}_{{\rm{th}}}^{2}({r}_{j})),\end{eqnarray} \tag{ 12 }$$

where w_e and w_l are the observed early- and late-forming projected correlation functions, and $D=O+T$. For notational simplicity, we have omitted the subscript p in w_p and r_p. For these calculations, we consider only r_p bins in the $5\mbox{--}35\;{h}^{-1}$ Mpc range. We find that the theory and our observation are inconsistent. Given 5 degrees of freedom, the probability to yield a ${\chi }^{2}$ as large as observed (27.3) if the data and theory were drawn from the same distribution is only $p=5.0\times {10}^{-5}$. The corresponding ${\chi }^{2}$ when z₅₀ is adopted is 64.1, corresponding to a probability well below 10⁻¹⁰.

We next turn to the comparison with the predictions from the Guo et al. (2011) semi-analytic model. Analogous to the approach described above, for a given central galaxy sample, we again start with simulated halo samples permitted by the observed ${\rm{\Delta }}{\rm{\Sigma }}$ profile, then select early-forming (late-forming) halos as those hosting central galaxies with ${z}_{{\rm{form}}}\gt {z}_{{\rm{div}}}$ (${z}_{{\rm{form}}}\lt {z}_{{\rm{div}}}$). Here ${z}_{{\rm{form}}}={z}_{{\rm{mah}}}$ is calculated from the stellar mass assembly history from the Guo et al. (2011) model. The uncertainties in the model predictions are also calculated with jackknife resampling of the simulation box. The projected correlation functions for the early- and late-forming models are shown as orange and purple long dashed curves in the top panel of Figure 5, respectively, and the ratio of the two is shown as the long dashed curve in the middle panel (the orange shaded region represents the $1\sigma$ uncertainties). Using Equation (12), we find that the probability that the model and data are consistent is p = 0.019 with ${\chi }^{2}=13.5$. We note, however, that the fraction of early-forming halos among the halos used in the models is only 0.04, far below that of the observed sample. Given that this model also produces a sSFR distribution that is not consistent with the observed sSFR, we will not further consider it in Section 4.2.

It is possible that at mass scales lower than ${M}_{{\rm{nl}}}$ (as is the case for our samples), the origin of the assembly bias is partially due to the so-called "backsplash" halos (Dalal et al. 2008), those that have been accreted onto massive halos but are on highly elongated orbits and thus would have spent a substantial amount of time outside of the parent halo's virial radius (Wetzel et al. 2014). The member galaxies of these backsplash halos would appear old because they have been affected by dense environments of the massive halos, and the bias of these small halos would be that of their massive parent halos. Our treatment of satellite removal may exclude some of such halos from our sample, and therefore inadvertently suppress the effect of assembly bias. Unless there are clever ways to observationally distinguish the backsplash halos from low-mass halos that are unrelated to nearby, massive halos, it seems to be a challenge to analyses like ours to only remove satellites that are bounded within massive halos, but not backsplash galaxies. Perhaps a better (and more feasible) approach is to mimic the effect of satellite removal in simulated data. Although a full-blown analysis is beyond the scope of the present paper, in Section 5.2 we will use a simple method to roughly estimate the effect of possible exclusion of backsplash halos to the non-detection of assembly bias.

Ideally, we would like to explore the assembly bias with halos over a wide mass range, similar to what is done in Yang et al. (2006). Unfortunately, at the low-mass end, we are limited by the inability of VESPA to yield sufficient number of early-forming galaxies (likely due to the quality of SDSS spectra), thus rendering the lensing measurement too noisy, while at higher mass end, both the fact that the number of central galaxies decreases precipitously and that most of them are quite old (beyond the temporal resolution offered by VESPA) make it more difficult to apply our method.

We next follow a procedure similar to that outlined in Section 4.1 for samples selected by sSFR. As we now need to consider high- and low-sSFR centrals separately, operationally it is much more convenient to use the mass estimates from the Y07 catalog, rather than using the M11 relations, in selecting samples 1 and 2. In constructing the galaxy samples, we again use the FOF algorithm to remove contaminations from satellite galaxies.

We have thus arrived at a pair of low- and high-sSFR samples with similar mean masses. They are selected by $\mathrm{log}({M}_{200c}/{M}_{\odot })=12.0\mbox{--}12.5$ and ${\rm{sSFR}}\lt {10}^{-11.8}$ yr⁻¹, and $\mathrm{log}({M}_{200c}/{M}_{\odot })=12.75\mbox{--}13.10$ and ${\rm{sSFR}}\gt {10}^{-11.8}$ yr⁻¹, respectively. After the FOF step that removes about 10% of the galaxies, the two samples contain 25,838 and 29,659 galaxies. The lensing measurements are shown in the bottom panel of Figure 6, giving the mean masses of ${M}_{200c}\;=({1.39}_{-0.21}^{+0.24})\times {10}^{12}\;{h}^{-1}\;{M}_{\odot }$ and $({1.26}_{-0.20}^{+0.24})\times {10}^{12}\;{h}^{-1}\;{M}_{\odot }$. The correlation functions and the ratio of low-to-high-sSFR correlation functions are shown in the top and middle panels. The square root of the mean ratio (over $5\mbox{--}35\;{h}^{-1}$ Mpc) is 1.07 ± 0.14. Given that the data points are all above unity, we use a procedure analogous to that described in Section 4.1 to examine the consistency between the two samples by calculating

$$\begin{eqnarray}&&{\chi }^{2}=\displaystyle \sum _{{ij}}({w}_{\mathrm{ls}}({r}_{i})-{w}_{{hs}}({r}_{i})){O}_{{ij}}^{-1}({w}_{\mathrm{ls}}({r}_{j})-{w}_{{hs}}({r}_{j})),\end{eqnarray} \tag{ 13 }$$

where w_ls and w_hs are the projected correlation function of the low-sSFR and high-sSFR samples, and O the covariance matrix built from the difference between w_ls and w_hs and their associated jackknife samples. With ${\chi }^{2}=5.5$ from 5 degrees of freedom, we find that the two samples have a 36% probability to be consistent.

Zoom In Zoom Out Reset image size

In the top and middle panels of the Figure, we compare our measurements with predictions from the age-matching model (Watson et al. 2015; see Section 2.4). Following the procedure outlined in Section 4.1, we construct halo samples that can represent the observed galaxy samples by approximating the halo mass distribution as log-normal, and considering halo samples that provide good fits to the observed ${\rm{\Delta }}{\rm{\Sigma }}$ profile. We adopt the same division criterion for low- and high-sSFR as in the real data (${\rm{sSFR}}={10}^{-11.8}$ yr⁻¹). From the curve in the middle panel (representing the ratio of the two model curves shown in the top panel), we see that the age-matching model prediction for ${b}_{{\rm{th}}}^{2}$ is ∼1.6. We find $({\chi }^{2},p)=(10.5,0.033)$, indicating the model is inconsistent with the observation.

As the age-matching model is calibrated against a SDSS galaxy sample that is slightly different from the one we use (Watson et al. 2015), it is possible the sSFR distribution in the model does not match perfectly with our data. After inspecting the cumulative sSFR distributions in both the real data and the model, we find it more appropriate to set the sSFR division to be ${10}^{-11.4}$ yr⁻¹ in the model for comparison. With this adjustment, the fraction of early-forming halos in the model is 0.41, which is not far from the observed value of 0.47. The results become $({\chi }^{2},p)=(20.5,4.0\times {10}^{-4})$. Therefore, our observations do not seem to be compatible with the age-matching model, in which the magnitude of assembly bias is expected to be maximal (Section 2.4). Such a result may not be surprising, given that the central galaxy stellar mass–halo mass relationships based on the age-matching model (particularly that for blue galaxies) do not match well with the real data, as indicated by lensing measurements (Mandelbaum et al. 2015).

A proper definition of both halo and galaxy formation time is critical, as it affects the expected amplitude of assembly bias as a function of halo mass and, more importantly, it determines whether one can faithfully link the observed galaxy population to the underlying dark matter halos. As we have seen before, the relative bias between early- and late-forming halos is a function of halo mass, and the detailed mass dependence actually depends on whether z₅₀ or ${z}_{{\rm{mah}}}$ is adopted for ${z}_{{\rm{form}}}$ (Figure 1). This in turn affects any comparison with observations, especially the inference of the statistical significance (see Section 4.1).

How can we best link the observed properties of galaxies to the formation history of the host dark matter halos? Various groups have adopted different ways to define the formation time for the central galaxies, including the current sSFR (or equivalently, broadband optical color) and luminosity-weighted mean age of the stellar populations (Section 1). One of the new aspects of our study is the use of resolved SFH from the VESPA algorithm. To gain insight into the answer to the above question, we make use of the results from the semi-analytic model of Guo et al. (2011). We have extracted information from ∼157,000 central galaxies whose present-day halo mass is ${M}_{200c}=(1\mbox{--}2)\times {10}^{12}\;{h}^{-1}\;{M}_{\odot }$, and computed z₅₀ and ${z}_{{\rm{mah}}}$ for both the total mass and the stellar mass assembly history (in the following denoted with a subscript "t" and "s," respectively).

After examining the correlations between the various galactic properties (including ${z}_{{\rm{50,s}}}$ and ${z}_{{\rm{mah,s}}}$) and the halo formation time (${z}_{{\rm{50,t}}}$, ${z}_{{\rm{mah,t}}}$), we have found that the best correlation is between ${z}_{{\rm{mah,s}}}$ and ${z}_{{\rm{mah,t}}}$ (with a Pearson correlation coefficient r = 0.83), followed by that between ${z}_{{\rm{50,s}}}$ and ${z}_{{\rm{50,t}}}$ (r = 0.54), age and ${z}_{{\rm{50,t}}}$ (r = 0.49), and ${z}_{{\rm{50,s}}}$ and ${z}_{{\rm{mah,t}}}$ (r = 0.36). Therefore, if the Guo et al. (2011) model is a good approximation to the real galaxy populations (see, e.g., Lin et al. 2013 for the agreement between the model prediction and observation for the stellar mass assembly history of brightest cluster galaxies at high redshift), with a suitable choice of formation time indicators, it is possible to infer the halo formation history from that of the central galaxies for the low-mass halos we consider here. Furthermore, if ${z}_{{\rm{50,s}}}$ derived from VESPA is representative of the true value on average, our results presented in Section 4.1 are on solid footing.

Ideally, we would like to calculate ${z}_{{\rm{mah,s}}}$ for the observed galaxy populations; however, we refrain from deriving it from the VESPA-based SFH here as the binning in lookback time for the public VESPA data is not optimal for this purpose. Furthermore, it remains to be seen if VESPA-based SFH would result in a unbiased ${z}_{{\rm{mah}}}$. In future work we plan to also investigate the use of other formation time indicators such as the strength of the 4000 Å break and the luminosity-weighted mean age for the study of assembly bias.

In principle, the signature of assembly bias would be stronger if one uses the extrema of the distribution. Due to the temporal resolution of VESPA SFH, in our analysis, we approximately designate galaxies that have ${z}_{{\rm{50,s}}}\gtrsim 1.8$ as early-forming, and those with ${z}_{{\rm{50,s}}}\lesssim 2.1$ as late-forming (Section 4.1). Should we have a higher resolution SFH, we could have examined the distribution of ${z}_{{\rm{50,s}}}$ (or ${z}_{{\rm{mah,s}}}$), and only used the earliest- and latest-forming 20% for the clustering and lensing measurements. Such an analysis would require both much better quality spectra and much larger sample size (or much deeper imaging data than SDSS), however, and is therefore currently not yet feasible.

One potential concern regarding our way of constructing pure samples of central galaxies is the FOF removal of satellites. While we believe this is a necessary operation, some of the galaxies thus removed may be those in the so-called backsplash halos, which may partially contribute to the signal of assembly bias at low-mass scales (Wang et al. 2007; Dalal et al. 2008; see Section 4.1). If true, it may complicate the interpretation of our results (in other words, the non-detection of assembly bias may be due in part to the removal of these galaxies).

Here we attempt to evaluate the effect of satellite removal on the magnitude of assembly bias using simulated data. We repeat the procedure of Section 4.1 and consider halo samples that have a mass distribution following the log-normal form. The acceptable combinations of (${M}_{{\rm{cen}}},{\sigma }_{\mathrm{log}M}$) are again those that give ${\chi }^{2}\leqslant {\chi }_{{\rm{min}}}^{2}+2.3$ for a given observed galaxy sample. For a halo sample constructed from a given parameter set of (${M}_{{\rm{cen}}},{\sigma }_{\mathrm{log}M}$), we remove halos that are located within 2r_200c from any of the more massive halos in the whole simulation box (that is, not restricted to those that satisfy the log-normal mass distribution). For the remaining "isolated" halos, which may correspond crudely to our satellite-trimmed galaxy sample, we apply the ${z}_{{\rm{form}}}\gt {z}_{{\rm{div}}}$ (${z}_{{\rm{form}}}\lt {z}_{{\rm{div}}}$) criterion as before to further filter the halos, depending on whether the observed sample in question is early- or late-forming. Finally, the theoretical expectation is obtained by averaging over the correlation functions from all acceptable models. For the L0250 simulation used in Section 4.1, this procedure typically removes about 8% of the halos, and reduces the ratio of the resulting correlation functions (early to late) by $\sim 10\%$ compared to that before the removal. Such a reduction would make theoretical predictions more compatible with our observations. We find that the probability that the theory and the observed data are drawn from the same distribution is p = 0.055 (with ${z}_{{\rm{form}}}={z}_{{\rm{mah}}}$). For the age-matching halo sample (Section 4.2), about 30% of the halos are removed, resulting in a larger reduction in the ratio of correlation functions, and the probability for the model to be consistent with the observations is p = 0.0041.

Admittedly, this procedure is rather crude, but it should capture the essence of the effect of backsplash halos. We see that indeed such halos may contribute partly to the assembly bias, which is consistent with the findings of Wang et al. (2009).

Another potential concern is the form and width of the mass distribution of the real galaxy samples. When matching an early-forming sample with a late-forming one, we only require their mean halo masses to be consistent (e.g., within $1\sigma$), but not the distribution in mass. Our assumption that the halo mass distribution follows the log-normal form is informed by the age-matching model (Section 4.1), and with this form the centroid and width of the distribution may be roughly determined (as constrained by the observed ${\rm{\Delta }}{\rm{\Sigma }}$ profile). Although we have limited the mass range during the initial sample selection to be $\approx 0.5$ dex or smaller (based on the central stellar mass–halo mass relation, or the mass given in the Y07 catalog), we note these proxies are derived with the presence of satellites, and thus may bias the actual mass spread in an unpredictable way. A possible way forward is to construct (satellite-trimmed) central–halo relations with SFH or sSFR-selected samples (rather than those selected by stellar mass or luminosity thresholds, as commonly used), and use them to better constrain $({M}_{{\rm{cen}}},{\sigma }_{\mathrm{log}M})$ of the galaxy samples in question.

Finally, when we fit the NFW profile to the lensing measurements, we have assumed that the concentration is only a function of halo mass, and ignored any dependence of concentration on the halo formation time. This is mainly because the signal-to-noise ratio of the galaxy–galaxy lensing profiles is insufficient to fit for concentration as well. Such an assumption has the effect of potentially biasing the masses of our galaxy samples, in the sense that the mass of the early-forming (late-forming) sample would be overestimated (underestimated) when the mean concentration is used. From the L0250 simulation, we have found that the concentration is $\sim 60\%$ higher in early-forming halos than in late-forming ones. Adopting such a difference in our lensing measurements leads to roughly a 10% reduction (increase) in halo mass for the early-forming (late-forming) galaxy sample. Given that the difference in mass as presented in Section 4 (i.e., when the mean concentration is used), as well as our mass measurement uncertainties, are both also at similar levels (10%–20%), and that observationally the difference in concentration of halos hosting red and blue central galaxies is much smaller than the 60% value used above (Mandelbaum et al. 2015), we conclude that assuming concentration is only a function of halo mass does not have an appreciable effect in our analysis.

In Sections 3 and 4 we have used the halo mass estimates from Y07 and M11, which are based on an approach similar to abundance matching, and on satellite kinematics, respectively, to guide our initial sample selection. Although we have not carried out a systematic comparison of these estimates with weak lensing (see Mandelbaum et al. 2015 for such an effort), during our two-step procedure of constructing galaxy samples, we have built three samples (selected only by color) as a by-product that can test the M11 relations. It appears that the red and blue central stellar mass–halo mass relations from M11 give rise to consistent halo masses, although these masses are $\sim 40\%$ higher than that indicated by lensing. This is consistent with the finding of Kravtsov et al. (2014).

As for the mass provided by Y07, which is obtained in a fashion similar to the abundance matching technique, our lensing measurements in Sections 3 and 4.2 imply that the scatter in their mass estimates somehow correlates with basic physical properties of galaxies such as SFH or sSFR, at least at the low halo mass regime we study here. Therefore, any study that assumes or requires the scatter in halo mass to be random should be cautious when adopting the mass estimates from this catalog.

As we have demonstrated in Section 2.4, together with numerous previous studies (e.g., Gao et al. 2005; Wechsler et al. 2006; Jing et al. 2007; Li et al. 2008), the halo assembly bias is a marked feature of the CDM model, especially at the low-mass scales we study ($\sim {10}^{12}\;{h}^{-1}\;{M}_{\odot }$). How can this be reconciled with our lack of detection in galaxy populations? It is possible that the baryonic processes of galaxy formation have rendered the signal small, the SFH derived by VESPA is too noisy so that the signal is washed out (for the case of Section 4.1), the sSFR measurements from SDSS are too noisy to be a good indicator of the halo formation epoch (for the case of Section 4.2), or we have not yet found a galactic property that is closely linked to the halo formation history.

To check the first possibility, we again make use of the galaxy catalog from the Guo et al. (2011) model. For central galaxies living in halos of mass ${M}_{200c}=(1\mbox{--}2)\times {10}^{12}\;{h}^{-1}\;{M}_{\odot }$, we compute the mean of ratio ${\xi }_{{\rm{early}}}/{\xi }_{{\rm{late}}}$ over $5\mbox{--}20\;{h}^{-1}$ Mpc, where the early- and late-forming samples are defined by several different conditions: (1) by the halo ${z}_{{\rm{mah,t}}}$, (2) by the stellar ${z}_{{\rm{mah,s}}}$, (3) by the stellar ${z}_{{\rm{50,s}}}$, and (4) by the stellar age. For each of the formation time indicators, we designate the galaxies with value higher (lower) than the mode of the distribution to be early-forming (late-forming). Comparing the ratio from cases (2)–(4) with that of case (1) would inform us whether or not the signature of assembly bias is erased with (this particular model of) galaxy formation. It is found that the ratio for cases when the formation history is inferred from galactic properties is similar to, and not smaller than, that derived purely from the dark matter halo assembly history, and thus it is probable that galaxy formation preserves assembly bias (see Wang et al. 2013).

It is certainly possible that the VESPA-based SFH is too noisy. Tests with mock galaxy spectra carried out in Tojeiro et al. (2009) indicate that for simple SFHs (exponential decay or dual-burst), the recovery of SFH is satisfactory. The SFH of our central galaxies inevitably would be much more complicated than the simple cases tested above, and whether VESPA can reliably decipher the SFH needs to be checked with other tools such as STARLIGHT (Cid Fernandes et al. 2005), FAST (Kriek et al. 2009), MAGPHYS (da Cunha et al. 2008), and FIREFLY (Wilkinson et al. 2015).

In Section 2.4 we have noted that in the Watson et al. (2015) model, the sSFR of a central galaxy is assumed to have a one-to-one correspondence with the formation epoch of the host halo. The assembly bias is thus maximally built into this model. The fact that the probability for the model and our observations to be consistent is at most only at the percent level (Section 5.2) suggests that either intrinsic scatters in the sSFR-formation epoch correspondence are much larger than assumed in the age-matching model (thus rendering the effect of assembly bias too small to be detectable), or the sSFR measurements from SDSS are not adequate for picking up the assembly bias signal, or a combination of both.

Finally, we entertain the possibility that we have yet to employ a new galactic property to better separate galactic systems into early- and late-forming ones before we can unambiguously identify the assembly bias signature in the universe. As discussed in Section 5.1, the best candidate appears to be ${z}_{{\rm{mah,s}}}$, followed by ${z}_{{\rm{50,s}}}$ and the mean age (see also the proxy discussed in Lim et al. 2016). These in principle could be obtained from high-quality spectra with VESPA, as well as the aforementioned codes, and will be subjects of our future investigation. Inspiration could also come from the analysis of Miyatake et al. (2015), who have recently claimed a detection of strong assembly bias using galaxy clusters. Their proxy for halo formation time is $\langle {R}_{{\rm{mem}}}\rangle$, the mean projected separation of member galaxies from the cluster center. In principle, it is possible to adopt a similar proxy for galaxy scale halos, although the number of satellites is much smaller, and the membership determination is less certain.

We conclude by noting that it is also imperative to better understand theoretically the origin of assembly bias across the mass spectrum of halos. If the cause at low-mass scales is the backsplash halos, then instead of identifying the best formation time indicator, one should look for ways to observationally distinguish such galactic systems, or to take such population fully into account when comparing with theoretical models. Regarding the high-mass end, since it is likely that the physical origin of assembly bias in high-mass halos is believed to be quite different than for low-mass halos (Dalal et al. 2008), the detection presented by Miyatake et al. (2015) is not necessarily inconsistent with the non-detection reported here. It is clear that the assembly bias phenomenon is far richer and complicated than expected, and warrants further investigations, such as exploring better proxies for halo formation time, and the scale dependence of assembly bias (Section 2.4.1).