Mon. Not. R. Astron. Soc. 000, 1–15 (2010)

Printed 22 March 2011

(MN LATEX style file v2.2)

arXiv:1012.4833v2 [astro-ph.CO] 21 Mar 2011

Modeling the clustering of dark-matter haloes in resummed perturbation theories A. Elia1⋆ , S. Kulkarni2, C. Porciani1, M. Pietroni3, S. Matarrese4,3 1

Argelander Institut f¨ ur Astronomie der Universit¨ at Bonn, Auf dem H¨ ugel 71, D-53121 Bonn, Germany Center for Theoretical Physics & Physikalisches Institut, Universit¨ at Bonn, Nußallee 12, D-53115 Bonn, Germany 3 INFN, Sezione di Padova, Via Marzolo 8, I-35131 Padova, Italy 4 Dipartimento di Fisica “G. Galilei”, Universit` a degli Studi di Padova, Via Marzolo 8, I-35131 Padova, Italy

2 Bethe

22 March 2011

ABSTRACT

We address the issue of the cosmological bias between matter and galaxy distributions, looking at dark-matter haloes as a first step to characterize galaxy clustering. Starting from the linear density field at high redshift, we follow the centre of mass trajectory of the material that will form each halo at late times (proto-halo). We adopt a fluid-like description for the evolution of perturbations in the proto-halo distribution, which is coupled to the matter density field via gravity. We present analytical solutions for the density and velocity fields, in the context of renormalized perturbation theory. We start from the linear solution, then compute one-loop corrections for the propagator and the power spectrum. Finally we analytically resum the propagator and we use a suitable extension of the time-renormalization-group method (Pietroni 2008) to resum the power spectrum. For halo masses M < 1014 h−1 M⊙ our results at z = 0 are in good agreement with N-body simulations. Our model is able to predict the halo-matter cross spectrum with an accuracy of 5 per cent up to k ≈ 0.1 h Mpc−1 approaching the requirements of future galaxy redshift surveys. Key words: cosmology: theory, large-scale structure of Universe– galaxies: haloes – methods: analytical, N-body simulations.

1

INTRODUCTION

Redshift surveys have shown that the clustering properties of galaxies strongly depend on their luminosity, color and morphological (or spectral) type (e.g. Norberg et al. 2002; Zehavi et al. 2004). This indicates that galaxies do not perfectly trace the distribution of the underlying dark matter, a phenomenon commonly referred to as ‘galaxy biasing’. Its origin lies in the details of the galaxy formation process which is shaped by the interplay between complex hydrodynamical and radiative processes together with the darkmatter driven formation of the large-scale structure. Attempts to infer cosmological parameters from galaxy clustering studies are severely hampered by galaxy biasing. A number of theoretical arguments and the outcome of numerical simulations both suggest that, on sufficiently large scales, the power spectra of galaxies and matter should be proportional to each other: Pg = b21 Pm where the linear bias factor b1 depends on galaxy type but is generally scale independent (e.g. Coles 1993; Mann et al. 1998). Similarly, to model higher-order statistics, such as the galaxy bispectrum,

⋆

E-mail: [email protected]

it is generally assumed that galaxy biasing is a local process 2 such that δg = b1 δm + b2 δm /2 + . . . where δg and δm are the (smoothed) galaxy and dark-matter density contrast, respectively (Fry & Gazta˜ naga 1993). However, the reliance of these phenomenological approximations limits cosmological studies to very large scales whereas data with better signal-to-noise ratio are already available on much smaller scales. Moreover, future studies of baryonic acoustic oscillations (e.g. Sugiyama 1995; Eisenstein & Hu 1998) will require measurements of the matter power-spectrum with percent or even sub-percent accuracy in order to shed new light on the source of cosmic acceleration. Understanding and controlling the effects of galaxy biasing with this precision will be challenging. All this provides a very strong motivation for developing more accurate (and physically driven) models of galaxy biasing. A number of authors have used the power spectrum statistics to explore the scale dependence of galaxy biasing based on numerical simulations (Cole et al. 2005; Seo & Eisenstein 2005; Huff et al. 2007; Manera et al. 2010; Manera & Gazta˜ naga 2009; Montesano et al. 2010) or analytical calculations (Seljak 2001; Schulz & White 2006; Guzik et al. 2007; Smith et al. 2007) stemming from either

2

A. Elia, S. Kulkarni, C. Porciani, M. Pietroni, S. Matarrese

perturbation theory or the halo model for the large-scale structure (see Cooray & Sheth 2002 for a review). The general picture is that galaxy biasing is expected to be scale dependent (i.e. Pg (k) = b(k)2 Pm (k)) and the functional form of b(k) can sensibly depend on the selected tracer of the large-scale structure. Since galaxies are expected to form within dark-matter haloes, understanding the clustering properties of the haloes is a key step to accurately model galaxy biasing. This is a much simpler problem, considering that dark-matter haloes form under the sole action of gravity. It is in fact expected that long-wavelength density fluctuations modulate halo formation by modifying the collapse time of localized short-wavelength density peaks (Bardeen et al. 1986; Cole & Kaiser 1989). This argument (known as the peakbackground split) predicts that, on large scales, the halo overdensity δh = b δm where the bias coefficient b varies with the halo mass (Mo & White 1996). The numerical value of the bias coefficient is determined by two different occurrences: first, haloes form out of highly biased regions in the linear density field (Kaiser 1984; Porciani et al. 1999) and, second, they move over time as they are accelerated towards the densest regions of the large-scale structure by gravity (Mo & White 1996). These two phenomena generally go under the name of “Lagrangian biasing” and “Lagrangian to Eulerian passage”, respectively. Mo & White (1996) dealt with the second problem by assuming that long-wavelength density perturbations evolve according to the spherical tophat model. A more sophisticated generalization of the peakbackground split has been presented by Catelan et al. (1998) who assumed that also the large-scale motion of the density “peaks” is fully determined by the long-wavelength component of the density field. Since the halo population and the matter feel the same large-scale gravitational potential, their density fluctuations are strongly coupled and their time evolution must be solved simultaneously. This makes the process of halo biasing non-linear and non-local even if one starts from a linear and local Lagrangian biasing scheme (Catelan et al. 1998; Matsubara 2008). The bispectrum can be used to test this model against the standard Eulerian local biasing scheme (Catelan et al. 2000). In this paper, we present a novel and very promising approach to model the clustering of dark matter haloes. Adopting the formalism by Catelan et al. (1998) combined with a non-local Lagrangian biasing scheme for the haloes (Matsubara 1999), we simultaneously follow the growth of perturbations in the matter and in the halo distribution over cosmic time. We present perturbative solutions for the corresponding overdensity and velocity fields and we are able to resum the perturbative series in the limit of large wavenumbers. Moreover, we write down a system of equations for the power spectra Pm and Ph using the time-renormalizationgroup (TRG) approach by Pietroni (2008) and numerically integrate them. Our results are in excellent agreement with the output of a high-resolution N-body simulation, showing an improvement over linear theory, and we are able to predict the matter-halo cross spectrum with a precision within 5 per cent for k < 0.15 h Mpc−1 . Related work has been very recently presented by Desjacques et al. (2010) who computed the two-point correlation function of linear density peaks and followed its time evolution assuming that peaks move according to the

Zel’dovich approximation. For massive haloes this results in a scale-dependent bias (with variations of ∼ 5 per cent) on the scales relevant for baryonic-oscillation studies. Contrary to their approach, we do not deal with a point process but describe large-scale fluctuations in the distribution of darkmatter haloes as perturbations in a continuous fluid. On the other hand, we account for the full gravitational motion of the haloes and do not rely on simplified dynamical models as like as the Zel’dovich approximation. The structure of our paper is as follows. In Section 2 we present our model for the joint evolution of the matter and halo power spectra. The initial conditions for our evolutionary equations are discussed in Section 3. The solution of the linearized equations is presented in Section 4 where we also quantify the importance of the halo velocity bias. Using a perturbative technique, in Section 5 we compute analytic solutions for the propagator of perturbations (the two-time correlator). We derive 1-loop corrections and, in the limit of large wavenumbers, the fully resummed propagator. The discussion in Sections 5.1 and 5.2 is very technical and the less experienced readers can safely skip it without compromising understanding of the remainder. In Section 6, we numerically integrate the full equations for the evolution of halo and matter power spectra in the TRG formalism. We then compare the results against the outcome of a high-resolution N-body simulation. Finally, in Section 7 we conclude.

2

THE MODEL

2.1

Dynamics of gravitational instability

The large-scale structure observed today in the universe is believed to be the result of gravitational amplification of primordial fluctuations, caused by the interaction among cold dark matter (CDM) particles. If we denote δm as the matter density contrast and v as the velocity, the Eulerian dynamics of a system of such particles, which interact only via gravity, is ruled by a set of three equations (continuity, Euler and Poisson) that in a ΛCDM model reads: ∂ δm + ∇ · [(1 + δm )v] = 0 , ∂τ ∂v + H v + (v · ∇)v = −∇φ , ∂τ 3 ∇2 φ = H2 Ωm δm , (1) 2 where τ is the conformal time. If we define the velocity divergence θ(x, τ ) = ∇ · v(x, τ ) and switch to Fourier space, the equations in (1) become: ∂ δm (k, τ ) + θ(k, τ ) Z∂ τ

+

d3 q d3 p δD (k − q − p)α(q, p)θ(q, τ )δm (p, τ ) = 0 ,

∂ θ(k, τ ) 3 + H θ(k, τ ) + H2 Ωm (τ )δm (k, τ ) ∂ τ 2 Z

+

d3 q d3 p δD (k − q − p)β(q, p)θ(q, τ )θ(p, τ ) = 0 , (2)

where

Halo clustering with resummed perturbation theories α(q, p) =

(p + q) · q , q2

β(q, p) =

(p + q)2 p · q , 2 p2 q 2

(3)

Equations (2) can be written in a compact form if we define a new time variable η ≡ ln(D+ /D+in ), being D+in the growth factor at an early epoch, and a doublet ϕa (a = 1, 2) δm (k, η) ϕ1 (k, η) , (4) ≡ e−η −θ(k, η)/(Hf+ ) ϕ2 (k, η) with f+ ≡ d ln D+ /d ln a. The velocity divergence is scaled such that it has the same dimension of the density contrast and in the linear regime −θ(k, η)/(Hf+ ) ≈ δm (k, η), i.e. ϕ1 (k, 0) = ϕ2 (k, 0). The system is therefore ∂η ϕa (k, η) = −Ωab (η)ϕb (k, η) +eη γabc (k, −p, −q)ϕb (p, η) ϕc (q, η),

(5)

where sum over repeated indices and integration over repeated momenta are understood. The vertex function γabc (k, p, q) (a, b, c, = 1, 2) has only three non-vanishing elements 1 γ121 (k, p, q) = δD (k + p + q) α(p, q) , 2 γ222 (k, p, q) = δD (k + p + q) β(p, q) , (6) and γ112 (k, q, p) = γ121 (k, p, q), with δD the Dirac-delta distribution. All the information about the cosmological model is contained in the matrix 1 −1 (7) Ω(η) = 3Ωm 3Ωm , − 2 2 2f+ 2f+ that, in the following, will be considered as a constant matrix, approximating Ωf m 2 ≈ 1. This ratio is indeed very close +

to unity for most of the history of the Universe. Making this approximation, one is modifying the behavior of the decaying mode, while the growing one is left unaltered. It has been shown in Pietroni (2008) that it affects the matter power spectrum at z = 0 at a less than percent level up to k ≈ 0.3 h Mpc−1 . The power spectrum, defined by an ensemble average, in this notation is a 2 × 2 matrix hϕa (k; η)ϕb (q; η)i ≡ δD (k + q)Pab (k; η) ,

(8)

and the bispectrum, defined by hϕa (k; η)ϕb (q; η)ϕc (p; η)i ≡ δD (k + q + p)Babc (k, q, p; η) ,

(9)

has 8 components. In the following, we will also consider a different-time two-point correlator, defined as hϕa (k, ηa )ϕb (q, ηb )i ≡ δD (k + q)Pab (k; ηa , ηb ) ,

(10)

which obviously coincides with (8) for ηa = ηb . 2.2

The distribution of dark-matter haloes

Let us consider a set of dark-matter haloes identified at a given redshift zid according to some predefined criterion. The material that forms the haloes can be traced back to its initial location in the linear overdensity field at z → ∞. We dub each of these regions as a proto-halo. In other words, a proto-halo is the Lagrangian region of space that

3

will collapse to form a halo at redshift zid . N-body simulations show that nearly all proto-haloes are simply connected (Porciani et al. 2002) even though this property is not key to our analysis. Let us now follow the evolution of a proto-halo over cosmic time in Eulerian space. Basically its shape and topology will be distorted (proto-haloes will first fragment into smaller substructures that will later merge to form the final halo) and its overall volume will be compressed while its centre of mass will move along a given trajectory determined by the mass density field via gravity. We focus our analysis onto this motion that connects the Lagrangian position of the proto-halo with the Eulerian location of the final halo. On scales much larger than the characteristic size of (and separation between) the proto-haloes, the density fluctuations traced by the centre-of-mass trajectories can be described with a continuous overdensity field δh (x, τ |zid ). Note that while τ labels conformal time along the trajectories, zid is just a tag that identifies the halo population. Unlike real haloes that undergo merging, by construction protohaloes always preserve their identity. Their total number is therefore conserved over time and we can write a continuity equation for their number density: ∂δh + ∇ · [(1 + δh )vh ] = 0 . (11) ∂τ Here the proto-halo density and velocity fields should be intended as coarse grained on a scale of a few times the mean inter-halo separation (so as to suppress discreteness effects as proto-haloes are individually separate units). Strictly speaking, the smoothed velocity field does not obey the Euler-Poisson system in equation (1) due to the presence of the non-linear term (v · ∇)v. In fact, the coarse graining procedure introduces new terms in the fluid equations generated by the degrees of freedom one has integrated out, namely: a velocity dispersion term and a correction to the mean- field gravitational acceleration due to density fluctuations on scales smaller than the smoothing radius (Buchert & Dominguez 2005). On the other hand, it is reasonable to assume that the large-scale motion of the protohaloes is generated by density fluctuations with wavelength larger than the characteristic halo size and is not influenced by perturbations with shorter wavelength. The very same assumption of neglecting the coupling to the small scales is routinely done when one writes equation (1) for the mass density field (see Section 3 in Buchert & Dominguez 2005) albeit adopting much narrower smoothing kernels than for the haloes. With the same spirit, in what follows, we will ignore the extra terms in the fluid equations generated by the coarse graining procedure. This is a working hypothesis which makes the problem mathematically treatable and whose accuracy will be tested by comparing our final results against high-resolution numerical simulations. We therefore write an Euler equation for the proto-halo fluid velocities ∂ vh + H vh + (vh · ∇)vh = −∇φ , (12) ∂τ where the gravitational potential is the same as in equation (1). Note that if vh matches v in the initial conditions then it will always do. On the contrary, any initial velocity bias between proto-halos and matter will be progressively erased by the gravitational acceleration.

4

A. Elia, S. Kulkarni, C. Porciani, M. Pietroni, S. Matarrese

Thus, given suitable initial conditions for δh and vh at τ → 0 (i.e. a prescription for the Lagrangian biasing of proto-haloes), we can in principle use equations (11) and (12) to follow the clustering evolution of the proto-halo population at all times. We are particularly interested in the solution of the fluid equations at the special time τ that corresponds to zid . In fact this solution has a particular physical meaning as it gives the density and velocity fields of the actual dark-matter haloes. 2.3

Growth of matter and halo perturbations

The system (1) is now extended by the inclusion of eq. (11) and eq. (12). We define a quadruplet ϕa (a = 1, 2, 3, 4) δm (k, η) ϕ1 (k, η) ϕ2 (k, η) −η −θ(k, η)/(Hf+ ) (13) , ϕ3 (k, η) ≡ e δh (k, η) −θh (k, η)/(Hf+ ) ϕ4 (k, η) in such a way that eq. (5) still holds, but with indices running from 1 to 4. There are three more non-vanishing elements of the vertex γ343 (k, p, q) = γ334 (k, q, p) = γ121 (k, p, q) and γ444 (k, q, p) = γ222 (k, p, q), and the 4 × 4 Ω matrix is 1 −1 0 0 3 3 0 0 − . 2 2 Ω= (14) 0 1 −1 0 3 3 0 0 − 2 2 From the definitions (8) and (9), with a = 1, 2, 3, 4, we get a 4 × 4 matrix for power spectrum; in the following, we will focus on the matter power spectrum P11 and the matter-halo cross spectrum P13 .

3

INITIAL CONDITIONS

In the previous section, we have presented a model that describes the non-linear evolution of the matter and halo density fields. Given suitable initial conditions, the formal equations we have introduced can be integrated numerically so that to compute the perturbative propagators and the TRG-evolved power spectra. The choice of the initial conditions therefore plays a very important role in our theory and will be the subject of this section. 3.1

N-body simulation

To gain insight into the properties of proto-haloes (and, later, to test our results at z = 0), we use one highresolution N-body simulation extracted from the suite presented by Pillepich et al. (2010). This contains 10243 darkmatter particles within a periodic cubic box with a side of Lbox = 1200h−1 Mpc and follows the formation of structure in a ΛCDM model with Gaussian initial conditions and cosmological parameters: h = 0.701, σ8 = 0.817, ns = 0.96, Ωm = 0.279, Ωb = 0.0462 and ΩΛ = 0.721. We identify dark-matter haloes at z = 0 using the friends-of-friends algorithm with a linking length equal to 0.2 times the mean interparticle distance. We only consider haloes containing more than 100 particles (i.e. with mass

M > 1.24 · 1013 h−1 M⊙ ) and we split them into four mass bins to keep track of their different clustering properties. The corresponding mass ranges and the total number of haloes in each bin are given in Table 1, along with an estimate of the highest wavevector up to which the fluid approximation for haloes holds. This value is determined by the number of haloes we require to be in a volume element to consider them as a fluid, and we set this number to 30. On smaller scales, our assumption breaks down, therefore we will look at results in the specified range, that, of course, decreases as the mass of the haloes increases. In the plots that will be shown in Section 6 the limit to which we can trust our model will be represented by vertical black dotted lines. Proto-haloes are identified by tracing the positions of the particles forming a halo at z = 0 back to the linear density field. The centre of mass of each proto-halo is used as a proxy for its spatial location. Similarly, the mass weighted linear velocity gives the proto-halo velocity. Halo and proto-halo density and momentum fields are computed with the cloud-in-cell grid assignement using a 5123 mesh. Velocity fields are obtained by taking the ratio of the momentum and density distributions (preventively smoothed to preclude the existence of empty cells) as shown in Scoccimarro (2004). Power spectra have been computed using FFT. In order to avoid uncertain shot-noise corrections for the haloes, we only consider their cross spectra with the matter density field. 3.2

Lagrangian halo bias

Concerning the matter density, the initial conditions are given by linear theory which directly follows from the adopted cosmological model (transfer function) and the statistics of primordial perturbations (spectral index, Gaussianity). On the other hand, for the dark-matter halos, we can follow two different approaches: (i) extract the relevant information directly from the simulation or (ii) use a model for the Lagrangian bias of the halos. The latter option offers a number of advantages. First, it allows us to make general predictions independently of the simulation specifics. Second, it allows us to include halo bispectra in our formalism (while it would be extremely demanding and time consuming to compute all possible triangular configurations from the simulation). For these reasons we will present below a model for the bias of the proto-haloes. Note, however, that any lack of accuracy of the adopted Lagrangian biasing scheme will propagate through the time evolution of our model and contribute to the imprecision of its final results. Therefore, in order to test the accuracy of our evolutionary equations alone, we will also extract initial conditions directly from the simulation and compare the corresponding outcome of the evolution model with the statistics of the simulated haloes at z = 0. Let us consider the overdensity of proto-haloes in Lagrangian space δh (q) and the corresponding mass-density fluctuation δ(q). We assume that their Fourier transforms are linked by the expression: δh (k) = (b1 + b2 · k2 ) δm (k),

(15)

which corresponds to a non-local relation in real space. This form was first proposed by Matsubara (1999) and describes

Halo clustering with resummed perturbation theories the clustering of linear density peaks (Desjacques 2008). In this case, the initial conditions for P33 and P13 are: P33 (k)

=

(b1 + b2 k2 )2 P11 (k)e−k

P13 (k)

=

(b1 + b2 k2 )P11 (k)e−k

2

2

R2

,

R2 /2

,

(16)

where the exponential functions accounts for the finite size of the density peaks corresponding to a given halo mass. We find that the expression above accurately describes the cross-spectrum of proto-haloes and matter in our N-body simulation when b1 , b2 and R are treated as fitting parameters (see Elia et al. 2011 for further details). 1 This is not surprising as Ludlow & Porciani (2010) have shown that most proto-haloes include a density peak of the corresponding mass scale within their Lagrangian volume. In Table 1 we quote, for each halo-mass bin, the parameters b1 , b2 and R that best fit the simulation data using eq. (16) where the linear matter power spectrum P11 is computed using the CAMB online tool Lewis et al. (2000). Whereas Gaussianity is a good approximation for the linear matter distribution, fluctuations in the halo counts are non-Gaussian even in the initial conditions. We can quantify their level of non-Gaussianity in terms of their auto and cross bispectra that, using equation (15), can all be reduced to one of the following forms: B333 (k1 , k2 , k3 )

=

(b1 + b2 k12 )(b1 + b2 k22 )(b1 + b2 k32 ) Bm (k1 , k2 , k3 ) ,

σn2 =

1 2π 2

Z

∞

dk k2(n+1) P11 (k) e−k

2 R2

.

5 (20)

0

In the equations above a Gaussian smoothing window of size R has been adopted to coarse grain the matter density and identify the peaks. Once the proto-halo mass is linked to R through M = (2π)3/2 ρ¯ R3 (with ρ¯ the mean comoving density of matter), this model has no free parameters, since the spectral moments are completely defined by P11 (k). It follows that the initial conditions for P24 and P44 are P44 (k)

=

bv (k)2 P22 (k) ,

P24 (k)

=

bv (k) P22 (k) ,

(21)

and these expressions are in very good agreement with the spectra computed from the N-body simulation (Elia et al. 2011). The corresponding values of σ02 /σ12 are listed in Table 1 as a function of halo mass. Note that the velocity bias becomes more and more important on large scales with increasing the halo mass.

4

LINEAR THEORY

The lowest order approximation to the perturbation equation (5) consists in setting γabc = 0. In this limit, the evolution of the field from the initial time η = 0 to a generic η is given by

B133 (k1 , k2 , k3 )

=

(b1 + b2 k22 )(b1 + b2 k32 )Bm (k1 , k2 , k3 ) ,

ϕa (k; η) = gab (η)ϕb (k; 0) ,

B113 (k1 , k2 , k3 )

=

(b1 + b2 k32 )Bm (k1 , k2 , k3 ) ,

where gab (η) is the linear propagator, defined by the equation

(17)

where the matter bispectrum Bm (k1 , k2 , k3 ) is computed using the tree-level expression of the standard perturbation theory, 1 2 k2 1 k1 1 µ12 + µ12 + + Bm (k1 , k2 , k3 ) = 2 2 2 k2 k1 2 P11 (k1 , zin )P11 (k2 , zin ) + cycl. , (18) with µ12 ≡

k1 ·k2 . k1 k2

In order to solve our evolutionary equations, we need to know also the linear velocity field of proto-haloes. In principle proto-haloes might not move with the same velocity as matter at the same location (at the very least they should match the mass velocity smoothed on the Lagrangian halo size). We model this effect assuming that proto-haloes are indeed related to linear density peaks which, as discussed above, gives a good description of their clustering properties. In particular we follow Desjacques & Sheth (2010) who proposed a model for the peak velocities, which in Fourier space assumes the form 2 2 σ2 (19) θh (k) = 1 − 02 k2 e−k R /2 θm (k) ≡ bv (k)θm (k) , σ1 with bV the scale-dependent “velocity bias” and σn (n = 0, 1) being the spectral moments of the matter power spectrum defined as

1

Note that, adopting a local-bias model, δh (q) = b1 · δm (q) + 2 (q), provides a worse fit to simulation results. In this case, b2 · δm the model P13 lacks power for all but the smallest wavenumbers.

(δab ∂η + Ωab ) gbc (η) = δac δD (η) ,

(22)

(23)

with δab the Kronecker delta. Solving eq. (23) with causal boundary conditions (gab (η) = 0 for η < 0, see e.g. Crocce & Scoccimarro 2006) one gets 3/5 2/5 0 0 3/5 2/5 0 0 gab (η) = 3/5 2/5 0 0 3/5 2/5 0 0 2/5 −2/5 0 0 −3/5 3/5 0 0 +e−5/2η 2/5 −2/5 0 0 −3/5 3/5 0 0 0 0 0 0 0 0 0 0 +e−3/2η 0 2 0 −2 0 −1 0 1 0 0 0 0 0 0 0 0 (24) +e−η −1 −2 1 2 θ(η) , 0 0 0 0 with θ(η) Heavyside’s step function. Notice that gab (η) → δab as η → 0+ . The first and second contributions represent the standard growing and decaying modes, respectively (Crocce & Scoccimarro 2006). The third and fourth contributions represent two new modes, decaying respectively as e−3η/2 and e−η compared with the growing one. To understand their physical effect we notice that an initial condition of the form

6

A. Elia, S. Kulkarni, C. Porciani, M. Pietroni, S. Matarrese Bin Bin Bin Bin Bin

1 2 3 4

Mass range (1013 M⊙ /h)

# haloes

kmax (Mpc−1 h)

b1

b2 (Mpc2 h−2 )

R (Mpc h−1 )

σ02 /σ12 (Mpc2 h−2 )

1.24 − 1.8 1.8 − 3.4 3.4 − 10 > 10

202948 211305 150105 48985

0.24 0.24 0.22 0.15

7.28 ± 0.38 14.2 ± 0.4 25.9 ± 0.4 66.2 ± 1.3

422 ± 102 356 ± 84 708 ± 103 1025 ± 401

2.7 ± 0.8 2.1 ± 0.7 2.9 ± 0.4 3.5 ± 0.8

9.1 12.3 20.5 50.9

Table 1. Mass range and number of the haloes in the four bins.

ϕ(k) ϕ(k) ϕa (k; 0) = ϕh (k) , ϕv (k)

(25)

evolves (using eq. (22)) into ϕa (k; η) given by ϕ ϕ ϕ + 2e−3η/2 (ϕ − ϕv ) + e−η (−3ϕ + ϕh + 2ϕv ) ϕ − e−3η/2 (ϕ − ϕv )

, (26)

i.e. both the halo density and velocity fields relax to the matter ones as η → ∞ (but at a different pace). Also note that in the absence of an initial density bias (i.e. φ3 = φ) but in the presence of an initial velocity bias (i.e. φ4 6= φ), the linear dynamics quickly generates a transient density bias that vanishes at late times as e−η −e−3η/2 . The initial power spectrum at η = 0, corresponding to the field configuration 0 (k), and it evolves forward in time as in eq. (25), is Pab 0 PL,ab (k; ηa , ηb ) = gac (ηa )gbd (ηb )Pcd (k) .

4.1

(27)

The importance of velocity bias

It is interesting to assess the role of the velocity bias in the linear solution previously discussed. Assuming ϕ4 = ϕ2 at all times, the linear propagator for the first three components ϕi with i = 1, 2, 3 becomes 3/5 2/5 0 gab (η) = 3/5 2/5 0 3/5 2/5 0 2/5 −2/5 0 −5/2η −3/5 3/5 0 +e 2/5 −2/5 0 0 0 0 0 0 θ(η) , (28) +e−η 0 −1 0 1 and the third component of (26) reduces to ϕ3 (k; η) = ϕ(k) + e−η (ϕh (k) − ϕ(k)) .

(29)

This expresses the well known linear debiasing between the halo and matter fields at late times, derived by Fry (1996) for tracers that do not undergo merging and move solely under the influence of gravity. The corresponding halo-matter cross spectrum is (3)

0 0 0 PL,13 = P11 + e−η (P13 − P11 ) ≡ PL,13 ,

(30)

while keeping ϕ4 6= ϕ2 one gets (3)

0 0 PL,13 = PL,13 + 2(P11 − P14 )(e−3η/2 − e−η ) .

(31)

Figure 1. A comparison between P14 (blue dashed line) and P11 (red solid line) in the four bins. The vertical black dotted lines represent the limit up to which we expect our fluid approximation for haloes to work. A smoothing scale of R = 7 Mpc/h has been used for P14 . For a fair comparison with P11 , which is not smoothed, P14 has been redivided by the smoothing function. 0 0 In Figure 1 we compare P11 against P14 extracted from the N-body simulation, for the different mass bins. While the spectra agree well on very large scales (k . 0.05 h Mpc−1 ), they progressively depart for smaller scales. This is in line with the model introduced in Section 3. Note that the last term in eq. (31) vanishes in the initial conditions, reaches a minimum for η ≈ 0.8 and it is suppressed at late times. We quantify its amplitude at z = 0 in Figure 2, where we (3) plot the ratio rL = PL,13 /PL,13 which ranges between 0 − 3 per cent, depending on halo mass and scale. This suggests that the effect of the velocity bias on the halo-matter cross spectrum is small for low redshifts.

5

ANALYTICAL TREATMENT OF NON-LINEARITIES

In this section we deal with the non-linear evolution of the matter-halo system. We first compute 1-loop corrections for the propagator and then perform the corresponding large-k resummation to all perturbative orders thus extending the

Halo clustering with resummed perturbation theories +(a ↔ b) , II ∆Pab (k; ηa , ηb )

=

Z

ηa

dsa

0

7 (35)

Z

ηb

dsb Φcd (k; sa , sb )

0

gac (k; ηa − sa ) gbd (k; ηb − sb ) , (36) with Φcd (k; sa , sb )

=

2esa +sb

Z

d3 q

γcei (k, −q, q − k)PL,ef (q, sa , sb ) PL,ih (|k − q|, sa , sb )γdf h (−k, q, k − q) . (37) The explicit expressions for

5.2

(3)

Figure 2. The ratio between PL,13 and PL,13 : going from bottom to top, Bin 1 (red), Bin 2 (blue), Bin 3 (green) and Bin 4 (black).

results presented by Crocce & Scoccimarro (2006) for the matter density field. Finally, we compute non-linear power spectra using the TRG approach. For simplicity we only consider the case with no velocity bias, which we have demonstrated to be accurate (at least in the linear regime) at low redshifts, where current observations are available. The 3×3 Ω matrix for this case is 1 −1 0 3 3 Ω= − (32) 0 . 2 2 0 −1 1 5.1

1-loop perturbation theory

The 1-loop correction to the linear propagator (Crocce & Scoccimarro 2006) is given by Z η Z s1 Z ∆gab (k; η) = 4 ds1 ds2 d3 q es1 +s2 0

0

gac (η − s1 )γcie (k, −q, q − k)gef (s1 − s2 ) γf hd (k − q, q, −k)gdb (s2 )PL,ih (q, s1 , s2 ) , (33) (see appendix A for its explicit expressions). The 1-loop correction to the two-point correlator (10) is given by the sum of two contributions I ∆Pab

(k; ηa , ηb ) +

II ∆Pab

(k; ηa , ηb ) ,

(34)

which are also known in the literature as “P13 ” and “P22 ”, respectively 2 . They are given by I ∆Pab

2

(k; ηa , ηb )

=

0 ∆gac (k; ηa ) gbd (k; ηb ) Pcd (k)

They should not be confused with the P13 and P22 of our notation!

II ∆Pab

are given in appendix A.

Large-k resummation for the propagator

In the large-k limit, the 1-loop correction for the propagator, eq. (33), grows as k2 , and eventually dominates over the (k-independent) linear propagator. Taking into account higher orders, the situation gets even worse. The 2-loop correction grows as k4 , the 3-loop as k6 , and so on. This a manifestation of the perturbative expansion breakdown in cosmological PT, which appears not only in the computation of the propagator, but also of the power spectrum, the bispectrum, and so on. However, for the case of the propagator, it was shown in Crocce & Scoccimarro (2006b) that the leading order corrections in the large k and large η limit can be resummed at all orders in perturbation theory, giving a well-behaved propagator. The propagator Gac (k; η) connects the initial correlators with the cross-correlations between final and initial field configurations, 0 NG Pab (k; η, 0) = Gac (k; η)Pcb (k) + ∆Pab (k; η, 0) ,

(38)

where the last term at the rhs comes from the initial nonGaussianity of the matter and halo fields. At leading order, it is given by (see appendix C) Z η NG ∆Pab (k; η, 0) = ds es gac (η − s)gdf (s)geg (s) 0 Z × d3 q γcde (k, −q, q − k)Bf gb (q, k − q, −k) , (39) where Babc is the initial bispectrum at η = 0 (z = zin ) (see also Bernardeau et al. 2010). In Section 6.1 we will use eq. (38) to assess the validity of different approximation schemes for the propagator. In the large-k limit, G decays as 2 2 2η k σ e Gab (k; η) = gab (η) exp − , (40) 2 with 1 σ = 3 2

Z

d3 q

P 0 (q) . q2

(41)

Therefore, at least in the case of the propagator, the bad ultraviolet behavior is just an artifact of the perturbative expansion, which, at any finite order, completely misses the nice –and physically motivated– Gaussian decay of eq. (40) (see Crocce & Scoccimarro 2006b and Matarrese & Pietroni 2007 for a detailed discussion).

8

A. Elia, S. Kulkarni, C. Porciani, M. Pietroni, S. Matarrese

Although the result (40) was obtained for the 2 × 2 propagator of the matter density-velocity system, it holds, taking into account proper modifications, also when halos are included, i.e. for the 3 × 3 propagator considered in this Section. As in Crocce & Scoccimarro (2006b), we obtain an improved propagator interpolating between the 1-loop result (eq. 33) at low k and the Gaussian decay (40) (with a modified pre-factor) at high k. The details of the derivation and the relevant formulae are given in appendix B.

5.3

TRG

Unlike the propagator, the power spectrum cannot be resummed analytically at large k. Different semi-analytical procedures (Crocce & Scoccimarro 2006, McDonald 2007, Taruya & Hiramatsu 2008) have been proposed to compute it in the mildly non-linear regime. In this paper we will consider the TRG technique introduced in Pietroni (2008). Starting from eq. (5), a hierarchy of differential equations for the power spectrum, the bispectrum and higher order correlations is obtained. We choose to truncate it at the level of the trispectrum Qabcd = 0, so that the equations for Pab and Babc form a closed system ∂η Pab (k; η) = − Ωac (η)Pcb (k; η) − Ωbc (η)Pac (k; η) Z +eη d3 q [γacd (k, −q, q − k) Bbcd (k, −q, q − k; η) +Bacd (k, −q, q − k; η) γbcd (k, −q, q − k)] ,

The vertical black dotted lines represent the limit up to which we expect our fluid approximation for the haloes to work.

∂η Babc (k, −q, q − k; η) = −Ωad (η)Bdbc (k, −q, q − k; η) −Ωbd (η)Badc (k, −q, q − k; η) −Ωcd (η)Babd (k, −q, q − k; η) +2eη [γade (k, −q, q − k)Pdb (q; η)Pec (|k − q|; η) +γbde (−q, q − k, k)Pdc (|k − q|; η)Pea (k; η) + γcde (q − k, k, −q)Pda (k; η)Peb (q; η)] ,

(42)

which integrated gives the power spectra at any redshift and for any momentum scale. The system (42) consists of coupled differential equations which are solved numerically, starting from given initial conditions, i.e. Pab (k; ηi ) and Babc (k, −q, q − k; ηi ). From eq. (32), we can observe that Ω13 and Ω23 are zero, which means that the evolution of δm and θ is not modified, with respect to the original TRG formulation, by the presence of δh , as it is expected.

6 6.1

Figure 3. Cross spectrum between the halo density at z = 0 and the matter density at z = 50. The outcome of the N-body simulation (black points with error bars) is compared against lineartheory PP L (blue dashed line) and resummed result PP R (green solid line). The red dotted line shows the effect of neglecting the NG . non-Gaussian term, i.e. PP R − ∆P31

RESULTS

do not have to deal with the shot-noise problem. We extract the cross spectra from the simulation and compare them against those obtained both using linear theory prop0 0 0 agators PP L ≡ g31 P11 + g32 P21 + g31 P31 and resummed 0 0 0 NG propagators PP R ≡ G31 P11 + G32 P21 + G31 P31 + ∆P31 ; the result is shown in Figure 3. We note that the linear model severely overpredicts the two-time cross spectrum for k > 0.05 h Mpc−1 . It is evident that the resummed theory improves considerably upon the linear one, and agrees with the simulation within 10 per cent accuracy up to the scale where the fluid approximation holds. We include in NG ∆P31 the effect of the initial non-Gaussianity of the halo field via its initial bispectra, computed as in eq. (17). It turns out that the components giving a non-vanishing effect are of the B113 type. Their contribution is suppressed by a D+ (zin )/D+ (z = 0) factor with respect to that of an hypothetical primordial non-Gaussianity in the matter field (which we do not consider here). Therefore, the effect is of modest entity but, nevertheless, it improves the agreement with the simulation with respect to the case in which it is neglected.

Propagator

In order to assess the validity of our analytical approach, we compare our results for the resummed propagator against the simulation; to this end, we consider the relation in eq. (38). In particular we choose the indices a = 3, b = 1, so that we can check the components related to the haloes that were not present in the original formalism by Crocce & Scoccimarro (2006b) and, at the same time, we

6.2

Power spectrum

The TRG equations presented above are integrated numerically starting from the initial conditions discussed in section 3. As a first step, we set all the initial bispectra to zero In Figure 4 we show a comparison between the halo-matter cross spectra extracted from the N-body simulation and

Halo clustering with resummed perturbation theories

9

investigate the case of P11 . In the one-loop computation, the initial P11 contributes to the final one with a “weight” of (D+ (z = 0)/D+ (zin ))2 = e2η , as it also happens in linear theory. The initial matter bispectrum B111 , instead, has a weight of eη , so it is suppressed by a factor of e−η 3 . If we now consider the haloes, we have showed in (28) that there is a new decaying mode, responsible for the linear debiasing. This new mode, that involves only the halo field, carries an extra e−η suppression factor. We can now rank the contributions to P13 according to their relevance: (i) P11 (ii) P13 , B111 (iii) P33 , B113 (iv) B133 (v) B333 .

Figure 4. The cross spectrum between matter and halo distribution at z = 0 is shown in the four bins; the black dots with error bars represent the simulation, the blue dashed line is linear theory, the violet dot-dashed line is 1-loop and the red solid line is TRG. The vertical black dotted lines represent the limit up to which we expect our fluid approximation for the haloes to work.

the results of TRG, one-loop and linear theory. In the first three bins, corresponding to lower halo masses, linear theory overpredicts the power on mildly non-linear scales; note that this departure arises on larger scales compared to the matter auto spectrum (not shown in the figure). The overprediction of linear theory is cured by the one-loop power spectrum only on very large scale, while the TRG manages to correct it up to a smaller scale, before starting to fail. The fourth bin, though, displays a different behaviour: linear theory lacks of power on small scales, and neither the 1-loop correction nor the TRG are much more accurate. This might originate from the fact that very massive haloes are large and rare in the initial conditions, therefore less suited for the fluid approximation. Moreover, they also display the strongest velocity bias, which we are neglecting in our current non-linear treatment. A more quantitative analysis is presented in Figure 5, where we plot the ratios between the spectra from the simulation and the theoretical results. The TRG gives a cross spectrum within 5% accuracy at least up to k = 0.1 h Mpc−1 (barring bin 4), while linear theory does so up to k = 0.05 h Mpc−1 . As previously pointed out, however, the halo-density field is not Gaussian at zin = 50. We use expressions (17) to account for initial bispectra in the TRG method (they are not present in the linear theory and at 1-loop level). Figure 5 illustrates that this indeed produces a slight improvement in the agreement with the simulation. The correction resulting from the introduction of the initial bispectra turns out to be quite small. The reason is that their contribution to the final cross spectrum P13 is suppressed. To understand why, we can use perturbation theory; first, let us

Each item is suppressed by a factor of e−η with respect to the previous one. We can see that only B111 has some relevance, while the other terms are highly suppressed. Even if the reasoning was based on perturbation theory, it is valid also for TRG, at a qualitative level. Incidentally, the fact that P11 is the most relevant contribution is another evidence of debiasing. We can now address the effect of the truncation we perform in the TRG, namely considering the trispectrum Q = 0. First of all, the matter trispectrum Q1111 can be neglected in the range of scales under consideration, as one can conclude from the comparison between TRG and simulations in Pietroni (2008). The contribution from initial mixed (i.e. matter-halo) or pure halo trispectra is furtherly suppressed with respect to Q1111 by extra e−η factors, for the same reason as above. However, the trispectrum has its own time evolution as well, and one might argue that it becomes more relevant for z < zin ; this seems to be not the case, because enlarging the TRG truncation scheme by including the running of the trispectrum gives a contribution to the power spectrum which is at least of two-loop order and it is certainly subdominant in the scales we are considering. From Section 5 onward, we neglected any velocity bias, since this approximation proved to be accurate enough at z = 0 (in linear theory). As a further check, it is interesting to observe from Figure 6 that the TRG is able to give a better prediction than full linear theory in eq. (31), even neglecting the velocity bias that the linear theory accounts for. We can now look at the model predictions for the halo bias, defined as the ratio P13 /P11 . This quantity is plotted as a function of wavenumber in Figure 7. While the linear-theory bias always increases with scale, irrespectively of halo mass, the TRG result closely follows the scale dependence of b(k) seen in the simulation for the first two bins. It also gives a nearly constant bias for the third bin, as the simulation does, even though with a slightly lower value. However, the linear model performs better in the last bin where the bias in the simulation increases with k. It is also interesting to look at the results displayed in a different way; we can investigate the evolution of the cross spectrum from the initial conditions to today and the

3 For the most experienced readers, this happens because one of the two vertices carrying the eη factor is replaced by the bispectrum.

10

A. Elia, S. Kulkarni, C. Porciani, M. Pietroni, S. Matarrese

Figure 5. The ratio between the spectrum from simulation and: linear theory (blue squares), TRG without bispectra (red triangles), TRG with bispectra (green circles). The vertical black dotted lines represent the limit up to which we expect our fluid approximation for haloes to work. The shaded area marks the 5 per cent accuracy interval.

Figure 7. The effective bias b ≡ Pmh /Pm as a function of the wavevector in the four bins: the black circles with error bars are from the simulation, the blue dashed line represents linear theory and the red solid line TRG. The vertical black dotted lines represent the limit up to which we expect our fluid approximation for the haloes to work.

evolution of the bias as well. In Figures 8 and 9 we plot, respectively, rp ≡

Pmh (z = 0) Pmh (z = 50)

and

rb ≡

b(z = 0) . b(z = 50)

(43)

Again, our model is able to match accurately the trend of the simulation and to improve upon linear theory, excluding bin 4. A key feature of the linear solution in eq. (26) is the debiasing between halo and matter distributions with time. It is worth noting from Figure 9 that this effect is stronger for high-mass haloes, but constant on all the scales, while for low-mass haloes it is weaker on large scales and presents a strong k-dependence.

7

Figure 6. The ratio between the cross spectrum from simulation and TRG with bispectra (green) and linear theory with the inclusion of velocity bias (magenta). The vertical black dotted lines represent the limit up to which we expect our fluid approximation for haloes to work.

CONCLUSIONS

We have presented a novel approach to modeling the clustering of dark-matter haloes on mildly non-linear scales. This follows the motion of the regions that will collapse to form haloes (that we dub proto-haloes). Since the number of proto-haloes is conserved over time, for sufficiently large scales (k < 0.2 h Mpc−1 ), we can write a set of fluid equations that govern their evolution under the effect of gravity, which couples perturbations in the halo and matter density fields. We provide analytical solutions for the linearized equations and 1-loop perturbative corrections for the halo and matter power spectra. For the propagator, quantifying the memory of the density and velocity fields to their initial conditions, we also perform a resummation of perturbative corrections. Finally, for the power spectrum we compute

Halo clustering with resummed perturbation theories

11

the non-linear evolution using a semi-analytical procedure, namely an extension of the time renormalization group. The initial conditions for our evolutionary equations are specified adopting a Lagrangian bias model, originally developed to describe the clustering and motion of linear density peaks. We fix the parameters of the model so that to reproduce the distribution of proto-haloes in a high-resolution N-body simulation at z = 50. We use the same simulation to test the predictions of our model at z = 0. Our main results can be summarized as follows:

Figure 8. The ratio between the halo-matter cross spectrum at z = 0 and z = 50 for the four mass bins: symbols represent the simulation, dashed lines the linear theory and solid lines TRG. From top to bottom, we have Bin 1 in red, Bin2 in blue, Bin 3 in green and Bin 4 in black.

• Independently of the initial conditions, in the linear solution the halo density and velocity fields asymptotically match the corresponding matter fields at late times. This ’debiasing’ develops at a different rate for the density and the velocity, being faster for the latter. • Even if there is no initial density bias, the presence of a velocity bias generates a transient density bias that vanishes at late times. • Neglecting any initial velocity bias alters the linear predictions for the halo-matter cross spectrum at redshift z = 0 only by less than 3 per cent, for k < 0.3 h Mpc−1 . This provides us with the motivation to ignore the velocity bias in the non-linear analysis. • Unlike its linear counterpart, the resummed propagator is in good agreement with the N-body simulation, independently of halo mass. • The halo-matter cross spectrum predicted by the TRG is accurate to 5 per cent up to k ≃ 0.1 h Mpc−1 for a broad range of halo masses. This does not hold for very massive haloes (M > 1014 h−1 M⊙ ), that have low number density and high initial velocity bias, for which discreteness effects are more important. • The TRG result improves upon both linear theory and 1-loop corrections. Its performance is slightly enhanced accounting for the initial non-Gaussianity of the halo distribution. • For low halo masses our model accurately describes the scale-dependent bias measured in the simulation at z = 0.

ACKNOWLEDGMENTS AE was supported through the SFB-Transregio 33 “The Dark Universe” by the Deutsche Forschungsgemeinschaft (DFG). AE and SK acknowledge BCGS support for part of the work. CP thanks the participants to the workshop “Modern Cosmology: Early Universe, CMB and LSS” held at the Centro de Ciencias de Benasque Pedro Pascual for stimulating discussions. Part of this work has been developed there in August 2010.

Figure 9. The ratio between the bias at z = 0 and z = 50 for the four mass bins: symbols represent the simulation, dashed lines the linear theory and solid lines TRG. From top to bottom, we have Bin 1 in red, Bin2 in blue, Bin 3 in green and Bin 4 in black.

REFERENCES Bardeen J.M., Bond J.R., Kaiser N., Szalay A.S., 1986, ApJ, 304, 15 Bartolo N., Beltran Almeida J.P., Matarrese S., Pietroni M., Riotto A., 2010, JCAP 1003, 011 Bernardeau F., Crocce M., Sefusatti E., 2010, Phys. Rev. D, 82, 083507 Buchert T., Dominguez A., 2005, A& A, 438, 443

12

A. Elia, S. Kulkarni, C. Porciani, M. Pietroni, S. Matarrese

Catelan P., Lucchin F., Matarrese S., Porciani C., 1998, MNRAS, 297, 692 Catelan P., Porciani C., Kamionkowski M., 2000, MNRAS, 318, L39 Cole S., Kaiser N., 1989, MNRAS, 237, 1127 Cole S. et al., 2005, MNRAS, 362, 505 Coles P., 1993, MNRAS, 262, 1065 Cooray A., Sheth R.K., 2002, PhR, 372, 1 Crocce M., Scoccimarro R., 2006, Phys. Rev. D, 73, 063519 Crocce M., Scoccimarro R., 2006, Phys. Rev. D, 73, 063520 Desjacques V., 2006, Phys. Rev. D, 78, 103503 Desjacques V., Sheth R.K., 2010, Phys. Rev. D, 81, 023526 Desjacques V., Crocce M., Scoccimarro R., Sheth R.K., 2010, Phys. Rev. D, 82, 103529 Eisenstein D.J., Hu W., 1998, ApJ, 496, 605 Elia A., Ludlow A., Porciani C., in preparation Fry J.N., 1996, ApJ, 461, L65 Fry J.N., Gazta˜ naga E., 1993, ApJ, 413, 447 Guzik J., Bernstein G., Smith R.E., 2007 MNRAS, 375, 1329 Huff E., Schulz A. E., White M., Schlegel D.J., Warren, M.S., 2007, APh, 26, 351 Kaiser N., 1984, ApJ, 284, L9 Lewis A., Challinor A., Lasenby A., 2000, ApJ, 538, 473 Ludlow A. D., Porciani C., 2010, arXiv 1011.2493 Manera M., Gazta˜ naga E., 2009, arXiv 0912.0446 Manera M., Sheth R.K., Scoccimarro R., 2010, MNRAS, 402, 589 Mann, R. G., Peacock, J. A., Heavens, A. F., 1998, MNRAS, 293, 209 Matarrese S., Pietroni M., 2007, JCAP, 6, 26 Matsubara T., 1999, ApJ, 525, 543 Matsubara T., 2008, Phys. Rev. D, 78, 3519 McDonald P., 2007, Phys. Rev. D, 75, 043514 Mo H.J., White S.D.M., 1996, MNRAS, 282, 347 Montesano F., Sanchez A.G., Phleps S., 2010, arXiv 1007.0755 Norberg P. et al., 2002, MNRAS, 332, 827 Pietroni M., 2008, JCAP, 10, 36 Pillepich A., Porciani C., Hahn O., 2010, MNRAS, 402, 191 Porciani C., Catelan P., Lacey C., 1999, ApJ, 513, L99 Porciani C., Dekel A., Hoffman Y., 2002, MNRAS, 332, 339 Schulz A. E., White, M., 2006, APH, 25, 172 Scoccimarro R., 2004, Phys. Rev. D, 70, 083007 Seljak U., 2001, MNRAS, 325, 1359 Seo H., Eisenstein D.J., 2005, ApJ, 633, 575 Smith R.E., Scoccimarro R., Sheth R.K., 2007, Phys. Rev. D, 75, 063512 Sugiyama N., 1995, ApJ, 100, 281 Taruya A., Hiramatsu T., 2008, ApJ, 674, 617 Zehavi I. et al., 2004, ApJ, 608, 16

Halo clustering with resummed perturbation theories

13

APPENDIX A: COMPLETE ANALYTIC EXPRESSIONS OF THE 1-LOOP PROPAGATOR AND POWER SPECTRUM Linear power spectra: PL,11 (k; η, η)

=

0 P11 (k) = PL,12 (k; η, η) = PL,22 (k; η, η) ,

PL,13 (k; η, η)

=

0 0 e−η ((eη − 1)P11 (k) + P13 (k)) ,

PL,33 (k; η, η)

=

0 0 0 e−2η ((eη − 1)2 P11 (k) + P33 (k) + 2(eη − 1)P13 (k)) .

A1

(A1)

Propagator

First we report the integrands of the one-loop corrections of the components we are interested in. The linear power spectrum can be split as PL,ci ,cj (qi ; si , sj ) = P 0 (qi )Uci ,cj + ∆Pci ,cj (qi ; si , sj ) ,

(A2)

where Uci ,cj = 1

for any ci , cj ,

(A3)

while ∆Pci ,cj (qi ; si , sj ) 6= 0

only if ci or cj = 3 .

(A4)

(1) (2) Looking at eq. (33) we can consistently split R the one-loop correction to the propagator as ∆gab = ∆ gab + ∆ gab and we denote them as δ (i) gab , so that ∆(i) gab = dq δ (i) gab , for i = 1, 2. h i 3 (k+q)2 0 π (eη − 1)2 P11 (q) 4 6k7 q − 79k5 q 3 + 50k3 q 5 − 21kq 7 − 3 k2 − q 2 2k2 + 7q 2 log (k−q) 2 δg11 = , 840k3 q 3

δ (1) g31 δ

(1)

=

g32

=

δ (1) g33

=

δ (2) g31

=

δ (2) g32

=

δ (2) g33

=

A2

δg11 , 2 δg11 , 3 2 0 − e−η (eη − 1)2 k2 πP11 (q) , 3 3 k+q 0 0 πeη P13 (q) − P11 (q) 4 9k5 q + 24k3 q 3 − 9kq 5 − 18 k2 − q 2 log k−q 280k3 q

,

2 (2) δ g31 , 3 0

(A5)

1-loop PS

From eq. (37), one can see that the one-loop corrections require integrals over ds1 , ds2 and d3 q = q 2 sin θ dq dθ dφ; the integration over q and over µ ≡ cos θ cannot be performed analytically. The following expressions are p therefore integral kernels, denoted by δqµ P (the factor of q 2 is already taken into account). If cos θ = µ = k·q and |k − q| = k2 − 2kµq + q 2 , kq 2 0 0 πk4 e2η P11 (q) 7kµ + 3 − 10µ2 q P11 (|k − q|) II δqµ P11 = , 49 (k2 − 2kµq + q 2 )2 (A6) II δqµ P13

=

πk3 eη 7kµ + 3 − 10µ2 q 49 (k2 − 2kµq + q 2 )2

0 0 0 0 7q(k − µq)P13 (q)P11 (|k − q|) + P11 (q)7µ k2 − 2kµq + q 2 P13 (|k − q|)

0 0 +kP11 (q)P11 (|k − q|) 7kµ (eη − 1) − q

II δqµ P33

=

10µ2 − 3 eη − 14µ2 + 7 .

0 πk2 2 2 0 P13 (q)P13 (|k − q|) 2 98µq(k − µq) k − 2kµq + q 2 − 2kµq + q ) 0 0 0 +7q(k − µq)P11 (|k − q|) 2kP13 (q) 7kµ (eη − 1) − q 10µ2 − 3 eη − 14µ2 + 7 + 7qP33 (q)(k − µq) 2 0 0 +k2 P11 (q)P11 (|k − q|) q 10µ2 − 3 eη − 14µ2 + 7 − 7kµ (eη − 1) 2 0 0 +49µ2 P11 (q) k2 − 2kµq + q 2 P33 (|k − q|) 0 0 2 2 +14kµP11 (q) k − 2kµq + q 7kµ (eη − 1) − q 10µ2 − 3 eη − 14µ2 + 7 P13 (|k − q|) .

(A7)

49 (k2

(A8)

14

A. Elia, S. Kulkarni, C. Porciani, M. Pietroni, S. Matarrese

APPENDIX B: THE RESUMMED PROPAGATOR As it was discussed in Crocce & Scoccimarro (2006b), the leading contribution in k2 exp(2η) to the propagator at a generic n-loop order contains a chain of propagators and vertices of the form ga b1 (η − s1 )γb1 c1 a1 (k, −q1 , q1 − k)ga1 b2 (s1 − s2 ) · · · ga2n−1 b2n (s2n−1 − s2n )γb2n c2n a2n (k + q2n , −q2n , −k)ga2n b (s2n ) . (B1) The ci indices have to be contracted in all possible pairings, by inserting n linear power spectra, each of the form δD (qi + qj )PL,ci ,cj (qi ; si , sj ) .

(B2)

The n-loop contribution is obtained by multiplying by exp( integrating over Z η Z dsi d3 qi . Π2n i=1

P2n

i=1

si ) and by the appropriate combinatoric factor, and then by

(B3)

0

Recall eqs. (A2-A4): if all the insertions are of the P 0 (qi )Uci ,cj type then the resummation goes exactly as in the standard case. Indeed 1q·k Uc1 ,cj gab1 (η − s1 )γb1 c1 a1 (k, −q1 , q1 − k)ga1 b2 (s1 − s2 ) → ucj ga b2 (η − s2 ) , (B4) 2 q2 in the k ≫ q limit. Besides the explicit form for the vertices (see Section 1), we have used the composition property of the propagators gab (η − s1 )gbc (s1 − s2 ) = gac (η − s2 ) ,

(B5)

and have defined u = (1, 1, 1). Since each vertex, contracted by a ucj vector becomes proportional to a delta function in its first and third index, the chain of propagators composes up to a single propagator gab (η) and the time integral can be easily performed. The momentum integrals factorize into n integrals of the type Z (k · q)2 − d3 qP 0 (q) = −k2 σ 2 , (B6) q4 where σ 2 has been defined in (41). Using the appropriate combinatoric factors, the leading contribution to the propagator at n-loop, in the large momentum limit, when all the power spectrum insertions are of the P 0 (qi )Uci ,cj type is n 1 (eη − 1)2 gab (η) , (B7) −k2 σ 2 n! 2 η

2

]. which resums to gab (η) exp[−k2 σ 2 (e −1) 2 As for the insertions including the ∆Pci ,cj (qi ; si , sj ) contribution to the linear power spectrum, the important point to realize is that, due to (A4) and to the structure of the linear propagator and the vertices, a chain like (B1) can be contracted by at most one ∆Pci ,cj (qi ; si , sj ) power spectrum, the remaining ones being of the type P 0 (qi )Uci ,cj . Moreover, these insertions only contribute to G31 and G32 . Therefore, at n-th order, we have only two types of contributions: those with all P 0 (qi )Uci ,cj insertions, giving (B7), and those with one ∆Pci ,cj insertions and n − 1 P 0 (qi )Uci ,cj ones, which can also be resummed in the large k limit. As a consequence, the complete resummed propagator in the large momentum limit is Gab (k; η) = (gab (η) + δa3 fb ∆(2) g31 (k; η)) exp[−k2 σ 2

(eη − 1)2 ], 2

(B8)

with fb = (1, 2/3, 0). In order to have an expression for the propagator valid at any k, one can proceed as in (Crocce & Scoccimarro (2006b)) and interpolate between the large k limit above and the 1-loop result at low k. This can be done, for instance for G11 , starting from its 1-loop expression 5 g11 (η) + ∆g11 (k; η) ≃ g11 (η) 1 + ∆g11 (k; η) , (B9) 3 where the above approximation is exact in the large η limit. Since (eη − 1)2 5 ∆g11 (k; η) → −k2 σ 2 , for large k , 3 2 the required interpolation is given by 5 G11 (k; η) = g11 (η) exp ∆g11 (k; η) . 3 Proceeding analogously for the other components, we have: 5 G12 (k; η) = g12 (η) exp ∆g11 (k; η) , 3

(B10)

(B11)

Halo clustering with resummed perturbation theories G22 (k; η)

=

G21 (k; η)

=

G31 (k; η)

=

G32 (k; η)

=

G33 (k; η)

=

5 ∆g22 (k; η) , 2 5 g21 (η) exp ∆g22 (k; η) , 2 5 g31 (η) + ∆(2) g31 (k; η) exp ∆g11 (k; η) , 3 2 5 g32 (η) + ∆(2) g31 (k; η) exp ∆g11 (k; η) , 3 3 h i η (1) g33 (η) exp e ∆ g33 (k; η) , g22 (η) exp

15

where we have used the identities 2 ∆g11 (k; η) , ∆g12 (k; η) = 3 3 ∆g21 (k; η) = ∆g22 (k; η) , 2 (1) ∆ g31 (k; η) = ∆g11 (k; η) 2 (1) ∆(1) g32 (k; η) = ∆ g31 (k; η) 3 2 (2) ∆(2) g32 (k; η) = ∆ g31 (k; η) . 3

(B12)

(B13)

APPENDIX C: TAKING INTO ACCOUNT THE INITIAL HALO NON-GAUSSIANITY. While at zin the matter field can be considered gaussian with high accuracy, the same not necessarily holds for haloes. An initial non-Gaussianity for the equal-time power spectrum can be easily incorporated in the TRG formalism as discussed in Bartolo et al. (2010), and has been done in sect. 6.2. In order to obtain the effect of an initial non-Gaussianity on the cross-correlator at different times considered in eq. (38), it suffices to compute the O(γ) correction to the linear evolution of the field, by inserting eq. (22) in (5), to get Z η Z ϕ(1) (C1) ds gac (η − s)es d3 q γcde (k, −q, q − k)gdf (s)ϕf (q; 0)geg (s)ϕg (k − q; 0) . a (k; η) = 0

The cross-correlator ′ hϕ(1) a (k; η)ϕb (k ; 0)i ,

(C2)

then includes the non-Gaussian expectation value hϕf (q; 0)ϕg (k − q; 0)ϕb (k′ ; 0)i = δD (k + k′ )Bf gb (q, k − q, −k) , and therefore gives eq. (39). This paper has been typeset from a TEX/ LATEX file prepared by the author.

(C3)

Printed 22 March 2011

(MN LATEX style file v2.2)

arXiv:1012.4833v2 [astro-ph.CO] 21 Mar 2011

Modeling the clustering of dark-matter haloes in resummed perturbation theories A. Elia1⋆ , S. Kulkarni2, C. Porciani1, M. Pietroni3, S. Matarrese4,3 1

Argelander Institut f¨ ur Astronomie der Universit¨ at Bonn, Auf dem H¨ ugel 71, D-53121 Bonn, Germany Center for Theoretical Physics & Physikalisches Institut, Universit¨ at Bonn, Nußallee 12, D-53115 Bonn, Germany 3 INFN, Sezione di Padova, Via Marzolo 8, I-35131 Padova, Italy 4 Dipartimento di Fisica “G. Galilei”, Universit` a degli Studi di Padova, Via Marzolo 8, I-35131 Padova, Italy

