Excursion set approach to the clustering of dark matter haloes in Lagrangian space

19 Pages·0.49 MB·English

by  Cristiano Porciani

Mon.Not.R.Astron.Soc.000,000–000 (0000) Printed1February2008 (MNLATEXstylefilev1.4)
Excursion set approach to the clustering of dark
matter haloes in Lagrangian space
Cristiano Porciani1, Sabino Matarrese2, Francesco Lucchin3 and Paolo Catelan4,5
1 SISSA, Scuola Internazionale di Studi Superiori Avanzati, via Beirut 2-4, I-34014 Trieste, Italy
2 Dipartimento di Fisica GalileoGalilei, Universit`adi Padova, via Marzolo 8, I-35131 Padova, Italy
8 3 Dipartimento di Astronomia, Universit`adi Padova, vicolo dell’Osservatorio 5, I-35122 Padova, Italy
9 4 Theoretical Astrophysics Center, Juliane Maries Vej30, DK-2100 Copenhagen Ø, Denmark
9
5 Department of Physics, Astrophysics, Nuclear Physics Laboratory, Keble Road OX1 3RH, Oxford, UK
1

n
a
J
8 1February2008
2

1
ABSTRACT
v
0 We present a stochastic approach to the spatial clustering of dark matter
9 haloes in Lagrangian space. Our formalism is based on a local formulation
2
of the ‘excursion set’ approach by Bond et al., which automatically accounts
1
0 for the ‘cloud-in-cloud’ problem in the identification of bound systems. Our
8
methodallowstocalculatecorrelationfunctionsofhaloesinLagrangianspace
9
/ usingeitheramulti-dimensionalFokker-Planckequationwithsuitablebound-
h
p ary conditions or an array of Langevin equations with spatially correlated
- random forces. We compare the results of our method with theoretical pre-
o
r dictionsforthe haloauto-correlationfunctionconsideredinthe literatureand
t
s find good agreement with the results recently obtained within a treatment of
a
: halo clustering in terms of ‘counting fields’ by Catelan et al.. The possible
v
i effect of spatial correlations on numerical simulations of halo merger trees is
X
finally discussed.
r
a
Key words: galaxies: clustering - cosmology: theory - large-scale structure
of Universe.
1 INTRODUCTION
The celebrated Press & Schechter (1974, hereafter PS) theory represents a fundamental tool to determine the mean
massdistributionofboundcondensations,nowadaysidentifiedasdarkmatter(DM)haloes,intheUniverse.Recently
a number of authors (e.g. Efstathiou et al. 1988; Cole & Kaiser 1989; Mo & White 1996; Mo, Jing & White 1996;
Catelanetal.1997,hereafterCLMP)wentbeyondthisapplication,extendingthePSmethodtostudyhaloclustering
both in Lagrangian and Eulerian space.
The clustering problem can be dealt with in two steps: i) identifying the preferential sites for halo formation in
Lagrangianspaceandii)providingatheoreticalframeworktofollowthenon-lineardynamicsofthehalodistribution.
The former step can be accurately modeled within linear theory by the PS scheme, which in fact suggests a simple
criterion to assign each Lagrangian point to a matter clump of some mass M which is going to collapse at some
formation redshift z . Thelatterstep,instead, requiressome knowledgeof thenon-linearmatterdynamics,which is
f
(cid:13)c 0000RAS
2 C. Porciani, S. Matarrese, F. Lucchin & P. Catelan
ultimately responsible for theactual Eulerian distribution of matterclumpsas theyform bythehierarchical process
of matter accretion and merging of subunits. A relevant progress in this direction has been made by Mo & White
(1996),whoshowedthatalinearlocalbiasrelationconnectstheEulerianhaloauto-correlationfunctiontothatofthe
mass,thusprovidingagood fittothetwo-pointfunctionofhaloesinN-bodysimulations.Amorerefined,non-linear
and non-local relation has been recently obtained by CLMP.
At theLagrangian level, however, one has to deal with theso-called cloud-in-cloud problem, inherent in thePS
scheme,namelythefactthattheirapproachselectsboundsystemsofgivenmassthatcanhavebeenalreadyincluded
in larger mass condensations of the same catalogue. A straightforward solution of this problem in connection with
halo correlations was sketched by CLMP, who adopted the ‘peak-background’ split (e.g. Bardeen et al. 1986; Cole
& Kaiser 1989) to solve the problem up to the resolution scale of the ‘background’ mass component, which is also
responsible for the motion of the haloes from their original Lagrangian positions.
The transition from the Lagrangian to the Eulerian halo distribution can be dealt with exactly thanks to the
continuityequation,asshownbyCLMP,whoobtainedananalyticformulafortheEulerian‘biasfield’connectingthe
halo number density fluctuations to the mass perturbations. This relation possesses the remarkable feature of being
non-linear and non-local, as it depends on the halo overdensity at the original Lagrangian positions. Therefore, the
statistical halo distribution in Lagrangian space, which can also be described in terms of a hierarchy of bias factors
(cf. CLMP), plays therole of the initial condition in generating theEulerian halo density field.
Aim of the present paper is to provide an exact treatment of the halo correlation properties in the Lagrangian
world, which extendsand putson sounder bases the results obtained by Mo & White (1996), as already sketched in
CLMP. Ourapproach hereisentirely based ontheexcursion set version ofthePStheory (Peacock & Heavens1990;
Cole 1991; Bond et al. 1991, hereafter BCEK), which allows to define halo populations free of the cloud-in-cloud
problem on all relevant scales.
Ourformalism hastobeconsideredineveryrespectasthenaturalextensionoftheBCEK modeltoincludethe
spatial clustering properties of the haloes. Once a cosmological model and a power-spectrum of density fluctuations
are selected the present formalism allows to construct halo correlations of any order in Lagrangian space. Here we
onlygiveexplicit resultsfor thetwo-pointfunction,butthegeneralization tohigherorderwould bestraightforward.
An obvious and important application of this type of study aims at modeling galaxy clustering at different
redshifts.ThekeyideabeingthatthespatialdistributionofDMhaloescanbeaclueforunderstandingtheclustering
propertiesofluminousobjectslikegalaxies withdifferentphysicalpropertiesandatvariousepochs(e.g.Kauffmann,
Nusser & Steinmetz 1997; Matarrese et al. 1997; Moscardini et al. 1997). The actual relation between the parent
DMhaloes andthegalaxies, however,dependson thewayhaloes accrete massandonhowtheygrowbytheprocess
of merging with the surrounding haloes. These processes are usually studied in terms of the so-called ‘halo merger
trees’,obtainedeitherbyMonteCarlomethods(e.g.Cole1991;Lacey&Cole1993;Kauffmann&White1993;Tozzi,
Governato & Cavaliere 1996; Somerville & Kolatt 1997; Salvador-Sol´e, Solanes & Manrique 1997 and references
therein) or by directly looking at the dynamics of particles in cosmological N-body simulations (e.g. Kauffmann et
al. 1997; Roukema et al. 1997; Governato et al. 1997). Realizations of halo merger trees can be in fact obtained by
Monte Carlo techniques in the frame of the excursion set theory, where halo accretion histories are associated to
random walks. The key idea is to average the linear density field over spheres of decreasing size centred on each
Lagrangian point, thus defining a set of trajectories of the smoothed field as a function of the resolution scale. At
each cosmic epoch the process of halo formation is here modeled as the upcrossing of a threshold by these random
trajectories.Therefore,everytrajectorygivesadetaileddescriptionofthemassaccretionhistoryofthecorresponding
halo. A merger tree is then obtained by suitably connecting the accretion histories of all the progenitors of a given
halo (which, in thetree-jargon, are usually called ‘branches’).
Inthelatterapproach,however,spatialcorrelationsbetweendifferentbranchesofthetreeaswellascorrelations
betweendifferenttreesarecompletelyignored.TheexistenceofthisshortcominghasbeenalsonoticedbyRodrigues&
Thomas(1996)inconnectionwiththeso-called‘blockmodel’(Cole&Kaiser1988;Cole1991),whichcanbeconsidered
asarudimentaryformofmergingtree,wherehaloescanonlygrowbydoublingtheirmass.Yano,Nagashima&Gouda
(cid:13)c 0000RAS,MNRAS000,000–000
Excursion set approach to halo clustering 3
Figure 1. Examples of trajectories obtained by smoothing the linear density fluctuation field with a series of sharp k-space
filters (with decreasing resolution scales Rf = 1/kf) at two points separated by the distance r. Here we consider r = 4Mpc
and a CDM power spectrum (with density parameter Ω = 1 and present Hubble constant H◦ = 50kms−1Mpc−1) linearly
extrapolated to σ8 =1. Theheavy continuous linerepresents thetrajectory associated to apoint xinLagrangian space. The
trajectoryassociatedtothepointysuchthat|y−x|=risplottedwithalightcontinuousline.Thedottedlineisobtainedby:
i)consideringthetrajectoryassociated to y,ii)artificiallyremovinganycorrelationwiththetrajectory atx andiii)suitably
rescalingthe result. Inpractice, the dotted linerepresents atrajectory whichis completely independent of the one associated
tox,butwhichhasthesamestatisticalone-pointproperties.Thisiswhatisgenerallyusedinbuildingupmergertreesofdark
matter haloes. However, the actual trajectories associated to neighbouring points are stronglycorrelated when the smoothing
length Rf is comparable to the physical separation between the points r. The same trajectories tend to become less and less
correlatedasRf decreases.Itisinterestingtonotethatamemoryofthestrongcorrelationexistingforr/Rf ≪1remainseven
whenr/Rf ≫1.Infact,whenRf issmallcomparedtor,therealtrajectoryassociated toy canbeaccurately approximated
byaddingaconstantvaluetotheindependenttrajectorysothatitmatchestherandomwalkatxwhenRf∼> r.Inotherwords,
since the trajectories are continuous functions of Rf, the value assumed at Rf ∼r (i.e. at the end of the strongly correlated
regime)represents the initialcondition fortherandom walkatsmallRf (i.e.inthe independent regime).Thisderives froma
sortofnaturalpeak-background split:byevaluatingthedensityfluctuationfieldintwopointsseparatedby r onesamplesthe
samebackground value(obtained bysummingupalltheFouriercomponents correspondingtowavelenghts largerthan r)but
differenthighfrequencymodes.
(1996)andNagashima&Gouda(1997)analysedthepossiblechangesinducedbyspatialcorrelationsonthePSmass
function.
TheproblemwithinthestandardMonteCarloapproachisthateachtrajectoryistreatedasbeingindependentof
alltheothers.InageneralGaussianfield,however,thissituationisnotrealizedand,forpower-spectraofcosmological
interest, the density fluctuations turn out to be strongly correlated on small scales. This statistical property can be
clearly detected by looking at a pair of trajectories associated to two neighbouring points. In Fig. 1 we show two
(cid:13)c 0000RAS,MNRAS000,000–000
4 C. Porciani, S. Matarrese, F. Lucchin & P. Catelan
trajectories obtainedbysmoothingacorrelated fieldwithasharpk-spacefilter.Whenthesmoothinglengthismuch
larger than the separation between the points in which the trajectories are computed the density contrast assumes
valuesthatarepractically coincident.Ontheotherhand,for filteringradiimuchsmaller thanthelag thatseparates
the points, the trajectories become indeed uncorrelated. The transition from the completely correlated situation to
the totally uncorrelated one is continuous, i.e. the trajectories become less and less correlated as the smoothing
length decreases. For comparison, in Fig. 1 we also show the same trajectories modified by artificially removing any
correlation between them. It is clear that the joint distribution of first upcrossings at different points is heavily
influenced bythe correlations of theunderlyingdensity fluctuation field.
Inthispaperwewilli)computethejointtwo-pointdistributionofthefirstupcrossingeventsandusethesolution
to obtain an approximate analytical expression for the halo-halo correlation function, ii) numerically solve suitable
setsofspatially correlated Langevinequationstodirectly computethehalotwo-point function atdifferentlags. Ina
forthcoming paper we will consider theeffects of correlations on themerging histories of dark matter haloes.
The plan of the paper is as follows. In Section 2 we present a straightforward extension of the excursion set
approach and review the standard application of this method which allows to determine the mean mass function
of DM haloes. In Section 3 we study the Lagrangian clustering of haloes identified by our scheme and obtain the
halo-halocorrelation functionbothfrom anapproximateanalytictreatmentandfrom anumericalintegration oftwo
spatially correlated Langevin equations.Section4containsadiscussion ofourresultsandhintsatsome applications
of our algorithm to generate spatially correlated halo merger trees.
2 EXCURSION SET APPROACH
In this section we review the formulation of the excursion set approach, mostly following the treatment given by
BCEK.
2.1 Langevin equation
Let us assume that thelinear density fluctuations form a homogeneous and isotropic Gaussian random field δ(x,z),
withxtheLagrangiancomovingpositionandztheredshift.Inlineartheoryonecanfactorouttheredshiftdependence
as D(z)δ(x), where D(z) is the linear growth factor (e.g. Peebles 1980). With the normalization D(z =0)=1, δ(x)
becomes themass density fluctuation linearly extrapolated to thepresent epoch.
ThestatisticalpropertiesoftheGaussian fieldδ(x)arecompletelyspecifiedbythetwo-pointfunctioninFourier
space, which is related to the power-spectrum P(k) by δ˜(k )δ˜(k ) = (2π)3δ (k +k )P(k ). The symbol δ
h 1 2 i D 1 2 1 D
representstheDiracdeltafunction,andthebrackets denoteensembleaveraging.OurFouriertransformconvention
h·i
is δ˜(k)= dxδ(x)eik·x.
Wewantnowtostudythestatistical propertiesofthedensityfluctuationfieldatsomeresolution scale R .This
R f
is introduced by convolving δ(x) bysome filterfunction W(x′ x,R ),
f
| − |
δ(x,R )= dx′W(x′ x,R )δ(x′)= 1 dkW(kR )δ˜(k)e−ik·x , (1)
f | − | f (2π)3 f
Z Z
where W is the Fourier transform of the filter. At eacfh point x the smoothed field represents the weighted average
of δ(x) over a spherical region of characteristic dimension R centred in x. The detailed properties of δ(x,R )
f f
f
clearly depend upon the specific choice of window function. The most commonly used smoothing kernels are the
top-hatfilterW (x,R )=3Θ(R x)/4πR3, whereΘ(x)is theHeaviside step function,and theGaussian one
TH | | f f −| | f
W (x,R )=(2πR2)−3/2exp( x2/2R2).Recently,forconvenienceofanalysis,top-hatfilteringhasbeenalsoapplied
G f f − f
inmomentumspaceW (k,R )=Θ(k k),wherek =1/R andk = k .Thiskernelisgenerallycalledsharp
SKS f f f f f f
− | |
k-spacefilter.Whileitiseasytoassociateamasstorealspacetop-hatfilteringM (R )=4πρ R3/3,thereisalways
TH f b f
f
abitofarbitrarinessinassigningamasstotheotherwindowfunctions.Themostcommonprocedureistomultiplythe
average density by the volume enclosed by thefilter, obtaining M (R )=(2π)3/2ρ R3 and M (R )=6π2ρ k−3
G f b f SKS f b f
(cid:13)c 0000RAS,MNRAS000,000–000
Excursion set approach to halo clustering 5
(Lacey & Cole 1993). An alternative procedure, originally introduced by Bardeen et al. (1986), corresponds to the
choice M (R ) = 4πρ R3 /3, where R is chosen by requiring σ2 (R ) = σ2 (R ), and similarly for the
SKS f b TH TH SKS f TH TH
Gaussian filter. In this way one obtains good agreement with numerical simulations of clustering growth (Lacey &
Cole 1994).
Inhierarchicalbottom-upmodelsofstructureformationateachepochmassfluctuationshavetypicallyundergone
non-linearcollapseuptosomescaleR ,andthefirstobjectsthatformarethoseoflowermass.Haloesoflargermass
f
arise by themerging and accretion of subunits.The mass distribution derivingfrom this hierarchical halo formation
process has been successfully modeled by Peacock & Heavens(1990), Cole (1991), BCEK and Lacey & Cole (1993).
In order to mimic the accretion of matter all these models consider a full hierarchy of decreasing resolution scales
R . The effect of varyingR can be obtained bydifferentiating eq. (1)
f f
∂δ(x,Rf) = 1 dkδ˜(k)∂W(kRf)e−ik·x η(x,R ). (2)
∂R (2π)3 ∂R ≡ f
f f
Z
f
Thishastheform ofaLangevinequation,whichshowshowaninfinitesimalchangeoftheresolution scale R affects
f
the value of the density fluctuation field δ(x,R ) in the given position x through the action of the stochastic force
f
η(x,R ). In the limit R one has δ(x;R ) 0, which can be adopted as initial condition for our first-order
f f f
→ ∞ →
stochastic differential equation. Thus, by solving it, we can associate to each point x a trajectory δ(x,R ) obtained
f
byvaryingtheresolution scale R .Noticethat,sinceeq.(2) islinear, η(x,R ) isalso azeromean Gaussian random
f f
field, uniquelyspecified by its auto-correlation function
∞
1 ∂W(k R ) ∂W(k R )
η(x ,R )η(x ,R ) = dk k2P(k ) 1 f1 1 f2 j (k r), (3)
h 1 f1 2 f2 i 2π2 1 1 1 ∂R ∂R 0 1
Z0 f1 f2
f f
where r = x x and j (x) is the zeroth order spherical Bessel function. Trajectories associated to different
1 2 0
| − |
neighbouring points will be statistically influenced by the correlation properties of the force η(x,R ), i.e. of the
f
underlying Gaussian field δ(x). On the other hand the coherence of each trajectory along the R direction depends
f
exclusivelyontheanalyticformofthefilterfunctionandvanishesforthesharpk-spacewindow(BCEK).Withsuch
afilter,bydecreasingthesmoothinglengthoneaddsupanewsetofFouriermodesoftheunsmootheddistributionto
δ(x,R ).ForaGaussianfieldthisiscompletelyindependentofthepreviousincrements,andeachtrajectoryδ(x,R )
f f
becomes a Brownian random walk.
Inthecaseofsharpk-spacefilteringthenotation greatly simplifiesifweuseastimevariablethevarianceofthe
filtered density field, Λ σ2(k )= δ(k )2 =(2π2)−1 kf dkk2P(k). In such a case the stochastic process reduces
≡ f h f i 0
to a Wiener one, namely
R
∂δ(x,Λ)
=ζ(x,Λ), (4)
∂Λ
with ζ(x,Λ) =0 and
h i
ζ(x,Λ )ζ(x,Λ ) =δ (Λ Λ ) (5)
h 1 2 i D 1− 2
[see eq. (17) for the spatial correlation]. In the following we will adopt Λ as time variable, unless explicitly stated.
The solution of the Langevin equation (4) in an arbitrary point of space (the position x is here understood), with
the initial condition δ(Λ=0)=0, is simply δ(Λ)= ΛdΛ′ζ(Λ′). By ensemble averaging this expression one obtains
0
δ(Λ) =0and δ(Λ )δ(Λ ) =min(Λ ,Λ ),whichuniquelydeterminetheGaussiandistributionofδ(Λ).Moredetails
h i h 1 2 i 1 2 R
can befound in Porciani et al. (1996).
2.2 Press-Schechter mass function
Press & Schechter (1974) proposed a simple model to compute the comoving number density of collapsed haloes
directlyfromthestatisticalpropertiesofthelinearGaussiandensityfield(seealsoDoroshkevich1967).Accordingto
the PS theory a patch of fluid is part of a collapsed region of scale larger than M(R ) if the value of the smoothed
f
lineardensitycontrastonthesamescaleexceedsathresholdt .Theideaistouseaglobalthresholdinordertomimic
f
non-linear dynamical effects ending up to halo collapse and virialization. An exact value for t can be obtained by
f
(cid:13)c 0000RAS,MNRAS000,000–000
6 C. Porciani, S. Matarrese, F. Lucchin & P. Catelan
describing theevolution of densityperturbationsaccording tothespherical top-hatmodel. Inthiscase afluctuation
of amplitude δ will collapse at a redshift z such that δ(x)=t δ /D(z ). In the Einstein de Sitter universe and
f f c f
≡
duringthematterdominatederathecriticalvalueδ doesnotdependonanycosmological parameterandisgivenby
c
δ =3(12π)2/3/20 1.686, whereas, forgeneral cosmologies, itshowsaweak dependenceonthevalueof thedensity
c
≃
parameter, thecosmological constant, theHubbleconstant, thuson redshift (e.g. Lacey & Cole 1993).
ThePS formulafor themass functionhasbeen thoroughlycompared with theoutcomes ofN-bodysimulations,
finding general good agreement (e.g. Efstathiou et al. 1988; Gelb & Bertschinger 1994; Lacey & Cole 1994). As
mentioned in theIntroduction,themodel is intrinsically flawed bythecloud-in-cloud problem, namely thefact that
afluctuationonagivenscalecancontainsubstructuresofsmallerscalesandthesamefluidelementscanbeassigned,
accordingtothePScollapsecriterion,tohaloesofdifferentmass.Moreover,inahierarchicalscenario,oneexpectsto
find all the mass collapsed in objects of some scale, while the PS model can account only for half of it: this problem
is intimately related to the fact that in a Gaussian field only half volume is overdense. Press and Schechter faced
the problem by simply multiplying their result by a fudge factor of 2. In this section we review how the Langevin
equations introduced above can be used to extend the PS theory in such a way to solve both problems. In the next
sections we will use thesame approach to computethe halo-halo correlation function.
Thesolutionofthecloud-in-cloudproblemhasbeengivenbyPeacock&Heavens(1990),Cole(1991)andBCEK.
Their approach consists in considering at any given point the trajectory δ(R ) as a function of the filtering radius,
f
and then determining the largest R at which δ(R ) upcrosses the threshold t (z ) corresponding to the formation
f f f f
redshift z . This determines the largest mass collapsed at that point, all sub-structures having been erased. So,
f
the computation of the mass function is equivalent to calculating the fraction of trajectories that first upcross the
thresholdt asthescaleM decreases.ThesolutionoftheproblemisenormouslysimplifiedforBrowniantrajectories,
f
that is for sharp k-space filtered density fields. In such a case one only has to solve the Fokker-Planck equation for
theprobability density (δ,Λ)dδ that the stochastic process at Λ assumes a valuein theinterval δ,δ+dδ,
W
∂ (δ,Λ) 1 ∂2 (δ,Λ)
W = W , (6)
∂Λ 2 ∂δ2
with the absorbing boundary condition (t ,Λ) = 0 and initial condition (δ,0) = δ (δ). The solution is well
f D
W W
known (Chandrasekhar 1943)
1 δ2 (δ 2t )2
W(δ,Λ;tf)dδ= √2πΛ exp −2Λ −exp − −2Λf dδ. (7)
(cid:20) (cid:18) (cid:19) (cid:18) (cid:19)(cid:21)
Defining S(Λ,t )= tf dδ (δ,Λ,t ) assurvivalprobabilityofthetrajectories, oneobtainsthedensityprobability
f −∞ W f
distribution of first-crossing variances by differentiation
R
∂S(Λ,t ) ∂ tf 1∂ (δ,Λ,t ) tf t t2
P1(Λ)=− ∂Λ f =−∂ΛZ−∞dδW(δ,Λ;tf)=(cid:20)−2 W ∂δ f (cid:21)−∞ = √2πfΛ3 exp(cid:18)−2Λf (cid:19) . (8)
The function (Λ)dΛ yields the probability that a realization of the random walk is absorbed by the barrier in
1
P
the interval (Λ,Λ+dΛ) or, by the ergodic theorem, the probability that a fluid element belongs to a structure with
mass in the range [M(Λ+dΛ),M(Λ)]. Finally, the comoving number density of haloes with mass in the interval
[M,M +dM] collapsed at redshift z is
f
ρ dΛ
n(M,z )dM = b (Λ) dM . (9)
f M P1 dM
(cid:12) (cid:12)
Inserting the expression of(cid:12) 1(Λ(cid:12)) of eq. (8) in the latter equation one obtains the well-known PS expression for the
(cid:12)P (cid:12)
mass function
ρ t (z ) 1 dlnΛ t (z )2
n(M,z )dM = b f f exp f f dM . (10)
f √2π M2 Λ(M) dlnM −2Λ(M)
(cid:12) (cid:12) (cid:18) (cid:19)
(cid:12) (cid:12)
The original fudge factor of 2pof the P(cid:12)S appr(cid:12)oach is now naturally justified. A generalization of this formalism to
simple cases of non-Gaussian initial conditions has been given by Porciani et al. (1996).
Previous investigations (e.g. Peacock & Heavens 1990; BCEK) have shown that only for sharp k-space filtering
itispossibletowriteananalyticformulaforthemassfunctionobtainedfrom theexcursionsetapproach.Numerical
(cid:13)c 0000RAS,MNRAS000,000–000
Excursion set approach to halo clustering 7
solutionsofthecloud-in-cloudproblemwithphysicallymoreacceptablesmoothingkernelslikeGaussian andtop-hat
result in mass functions that are a factor of two lower in the high-mass tail and have different small-mass slopes
compared with the PS formula. The standard interpretation of this result is that the excursion set method is not
reliable for M ≪M∗, where M∗, definedby Λ(M∗)=t2f, is thetypical mass collapsing at zf.
3 LAGRANGIAN CLUSTERING
In this section we show how the excursion set approach can be extended to derive the clustering properties of
dark matter haloes. Specifically, we study the evolution of a set of density fluctuation processes at different spatial
locations, asthesmoothing scale progressively shrinks.Thetrajectories turnouttobeintimatelycorrelated andthe
jointdistributionofthefirstupcrossingfilteringradiiisusedtoextracttheLagrangianhalo-halocorrelationfunction.
3.1 Two-point correlation function from joint upcrossing distribution
LetusselecttwopointsinLagrangianspacex ex =x +r.Wewanttostudytheevolutionofthestochasticprocesses
1 2 1
δ (Λ) = δ(x ,Λ) and δ (Λ) = δ(x ,Λ) as Λ varies. In particular, let us suppose that we know the joint probability
1 1 2 2
distribution (Λ ,Λ ;x ,r) of those pairsof variances (Λ ,Λ ) that correspond to thefirst upcrossing scales of the
2 1 2 1 1 2
P
thresholdt bythetwoprocessesδ (Λ)andδ (Λ).Becauseoftheunderlyinghomogeneityandisotropy,theprobability
f 1 2
density cannotdependonthevectorx andontheorientationofr,i.e. (Λ ,Λ ;x ,r)= (Λ ,Λ ;r).Moreover,
2 1 2 1 2 1 2 1 2
P P P
bytheergodictheoremonecanidentify (Λ ,Λ ;r)withtheprobabilitydistributionofthepairs(Λ ,Λ ),obtained
2 1 2 1 2
P
by randomly selecting points spatially separated by the distance r, within a given realization of the density field.
Finally,followingtheargumentsgivenintheprevioussection,wecaninterpret (Λ ,Λ ;r)dΛ dΛ astheprobability
2 1 2 1 2
P
offindingtwopointsseparatedbyrwithintwohaloeswithmassintheintervals(M dM ,M )and(M dM ,M ),
1 1 1 2 2 2
− −
as fixed by the corresponding variance ranges (Λ ,Λ +dΛ ) and (Λ ,Λ +dΛ ). As we will discuss in detail in the
1 1 1 2 2 2
nextsections, theprobabilitydensity (Λ ,Λ ;r)canbeobtained byintegratingthesystem ofcorrelated Langevin
2 1 2
P
equations that describe the evolution of theprocesses δ (Λ) and δ (Λ).
1 2
In order to compute the halo-halo correlation function we adopt the following procedure. A class of objects is
selected by the mass interval corresponding to the Λ-range I [Λ ,Λ ]. The probability of determining two
min max
≡
points separated by r contained within collapsed objects of class I is
(r)= dΛ dΛ (Λ ,Λ ;r). (11)
II 1 2 2 1 2
P P
ZIZI
Similarly,theprobabilityoffindingapointcontainedinanobjectoftypeI is = dΛ (Λ).Fromthedefinition
PI I P1
of correlation function we then obtain
R
dΛ dΛ (Λ ,Λ ;r)
ξpts(r)= PII 1= I I 1 2P2 1 2 1. (12)
II PI2 − R R IdΛP1(Λ) 2 −
In a similar fashion, consid(cid:2)eRring disjoint(cid:3)classes I1 and I2,
dΛ dΛ (Λ ,Λ ;r)
ξpts (r)= I1 I2 1 2P2 1 2 1. (13)
I1I2 R I1RdΛP1(Λ) I2dΛP1(Λ) −
We stress thRat the quantiRties ξpts(r) and ξpts (r) are the Lagrangian correlations of the mass elements contained in
II I1I2
thecollapsedhaloes.ThisquantityhasbeenrecentlyusedbyBagla(1997)tostudytheevolutionofgalaxyclustering.
However, we are ultimately concerned with the calculation of the halo correlations ξhh(r), so we have to properly
II
weigh the statistical contribution for each extended halo. The problem shows particularly simple if we adopt the
PS theory. Indeed, in this case, the sets of points where the first upcrossings occur at the same Λ are point-like
disconnected regions (see Bond & Myers 1996 for a different approach). Only in a statistical sense they originate
collapsed haloes, eachcontributingby1/V(Λ),whereV(Λ) M/ρ isthetypicalLagrangian volumeofanobject of
b
≡
mass M associated to the variance Λ. Therefore, adopting the same procedure used to calculate the mass function,
(cid:13)c 0000RAS,MNRAS000,000–000
8 C. Porciani, S. Matarrese, F. Lucchin & P. Catelan
themeannumberdensityofcollapsed objects ofscale Λbecomesn(Λ)= (Λ)/V(Λ).Similarly,forthedistribution
1
P
of pairs at distance r,we define
(Λ ,Λ ;r)
n2(Λ1,Λ2;r)= PV2(Λ1)V(2Λ ) . (14)
1 2
The idea is, once again, to allow for the finitesize of thehaloes. Thehalo-halo correlations become, respectively,
ξhh(r)= nII 1 I IdΛ1dΛ2n2(Λ1,Λ2;r) 1, (15)
II n2I − ≡ R R dΛn(Λ) 2 −
I
and (cid:2)R (cid:3)
dΛ dΛ n (Λ ,Λ ,r)
ξhh (r)= I1 I2 1 2 2 1 2 1. (16)
I1I2 R IR1dΛn(Λ) I2dΛn(Λ) −
R R
3.2 Correlated Langevin equations: sharp k-space filtering
In Section 3.1 we showed that, in order to obtain the halo correlation function, it is crucial to know the joint
distribution (Λ ,Λ ;r).This quantitycan beobtained bysolving thesystem ofequations governingtheevolution
2 1 2
P
of the pair of correlated processes δ (Λ) and δ (Λ):
1 2
∂δ (Λ)
1 =ζ (Λ), δ (0)=0,
∂Λ 1 1

 h∂ζδ1∂2(Λ(ΛΛ))i==hζζ22((ΛΛ)),i=0, δζ21(0a)nd=ζ02,Gaussian processes, (17)
