Halo Assembly Bias in Hierarchical Structure Formation
Neal Dalal, Martin White, J. Richard Bond, Alexander Shirokov
aa r X i v : . [ a s t r o - ph ] M a r Draft version October 22, 2018
Preprint typeset using L A TEX style emulateapj v. 08/22/09
HALO ASSEMBLY BIAS IN HIERARCHICAL STRUCTURE FORMATION
Neal Dalal , Martin White , J. Richard Bond and Alexander Shirokov Draft version October 22, 2018
ABSTRACTWe investigate the origin of halo assembly bias, the dependence of halo clustering on assemblyhistory. We relate halo assembly to peak properties measured in the Lagrangian space of the initiallinear Gaussian random density field, and show how these same Lagrangian properties determinelarge-scale bias. We focus on the two regimes where assembly bias has been observed to be significant:at masses very large and very small compared to the nonlinear mass scale. At high masses, we showthat assembly bias is expected from the statistics of the peaks of Gaussian random fluctuations, andwe show that the extent of assembly bias found in N-body simulations of rare halos is in excellentagreement with our theoretical prediction. At low masses, we argue that assembly bias largely arisesfrom a sub-population of low mass halos whose mass accretion has ceased. Due to their arresteddevelopment, these halos naturally become unbiased, in contrast to their anti-biased peers. We showthat a simple toy model incorporating these effects can roughly reproduce the bias trends found inN-body simulations.
Subject headings: cosmology:theory – methods:numerical – dark matter: merging histories – galaxies:clusters INTRODUCTION
Hierarchical structure formation is a basic predictionof the highly successful cold dark matter cosmologicalmodel, and a key element of hierarchical growth is theformation of bound virialized objects called dark mat-ter halos. Dark matter halos are not distributed uni-formly throughout the universe, but rather compose abeaded filamentary network, sometimes referred to asthe cosmic web (Bond et al. 1996). The clustering ofthe rare, massive dark matter halos is enhanced rel-ative to the general mass distribution (Kaiser 1984;Efstathiou et al. 1988; Cole & Kaiser 1989; Bond et al.1991; Mo & White 1996; Sheth & Tormen 1999), an ef-fect known as bias. The more massive the halo, thelarger the bias. As a result, the mass of halos hostinga given population of objects is sometimes inferred bymeasuring their degree of clustering – allowing a statis-tical route to the notoriously difficult problem of mea-suring masses of cosmological objects. There are manyapplications for statistical mass determination, includ-ing (for example) cluster mass-observable self-calibration(e.g. Lima & Hu 2004; Wu et al. 2008) or reconstructionof the halo occupation distribution of galaxies or quasars(e.g. Cooray & Sheth 2002).Since halos of a given mass can differ in their forma-tion history, properties and large-scale environment , anatural question arises: do these details affect halo clus-tering? Numerical simulations are now able to producesamples with sufficient statistics to test for the depen-dence of clustering on history, properties and environ-ment and recent numerical simulations reveal subtle de-pendences of halo clustering on formation redshift, in- Canadian Institute for Theoretical Astrophysics, University ofToronto, 60 St. George St., Toronto, ON, M5S 3H8, Canada Department of Astronomy, University of California, Berkeley,CA 94720 The large-scale environment of a halo refers to the density,smoothed on some suitably large scale, e.g. 8 h − Mpc. ternal structure and spin (Gao, Springel & White 2005;Wechsler et al. 2006; Harker et al. 2006; Bett et al. 2007;Wetzel et al. 2007; Jing, Suto & Mo 2007; Gao & White2007; Angulo, Baugh & Lacey 2008)Using a marked correlation function, Sheth & Tormen(2004) found that close halo pairs tend to have earlierformation times than more distant pairs, work whichwas extended and confirmed by Harker et al. (2006).Gao, Springel & White (2005) found that later forming,low-mass halos are less clustered than typical halos of thesame mass at the present. Wechsler et al. (2006) andWetzel et al. (2007) found a similar dependence uponhalo formation time (defined in a slightly different way),showing that the trend reversed for more massive halosand that the clustering depended on halo concentration.Recently, Gao & White (2007) showed that a wide vari-ety of halo properties correlate with bias, including spinand substructure content, while Li et al. (2008) showedthat different definitions of halo age, when applied to thesame halo samples, give different amounts of assemblybias.What remains unclear from this numerical work isthe physical origin of this effect. Gao, Springel & White(2005) argued that assembly bias is inconsistent with thePress & Schechter (1974) model, which explicitly pre-dicts that halo abundance (and hence clustering) de-pends only upon halo mass and no other internal prop-erties related to assembly history. Attempts to explainassembly bias based on the Press-Schechter approach(i.e. the excursion set using a sharp k -space filter) wouldtherefore appear doomed from the outset, however sev-eral workers have explored modifications to the Press-Schechter model in the hopes of accounting for assem-bly bias. Sandvik et al. (2007) considered excursion setmass functions using barrier heights determined from theellipsoidal collapse model of Bond & Myers (1996), butwere unable to account for the magnitude of assemblybias seen in simulations. Mo et al. (2005) suggested that“pre-heating”, arising from the partial collapse of struc- Dalal et al.ture into sheets and filaments, could suppress mass accre-tion onto halos within those sheets and filaments. Thismechanism would clearly be incapable of explaining as-sembly bias in high-mass halos, which collapse before theformation of sheets and filaments (Bond et al. 1996), andSandvik et al. (2007) found that this mechanism failsin low-mass halos as well. Desjacques (2007) investi-gated the effects of environment on the Bond & Myers(1996) ellipsoidal collapse model, but could not repro-duce the sign reversal of assembly bias between high andlow masses, or the magnitude of assembly bias seen inN-body simulations.In this paper, we will argue that many aspects of as-sembly bias may be understood quite simply. In par-ticular, we show that the statistics of Gaussian randomfields require significant assembly bias in high-mass ha-los, and we show using N-body simulations that the clus-tering of massive halos matches precisely the assemblybias that we predict. At low halo masses, we argue thatassembly bias arises largely due to the cessation of massaccretion onto a sub-population of small halos. Thesenon-accreting halos then naturally become unbiased overtime, and therefore cluster more strongly than the ma-jority of halos of similar mass, which are anti-biased. BIAS
In this section, we briefly review previous results onbiasing and assembly bias.We start by assuming that massive dark matter ha-los arise from the collapse of peaks in the initial (lin-ear) density field; this allows us to infer many of theproperties of halo statistics using the statistics of gaus-sian fields (Press & Schechter 1974; Bardeen et al. 1986;Bond et al. 1991; Bond & Myers 1996). High peaks ofthe density field will in general cluster more strongly thanthe average matter clustering, which leads to halo bias.In the linear bias model, we directly relate the halo fluc-tuations to the local matter density fluctuations, δ h = b δ, (1)where δ ≡ δρ/ ¯ ρ is the overdensity and b is the bias pa-rameter. On scales over which this expression is valid ,we may express halo two-point functions in terms of thematter two-point function, scaled by appropriate powersof the bias. For example, the halo power spectrum wouldbecome P h ( k ) = b P m ( k ).The clustering of halos reflects the clustering of ini-tial peaks only in part, however. Another contributionto halo clustering comes from the motion of peaks alonglarge-scale matter flows. At late times, the halo overden-sity becomes δ h ( a ) = δ L + δ m ( a ) (2)where the Lagrangian overdensity δ L corresponds to theclustering of the peaks in the initial density field, andthe second term δ m reflects the motion of the peaks dueto advection along matter flows. On sufficiently largescales, in which we can treat the motion of the peaksas tracing the average bulk motion, the second term is This approximation should be valid on scales where the mattercorrelation is small, ξ ( r ) ≪ ( k ) = k P ( k ) / π ≪ simply the matter overdensity. Then the (Eulerian) biasbecomes b ( a ) = δ h ( a ) /δ m ( a ) = 1 + δ L /δ m ( a )= 1 + b L ( a ) . (3)Note that we explicitly include the time dependence ofthe matter overdensity, but that the Lagrangian over-density of peaks is not time dependent: it simply reflectsthe clustering of peaks in the initial conditions. There-fore, over time, the clustering associated with motion willdominate over the clustering from the initial conditions.That is, if we can treat the halos as test particles, we ex-pect their Lagrangian bias to decay with time inverselywith the linear growth factor, b L ( a ) ∝ /D ( a ) (e.g. Fry1996; Tegmark & Peebles 1998).Henceforth, we will focus on the Lagrangian bias b L ,noting that the conversion to Eulerian bias b is simplygiven by Eqn. (3) on linear scales. The next step is todetermine b L . From Eqn. (1), we can think of bias asthe rate at which the number of halos or peaks changesas the background density is varied. If we associate darkmatter halos with rare peaks above some high thresholddensity, then it becomes easy to see why massive halosare biased. In the Press-Schechter model, for example,the (Lagrangian) number density of halos is n PS ∝ δ c σ exp (cid:18) − δ c σ (cid:19) , (4)and changing the background density level by δ is equiv-alent to changing the threshold density to δ c − δ . Thenthe bias becomes b PS = n − dndδ = − n − dndδ c = δ c σ − δ c (5)(Cole & Kaiser 1989). This precise expression appliesonly for the Press-Schechter mass function, but its qual-itative features are generally valid. At high masses( σ ≪ δ c ), peaks above threshold are much more clusteredthan matter, because changing the background densitylevel dramatically changes the probability out on the tailof the Gaussian distribution. In contrast, at low masses( σ ≫ δ c ), halos are anti-biased ( b L <
0) because rais-ing the background density level reduces the number ofpeaks at low mass, by converting them to higher mass.Clearly, many aspects of the mass dependence of biasmay be inferred by consideration of the initial densitypeaks from which halos form. We can largely understandthe clustering of Eulerian halos found in complex nonlin-ear simulations by adopting a Lagrangian viewpoint, andexamining the properties of peaks of the initial lineardensity field. As mentioned in the introduction, simu-lations have shown that bias depends upon additionalhalo properties besides mass, such as assembly historyor structural parameters like concentration. The philos-ophy of this paper will be to assume that this additionaldependence may also be understood from the statistics ofthe initial Gaussian peaks. Accordingly, we will study N-body halos found in simulations described below, usinga Lagrangian viewpoint. SIMULATION AND POST-PROCESSING
We have performed N-body simulations of structureformation in scale-free cosmologies, with Ω m = 1 andalo assembly bias 3power-law power spectra. We have chosen to run scale-free models rather than more realistic ΛCDM models inorder to simplify the problem as much as possible, andremove possible effects from late-time cosmic accelera-tion or a varying spectral index. On scales far from thebox size and far from the Nyquist frequency, the onlyrelevant scale in the problem is M ⋆ ( a ), the characteristicnonlinear mass scale satisfying σ ( M ⋆ , a ) = δ c = 1 .
68. Ifthe power spectrum slope is n , such that P ( k ) ∝ k n , thenthe rms mass fluctuations scale like σ ( M ) ∝ M − (3+ n ) / .Because the linear growth factor behaves as D ( a ) = a for this critical-density cosmology, the nonlinear massscales like M ⋆ ( a ) ∝ a / (3+ n ) . Since the nonlinear massvaries rapidly with expansion factor a , by examining haloproperties at multiple redshifts we can study halo prop-erties over a wide range of M/M ⋆ (Efstathiou et al. 1988;Wechsler et al. 2006).In this paper, we will focus on the largest simulation weran, with 1024 particles, using a power spectral index n = −
2. We normalized the power spectrum such that σ ( M, a ) = 10 aM − / , where the mass M is measured innumber of particles. The simulation was evolved usingthe adaptive P M code gracos (Shirokov & Bertschinger2005), using a Plummer softening length ε P = 0 . a = 0 . a from a = 0 . a = 1, con-taining phase-space coordinates (comoving positions andpeculiar velocities) and ID’s for each particle.The phase space data for each output is used to findthe (sub-)halos in a two-step process. First we generate acatalog of halos using the Friends-of-Friends (FoF) algo-rithm (Davis et al. 1985) with a linking length of b = 0 . b ,with a density of roughly ρ > / (2 πb ) ≃
100 times thebackground density. We keep all halos with more than32 particles.We identify subhalos of the FOF groups using a newimplementation of the
Subfind method of Springel et al.(2001) which defines subhalos as gravitationally self-bound aggregations of particles bounded by a densitysaddle point. After some experimentation with differenttechniques we found this method gave a good match towhat would be selected “by eye” as subhalos. We usea spline kernel with 16 neighbors to estimate the den-sity and keep all subhalos with more than 20 particles.We call the most massive subhalo in any FoF group thecentral halo. The other subhalos we refer to as satellites.For each subhalo we compute and store a number ofproperties including the bound mass, velocity dispersion,peak circular velocity ( v c ), total potential energy, posi-tion and velocity. Since mass loss can be quite extremefor some subhalos we will use v c rather than mass toquote subhalo ‘size’. The mass, peak circular velocityand potential are highly correlated, however since oursubhalos are not spherical, there is a non-trivial scatterbetween v c and mass at fixed time.A merger tree is computed from the set of subhalocatalogs by identifying for each subhalo a “child” at alater time. We process 4 consecutive simulation outputsat a time and for each subhalo in the earliest output we define its child as the subhalo at a later time whichmaximizes α = ln − (cid:18) a a (cid:19) (cid:20) − | M − M | M + M (cid:21) X i ∈ φ i (6)where a and M are the scale-factor and mass of the pro-genitor, a and M are the candidate halo at the latertime, φ i is the potential of particle i computed usingthe particles in subhalo 1 and the sum is over all par-ticles in the progenitor that also lie in the candidate.The weighting by φ ensures that the “core” particlesin the progenitor are given more weight than the lessbound particles. We found that weighting by φ gavebetter results than weighting by φ , but higher powers of φ did not perform appreciably better than φ . We find ∼
95% of the subhalos have a child in the next time step,with one or two percent skipping one or more outputtimes. A similar method for subhalo tracking has beenemployed by Springel et al. (2005), Faltenbacher et al.(2005), Kravtsov et al. (2004), Allgood et al. (2006), andHarker et al. (2006).In addition to the Eulerian halo properties mentionedabove, we also require Lagrangian properties of halos,such as their centroids or mean initial velocities. Fromthe initial conditions for the simulation, we have the ini-tial density field δ ( x ), the initial comoving displacementfield d ( x ) satisfying δ + ∇ · d = 0, the peculiar velocityfield, which in the Zeldovich approximation for Ω m = 1is v = aH d , and the strain field S = ∇ d . Note thatTr S = − δ . These fields are all computed on the samegrid as the initial (undisplaced) particle positions, byFourier transform of the initial density field. For exam-ple, the strain field S ij is computed by FFT of δ k k i k j /k .We would like to average these quantities over the La-grangian volumes occupied by the halos. We may do sousing the unique ID’s for each particle. For each field F ,we compute the average h F i over every halo by a simpleaverage over the halo’s particles: h F i i = N − i N i X j =1 F ( x i,j ) , (7)where x i,j is the Lagrangian position of the j th parti-cle of halo i . Lastly, we will also require derivativesof these smoothed properties with respect to mass, e.g. d h δ i /d (log M ). We use the merger tree described aboveto compute derivatives with respect to mass. For ev-ery halo, we move along the most massive branch of itsmerger tree for a factor of 2 in mass, tabulating mass M and the average h F i . We then fit a linear expansion, F = F + F ′ ∆(log M ) in order to measure the derivative F ′ ≡ dF/d log M .Lastly, we will determine bias parameters in Fourierspace. We measure the matter power spectrum P ( k ),the halo auto-power spectrum P h ( k ), and the halo-matter cross spectrum P c ( k ), and define the halo bias as b = P c ( k ) /P ( k ). We use the ratio of the cross-spectrumrather than the halo auto-spectrum, because the formeris less sensitive to shot noise in the halo counts. HIGH MASS HALOS: M ≫ M ⋆ We first focus on halos well above the nonlinear massscale, M ≫ M ⋆ , where σ ( M ) ≪ δ c . These halos are par- Dalal et al. Fig. 1.—
Average Lagrangian overdensity h δ i for halos as a func-tion of σ ( M ). In the high mass limit, σ →
0, collapse becomesspherical and h δ i → δ c = 1 .
68. Point colors correspond to peakcurvature, d h δ i /d log M , as denoted in the colorbar on the right.The solid black curve and error bars depict the mean and standarddeviation of δ in bins of σ . ticularly simple : because they are so massive, they dom-inate their environments, and because they are far out onthe tail of the probability distribution of δ , the mass de-pendence of their bias scales like b ≃ δ c /σ . Anti-biasingfrom larger-scale halos should be negligible. Another keysimplification in this regime is that the collapse of theserare peaks is nearly spherical (Bardeen et al. 1986, Ap-pendix C). Because the collapse is roughly spherical, weexpect the collapse threshold predicted by the sphericalcollapse model (Gunn & Gott 1972) to be a good approx-imation: δ c ≈ .
68. This expectation is borne out by ourN-body results. In figure 1, we plot the average linearlyevolved overdensity for FOF halos well above M ⋆ , andindeed we find h δ i ≃ . & σ peaks. Note thatmultiple different redshifts and halo masses are plottedin this figure.The validity of δ c ≈ .
68 for halos in the regime σ ≪ b ∼ δ c /σ as mentioned above. This is the averageover all peaks of a given mass, which marginalizes overany other peak parameters which could affect bias. Oneobvious additional parameter which the bias must de-pend upon is the curvature of the peak, which we willparametrize by s = d h δ i /d (log M ). To see this, notethat we can estimate the large-scale background den-sity δ b in which the halos reside by a Taylor expansion: δ b ≈ δ + ∆(log M ) × dδ/d (log M ) + . . . . Two peaks ofthe same mass M and peak height δ but different peakcurvatures s > s (i.e. | s | < | s | ) will tend to live indifferent environments: peak 1 will tend to have a largerbackground density than peak 2, and therefore will havea larger bias. It is straightforward to estimate the depen-dence of bias on curvature. Writing ν = δ/σ , x = s/σ s ,and the cross-correlation coefficient γ = h νx i , the biasbecomes (Bardeen et al. 1986) b L ≈ σ − ν − γx − γ . (8)We can easily sketch a derivation of this expression fol-lowing the argument of Kaiser (1984). Let us com- Bardeen et al. (1986) parametrized curvature not with dδ/d log M as we do, but rather with the closely related ∇ δ . −0.5 −0.4 −0.3 −0.2 −0.15.05.56.06.57.07.58.08.59.0 d(log δ )/d(log M ) b E Fig. 2.—
Eulerian bias b E = 1+ b L for halos with 110-130 particlesat a = 0 . σ = 0 .
49, and are roughly ∼ M ⋆ at this redshift. Halos arebinned into quartiles of d log δ/d log M , and halo bias in each binis measured by the ratio of the halo-matter cross-spectrum to thematter power spectrum. The solid black curve shows Eqn. (8). pute the correlation function for regions above a highthreshold ν t . Consider regions 1 and 2, each of size R and separated by distance r , with peak heights ν , ν and curvatures x , x . These Gaussian fields are de-scribed by a covariance matrix with diagonal elements h ν i = h x i = 1. Of the off-diagonal elements, only h ν ν i = ψ ( r ) and h ν x i = h ν x i = γ are important;long-range correlations involving derivatives of the den-sity field fall off more quickly than these terms by powersof R/r . We compute the two-point correlation of re-gions with ν > ν t by integrating the Gaussian probabilityover ν > ν t , ν > ν t . If we are uninterested in the cur-vature x , then we first marginalize over x and x , leav-ing a 2 × ψ , and we find that the two-point function for regionsabove threshold becomes ξ pk = ν t ψ , giving a Lagrangianbias b L = ν t /σ (Kaiser 1984). We can similarly computethe two-point function for regions above threshold witha specified curvature x by repeating this calculation, butnot marginalizing over x and x . Recalling that theone-point probability distribution for ν at a given x is aGaussian centered at h ν | x i = γx with variance 1 − γ ,we change variables from ν to ν ′ = ( ν − γx ) / (1 − γ ) / .Noting that h ν ′ ν ′ i ≈ ψ/ (1 − γ ), we immediately seethat the two-point correlation function for regions abovethreshold becomes ξ pk = ( ν ′ t ) ψ/ (1 − γ ), giving a bias b L = ν ′ t / ( σ p − γ ), i.e. Eqn. (8).The bias of halos found in our simulation matchesEqn. (8) quite well. For a power spectrum index n = − σ s = σ/ √ γ = − / √ ∼ b corresponds to a factor of 2 variation in thecorrelation function for halos of this mass.One subtlety which we have glossed over in computingthe statistics of high-mass peaks is the issue of the ap-alo assembly bias 5 r / r sm P ( r ) Fig. 3.—
Probability for a particle to be linked into a high massFOF group as a function of initial Lagrangian radius. The thicksolid black curve shows the stacked Lagrangian volumes for halosin the mass range 200-300 particles at a = 0 . ∼ M ⋆ . Note that mass in subhalos is excluded fromthis average. The blue dotted curve shows the expected averageprofile for a tophat window, while the red dashed curve correspondsto a Gaussian of this mass range. The thin solid black curves showthe average Lagrangian profiles for halos of mass between 200-300particles at later outputs, a = 0 . , . , . M/M ⋆ ≃ , , . propriate smoothing to use. We have assumed a top-hatsmoothing, however there is no compelling theoreticalreason to favor a top-hat kernel versus other possibili-ties, such as a Gaussian (although the sharp k -space filterimplicit in the Press-Schechter model appears unphysi-cal). As a sanity check on our top-hat smoothing, wehave computed the average Lagrangian volumes corre-sponding to high peaks, by stacking roughly 1000 halosof mass ∼ M ⋆ . The results are shown in figure 3.The average Lagrangian volume occupied by FOF halosis reasonably compact, with the probability for halo in-clusion steeply varying from near 1 to near 0 over a range0 . < r TH < .
2, where r TH is the top-hat smoothing ra-dius for mass M . In comparison, a Gaussian smoothingwindow is clearly more diffuse. We conclude that ouruse of top-hat smoothing is reasonable, though we cau-tion that these results may depend upon the slope of thepower spectrum. Finally, we note that there is a nontriv-ial amount of mass inside r < r TH that does not end upin the halo; we will return to this point below in section5. Halo observables
We have shown that there is a strong dependence ofhalo clustering upon peak curvature. How does this re-late to the assembly bias observed in previous simula-tions?In the high mass regime, there is a direct relationshipbetween peak curvature and halo assembly history. Re-call that in the spherical collapse model, regions com-plete their collapse once the average linearly evolvedoverdensity δ reaches δ c . This implies that we can es-timate the assembly history of halos simply by mea-suring how the smoothed overdensity h δ i varies withsmoothing scale. We would expect the accretion rate a N f o f Fig. 4.—
Mass accretion histories as a function of peak cur-vature. The upper solid (blue) curve shows the average accre-tion history for the halos comprising the leftmost bin in Fig. 2,while the lower solid (red) curve shows the average history for therightmost bin. The dotted curves depict the behavior M ∝ a α for α = − / ( d log δ/d log M ) expected from the spherical collapsemodel. and the curvature to be related by d (log M ) /d (log a ) = − [ d (log δ ) /d (log M )] − . The average mass accretion his-tories for high mass halos appears consistent with this ex-pectation. In figure 4 we plot the stacked mass accretionhistory profiles for halos of 120 particles at a = 0 . a = 0 . r < r s ) of curvature-selected halosmatch well with profiles of concentration-selected halos,the outer profiles (at r > r s ) begin to diverge.Within the statistical errors of our finite sample of rare Dalal et al. r / r ( r ρ ) / ( r ρ m ) Fig. 5.—
Radial profiles for halos as a function of peak curvature.The various curves show the stacked radial profiles for subsets ofhalos of 4000-5000 particles at a = 0 . M ∼ M ⋆ .The solid curves correspond to the upper and lower quartiles of d log δ/d log M , while the dashed curves correspond to the upperand lower quartiles of fitted concentration c . For comparison,the vertical dotted black line shows r = 2 . ε P , the softening length. halos, we have found no residual dependence of bias onconcentration, once we account for the peak curvaturedependence described above. Our simulation had rela-tively low force resolution, and so our fitted concentra-tions may be sufficiently noisy to mask any underlyingresidual concentration dependence of bias, but it appearsthat assembly bias in rare, high mass halos is in agree-ment with theoretical expectations. The interpretationof this effect is straightforward: halos residing in highdensity environments have an enhanced mass supply, andtherefore grow faster, than halos in low density environ-ments. LOW MASS HALOS: M ≪ M ⋆ In this section we consider assembly bias for low masshalos with M ≪ M ⋆ , for which σ ( M ) ≫ δ c . The age-and concentration-dependence of bias at low masses isconsiderably stronger than in the high mass regime, andof the opposite sign: the oldest, highly concentratedhalos have a bias more than twice as large as that ofthe youngest, low-concentration halos of similar mass M ≪ M ⋆ (Gao, Springel & White 2005; Wechsler et al.2006; Gao & White 2007).In agreement with previous workers, we also find thisreversed sense of assembly bias. In keeping with our La-grangian viewpoint, we use mean Lagrangian overdensity h δ i as a proxy for age. In figure 6, we plot average ac-cretion histories for halos of 200 particles at a = 1, splitinto quartiles of h δ i . The oldest halos show little growthat late times, accreting only ∼
20% of their mass af-ter a = 0 .
5. In comparison, the nonlinear mass scaleincreases by a factor ≈
64 over this same time.Because the oldest halos do not grow appreciably inmass, we can think of their dynamics as those of testparticles following large-scale bulk flows. As discussed insection 2, we therefore expect the bias of this populationto relax over time to b →
1. This is indeed the behav-ior observed for our oldest halos, as shown in figure 7, a N f o f Fig. 6.—
Mass accretion histories for halos of 190-210 particlesat a = 1, corresponding to M = 0 . M ⋆ . The halos are split intoquartiles of h δ i , which correlates well with mass accretion history. δ b E Fig. 7.—
Eulerian bias for the same halos of ∼
200 particlesshown in Figure 6. The oldest halos, at high h δ i , are nearly unbi-ased while the rest of the population is strongly antibiased. Remov-ing all halos which previously were subhalos, and all halos bound tomore massive halos, removes the old unbiased subpopulation (reddashed curve). which plots the (Eulerian) bias of the same halos shownin figure 6. The oldest ∼
20% of halos are nearly un-biased, while the remainder are strongly anti-biased, toa degree consistent with predictions of even the simplestPress-Schechter model, b → − δ − c ≈ . ∼
200 particle halos, but excluding thosewhich in earlier outputs were subhalos. This removesalo assembly bias 7 σ / v d P / d l og σ Fig. 8.—
Probability distribution for ambient 3-D velocity dis-persion σ within a radius 5 r of halos with 200-300 particles at a = 1, roughly 0.005 M ⋆ . Halos are split into 5 bins of h δ i , withred, magenta, green, cyan, and blue from lowest to highest respec-tively. The oldest halos tend to live in hot environments, whereasthe material in the vicinity of most low-mass halos is sufficientlycold to permit accretion. ∼ −
25% of the total population. Because the timesampling of our simulation outputs do not finely coverthe dynamical times of all halos, this cut could miss sub-halos with orbital pericenters inside of larger halos, butapocenters well outside. Ludlow et al. (2008) have re-cently argued that a significant fraction of subhalos haveorbital apocenters extending out to many times the virialradius of their host halo. To account for this, we alsoconservatively exclude all halos which are within 3 r vir and gravitationally bound to larger halos. This will ex-clude the population discussed by Ludlow et al. (2008),as well as any halos entering the infall region of massiveobjects. This cut removes a further ∼
15% of low masshalos, leaving &
60% of the original low-mass popula-tion. The effect of these cuts is effectively to exclude theold subpopulation, even though the cuts do not explic-itly refer to halo age or accretion history. It thereforedoes seem likely that, as suggested by previous authors,the old subpopulation of low-mass halos corresponds toearly-formers whose growth was stunted by environmen-tal effects (Wang et al. 2007).One plausible environmental effect that would accountfor the arrested development of the oldest low-mass halosis simply the enhanced velocity dispersion in the vicinityof massive halos. In environments where the ambientvelocity dispersion greatly exceeds the escape velocityof low-mass halos, those halos will be unable to accretesignificant material. We indeed find that the oldest halospreferentially reside within hot environments, as shownin Fig. 8.Given that the oldest low-mass halos tend to exist inproximity to larger halos, we would expect them to orig-inate (in Lagrangian space) from the vicinities of highmass halos as well, c.f. Fig. 3. As discussed in the previ-ous section, much of the material within the Lagrangianvolumes of high mass halos ends up within the final col-lapsed halo, so how do we understand the survival ofthe old low-mass halos in the presence of their rapidly r / r sm δ E / ( a H r s m ) Fig. 9.—
Binding energy versus initial radius for peaks that areaccreted (red), that become subhalos (green), or that escape (blue). accreting neighbors? One possible explanation noted byLudlow et al. (2008) is that these halos are ejected out oflarger halos by 3-body interactions. Because our goal isto understand the clustering of halos in terms of their La-grangian properties of known (Gaussian) statistics, andsince complicated multi-body interactions are clearly dif-ficult to model using only Lagrangian quantities, we havenot attempted to address this mechanism using our sim-ulations. Rather, we have tried to find Lagrangian prop-erties of halos which correlate with their survival aroundlarger structures.One Lagrangian quantity which correlates with sur-vival is plotted in Figure 9. We have taken halos ofmass 40,000-60,000 particles at a = 1, nearly M ≃ M ⋆ ,and identified smaller-scale peaks within the sphericalLagrangian volumes centered on the larger halos. Asshown in the figure, the relative binding energy of thesubpeaks roughly correlates with their escape likelihood.We have defined the binding energy δE in the followingmanner. The gravitational potential near a peak may bewritten φ ( x ) ≈ φ + x · ∇ φ + 12 x · ∇∇ φ · x + . . . (9)The first two terms do not affect the collapse of the peak,and so we define the perturbed potential energy as thequadratic term, δφ ( x ) = 12 x · ∇∇ φ · x = 3Ω H a x · S · x (10)where we have used the Poisson equation to relate thestrain tensor S = ∇∇ ( ∇ − δ ) to the potential’s curvaturetensor. We can similarly write the kinetic energy near thepeak as T ( x ) = 12 | v | = 12 | aH ( x + d ) + v p | ≈
12 ( aH | x | ) + 2( aH ) x · d + . . . (11)where we have used the Zeldovich approximation forΩ m = 1 to relate the peculiar velocity and comovingdisplacement, v p = aH d . The unperturbed piece againdoes not affect the peak dynamics, and so we write the Dalal et al.perturbed kinetic energy as δT = 2( aH ) x · d , (12)and the total perturbed energy becomes δE = δφ + δT .In computing δφ , we smooth the strain tensor S over thelarge halo, whereas in δT we use the difference in themean Lagrangian peculiar velocities, ∆ v L = v L2 − v L1 .As can be seen from Figure 9, the initial Lagrangianbinding energy does not solely determine the survival ofsubpeaks; the escape threshold is radially dependent aswell. Without any physical justification, we find that anapproximate threshold E ′ ≡ δE + ( aH ) [2 r − r ] = 0 (13)roughly separates accreted subpeaks from escaping sub-peaks.Our assumption, therefore, is that subpeaks inside oflarger peaks can escape accretion onto their hosts if therelative energy is large enough, E ′ >
0. Due to environ-mental effects, however, they are unable to grow, andthereby produce the population of old, unbiased low-mass halos. In the next section, we test how well thissimple picture can reproduce our numerical bias results. A TOY MODEL FOR ASSEMBLY BIAS
In the previous two sections, we have discussed assem-bly bias in two asymptotic regimes, σ ≪ σ ≫ δ c .We found that two distinct effects generated assemblybias in these two regimes. At intermediate masses, σ ∼ h δ i and on peak curvature d log h δ i /d log M . Recall that these peak properties bothaffect halo properties like age and concentration. At highmasses, the scatter in peak height among peaks abovethreshold is small, and so the scatter in peak curvaturedominates the assembly bias. At low mass, the peakheight is the dominant factor determining halo age orconcentration, causing a reversal in the sense of assemblybias. To illustrate, we plot in the right panel of Fig. 10the bias for halos selected on h δ i /δ c − d log h δ i /d log M . Athigh masses, the curvature term dominates, while at lowmasses, the first term dominates. Halo samples selectedon formation time or concentration show similar behav-ior. At intermediate masses, both of these effects aresignificant and compete with each other, and the magni-tude and sign of assembly bias will depend in detail uponthe precise age indicator that is employed. For example,an age indicator that correlates more strongly with peakheight rather than peak curvature will show a crossoverat higher masses than the example shown in Fig. 10c.This may explain why different measures of halo age givedifferent amounts of assembly bias when applied to thesame samples of halos (Li et al. 2008).Using our results on the assembly bias of low- andhigh-mass halos, we have constructed a simple toy modelfor halo biasing. We generate a Gaussian random den-sity field, smooth the field on multiple scales, and start-ing from the largest smoothing scales, search for den-sity peaks exceeding the ellipsoidal collapse threshold (Bond & Myers 1996, hereafter BM96). Isolated peaksabove threshold are labeled as collapsed halos, whereaspeaks falling within 1.2 smoothing radii (c.f. Fig. 9) oflarger scale collapsed regions are only labeled as collapsedif their binding energies relative to their larger hosts ex-ceeds the threshold found empirically from our N-bodyhalos, i.e. Eqn. (13).A key component of this calculation is the collapsethreshold. A comparison of the simple BM96 modelwith our N-body results shows that the model workswell at predicting collapse thresholds for reasonably highpeaks, σ .
1, however at low peaks σ ≫
1, the modelbegins to break down. Very roughly, the ellipsoidalcollapse threshold is larger than the spherical collapsethreshold by a ratio δ ec /δ sc ∼ (1 − δ sc e v ) − , wherethe ellipticity e v of the strain tensor is given by thedifference between its largest and smallest eigenvalues, e v = ( S − S ) / δ . The ellipticity and prolaticity arerandom variables whose probability distribution is re-lated to peak height (Bardeen et al. 1986, BM96), butvery roughly the mean ellipticity scales as h e v i ∼ / (2 ν ).For ν <
1, the BM96 collapse threshold quickly grows,and for sufficiently large e v and − p v the model predictsthat collapse never occurs. However we find many lowmass halos with average peak heights h δ i falling far be-low the collapse thresholds expected given their averageLagrangian strain tensors h S i . This is not surprising. Asnoted explicitly by BM96, their homogeneous ellipsoidmodel assumes that the evolution of the external tidalfield acting upon the collapsing peak may be determinedfrom the strain tensor h S i averaged over the local volumeof the peak. While it is certainly true that the exter-nal tidal field and local strain tensor are closely relatedat early times, in the nonlinear regime the one-to-onecorrespondence between them breaks down. On scaleswhere σ ≫
1, the external tidal field evolves nonlinearlyin a manner that cannot be predicted from purely lo-cal measurements over the peak volume. Therefore, wewould not expect the simple BM96 model to apply inthis regime, as those authors note.Nevertheless, we require some choice of collapse thresh-old for our toy model. At low masses, the lowest h δ i ’sare roughly ∼ . δ c = min( δ BM96 , . z c by1 + z c = h δ i /δ c . Writing z s and z b as the collapse red-shifts for the sub-peak and large peak, respectively, werequire (1 + z s ) ≥ z b ) for a subpeak to be labeledas a collapsed halo.Figure 11 shows results of this simple calculation. Inorder to compare with Fig. 10, we plot bias as a functionof the same peak parameters. At high masses, σ ≪ σ ≫
1) the bias of isolated peaks asymptotes to b E ≈ . b E →
1. Thesefeatures appear fairly generic. Behavior in the transitionregion, σ ∼ −
2, depends upon details of the model. Inparticular, the bias of the old halos in the intermediatemass regime depends upon the fraction of halos that aresubpeaks, and varying parameters in the model changesthis fraction.Overall, however, we find reasonable agreement be-alo assembly bias 9 σ b E σ b E σ b E Fig. 10.—
Bias as a function of peak parameters. Bias is plotted for halos in the range 200-300 particles at various redshifts. In the leftpanel, the red and blue curves show bias for the upper and lower 20% of peak height h δ i /δ c , where δ c is the ellipsoidal collapse thresholdcalculated from each halo’s local strain tensor. In the middle panel, we split halos on peak curvature, d log h δ i /d log M . In the right panel,we split halos on an average of these two quantities, h δ i /δ c − d log h δ i /d log M . At high masses, the scatter in curvature dominates, whileat low masses, the scatter in peak height dominates, leading to a transition in assembly bias at intermediate masses. Formation time orconcentration shows similar behavior. σ b E σ b E σ b E Fig. 11.—
Similar to Fig. 10, but for peaks found in our toy model calculation rather than N-body halos. tween the toy model and our N-body simulation in theregimes where assembly bias is significant, for M ≫ M ⋆ and M ≪ M ⋆ . The point of this exercise is not toclaim that we have solved the halo formation problem,but merely to illustrate that the assembly bias in theseregimes is generic. The assembly bias at high masses isan unavoidable consequence of Gaussian statistics, andtherefore any halo formation model which relates ini-tial density peaks to halos will also find similar levelsof assembly bias, depending somewhat on the peak def-inition. The behavior at low masses may appear moremodel dependent, however we do not believe our resultsare specific to our particular mechanism for production ofnon-accreting halos. Rather, we expect that any mecha-nism which produces roughly the correct fraction of non-growing, low-mass halos will also match the assemblybias found in N-body simulations. DISCUSSION & CONCLUSIONS
We have investigated halo assembly bias in hierarchicalstructure formation. We find that assembly bias is sig-nificant both at high masses ( M ≫ M ⋆ ) and low masses( M ≪ M ⋆ ). Two different competing effects generate as-sembly bias in these two regimes. At high masses, peaksof low curvature are more strongly clustered than peaksof large curvature, in a manner entirely consistent with the statistics of peaks of random Gaussian fields. Weprovide a simple expression for the bias in the high-massregime which accounts for assembly bias and reduces tothe Press-Schechter bias when marginalized over all pa-rameters excluding mass (note that fitting formulae forbias calibrated at lower mass, such as those proposedby Sheth & Tormen (1999) or Seljak & Warren (2004),fail in the high mass regime (Cohn & White 2007)). Atlow masses, assembly bias appears largely driven by theformation of a non-accreting sub-population of low-masshalos in the vicinity of larger halos. Because the non-accreting halos become unbiased, they are more stronglyclustered than accreting halos of similar mass, which areanti-biased. We found that a simple toy model is ableto reproduce our numerical bias results in the asymp-totically large and small mass regimes. At intermediatemasses ( M ∼ M ⋆ ), both of the above effects become im-portant, and the assembly bias for a sample of halos willinvolve a weighted average of both effects.The application of these results to observed objectslike galaxies or quasars is not straightforward, becauseit will depend upon the halo occupation distribution(HOD) of those objects. The very simplest HOD mod-els assume that the galaxy content depends only uponhalo mass, and not upon assembly history, whereas semi-analytic models used in conjunction with N-body simu-0 Dalal et al.lations tend to find a strong dependence of galaxy prop-erties on assembly history at fixed mass (Springel et al.2005; Croton et al. 2007, e.g.). If the former assump-tion is valid, then assembly bias will not be importantfor galaxy clustering. Perhaps surprisingly, Tinker et al.(2007) have recently found that properties of observedSDSS galaxies do not appear to correlate significantlywith local environment, once the dependence on hosthalo mass has been accounted for, suggesting that as-sembly bias may indeed be unimportant for modelinggalaxy clustering.Besides the clustering of halos, our results also haveimplications for our understanding of halo assembly. Inparticular, we have found that massive halos are verywell described by spherical collapse. The average pre-collapse Lagrangian volumes occupied by the most mas-sive halos are approximately spherical regions of size r ∼ r TH = (3 M/ πρ ) / . However, r TH is not a sharpboundary in Lagrangian space: there is significant massoutside r > r TH at early times which is found inside thevirial radius at late times, and similarly there is consid- erable mass at r < r TH which avoids accretion at latetimes. Indeed, we might speculate that the post-collapseradial distribution of the matter inside and around a darkmatter halo is simply related to the pre-collapse profile oflocal density and potential. Since both of these quanti-ties are described by well-understood Gaussian statistics,this raises the intriguing possibility of understanding theinternal structure of virialized objects from first princi-ples. This topic is the subject of further study.We thank Andrey Kravtsov and Simon White for sev-eral helpful discussions. We also thank Joanne Cohn,Andrey Kravtsov and Andrew Wetzel for comments onan earlier version of this paper. All simulations wereperformed on CITA’s Sunnyvale cluster, funded by theCanada Foundation for Innovation and the Ontario Re-search Fund for Research Infrastructure. N.D. was sup-ported by the Natural Sciences and Engineering ResearchCouncil of Canada (NSERC). M.W. was supported inpart by NASA. REFERENCESAllgood, B., Flores, R. A., Primack, J. R., Kravtsov, A. V.,Wechsler, R. H., Faltenbacher, A., & Bullock, J. S. 2006,MNRAS, 367, 1781Angulo R.E., Baugh C.M., Lacey C.G., 2008, preprint[arXiv:0712.2280]Avila-Reese, V., Col´ın, P., Gottl¨ober, S., Firmani, C., &Maulbetsch, C. 2005, ApJ, 634, 51Bett P., Eke V., Frenk C.S., Jenkins A., Helly J., Navarro J.,2007, MNRAS, 376, 215Bardeen J.M., Bond J.R., Kaiser N., Szalay A.S., 1986, ApJ, 304,15Bond J.R., Cole S., Efstathiou G., & Kaiser N. 1991, ApJ, 379,440Bond, J. R., Kofman, L., & Pogosyan, D. 1996, Nature, 380, 603Bond J.R., Myers S.T., 1996, ApJS, 103, 1Cohn, J. D., & White, M. 2007, arXiv:0706.0208Cole, S. & Kaiser, N. 1989, MNRAS, 237, 1127Conroy, C., Wechsler, R.H., & Kravstov, A.V. 2006, ApJ, 647, 201Cooray, A. & Sheth, R. 2002, Phys. Rep., 372, 1Croton, D. J., Gao, L., & White, S. D. M. 2007, MNRAS, 374,1303Cuesta, A. J., Prada, F., Klypin, A., & Moles, M. 2007,arXiv:0710.5520Davis, M., Efstathiou, G., Frenk, C.S., & White, S.D.M. 1985,ApJ, 292, 371Desjacques, V. 2007, arXiv:0707.4670Efstathiou, G., et al. 1988, MNRAS, 235, 715Faltenbacher A., Allgood B., Gottloeber S., Yepes G., HoffmanY., 2005, MNRAS, 362, 1099Fry, J. N. 1996, ApJ, 461, L65Gao, L., Springel, V., & White, S.D.M 2005, MNRAS, 363, L66Gao, L., & White, S. D. M. 2007, MNRAS, 377, L5Gunn, J. E., & Gott, J. R. I. 1972, ApJ, 176, 1Harker, G., Cole, S., Helly, J., Frenk, C., & Jenkins, A. 2006,MNRAS, 367, 1039Jing Y.-P., Suto Y., Mo H.-J., 2007, ApJ, 657, 664Kaiser, N. 1984, ApJ, 284, L9Kauffmann, G. & Haehnelt, M. 2000, MNRAS, 311, 576Kravtsov, A. V., Gnedin, O. Y., & Klypin, A. A. 2004, ApJ, 609,482Lemson, G. & Kauffmann, G. 1999, MNRAS, 302, 111 Li, Y., Mo, H. J., & Gao, L. 2008, arXiv:0803.2250Lima, M., & Hu, W. 2004, Phys. Rev. D, 70, 043504Ludlow, A. D., Navarro, J. F., Springel, V., Jenkins, A., Frenk,C. S., & Helmi, A. 2008, arXiv:0801.1127Mo, H. J., Yang, X., van den Bosch, F. C., & Katz, N. 2005,MNRAS, 363, 1155Mo, H.J., & White, S.D.M. 1996, MNRAS, 282, 347Navarro, J. F., Frenk, C.S., & White, S.D.M. 1997, ApJ, 490, 493Peacock, J. A. 1999,