2 Bethe

22 March 2011

ABSTRACT

We address the issue of the cosmological bias between matter and galaxy distributions, looking at dark-matter haloes as a first step to characterize galaxy clustering. Starting from the linear density field at high redshift, we follow the centre of mass trajectory of the material that will form each halo at late times (proto-halo). We adopt a fluid-like description for the evolution of perturbations in the proto-halo distribution, which is coupled to the matter density field via gravity. We present analytical solutions for the density and velocity fields, in the context of renormalized perturbation theory. We start from the linear solution, then compute one-loop corrections for the propagator and the power spectrum. Finally we analytically resum the propagator and we use a suitable extension of the time-renormalization-group method (Pietroni 2008) to resum the power spectrum. For halo masses M < 1014 h−1 M⊙ our results at z = 0 are in good agreement with N-body simulations. Our model is able to predict the halo-matter cross spectrum with an accuracy of 5 per cent up to k ≈ 0.1 h Mpc−1 approaching the requirements of future galaxy redshift surveys. Key words: cosmology: theory, large-scale structure of Universe– galaxies: haloes – methods: analytical, N-body simulations.

1

INTRODUCTION

Redshift surveys have shown that the clustering properties of galaxies strongly depend on their luminosity, color and morphological (or spectral) type (e.g. Norberg et al. 2002; Zehavi et al. 2004). This indicates that galaxies do not perfectly trace the distribution of the underlying dark matter, a phenomenon commonly referred to as ‘galaxy biasing’. Its origin lies in the details of the galaxy formation process which is shaped by the interplay between complex hydrodynamical and radiative processes together with the darkmatter driven formation of the large-scale structure. Attempts to infer cosmological parameters from galaxy clustering studies are severely hampered by galaxy biasing. A number of theoretical arguments and the outcome of numerical simulations both suggest that, on sufficiently large scales, the power spectra of galaxies and matter should be proportional to each other: Pg = b21 Pm where the linear bias factor b1 depends on galaxy type but is generally scale independent (e.g. Coles 1993; Mann et al. 1998). Similarly, to model higher-order statistics, such as the galaxy bispectrum,