′ ′ ′
Thehhζζla11t((tΛΛe))rζζe21q((uΛΛa′t))iiio==n,h∂oζξb2(∂(trΛaΛ;i)nΛζe)2d(δΛDby()Λiin=−trδΛoDd′)(uΛ.cin−gΛsh)a,rp k-space filtering in equation (3), completes the definition of the
stochasticfieldζ(x,Λ)giveninequations(4)and(5).Hereξ(r;Λ)isthelineartwo-pointcorrelation functionforthe
mass density fluctuationssmoothed on thescale R 1/k associated to thevariance Λ
f f
≡
1 kf(Λ)
ξ(r;Λ) dkk2P(k)j (kr) (18)
≡ 2π2 0
Z0
and one has
∂ξ(r;Λ)
=j [k (Λ)r]. (19)
∂Λ 0 f
Byintegratingtheabovedifferentialequationsandaveragingovertheensembleoneobtainstheunconstrainedprob-
ability distribution
1 Λ(δ2+δ2) 2ξ(r;Λ)δ δ
Wr(δ1,δ2;Λ)= 2π Λ2 ξ(r;Λ)2 exp − 12[Λ22 −ξ(r;Λ)2] 1 2 . (20)
− (cid:20) − (cid:21)
This function solveps thetwo-dimensional Fokker-Planckequation associated to theLangevin equation (17), namely
∂ (δ ,δ ;Λ) 1 ∂2 ∂2 ∂ξ(r;Λ) ∂2
Wr ∂1Λ 2 = 2(cid:20)∂δ12 + ∂δ22 + 2 ∂Λ ∂δ1∂δ2(cid:21) Wr(δ1,δ2;Λ), (21)
with initial condition (δ ,δ ,Λ=0)=δ (δ )δ (δ ).The problem of findingthedistribution of thefirst upcross-
Wr 1 2 D 1 D 2
ingsofthethresholdt bythebinaryprocess δ ,δ reducestothatofimposingproperboundaryconditionstothe
f 1 2
{ }
equation (21). As done for the one-dimensional Fokker-Planck equation, we adopt the absorbing barrier approach.
However,inthetwo-dimensionalcaseweareconsidering,thedistribution (Λ ,Λ ;r)cannotbeeventuallyobtained
2 1 2
P
from (δ ,δ ;Λ) simply by differentiation. This is because the whole binary system automatically disappears as
r 1 2
W
(cid:13)c 0000RAS,MNRAS000,000–000
Excursion set approach to halo clustering 9
soon as when one ‘Brownian particle’ is first absorbed. Nonetheless, the joint distribution (Λ ,Λ ;r) can be in
2 1 2
P
principle calculated bya two-step procedure as follows.
Assuming the same initial condition, one solves the Fokker-Planck equation with absorbing barriers at δ = t
1 f
and δ = t , thus finding the survival probability density for the pairs which have never crossed the thresholds.
2 f
Havingfound this quantityone can compute theprobability current through each point,
1 ∂ ∂ξ(r;Λ) ∂ ∂ξ(r;Λ) ∂ ∂
J(δ1,δ2;Λ)=−2 ∂Wδ1r + ∂Λ ∂Wδ2r , ∂Λ ∂Wδ1r + ∂Wδ2r . (22)
(cid:18) (cid:19)
On a boundary wall, e.g. δ =t , where (t ,δ ;Λ)=0 (implying ∂ /∂δ =0), this reducesto
1 f r f 2 r 2
W W
1 ∂ ∂ξ(r;Λ) ∂
J(tf,δ2;Λ)=−2(cid:18) ∂Wδ1r(cid:12)tf , ∂Λ ∂Wδ1r(cid:12)tf(cid:19) . (23)
The flux through any poi(cid:12)nt of the barrier δ =(cid:12) t is then given by the scalar product J n, where n (1,0) is the
(cid:12) 1 (cid:12) f
· ≡
unit vector perpendicularto the absorbing wall,
1 ∂
Fr(tf,δ2;Λ)=−2 ∂Wδ1r tf . (24)
(cid:12)
Thequantity (t ,δ ;Λ(cid:12))dδ representstheprobability that thepair ofprocesses (δ ,δ ) leave thepermitted region
r f 2 (cid:12) 2 1 2
F
passingthroughthe‘gate’[(t ,δ ),(t ,δ +dδ )]atthetimeΛ.Thisfluxcontainsalltheinformationweneedforthe
f 2 f 2 2
computation of (Λ ,Λ ;r).Infact, onceδ hascrossed thebarrieratΛ ,forΛ>Λ weareinterestedin studying
2 1 2 1 1 1
P
only the evolution of the surviving process δ up to its first upcrossing through the boundary δ = t . Therefore,
2 2 f
since we are considering Brownian trajectories, free of correlations along the Λ axis, for Λ > Λ the evolution of
1
δ is simply governed by its own Langevin equation, and its probability distribution can be calculated from the
2
one-dimensionalFokker-Planckequation(6),withabsorbingboundaryδ =t ,butwithinitialcondition(atΛ=Λ )
2 f 1
δ2∗ δ2(Λ1δ1 =tf).Thus, bysimply modifying eq.(8),wefindthat thedistribution of thevariances Λ2 associated
≡ |
to first upcrossing eventsof the threshold by theprocess δ , given that δ upcrossed thecritical level at Λ , is:
2 1 1
P1(Λ2−Λ1,tf −δ2∗)= √2π(t(fΛ−2−δ2Λ∗)1)3/2 exp(cid:20)−2(t(fΛ−2−δ2Λ∗)12)(cid:21) . (25)
The joint distribution (Λ ,Λ ;r) is eventually obtained by a convolution
2 1 2
P
tf tf
(Λ ,Λ ;r)= dδ (t ,δ ;Λ ) (Λ Λ ,t δ )+ dδ (δ ,t ;Λ ) (Λ Λ ,t δ ), (26)
2 1 2 2 r f 2 1 1 2 1 f 2 1 r 1 f 2 1 1 2 f 1
P F P − − F P − −
Z−∞ Z−∞
where the first and second integrals on the r.h.s. represent, respectively, the contributions of those pairs for which
Λ Λ and Λ <Λ .
2 1 2 1
≥
However, thisformal expression isuseless unlesswe solvefor theprobability density . Thereare twocases in
r
W
whichthecalculation of istrivial: onefor r andtheotherforr 0.Thegeneralcaseof finitenon-zerolag
r
W →∞ →
r is far more complex, as we will discuss.
3.2.1 Perfectly uncorrelated processes
At infinitelag (r ) the two processes become independentand thesolution of theFokker-Planck equation is
→∞
∞(δ1,δ2;Λ)= (δ1,Λ;tf) (δ2,Λ;tf), (27)
W W W
where (δ,Λ;t ) denotes the probability distribution for a one-dimensional process with absorbing boundary at t
f f
W
givenineq.(7).Thissolutionconsistsofalinearsuperpositionoffourunconstrainedindependentdensitydistributions
deriving from different initial conditions. In practice one has to consider the ‘real’ initial distribution δ (0,0), an
D
image source δ (t ,t ) and two image sinksδ (t ,0) and δ (0,t ), that is to say
D f f D f D f
∞(δ1,δ2;Λ)= ∞(δ1,δ2;Λ) ∞(δ1 2tf,δ2;Λ) ∞(δ1,δ2 2tf;Λ)+ ∞(δ1 2tf,δ2 2tf;Λ), (28)
W G −G − −G − G − −
where (δ ,δ ;Λ) [2πΛ]−1exp[ (δ2+δ2)/2Λ] is the solution of the two-dimensional Fokker-Planck equation with
G 1 2 ≡ − 1 2
natural boundary conditions: limδi→∞G(δ1,δ2;Λ)=0,i=1,2. Obviously,using eq. (26) we obtain
(cid:13)c 0000RAS,MNRAS000,000–000
10 C. Porciani, S. Matarrese, F. Lucchin & P. Catelan
(Λ ,Λ ;r )= (Λ ,t ) (Λ ,t ), (29)
2 1 2 1 1 f 1 2 f
P →∞ P P
that,insertedinequations(12)and(13)orinequations(15)and(16),gives,asexpectedforinfinitelag,ξpts=ξhh=0.
3.2.2 Perfectly correlated processes
When r 0 the two processes become more and more correlated so that, in the limit, we end up with a one-
→
dimensional problem. In this case the solution of theFokker-Planck equation is
(δ ,δ ;Λ)= (δ +δ ,4Λ;t )δ (δ δ ), (30)
W0 1 2 W 1 2 f D 1− 2
where (δ,Λ;t ) denotestheprobability distribution given in eq.(7). Expandingthisexpression as a superposition
f
W
of Green’s functions, by using themethod of image sources, we obtain
(δ ,δ ;Λ)= (δ ,δ ;Λ) (δ 2t ,δ 2t ;Λ), (31)
0 1 2 0 1 2 0 1 f 2 f
W G −G − −
where (x,y;Λ) (x,Λ;t )δ (x y).Inthiscase, for thejoint distribution offirst upcrossing variances, wefind
G0 ≡W f D −
(Λ ,Λ ;r 0)= (Λ ,t )δ (Λ Λ ). (32)
P2 1 2 → P1 1 f D 1− 2
Consequently,we obtain ξpts [1/ dΛ (Λ)] 1 and ξhh dΛP(Λ)/V(Λ)2/[ dΛP(Λ)/V(Λ)]2 1. Simi-
II → I P − II →{ I I }−
larly, ξpts 0 and ξhh =0.
I1I2 → I1I2 R R R
3.2.3 Approximate general solution
In the limiting cases just discussed we were able to account for the boundary conditions by writing the solution in
termsofasuperpositionofimage-distributions,eachonesolvingtheFokker-Planckequation.Regrettably,theposition
and the sign of the image sources of probability come out dependent on the correlation between the processes (i.e.
on r). This fact suggests that we cannot write a general analytical solution by simply applying the image method.
However,wewillbuildhereasimplefunctionthatsatisfiestheabsorbingboundaryconditionsbeingalsoanaccurate
approximation for the solution of the correlated diffusion equation. In the next section we will test the accuracy of
this solution against the numerical one.
It is evident that, for small separation r R (we remind the reader that we change the smoothing radius R
f f
≪
at fixed separation r), the perfectly correlated solution will represent a very good approximation to the true one.
This can beeasily deducedbythefollowing argument. Thetwo Gaussian stochastic processes Σ(Λ)=δ (Λ)+δ (Λ)
1 2
and ∆(Λ) = δ (Λ) δ (Λ) are statistically independent, i.e. Σ(Λ)∆(Λ) = 0. The variances of their unconstrained
1 2
− h i
probability distributions are, respectively, σ2 = 2[Λ+ξ(r;Λ)] and σ2 = 2[Λ ξ(r;Λ)]. Therefore, for r R
Σ ∆ − ≪ f
[i.e. for Λ σ2(r), where σ2(r) denotes the variance of the mass density fluctuations smoothed on the scale r]
≪
where ξ(r;Λ) Λ, we have σ2 0 and the probability distribution of the variable ∆ is practically a Dirac delta
≃ ∆ ≃
functioncentredonthezerovalue.Thiscorrespondstotheperfectlycorrelatedsituation.However,thisregimeisnot
interesting for computing the halo correlation function since the condition r R implies that the points in which
f
≪
we follow thetrajectories that upcross thethreshold are involved in the collapse of the same halo.
Ontheotherhand,forr R (i.e.forΛ σ2(r)),wecanreplaceineq.(21)ξ(r;Λ)withtheunsmoothedlinear
f
≫ ≫
mass-correlation ξ (r).Therefore, eq.(21)simplybecomestheuncorrelatedtwo-dimensionaldiffusionequationthat
m
can be easily solved using the image method. Our ansatz for the full (δ ,δ ;Λ) is then obtained by keeping the
r 1 2
W
analytic form of the solution just obtained for r R but inserting in it the correlation function ξ(r;Λ) to replace
f
≫
its large lag limit ξ (r). Therefore, we have
m
(δ ,δ ;Λ)= +(δ ,δ ;Λ) −(δ 2t ,δ ;Λ) −(δ ,δ 2t ;Λ)+ +(δ 2t ,δ 2t ;Λ), (33)
Wr 1 2 Gr 1 2 −Gr 1− f 2 −Gr 1 2− f Gr 1− f 2− f
where
± = 1 exp Λ(δ12+δ22)−2[±ξ(r;Λ)]δ1δ2 . (34)
Gr 2π Λ2 ξ(r;Λ)2 − 2[Λ2 ξ(r)2]
− (cid:20) − (cid:21)
p
(cid:13)c 0000RAS,MNRAS000,000–000