⋆

E-mail: [email protected]

it is generally assumed that galaxy biasing is a local process 2 such that δg = b1 δm + b2 δm /2 + . . . where δg and δm are the (smoothed) galaxy and dark-matter density contrast, respectively (Fry & Gazta˜ naga 1993). However, the reliance of these phenomenological approximations limits cosmological studies to very large scales whereas data with better signal-to-noise ratio are already available on much smaller scales. Moreover, future studies of baryonic acoustic oscillations (e.g. Sugiyama 1995; Eisenstein & Hu 1998) will require measurements of the matter power-spectrum with percent or even sub-percent accuracy in order to shed new light on the source of cosmic acceleration. Understanding and controlling the effects of galaxy biasing with this precision will be challenging. All this provides a very strong motivation for developing more accurate (and physically driven) models of galaxy biasing. A number of authors have used the power spectrum statistics to explore the scale dependence of galaxy biasing based on numerical simulations (Cole et al. 2005; Seo & Eisenstein 2005; Huff et al. 2007; Manera et al. 2010; Manera & Gazta˜ naga 2009; Montesano et al. 2010) or analytical calculations (Seljak 2001; Schulz & White 2006; Guzik et al. 2007; Smith et al. 2007) stemming from either

2

A. Elia, S. Kulkarni, C. Porciani, M. Pietroni, S. Matarrese

perturbation theory or the halo model for the large-scale structure (see Cooray & Sheth 2002 for a review). The general picture is that galaxy biasing is expected to be scale dependent (i.e. Pg (k) = b(k)2 Pm (k)) and the functional form of b(k) can sensibly depend on the selected tracer of the large-scale structure. Since galaxies are expected to form within dark-matter haloes, understanding the clustering properties of the haloes is a key step to accurately model galaxy biasing. This is a much simpler problem, considering that dark-matter haloes form under the sole action of gravity. It is in fact expected that long-wavelength density fluctuations modulate halo formation by modifying the collapse time of localized short-wavelength density peaks (Bardeen et al. 1986; Cole & Kaiser 1989). This argument (known as the peakbackground split) predicts that, on large scales, the halo overdensity δh = b δm where the bias coefficient b varies with the halo mass (Mo & White 1996). The numerical value of the bias coefficient is determined by two different occurrences: first, haloes form out of highly biased regions in the linear density field (Kaiser 1984; Porciani et al. 1999) and, second, they move over time as they are accelerated towards the densest regions of the large-scale structure by gravity (Mo & White 1996). These two phenomena generally go under the name of “Lagrangian biasing” and “Lagrangian to Eulerian passage”, respectively. Mo & White (1996) dealt with the second problem by assuming that long-wavelength density perturbations evolve according to the spherical tophat model. A more sophisticated generalization of the peakbackground split has been presented by Catelan et al. (1998) who assumed that also the large-scale motion of the density “peaks” is fully determined by the long-wavelength component of the density field. Since the halo population and the matter feel the same large-scale gravitational potential, their density fluctuations are strongly coupled and their time evolution must be solved simultaneously. This makes the process of halo biasing non-linear and non-local even if one starts from a linear and local Lagrangian biasing scheme (Catelan et al. 1998; Matsubara 2008). The bispectrum can be used to test this model against the standard Eulerian local biasing scheme (Catelan et al. 2000). In this paper, we present a novel and very promising approach to model the clustering of dark matter haloes. Adopting the formalism by Catelan et al. (1998) combined with a non-local Lagrangian biasing scheme for the haloes (Matsubara 1999), we simultaneously follow the growth of perturbations in the matter and in the halo distribution over cosmic time. We present perturbative solutions for the corresponding overdensity and velocity fields and we are able to resum the perturbative series in the limit of large wavenumbers. Moreover, we write down a system of equations for the power spectra Pm and Ph using the time-renormalizationgroup (TRG) approach by Pietroni (2008) and numerically integrate them. Our results are in excellent agreement with the output of a high-resolution N-body simulation, showing an improvement over linear theory, and we are able to predict the matter-halo cross spectrum with a precision within 5 per cent for k < 0.15 h Mpc−1 . Related work has been very recently presented by Desjacques et al. (2010) who computed the two-point correlation function of linear density peaks and followed its time evolution assuming that peaks move according to the

Zel’dovich approximation. For massive haloes this results in a scale-dependent bias (with variations of ∼ 5 per cent) on the scales relevant for baryonic-oscillation studies. Contrary to their approach, we do not deal with a point process but describe large-scale fluctuations in the distribution of darkmatter haloes as perturbations in a continuous fluid. On the other hand, we account for the full gravitational motion of the haloes and do not rely on simplified dynamical models as like as the Zel’dovich approximation. The structure of our paper is as follows. In Section 2 we present our model for the joint evolution of the matter and halo power spectra. The initial conditions for our evolutionary equations are discussed in Section 3. The solution of the linearized equations is presented in Section 4 where we also quantify the importance of the halo velocity bias. Using a perturbative technique, in Section 5 we compute analytic solutions for the propagator of perturbations (the two-time correlator). We derive 1-loop corrections and, in the limit of large wavenumbers, the fully resummed propagator. The discussion in Sections 5.1 and 5.2 is very technical and the less experienced readers can safely skip it without compromising understanding of the remainder. In Section 6, we numerically integrate the full equations for the evolution of halo and matter power spectra in the TRG formalism. We then compare the results against the outcome of a high-resolution N-body simulation. Finally, in Section 7 we conclude.

2

THE MODEL

2.1

Dynamics of gravitational instability

The large-scale structure observed today in the universe is believed to be the result of gravitational amplification of primordial fluctuations, caused by the interaction among cold dark matter (CDM) particles. If we denote δm as the matter density contrast and v as the velocity, the Eulerian dynamics of a system of such particles, which interact only via gravity, is ruled by a set of three equations (continuity, Euler and Poisson) that in a ΛCDM model reads: ∂ δm + ∇ · [(1 + δm )v] = 0 , ∂τ ∂v + H v + (v · ∇)v = −∇φ , ∂τ 3 ∇2 φ = H2 Ωm δm , (1) 2 where τ is the conformal time. If we define the velocity divergence θ(x, τ ) = ∇ · v(x, τ ) and switch to Fourier space, the equations in (1) become: ∂ δm (k, τ ) + θ(k, τ ) Z∂ τ

+

d3 q d3 p δD (k − q − p)α(q, p)θ(q, τ )δm (p, τ ) = 0 ,

∂ θ(k, τ ) 3 + H θ(k, τ ) + H2 Ωm (τ )δm (k, τ ) ∂ τ 2 Z

+

d3 q d3 p δD (k − q − p)β(q, p)θ(q, τ )θ(p, τ ) = 0 , (2)

where

Halo clustering with resummed perturbation theories α(q, p) =

(p + q) · q , q2

β(q, p) =

(p + q)2 p · q , 2 p2 q 2

(3)

Equations (2) can be written in a compact form if we define a new time variable η ≡ ln(D+ /D+in ), being D+in the growth factor at an early epoch, and a doublet ϕa (a = 1, 2) δm (k, η) ϕ1 (k, η) , (4) ≡ e−η −θ(k, η)/(Hf+ ) ϕ2 (k, η) with f+ ≡ d ln D+ /d ln a. The velocity divergence is scaled such that it has the same dimension of the density contrast and in the linear regime −θ(k, η)/(Hf+ ) ≈ δm (k, η), i.e. ϕ1 (k, 0) = ϕ2 (k, 0). The system is therefore ∂η ϕa (k, η) = −Ωab (η)ϕb (k, η) +eη γabc (k, −p, −q)ϕb (p, η) ϕc (q, η),

(5)

where sum over repeated indices and integration over repeated momenta are understood. The vertex function γabc (k, p, q) (a, b, c, = 1, 2) has only three non-vanishing elements 1 γ121 (k, p, q) = δD (k + p + q) α(p, q) , 2 γ222 (k, p, q) = δD (k + p + q) β(p, q) , (6) and γ112 (k, q, p) = γ121 (k, p, q), with δD the Dirac-delta distribution. All the information about the cosmological model is contained in the matrix 1 −1 (7) Ω(η) = 3Ωm 3Ωm , − 2 2 2f+ 2f+ that, in the following, will be considered as a constant matrix, approximating Ωf m 2 ≈ 1. This ratio is indeed very close +

to unity for most of the history of the Universe. Making this approximation, one is modifying the behavior of the decaying mode, while the growing one is left unaltered. It has been shown in Pietroni (2008) that it affects the matter power spectrum at z = 0 at a less than percent level up to k ≈ 0.3 h Mpc−1 . The power spectrum, defined by an ensemble average, in this notation is a 2 × 2 matrix hϕa (k; η)ϕb (q; η)i ≡ δD (k + q)Pab (k; η) ,

(8)

and the bispectrum, defined by hϕa (k; η)ϕb (q; η)ϕc (p; η)i ≡ δD (k + q + p)Babc (k, q, p; η) ,

(9)

has 8 components. In the following, we will also consider a different-time two-point correlator, defined as hϕa (k, ηa )ϕb (q, ηb )i ≡ δD (k + q)Pab (k; ηa , ηb ) ,

(10)

which obviously coincides with (8) for ηa = ηb . 2.2

The distribution of dark-matter haloes

Let us consider a set of dark-matter haloes identified at a given redshift zid according to some predefined criterion. The material that forms the haloes can be traced back to its initial location in the linear overdensity field at z → ∞. We dub each of these regions as a proto-halo. In other words, a proto-halo is the Lagrangian region of space that

3

will collapse to form a halo at redshift zid . N-body simulations show that nearly all proto-haloes are simply connected (Porciani et al. 2002) even though this property is not key to our analysis. Let us now follow the evolution of a proto-halo over cosmic time in Eulerian space. Basically its shape and topology will be distorted (proto-haloes will first fragment into smaller substructures that will later merge to form the final halo) and its overall volume will be compressed while its centre of mass will move along a given trajectory determined by the mass density field via gravity. We focus our analysis onto this motion that connects the Lagrangian position of the proto-halo with the Eulerian location of the final halo. On scales much larger than the characteristic size of (and separation between) the proto-haloes, the density fluctuations traced by the centre-of-mass trajectories can be described with a continuous overdensity field δh (x, τ |zid ). Note that while τ labels conformal time along the trajectories, zid is just a tag that identifies the halo population. Unlike real haloes that undergo merging, by construction protohaloes always preserve their identity. Their total number is therefore conserved over time and we can write a continuity equation for their number density: ∂δh + ∇ · [(1 + δh )vh ] = 0 . (11) ∂τ Here the proto-halo density and velocity fields should be intended as coarse grained on a scale of a few times the mean inter-halo separation (so as to suppress discreteness effects as proto-haloes are individually separate units). Strictly speaking, the smoothed velocity field does not obey the Euler-Poisson system in equation (1) due to the presence of the non-linear term (v · ∇)v. In fact, the coarse graining procedure introduces new terms in the fluid equations generated by the degrees of freedom one has integrated out, namely: a velocity dispersion term and a correction to the mean- field gravitational acceleration due to density fluctuations on scales smaller than the smoothing radius (Buchert & Dominguez 2005). On the other hand, it is reasonable to assume that the large-scale motion of the protohaloes is generated by density fluctuations with wavelength larger than the characteristic halo size and is not influenced by perturbations with shorter wavelength. The very same assumption of neglecting the coupling to the small scales is routinely done when one writes equation (1) for the mass density field (see Section 3 in Buchert & Dominguez 2005) albeit adopting much narrower smoothing kernels than for the haloes. With the same spirit, in what follows, we will ignore the extra terms in the fluid equations generated by the coarse graining procedure. This is a working hypothesis which makes the problem mathematically treatable and whose accuracy will be tested by comparing our final results against high-resolution numerical simulations. We therefore write an Euler equation for the proto-halo fluid velocities ∂ vh + H vh + (vh · ∇)vh = −∇φ , (12) ∂τ where the gravitational potential is the same as in equation (1). Note that if vh matches v in the initial conditions then it will always do. On the contrary, any initial velocity bias between proto-halos and matter will be progressively erased by the gravitational acceleration.

4

A. Elia, S. Kulkarni, C. Porciani, M. Pietroni, S. Matarrese

Thus, given suitable initial conditions for δh and vh at τ → 0 (i.e. a prescription for the Lagrangian biasing of proto-haloes), we can in principle use equations (11) and (12) to follow the clustering evolution of the proto-halo population at all times. We are particularly interested in the solution of the fluid equations at the special time τ that corresponds to zid . In fact this solution has a particular physical meaning as it gives the density and velocity fields of the actual dark-matter haloes. 2.3

Growth of matter and halo perturbations

The system (1) is now extended by the inclusion of eq. (11) and eq. (12). We define a quadruplet ϕa (a = 1, 2, 3, 4) δm (k, η) ϕ1 (k, η) ϕ2 (k, η) −η −θ(k, η)/(Hf+ ) (13) , ϕ3 (k, η) ≡ e δh (k, η) −θh (k, η)/(Hf+ ) ϕ4 (k, η) in such a way that eq. (5) still holds, but with indices running from 1 to 4. There are three more non-vanishing elements of the vertex γ343 (k, p, q) = γ334 (k, q, p) = γ121 (k, p, q) and γ444 (k, q, p) = γ222 (k, p, q), and the 4 × 4 Ω matrix is 1 −1 0 0 3 3 0 0 − . 2 2 Ω= (14) 0 1 −1 0 3 3 0 0 − 2 2 From the definitions (8) and (9), with a = 1, 2, 3, 4, we get a 4 × 4 matrix for power spectrum; in the following, we will focus on the matter power spectrum P11 and the matter-halo cross spectrum P13 .

3

INITIAL CONDITIONS

In the previous section, we have presented a model that describes the non-linear evolution of the matter and halo density fields. Given suitable initial conditions, the formal equations we have introduced can be integrated numerically so that to compute the perturbative propagators and the TRG-evolved power spectra. The choice of the initial conditions therefore plays a very important role in our theory and will be the subject of this section. 3.1

N-body simulation

To gain insight into the properties of proto-haloes (and, later, to test our results at z = 0), we use one highresolution N-body simulation extracted from the suite presented by Pillepich et al. (2010). This contains 10243 darkmatter particles within a periodic cubic box with a side of Lbox = 1200h−1 Mpc and follows the formation of structure in a ΛCDM model with Gaussian initial conditions and cosmological parameters: h = 0.701, σ8 = 0.817, ns = 0.96, Ωm = 0.279, Ωb = 0.0462 and ΩΛ = 0.721. We identify dark-matter haloes at z = 0 using the friends-of-friends algorithm with a linking length equal to 0.2 times the mean interparticle distance. We only consider haloes containing more than 100 particles (i.e. with mass

M > 1.24 · 1013 h−1 M⊙ ) and we split them into four mass bins to keep track of their different clustering properties. The corresponding mass ranges and the total number of haloes in each bin are given in Table 1, along with an estimate of the highest wavevector up to which the fluid approximation for haloes holds. This value is determined by the number of haloes we require to be in a volume element to consider them as a fluid, and we set this number to 30. On smaller scales, our assumption breaks down, therefore we will look at results in the specified range, that, of course, decreases as the mass of the haloes increases. In the plots that will be shown in Section 6 the limit to which we can trust our model will be represented by vertical black dotted lines. Proto-haloes are identified by tracing the positions of the particles forming a halo at z = 0 back to the linear density field. The centre of mass of each proto-halo is used as a proxy for its spatial location. Similarly, the mass weighted linear velocity gives the proto-halo velocity. Halo and proto-halo density and momentum fields are computed with the cloud-in-cell grid assignement using a 5123 mesh. Velocity fields are obtained by taking the ratio of the momentum and density distributions (preventively smoothed to preclude the existence of empty cells) as shown in Scoccimarro (2004). Power spectra have been computed using FFT. In order to avoid uncertain shot-noise corrections for the haloes, we only consider their cross spectra with the matter density field. 3.2

Lagrangian halo bias

Concerning the matter density, the initial conditions are given by linear theory which directly follows from the adopted cosmological model (transfer function) and the statistics of primordial perturbations (spectral index, Gaussianity). On the other hand, for the dark-matter halos, we can follow two different approaches: (i) extract the relevant information directly from the simulation or (ii) use a model for the Lagrangian bias of the halos. The latter option offers a number of advantages. First, it allows us to make general predictions independently of the simulation specifics. Second, it allows us to include halo bispectra in our formalism (while it would be extremely demanding and time consuming to compute all possible triangular configurations from the simulation). For these reasons we will present below a model for the bias of the proto-haloes. Note, however, that any lack of accuracy of the adopted Lagrangian biasing scheme will propagate through the time evolution of our model and contribute to the imprecision of its final results. Therefore, in order to test the accuracy of our evolutionary equations alone, we will also extract initial conditions directly from the simulation and compare the corresponding outcome of the evolution model with the statistics of the simulated haloes at z = 0. Let us consider the overdensity of proto-haloes in Lagrangian space δh (q) and the corresponding mass-density fluctuation δ(q). We assume that their Fourier transforms are linked by the expression: δh (k) = (b1 + b2 · k2 ) δm (k),

(15)

which corresponds to a non-local relation in real space. This form was first proposed by Matsubara (1999) and describes

Halo clustering with resummed perturbation theories the clustering of linear density peaks (Desjacques 2008). In this case, the initial conditions for P33 and P13 are: P33 (k)

=

(b1 + b2 k2 )2 P11 (k)e−k

P13 (k)

=

(b1 + b2 k2 )P11 (k)e−k

2

2

R2

,

R2 /2

,

(16)

where the exponential functions accounts for the finite size of the density peaks corresponding to a given halo mass. We find that the expression above accurately describes the cross-spectrum of proto-haloes and matter in our N-body simulation when b1 , b2 and R are treated as fitting parameters (see Elia et al. 2011 for further details). 1 This is not surprising as Ludlow & Porciani (2010) have shown that most proto-haloes include a density peak of the corresponding mass scale within their Lagrangian volume. In Table 1 we quote, for each halo-mass bin, the parameters b1 , b2 and R that best fit the simulation data using eq. (16) where the linear matter power spectrum P11 is computed using the CAMB online tool Lewis et al. (2000). Whereas Gaussianity is a good approximation for the linear matter distribution, fluctuations in the halo counts are non-Gaussian even in the initial conditions. We can quantify their level of non-Gaussianity in terms of their auto and cross bispectra that, using equation (15), can all be reduced to one of the following forms: B333 (k1 , k2 , k3 )

=

(b1 + b2 k12 )(b1 + b2 k22 )(b1 + b2 k32 ) Bm (k1 , k2 , k3 ) ,

σn2 =

1 2π 2

Z

∞

dk k2(n+1) P11 (k) e−k

2 R2

.

5 (20)

0

In the equations above a Gaussian smoothing window of size R has been adopted to coarse grain the matter density and identify the peaks. Once the proto-halo mass is linked to R through M = (2π)3/2 ρ¯ R3 (with ρ¯ the mean comoving density of matter), this model has no free parameters, since the spectral moments are completely defined by P11 (k). It follows that the initial conditions for P24 and P44 are P44 (k)

=

bv (k)2 P22 (k) ,

P24 (k)

=

bv (k) P22 (k) ,

(21)

and these expressions are in very good agreement with the spectra computed from the N-body simulation (Elia et al. 2011). The corresponding values of σ02 /σ12 are listed in Table 1 as a function of halo mass. Note that the velocity bias becomes more and more important on large scales with increasing the halo mass.

4

LINEAR THEORY

The lowest order approximation to the perturbation equation (5) consists in setting γabc = 0. In this limit, the evolution of the field from the initial time η = 0 to a generic η is given by

B133 (k1 , k2 , k3 )

=

(b1 + b2 k22 )(b1 + b2 k32 )Bm (k1 , k2 , k3 ) ,

ϕa (k; η) = gab (η)ϕb (k; 0) ,

B113 (k1 , k2 , k3 )

=

(b1 + b2 k32 )Bm (k1 , k2 , k3 ) ,

where gab (η) is the linear propagator, defined by the equation

(17)

where the matter bispectrum Bm (k1 , k2 , k3 ) is computed using the tree-level expression of the standard perturbation theory, 1 2 k2 1 k1 1 µ12 + µ12 + + Bm (k1 , k2 , k3 ) = 2 2 2 k2 k1 2 P11 (k1 , zin )P11 (k2 , zin ) + cycl. , (18) with µ12 ≡

k1 ·k2 . k1 k2

In order to solve our evolutionary equations, we need to know also the linear velocity field of proto-haloes. In principle proto-haloes might not move with the same velocity as matter at the same location (at the very least they should match the mass velocity smoothed on the Lagrangian halo size). We model this effect assuming that proto-haloes are indeed related to linear density peaks which, as discussed above, gives a good description of their clustering properties. In particular we follow Desjacques & Sheth (2010) who proposed a model for the peak velocities, which in Fourier space assumes the form 2 2 σ2 (19) θh (k) = 1 − 02 k2 e−k R /2 θm (k) ≡ bv (k)θm (k) , σ1 with bV the scale-dependent “velocity bias” and σn (n = 0, 1) being the spectral moments of the matter power spectrum defined as

1

Note that, adopting a local-bias model, δh (q) = b1 · δm (q) + 2 (q), provides a worse fit to simulation results. In this case, b2 · δm the model P13 lacks power for all but the smallest wavenumbers.

(δab ∂η + Ωab ) gbc (η) = δac δD (η) ,

(22)

(23)

with δab the Kronecker delta. Solving eq. (23) with causal boundary conditions (gab (η) = 0 for η < 0, see e.g. Crocce & Scoccimarro 2006) one gets 3/5 2/5 0 0 3/5 2/5 0 0 gab (η) = 3/5 2/5 0 0 3/5 2/5 0 0 2/5 −2/5 0 0 −3/5 3/5 0 0 +e−5/2η 2/5 −2/5 0 0 −3/5 3/5 0 0 0 0 0 0 0 0 0 0 +e−3/2η 0 2 0 −2 0 −1 0 1 0 0 0 0 0 0 0 0 (24) +e−η −1 −2 1 2 θ(η) , 0 0 0 0 with θ(η) Heavyside’s step function. Notice that gab (η) → δab as η → 0+ . The first and second contributions represent the standard growing and decaying modes, respectively (Crocce & Scoccimarro 2006). The third and fourth contributions represent two new modes, decaying respectively as e−3η/2 and e−η compared with the growing one. To understand their physical effect we notice that an initial condition of the form

6

A. Elia, S. Kulkarni, C. Porciani, M. Pietroni, S. Matarrese Bin Bin Bin Bin Bin

1 2 3 4

Mass range (1013 M⊙ /h)

# haloes

kmax (Mpc−1 h)

b1

b2 (Mpc2 h−2 )

R (Mpc h−1 )

σ02 /σ12 (Mpc2 h−2 )

1.24 − 1.8 1.8 − 3.4 3.4 − 10 > 10

202948 211305 150105 48985

0.24 0.24 0.22 0.15

7.28 ± 0.38 14.2 ± 0.4 25.9 ± 0.4 66.2 ± 1.3

422 ± 102 356 ± 84 708 ± 103 1025 ± 401

2.7 ± 0.8 2.1 ± 0.7 2.9 ± 0.4 3.5 ± 0.8

9.1 12.3 20.5 50.9

Table 1. Mass range and number of the haloes in the four bins.

ϕ(k) ϕ(k) ϕa (k; 0) = ϕh (k) , ϕv (k)

(25)

evolves (using eq. (22)) into ϕa (k; η) given by ϕ ϕ ϕ + 2e−3η/2 (ϕ − ϕv ) + e−η (−3ϕ + ϕh + 2ϕv ) ϕ − e−3η/2 (ϕ − ϕv )

, (26)

i.e. both the halo density and velocity fields relax to the matter ones as η → ∞ (but at a different pace). Also note that in the absence of an initial density bias (i.e. φ3 = φ) but in the presence of an initial velocity bias (i.e. φ4 6= φ), the linear dynamics quickly generates a transient density bias that vanishes at late times as e−η −e−3η/2 . The initial power spectrum at η = 0, corresponding to the field configuration 0 (k), and it evolves forward in time as in eq. (25), is Pab 0 PL,ab (k; ηa , ηb ) = gac (ηa )gbd (ηb )Pcd (k) .

4.1

(27)

The importance of velocity bias

It is interesting to assess the role of the velocity bias in the linear solution previously discussed. Assuming ϕ4 = ϕ2 at all times, the linear propagator for the first three components ϕi with i = 1, 2, 3 becomes 3/5 2/5 0 gab (η) = 3/5 2/5 0 3/5 2/5 0 2/5 −2/5 0 −5/2η −3/5 3/5 0 +e 2/5 −2/5 0 0 0 0 0 0 θ(η) , (28) +e−η 0 −1 0 1 and the third component of (26) reduces to ϕ3 (k; η) = ϕ(k) + e−η (ϕh (k) − ϕ(k)) .

(29)

This expresses the well known linear debiasing between the halo and matter fields at late times, derived by Fry (1996) for tracers that do not undergo merging and move solely under the influence of gravity. The corresponding halo-matter cross spectrum is (3)

0 0 0 PL,13 = P11 + e−η (P13 − P11 ) ≡ PL,13 ,

(30)

while keeping ϕ4 6= ϕ2 one gets (3)

0 0 PL,13 = PL,13 + 2(P11 − P14 )(e−3η/2 − e−η ) .

(31)

Figure 1. A comparison between P14 (blue dashed line) and P11 (red solid line) in the four bins. The vertical black dotted lines represent the limit up to which we expect our fluid approximation for haloes to work. A smoothing scale of R = 7 Mpc/h has been used for P14 . For a fair comparison with P11 , which is not smoothed, P14 has been redivided by the smoothing function. 0 0 In Figure 1 we compare P11 against P14 extracted from the N-body simulation, for the different mass bins. While the spectra agree well on very large scales (k . 0.05 h Mpc−1 ), they progressively depart for smaller scales. This is in line with the model introduced in Section 3. Note that the last term in eq. (31) vanishes in the initial conditions, reaches a minimum for η ≈ 0.8 and it is suppressed at late times. We quantify its amplitude at z = 0 in Figure 2, where we (3) plot the ratio rL = PL,13 /PL,13 which ranges between 0 − 3 per cent, depending on halo mass and scale. This suggests that the effect of the velocity bias on the halo-matter cross spectrum is small for low redshifts.

5

ANALYTICAL TREATMENT OF NON-LINEARITIES

In this section we deal with the non-linear evolution of the matter-halo system. We first compute 1-loop corrections for the propagator and then perform the corresponding large-k resummation to all perturbative orders thus extending the

Halo clustering with resummed perturbation theories +(a ↔ b) , II ∆Pab (k; ηa , ηb )

=

Z

ηa

dsa

0

7 (35)

Z

ηb

dsb Φcd (k; sa , sb )

0

gac (k; ηa − sa ) gbd (k; ηb − sb ) , (36) with Φcd (k; sa , sb )

=

2esa +sb

Z

d3 q

γcei (k, −q, q − k)PL,ef (q, sa , sb ) PL,ih (|k − q|, sa , sb )γdf h (−k, q, k − q) . (37) The explicit expressions for

5.2

(3)

Figure 2. The ratio between PL,13 and PL,13 : going from bottom to top, Bin 1 (red), Bin 2 (blue), Bin 3 (green) and Bin 4 (black).

results presented by Crocce & Scoccimarro (2006) for the matter density field. Finally, we compute non-linear power spectra using the TRG approach. For simplicity we only consider the case with no velocity bias, which we have demonstrated to be accurate (at least in the linear regime) at low redshifts, where current observations are available. The 3×3 Ω matrix for this case is 1 −1 0 3 3 Ω= − (32) 0 . 2 2 0 −1 1 5.1

1-loop perturbation theory

The 1-loop correction to the linear propagator (Crocce & Scoccimarro 2006) is given by Z η Z s1 Z ∆gab (k; η) = 4 ds1 ds2 d3 q es1 +s2 0

0

gac (η − s1 )γcie (k, −q, q − k)gef (s1 − s2 ) γf hd (k − q, q, −k)gdb (s2 )PL,ih (q, s1 , s2 ) , (33) (see appendix A for its explicit expressions). The 1-loop correction to the two-point correlator (10) is given by the sum of two contributions I ∆Pab

(k; ηa , ηb ) +

II ∆Pab

(k; ηa , ηb ) ,

(34)

which are also known in the literature as “P13 ” and “P22 ”, respectively 2 . They are given by I ∆Pab

2

(k; ηa , ηb )

=

0 ∆gac (k; ηa ) gbd (k; ηb ) Pcd (k)

They should not be confused with the P13 and P22 of our notation!

II ∆Pab

are given in appendix A.

Large-k resummation for the propagator

In the large-k limit, the 1-loop correction for the propagator, eq. (33), grows as k2 , and eventually dominates over the (k-independent) linear propagator. Taking into account higher orders, the situation gets even worse. The 2-loop correction grows as k4 , the 3-loop as k6 , and so on. This a manifestation of the perturbative expansion breakdown in cosmological PT, which appears not only in the computation of the propagator, but also of the power spectrum, the bispectrum, and so on. However, for the case of the propagator, it was shown in Crocce & Scoccimarro (2006b) that the leading order corrections in the large k and large η limit can be resummed at all orders in perturbation theory, giving a well-behaved propagator. The propagator Gac (k; η) connects the initial correlators with the cross-correlations between final and initial field configurations, 0 NG Pab (k; η, 0) = Gac (k; η)Pcb (k) + ∆Pab (k; η, 0) ,

(38)

where the last term at the rhs comes from the initial nonGaussianity of the matter and halo fields. At leading order, it is given by (see appendix C) Z η NG ∆Pab (k; η, 0) = ds es gac (η − s)gdf (s)geg (s) 0 Z × d3 q γcde (k, −q, q − k)Bf gb (q, k − q, −k) , (39) where Babc is the initial bispectrum at η = 0 (z = zin ) (see also Bernardeau et al. 2010). In Section 6.1 we will use eq. (38) to assess the validity of different approximation schemes for the propagator. In the large-k limit, G decays as 2 2 2η k σ e Gab (k; η) = gab (η) exp − , (40) 2 with 1 σ = 3 2

Z

d3 q

P 0 (q) . q2

(41)

Therefore, at least in the case of the propagator, the bad ultraviolet behavior is just an artifact of the perturbative expansion, which, at any finite order, completely misses the nice –and physically motivated– Gaussian decay of eq. (40) (see Crocce & Scoccimarro 2006b and Matarrese & Pietroni 2007 for a detailed discussion).

8

A. Elia, S. Kulkarni, C. Porciani, M. Pietroni, S. Matarrese

Although the result (40) was obtained for the 2 × 2 propagator of the matter density-velocity system, it holds, taking into account proper modifications, also when halos are included, i.e. for the 3 × 3 propagator considered in this Section. As in Crocce & Scoccimarro (2006b), we obtain an improved propagator interpolating between the 1-loop result (eq. 33) at low k and the Gaussian decay (40) (with a modified pre-factor) at high k. The details of the derivation and the relevant formulae are given in appendix B.

5.3

TRG

Unlike the propagator, the power spectrum cannot be resummed analytically at large k. Different semi-analytical procedures (Crocce & Scoccimarro 2006, McDonald 2007, Taruya & Hiramatsu 2008) have been proposed to compute it in the mildly non-linear regime. In this paper we will consider the TRG technique introduced in Pietroni (2008). Starting from eq. (5), a hierarchy of differential equations for the power spectrum, the bispectrum and higher order correlations is obtained. We choose to truncate it at the level of the trispectrum Qabcd = 0, so that the equations for Pab and Babc form a closed system ∂η Pab (k; η) = − Ωac (η)Pcb (k; η) − Ωbc (η)Pac (k; η) Z +eη d3 q [γacd (k, −q, q − k) Bbcd (k, −q, q − k; η) +Bacd (k, −q, q − k; η) γbcd (k, −q, q − k)] ,

The vertical black dotted lines represent the limit up to which we expect our fluid approximation for the haloes to work.

∂η Babc (k, −q, q − k; η) = −Ωad (η)Bdbc (k, −q, q − k; η) −Ωbd (η)Badc (k, −q, q − k; η) −Ωcd (η)Babd (k, −q, q − k; η) +2eη [γade (k, −q, q − k)Pdb (q; η)Pec (|k − q|; η) +γbde (−q, q − k, k)Pdc (|k − q|; η)Pea (k; η) + γcde (q − k, k, −q)Pda (k; η)Peb (q; η)] ,

(42)

which integrated gives the power spectra at any redshift and for any momentum scale. The system (42) consists of coupled differential equations which are solved numerically, starting from given initial conditions, i.e. Pab (k; ηi ) and Babc (k, −q, q − k; ηi ). From eq. (32), we can observe that Ω13 and Ω23 are zero, which means that the evolution of δm and θ is not modified, with respect to the original TRG formulation, by the presence of δh , as it is expected.

6 6.1

Figure 3. Cross spectrum between the halo density at z = 0 and the matter density at z = 50. The outcome of the N-body simulation (black points with error bars) is compared against lineartheory PP L (blue dashed line) and resummed result PP R (green solid line). The red dotted line shows the effect of neglecting the NG . non-Gaussian term, i.e. PP R − ∆P31

RESULTS

do not have to deal with the shot-noise problem. We extract the cross spectra from the simulation and compare them against those obtained both using linear theory prop0 0 0 agators PP L ≡ g31 P11 + g32 P21 + g31 P31 and resummed 0 0 0 NG propagators PP R ≡ G31 P11 + G32 P21 + G31 P31 + ∆P31 ; the result is shown in Figure 3. We note that the linear model severely overpredicts the two-time cross spectrum for k > 0.05 h Mpc−1 . It is evident that the resummed theory improves considerably upon the linear one, and agrees with the simulation within 10 per cent accuracy up to the scale where the fluid approximation holds. We include in NG ∆P31 the effect of the initial non-Gaussianity of the halo field via its initial bispectra, computed as in eq. (17). It turns out that the components giving a non-vanishing effect are of the B113 type. Their contribution is suppressed by a D+ (zin )/D+ (z = 0) factor with respect to that of an hypothetical primordial non-Gaussianity in the matter field (which we do not consider here). Therefore, the effect is of modest entity but, nevertheless, it improves the agreement with the simulation with respect to the case in which it is neglected.

Propagator

In order to assess the validity of our analytical approach, we compare our results for the resummed propagator against the simulation; to this end, we consider the relation in eq. (38). In particular we choose the indices a = 3, b = 1, so that we can check the components related to the haloes that were not present in the original formalism by Crocce & Scoccimarro (2006b) and, at the same time, we

6.2

Power spectrum

The TRG equations presented above are integrated numerically starting from the initial conditions discussed in section 3. As a first step, we set all the initial bispectra to zero In Figure 4 we show a comparison between the halo-matter cross spectra extracted from the N-body simulation and

Halo clustering with resummed perturbation theories

9

investigate the case of P11 . In the one-loop computation, the initial P11 contributes to the final one with a “weight” of (D+ (z = 0)/D+ (zin ))2 = e2η , as it also happens in linear theory. The initial matter bispectrum B111 , instead, has a weight of eη , so it is suppressed by a factor of e−η 3 . If we now consider the haloes, we have showed in (28) that there is a new decaying mode, responsible for the linear debiasing. This new mode, that involves only the halo field, carries an extra e−η suppression factor. We can now rank the contributions to P13 according to their relevance: (i) P11 (ii) P13 , B111 (iii) P33 , B113 (iv) B133 (v) B333 .

Figure 4. The cross spectrum between matter and halo distribution at z = 0 is shown in the four bins; the black dots with error bars represent the simulation, the blue dashed line is linear theory, the violet dot-dashed line is 1-loop and the red solid line is TRG. The vertical black dotted lines represent the limit up to which we expect our fluid approximation for the haloes to work.

the results of TRG, one-loop and linear theory. In the first three bins, corresponding to lower halo masses, linear theory overpredicts the power on mildly non-linear scales; note that this departure arises on larger scales compared to the matter auto spectrum (not shown in the figure). The overprediction of linear theory is cured by the one-loop power spectrum only on very large scale, while the TRG manages to correct it up to a smaller scale, before starting to fail. The fourth bin, though, displays a different behaviour: linear theory lacks of power on small scales, and neither the 1-loop correction nor the TRG are much more accurate. This might originate from the fact that very massive haloes are large and rare in the initial conditions, therefore less suited for the fluid approximation. Moreover, they also display the strongest velocity bias, which we are neglecting in our current non-linear treatment. A more quantitative analysis is presented in Figure 5, where we plot the ratios between the spectra from the simulation and the theoretical results. The TRG gives a cross spectrum within 5% accuracy at least up to k = 0.1 h Mpc−1 (barring bin 4), while linear theory does so up to k = 0.05 h Mpc−1 . As previously pointed out, however, the halo-density field is not Gaussian at zin = 50. We use expressions (17) to account for initial bispectra in the TRG method (they are not present in the linear theory and at 1-loop level). Figure 5 illustrates that this indeed produces a slight improvement in the agreement with the simulation. The correction resulting from the introduction of the initial bispectra turns out to be quite small. The reason is that their contribution to the final cross spectrum P13 is suppressed. To understand why, we can use perturbation theory; first, let us

Each item is suppressed by a factor of e−η with respect to the previous one. We can see that only B111 has some relevance, while the other terms are highly suppressed. Even if the reasoning was based on perturbation theory, it is valid also for TRG, at a qualitative level. Incidentally, the fact that P11 is the most relevant contribution is another evidence of debiasing. We can now address the effect of the truncation we perform in the TRG, namely considering the trispectrum Q = 0. First of all, the matter trispectrum Q1111 can be neglected in the range of scales under consideration, as one can conclude from the comparison between TRG and simulations in Pietroni (2008). The contribution from initial mixed (i.e. matter-halo) or pure halo trispectra is furtherly suppressed with respect to Q1111 by extra e−η factors, for the same reason as above. However, the trispectrum has its own time evolution as well, and one might argue that it becomes more relevant for z < zin ; this seems to be not the case, because enlarging the TRG truncation scheme by including the running of the trispectrum gives a contribution to the power spectrum which is at least of two-loop order and it is certainly subdominant in the scales we are considering. From Section 5 onward, we neglected any velocity bias, since this approximation proved to be accurate enough at z = 0 (in linear theory). As a further check, it is interesting to observe from Figure 6 that the TRG is able to give a better prediction than full linear theory in eq. (31), even neglecting the velocity bias that the linear theory accounts for. We can now look at the model predictions for the halo bias, defined as the ratio P13 /P11 . This quantity is plotted as a function of wavenumber in Figure 7. While the linear-theory bias always increases with scale, irrespectively of halo mass, the TRG result closely follows the scale dependence of b(k) seen in the simulation for the first two bins. It also gives a nearly constant bias for the third bin, as the simulation does, even though with a slightly lower value. However, the linear model performs better in the last bin where the bias in the simulation increases with k. It is also interesting to look at the results displayed in a different way; we can investigate the evolution of the cross spectrum from the initial conditions to today and the

3 For the most experienced readers, this happens because one of the two vertices carrying the eη factor is replaced by the bispectrum.

10

A. Elia, S. Kulkarni, C. Porciani, M. Pietroni, S. Matarrese

Figure 5. The ratio between the spectrum from simulation and: linear theory (blue squares), TRG without bispectra (red triangles), TRG with bispectra (green circles). The vertical black dotted lines represent the limit up to which we expect our fluid approximation for haloes to work. The shaded area marks the 5 per cent accuracy interval.

Figure 7. The effective bias b ≡ Pmh /Pm as a function of the wavevector in the four bins: the black circles with error bars are from the simulation, the blue dashed line represents linear theory and the red solid line TRG. The vertical black dotted lines represent the limit up to which we expect our fluid approximation for the haloes to work.

evolution of the bias as well. In Figures 8 and 9 we plot, respectively, rp ≡

Pmh (z = 0) Pmh (z = 50)

and

rb ≡

b(z = 0) . b(z = 50)

(43)

Again, our model is able to match accurately the trend of the simulation and to improve upon linear theory, excluding bin 4. A key feature of the linear solution in eq. (26) is the debiasing between halo and matter distributions with time. It is worth noting from Figure 9 that this effect is stronger for high-mass haloes, but constant on all the scales, while for low-mass haloes it is weaker on large scales and presents a strong k-dependence.

7

Figure 6. The ratio between the cross spectrum from simulation and TRG with bispectra (green) and linear theory with the inclusion of velocity bias (magenta). The vertical black dotted lines represent the limit up to which we expect our fluid approximation for haloes to work.

CONCLUSIONS

We have presented a novel approach to modeling the clustering of dark-matter haloes on mildly non-linear scales. This follows the motion of the regions that will collapse to form haloes (that we dub proto-haloes). Since the number of proto-haloes is conserved over time, for sufficiently large scales (k < 0.2 h Mpc−1 ), we can write a set of fluid equations that govern their evolution under the effect of gravity, which couples perturbations in the halo and matter density fields. We provide analytical solutions for the linearized equations and 1-loop perturbative corrections for the halo and matter power spectra. For the propagator, quantifying the memory of the density and velocity fields to their initial conditions, we also perform a resummation of perturbative corrections. Finally, for the power spectrum we compute

Halo clustering with resummed perturbation theories

11

the non-linear evolution using a semi-analytical procedure, namely an extension of the time renormalization group. The initial conditions for our evolutionary equations are specified adopting a Lagrangian bias model, originally developed to describe the clustering and motion of linear density peaks. We fix the parameters of the model so that to reproduce the distribution of proto-haloes in a high-resolution N-body simulation at z = 50. We use the same simulation to test the predictions of our model at z = 0. Our main results can be summarized as follows:

Figure 8. The ratio between the halo-matter cross spectrum at z = 0 and z = 50 for the four mass bins: symbols represent the simulation, dashed lines the linear theory and solid lines TRG. From top to bottom, we have Bin 1 in red, Bin2 in blue, Bin 3 in green and Bin 4 in black.

• Independently of the initial conditions, in the linear solution the halo density and velocity fields asymptotically match the corresponding matter fields at late times. This ’debiasing’ develops at a different rate for the density and the velocity, being faster for the latter. • Even if there is no initial density bias, the presence of a velocity bias generates a transient density bias that vanishes at late times. • Neglecting any initial velocity bias alters the linear predictions for the halo-matter cross spectrum at redshift z = 0 only by less than 3 per cent, for k < 0.3 h Mpc−1 . This provides us with the motivation to ignore the velocity bias in the non-linear analysis. • Unlike its linear counterpart, the resummed propagator is in good agreement with the N-body simulation, independently of halo mass. • The halo-matter cross spectrum predicted by the TRG is accurate to 5 per cent up to k ≃ 0.1 h Mpc−1 for a broad range of halo masses. This does not hold for very massive haloes (M > 1014 h−1 M⊙ ), that have low number density and high initial velocity bias, for which discreteness effects are more important. • The TRG result improves upon both linear theory and 1-loop corrections. Its performance is slightly enhanced accounting for the initial non-Gaussianity of the halo distribution. • For low halo masses our model accurately describes the scale-dependent bias measured in the simulation at z = 0.

ACKNOWLEDGMENTS AE was supported through the SFB-Transregio 33 “The Dark Universe” by the Deutsche Forschungsgemeinschaft (DFG). AE and SK acknowledge BCGS support for part of the work. CP thanks the participants to the workshop “Modern Cosmology: Early Universe, CMB and LSS” held at the Centro de Ciencias de Benasque Pedro Pascual for stimulating discussions. Part of this work has been developed there in August 2010.

Figure 9. The ratio between the bias at z = 0 and z = 50 for the four mass bins: symbols represent the simulation, dashed lines the linear theory and solid lines TRG. From top to bottom, we have Bin 1 in red, Bin2 in blue, Bin 3 in green and Bin 4 in black.

REFERENCES Bardeen J.M., Bond J.R., Kaiser N., Szalay A.S., 1986, ApJ, 304, 15 Bartolo N., Beltran Almeida J.P., Matarrese S., Pietroni M., Riotto A., 2010, JCAP 1003, 011 Bernardeau F., Crocce M., Sefusatti E., 2010, Phys. Rev. D, 82, 083507 Buchert T., Dominguez A., 2005, A& A, 438, 443

12

A. Elia, S. Kulkarni, C. Porciani, M. Pietroni, S. Matarrese

Catelan P., Lucchin F., Matarrese S., Porciani C., 1998, MNRAS, 297, 692 Catelan P., Porciani C., Kamionkowski M., 2000, MNRAS, 318, L39 Cole S., Kaiser N., 1989, MNRAS, 237, 1127 Cole S. et al., 2005, MNRAS, 362, 505 Coles P., 1993, MNRAS, 262, 1065 Cooray A., Sheth R.K., 2002, PhR, 372, 1 Crocce M., Scoccimarro R., 2006, Phys. Rev. D, 73, 063519 Crocce M., Scoccimarro R., 2006, Phys. Rev. D, 73, 063520 Desjacques V., 2006, Phys. Rev. D, 78, 103503 Desjacques V., Sheth R.K., 2010, Phys. Rev. D, 81, 023526 Desjacques V., Crocce M., Scoccimarro R., Sheth R.K., 2010, Phys. Rev. D, 82, 103529 Eisenstein D.J., Hu W., 1998, ApJ, 496, 605 Elia A., Ludlow A., Porciani C., in preparation Fry J.N., 1996, ApJ, 461, L65 Fry J.N., Gazta˜ naga E., 1993, ApJ, 413, 447 Guzik J., Bernstein G., Smith R.E., 2007 MNRAS, 375, 1329 Huff E., Schulz A. E., White M., Schlegel D.J., Warren, M.S., 2007, APh, 26, 351 Kaiser N., 1984, ApJ, 284, L9 Lewis A., Challinor A., Lasenby A., 2000, ApJ, 538, 473 Ludlow A. D., Porciani C., 2010, arXiv 1011.2493 Manera M., Gazta˜ naga E., 2009, arXiv 0912.0446 Manera M., Sheth R.K., Scoccimarro R., 2010, MNRAS, 402, 589 Mann, R. G., Peacock, J. A., Heavens, A. F., 1998, MNRAS, 293, 209 Matarrese S., Pietroni M., 2007, JCAP, 6, 26 Matsubara T., 1999, ApJ, 525, 543 Matsubara T., 2008, Phys. Rev. D, 78, 3519 McDonald P., 2007, Phys. Rev. D, 75, 043514 Mo H.J., White S.D.M., 1996, MNRAS, 282, 347 Montesano F., Sanchez A.G., Phleps S., 2010, arXiv 1007.0755 Norberg P. et al., 2002, MNRAS, 332, 827 Pietroni M., 2008, JCAP, 10, 36 Pillepich A., Porciani C., Hahn O., 2010, MNRAS, 402, 191 Porciani C., Catelan P., Lacey C., 1999, ApJ, 513, L99 Porciani C., Dekel A., Hoffman Y., 2002, MNRAS, 332, 339 Schulz A. E., White, M., 2006, APH, 25, 172 Scoccimarro R., 2004, Phys. Rev. D, 70, 083007 Seljak U., 2001, MNRAS, 325, 1359 Seo H., Eisenstein D.J., 2005, ApJ, 633, 575 Smith R.E., Scoccimarro R., Sheth R.K., 2007, Phys. Rev. D, 75, 063512 Sugiyama N., 1995, ApJ, 100, 281 Taruya A., Hiramatsu T., 2008, ApJ, 674, 617 Zehavi I. et al., 2004, ApJ, 608, 16

Halo clustering with resummed perturbation theories

13

APPENDIX A: COMPLETE ANALYTIC EXPRESSIONS OF THE 1-LOOP PROPAGATOR AND POWER SPECTRUM Linear power spectra: PL,11 (k; η, η)

=

0 P11 (k) = PL,12 (k; η, η) = PL,22 (k; η, η) ,

PL,13 (k; η, η)

=

0 0 e−η ((eη − 1)P11 (k) + P13 (k)) ,

PL,33 (k; η, η)

=

0 0 0 e−2η ((eη − 1)2 P11 (k) + P33 (k) + 2(eη − 1)P13 (k)) .

A1

(A1)

Propagator

First we report the integrands of the one-loop corrections of the components we are interested in. The linear power spectrum can be split as PL,ci ,cj (qi ; si , sj ) = P 0 (qi )Uci ,cj + ∆Pci ,cj (qi ; si , sj ) ,

(A2)

where Uci ,cj = 1

for any ci , cj ,

(A3)

while ∆Pci ,cj (qi ; si , sj ) 6= 0

only if ci or cj = 3 .

(A4)

(1) (2) Looking at eq. (33) we can consistently split R the one-loop correction to the propagator as ∆gab = ∆ gab + ∆ gab and we denote them as δ (i) gab , so that ∆(i) gab = dq δ (i) gab , for i = 1, 2. h i 3 (k+q)2 0 π (eη − 1)2 P11 (q) 4 6k7 q − 79k5 q 3 + 50k3 q 5 − 21kq 7 − 3 k2 − q 2 2k2 + 7q 2 log (k−q) 2 δg11 = , 840k3 q 3

δ (1) g31 δ

(1)

=

g32

=

δ (1) g33

=

δ (2) g31

=

δ (2) g32

=

δ (2) g33

=

A2

δg11 , 2 δg11 , 3 2 0 − e−η (eη − 1)2 k2 πP11 (q) , 3 3 k+q 0 0 πeη P13 (q) − P11 (q) 4 9k5 q + 24k3 q 3 − 9kq 5 − 18 k2 − q 2 log k−q 280k3 q

,

2 (2) δ g31 , 3 0

(A5)

1-loop PS

From eq. (37), one can see that the one-loop corrections require integrals over ds1 , ds2 and d3 q = q 2 sin θ dq dθ dφ; the integration over q and over µ ≡ cos θ cannot be performed analytically. The following expressions are p therefore integral kernels, denoted by δqµ P (the factor of q 2 is already taken into account). If cos θ = µ = k·q and |k − q| = k2 − 2kµq + q 2 , kq 2 0 0 πk4 e2η P11 (q) 7kµ + 3 − 10µ2 q P11 (|k − q|) II δqµ P11 = , 49 (k2 − 2kµq + q 2 )2 (A6) II δqµ P13

=

πk3 eη 7kµ + 3 − 10µ2 q 49 (k2 − 2kµq + q 2 )2

0 0 0 0 7q(k − µq)P13 (q)P11 (|k − q|) + P11 (q)7µ k2 − 2kµq + q 2 P13 (|k − q|)

0 0 +kP11 (q)P11 (|k − q|) 7kµ (eη − 1) − q

II δqµ P33

=

10µ2 − 3 eη − 14µ2 + 7 .

0 πk2 2 2 0 P13 (q)P13 (|k − q|) 2 98µq(k − µq) k − 2kµq + q 2 − 2kµq + q ) 0 0 0 +7q(k − µq)P11 (|k − q|) 2kP13 (q) 7kµ (eη − 1) − q 10µ2 − 3 eη − 14µ2 + 7 + 7qP33 (q)(k − µq) 2 0 0 +k2 P11 (q)P11 (|k − q|) q 10µ2 − 3 eη − 14µ2 + 7 − 7kµ (eη − 1) 2 0 0 +49µ2 P11 (q) k2 − 2kµq + q 2 P33 (|k − q|) 0 0 2 2 +14kµP11 (q) k − 2kµq + q 7kµ (eη − 1) − q 10µ2 − 3 eη − 14µ2 + 7 P13 (|k − q|) .

(A7)

49 (k2

(A8)

14

A. Elia, S. Kulkarni, C. Porciani, M. Pietroni, S. Matarrese

APPENDIX B: THE RESUMMED PROPAGATOR As it was discussed in Crocce & Scoccimarro (2006b), the leading contribution in k2 exp(2η) to the propagator at a generic n-loop order contains a chain of propagators and vertices of the form ga b1 (η − s1 )γb1 c1 a1 (k, −q1 , q1 − k)ga1 b2 (s1 − s2 ) · · · ga2n−1 b2n (s2n−1 − s2n )γb2n c2n a2n (k + q2n , −q2n , −k)ga2n b (s2n ) . (B1) The ci indices have to be contracted in all possible pairings, by inserting n linear power spectra, each of the form δD (qi + qj )PL,ci ,cj (qi ; si , sj ) .

(B2)

The n-loop contribution is obtained by multiplying by exp( integrating over Z η Z dsi d3 qi . Π2n i=1

P2n

i=1

si ) and by the appropriate combinatoric factor, and then by

(B3)

0

Recall eqs. (A2-A4): if all the insertions are of the P 0 (qi )Uci ,cj type then the resummation goes exactly as in the standard case. Indeed 1q·k Uc1 ,cj gab1 (η − s1 )γb1 c1 a1 (k, −q1 , q1 − k)ga1 b2 (s1 − s2 ) → ucj ga b2 (η − s2 ) , (B4) 2 q2 in the k ≫ q limit. Besides the explicit form for the vertices (see Section 1), we have used the composition property of the propagators gab (η − s1 )gbc (s1 − s2 ) = gac (η − s2 ) ,

(B5)

and have defined u = (1, 1, 1). Since each vertex, contracted by a ucj vector becomes proportional to a delta function in its first and third index, the chain of propagators composes up to a single propagator gab (η) and the time integral can be easily performed. The momentum integrals factorize into n integrals of the type Z (k · q)2 − d3 qP 0 (q) = −k2 σ 2 , (B6) q4 where σ 2 has been defined in (41). Using the appropriate combinatoric factors, the leading contribution to the propagator at n-loop, in the large momentum limit, when all the power spectrum insertions are of the P 0 (qi )Uci ,cj type is n 1 (eη − 1)2 gab (η) , (B7) −k2 σ 2 n! 2 η

2

]. which resums to gab (η) exp[−k2 σ 2 (e −1) 2 As for the insertions including the ∆Pci ,cj (qi ; si , sj ) contribution to the linear power spectrum, the important point to realize is that, due to (A4) and to the structure of the linear propagator and the vertices, a chain like (B1) can be contracted by at most one ∆Pci ,cj (qi ; si , sj ) power spectrum, the remaining ones being of the type P 0 (qi )Uci ,cj . Moreover, these insertions only contribute to G31 and G32 . Therefore, at n-th order, we have only two types of contributions: those with all P 0 (qi )Uci ,cj insertions, giving (B7), and those with one ∆Pci ,cj insertions and n − 1 P 0 (qi )Uci ,cj ones, which can also be resummed in the large k limit. As a consequence, the complete resummed propagator in the large momentum limit is Gab (k; η) = (gab (η) + δa3 fb ∆(2) g31 (k; η)) exp[−k2 σ 2

(eη − 1)2 ], 2

(B8)

with fb = (1, 2/3, 0). In order to have an expression for the propagator valid at any k, one can proceed as in (Crocce & Scoccimarro (2006b)) and interpolate between the large k limit above and the 1-loop result at low k. This can be done, for instance for G11 , starting from its 1-loop expression 5 g11 (η) + ∆g11 (k; η) ≃ g11 (η) 1 + ∆g11 (k; η) , (B9) 3 where the above approximation is exact in the large η limit. Since (eη − 1)2 5 ∆g11 (k; η) → −k2 σ 2 , for large k , 3 2 the required interpolation is given by 5 G11 (k; η) = g11 (η) exp ∆g11 (k; η) . 3 Proceeding analogously for the other components, we have: 5 G12 (k; η) = g12 (η) exp ∆g11 (k; η) , 3

(B10)

(B11)

Halo clustering with resummed perturbation theories G22 (k; η)

=

G21 (k; η)

=

G31 (k; η)

=

G32 (k; η)

=

G33 (k; η)

=

5 ∆g22 (k; η) , 2 5 g21 (η) exp ∆g22 (k; η) , 2 5 g31 (η) + ∆(2) g31 (k; η) exp ∆g11 (k; η) , 3 2 5 g32 (η) + ∆(2) g31 (k; η) exp ∆g11 (k; η) , 3 3 h i η (1) g33 (η) exp e ∆ g33 (k; η) , g22 (η) exp

15

where we have used the identities 2 ∆g11 (k; η) , ∆g12 (k; η) = 3 3 ∆g21 (k; η) = ∆g22 (k; η) , 2 (1) ∆ g31 (k; η) = ∆g11 (k; η) 2 (1) ∆(1) g32 (k; η) = ∆ g31 (k; η) 3 2 (2) ∆(2) g32 (k; η) = ∆ g31 (k; η) . 3

(B12)

(B13)

APPENDIX C: TAKING INTO ACCOUNT THE INITIAL HALO NON-GAUSSIANITY. While at zin the matter field can be considered gaussian with high accuracy, the same not necessarily holds for haloes. An initial non-Gaussianity for the equal-time power spectrum can be easily incorporated in the TRG formalism as discussed in Bartolo et al. (2010), and has been done in sect. 6.2. In order to obtain the effect of an initial non-Gaussianity on the cross-correlator at different times considered in eq. (38), it suffices to compute the O(γ) correction to the linear evolution of the field, by inserting eq. (22) in (5), to get Z η Z ϕ(1) (C1) ds gac (η − s)es d3 q γcde (k, −q, q − k)gdf (s)ϕf (q; 0)geg (s)ϕg (k − q; 0) . a (k; η) = 0

The cross-correlator ′ hϕ(1) a (k; η)ϕb (k ; 0)i ,

(C2)

then includes the non-Gaussian expectation value hϕf (q; 0)ϕg (k − q; 0)ϕb (k′ ; 0)i = δD (k + k′ )Bf gb (q, k − q, −k) , and therefore gives eq. (39). This paper has been typeset from a TEX/ LATEX file prepared by the author.

(C3)