Modelling the clustering of dark matter haloes in resummed perturbation theories

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 (protohalo). We adopt a fluid-like description for the evolution of perturbations in the protohalo 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 to resum the power spectrum. For halo masses M < 10¹⁴h⁻¹ 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⁻¹ approaching the requirements of future galaxy redshift surveys.

methods: analytical, galaxies: haloes, cosmology: theory, large-scale structure of Universe

1 INTRODUCTION

Redshift surveys have shown that the clustering properties of galaxies strongly depend on their luminosity, colour 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 dark-matter-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: P_g=b²₁P_m, where the linear bias factor b₁ depends on galaxy type but is generally scale independent (e.g. Coles 1993; Mann, Peacock & Heavens 1998). Similarly, to model higher order statistics, such as the galaxy bispectrum, it is generally assumed that galaxy biasing is a local process such that δ_g=b₁δ_m+b₂δ²_m/2 +⋯, where δ_g and δ_m are the (smoothed) galaxy and dark matter density contrast, respectively (Fry & Gaztañaga 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 per cent or even sub-per cent 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 & Gaztañaga 2010; Manera, Sheth & Scoccimarro 2010; Montesano, Sanchez & Phleps 2010) or analytical calculations (Seljak 2001; Schulz & White 2006; Guzik, Bernstein & Smith 2007; Smith, Scoccimarro & Sheth 2007) stemming from either perturbation theory (PT) 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. P_g(k) =b(k)²P_m(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 peak-background 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, Catelan & Lacey 1999) and, secondly, 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 top-hat model. A more sophisticated generalization of the peak-background split has been presented by Catelan et al. (1998) who assumed that 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, Porciani & Kamionkowski 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 P_m and P_h using the time-renormalization-group (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⁻¹.

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 dark matter 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 such 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 one-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 the 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

1

where τ is the conformal time. If we define the velocity divergence θ (x, τ) =∇·v (x, τ) and switch to Fourier space, the equations in (1) become

2

where

3

Equations (2) can be written in a compact form if we define a new time variable η≡ ln (D₊/D_+in), D_+in being the growth factor at an early epoch, and a doublet φ_a (a = 1, 2):

4

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 ⁠, i.e. φ₁(k, 0) =φ₂(k, 0). The system is therefore

5

where the 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,

6

and γ₁₁₂(k, q, p) =γ₁₂₁(k, p, q), with δ_D the Dirac delta distribution. All the information about the cosmological model is contained in the matrix

7

that, in the following, will be considered as a constant matrix, approximating Ω_m/f²₊≈ 1. This ratio is indeed very close to unity for most of the history of the Universe. Making this approximation, one is modifying the behaviour 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 per cent level up to k≈ 0.3 h Mpc⁻¹.

The power spectrum, defined by an ensemble average, in this notation is a 2 × 2 matrix

8

and the bispectrum, defined by

9

has eight components. In the following, we will also consider a different-time two-point correlator, defined as

10

which obviously coincides with equation (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 z_id according to some pre-defined 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 protohalo. In other words, a protohalo is the Lagrangian region of space that will collapse to form a halo at redshift z_id. N-body simulations show that nearly all protohaloes are simply connected (Porciani, Dekel & Hoffman 2002) even though this property is not key to our analysis.

Let us now follow the evolution of a protohalo over cosmic time in Eulerian space. Basically, its shape and topology will be distorted (protohaloes 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 on to this motion that connects the Lagrangian position of the protohalo with the Eulerian location of the final halo.

On scales much larger than the characteristic size of (and separation between) the protohaloes, the density fluctuations traced by the centre-of-mass trajectories can be described with a continuous overdensity field δ_h(x, τ | z_id). Note that while τ labels conformal time along the trajectories, z_id 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:

11

Here the protohalo density and velocity fields should be intended as coarse grained on a scale of a few times the mean interhalo separation (so as to suppress discreteness effects as protohaloes 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 protohalo fluid velocities as

12

where the gravitational potential is the same as in equation (1). Note that if v_h matches v in the initial conditions, then it will always do. On the contrary, any initial velocity bias between protohaloes and matter will be progressively erased by the gravitational acceleration.

Thus, given the suitable initial conditions for δ_h and v_h at τ→ 0 (i.e. a prescription for the Lagrangian biasing of protohaloes), we can in principle use equations (11) and (12) to follow the clustering evolution of the protohalo population at all times. We are particularly interested in the solution of the fluid equations at the special time τ that corresponds to z_id. 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 equations (11) and (12). We define a quadruplet φ_a (a = 1, 2, 3, 4) as

13

in such a way that equation (5) still holds, but with indices running from 1 to 4. There are three more non-vanishing elements of the vertex γ₃₄₃(k, p, q) =γ₃₃₄(k, q, p) =γ₁₂₁(k, p, q) and γ₄₄₄(k, q, p) =γ₂₂₂(k, p, q), and the 4 × 4 matrix is

14

From 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 P₁₁ and the matter–halo cross-spectrum P₁₃.

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 as 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 protohaloes (and, later, to test our results at z = 0), we use one high-resolution N-body simulation extracted from the suite presented by Pillepich et al. (2010). This contains 1024³ dark matter particles within a periodic cubic box with a side of L_box = 1200 h⁻¹ Mpc and follows the formation of structure in a ΛCDM model with Gaussian initial conditions and cosmological parameters: h = 0.701, σ₈ = 0.817, n_s = 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 × 10¹³h⁻¹ 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, which, 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.

Table 1

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

Bin	Mass range (1013 M⊙/h)	\# haloes	kmax (Mpc−1h)	b1	b2 (Mpc2h−2)	R (Mpc h−1)	σ20/σ21 (Mpc2h−2)
Bin 1	1.24–1.8	202 948	0.24	7.28 ± 0.38	422 ± 102	2.7 ± 0.8	9.1
Bin 2	1.8–3.4	211 305	0.24	14.2 ± 0.4	356 ± 84	2.1 ± 0.7	12.3
Bin 3	3.4–10	150 105	0.22	25.9 ± 0.4	708 ± 103	2.9 ± 0.4	20.5
Bin 4	>10	48 985	0.15	66.2 ± 1.3	1025 ± 401	3.5 ± 0.8	50.9

Protohaloes 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 protohalo is used as a proxy for its spatial location. Similarly, the mass-weighted linear velocity gives the protohalo velocity.

Halo and protohalo density and momentum fields are computed with the cloud-in-cell grid assignment using a 512³ 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 Fast Fourier Transform (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 haloes, 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 haloes. The latter option offers a number of advantages. First, it allows us to make general predictions independently of the simulation specifics. Secondly, 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 protohaloes. 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 protohaloes in Lagrangian space δ_h(q) and the corresponding mass-density fluctuation δ (q). We assume that their Fourier transforms are linked by the expression

15

which corresponds to a non-local relation in real space. This form was first proposed by Matsubara (1999) and describes the clustering of linear density peaks (Desjacques 2008). In this case, the initial conditions for P₃₃ and P₁₃ are

16

where the exponential functions account 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 protohaloes and matter in our N-body simulation when b₁, b₂ and R are treated as fitting parameters (see Elia, Ludlow & Porciani, in preparation, for further details).1 This is not surprising as Ludlow & Porciani (2011) have shown that most protohaloes 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 b₁, b₂ and R that best fit the simulation data using equation (16), where the linear matter power spectrum P₁₁ is computed using the camb online tool (Lewis, Challinor & Lasenby 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:

17

where the matter bispectrum B_m(k₁, k₂, k₃) is computed using the tree-level expression of the standard PT,

18

with μ₁₂≡k₁·k₂/(k₁k₂).

In order to solve our evolutionary equations, we need to know also the linear velocity field of protohaloes. In principle, protohaloes 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 protohaloes 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

19

with b_v the scale-dependent ‘velocity bias’ and σ_n (n = 0, 1) being the spectral moments of the matter power spectrum defined as

20

In the above equations, a Gaussian smoothing window of size R has been adopted to coarse grain the matter density and identify the peaks. Once the protohalo mass is linked to R through (with the mean comoving density of matter), this model has no free parameters, since the spectral moments are completely defined by P₁₁(k). It follows that the initial conditions for P₂₄ and P₄₄ are

21

and these expressions are in very good agreement with the spectra computed from the N-body simulation (Elia et al., in preparation). The corresponding values of σ²₀/σ²₁ 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 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

22

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

23

with δ_ab the Kronecker delta. Solving equation (23) with causal boundary conditions [g_ab(η) = 0 for η < 0; see e.g. Crocce & Scoccimarro 2006a] one gets

24

with θ(η) Heaviside’s step function. Note that g_ab(η) →δ_ab as η→ 0⁺. The first and second contributions represent the standard growing and decaying modes, respectively (Crocce & Scoccimarro 2006a). The third and fourth contributions represent two new modes, decaying, respectively, as e^−3η/2 and e^−η are compared with the growing one. To understand their physical effect, we note that an initial condition of the form

25

evolves (using equation 22) into φ_a(k; η) given by

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. ϕ₃=ϕ) but in the presence of an initial velocity bias (i.e. ϕ₄≠ϕ), 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 in equation (25), is P⁰_ab(k), and it evolves forward in time as

27

4.1 The importance of velocity bias

It is interesting to assess the role of the velocity bias in the linear solution previously discussed. Assuming φ₄=φ₂ at all times, the linear propagator for the first three components φ_i with i = 1, 2, 3 becomes

28

and the third component of equation (26) reduces to

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

30

while keeping φ₄≠φ₂ one gets

31

In Fig. 1 we compare P⁰₁₁ against P⁰₁₄ 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⁻¹), they progressively depart for smaller scales. This is in line with the model introduced in Section 3. Note that the last term in equation (31) vanishes in the initial conditions, reaches a minimum for η≈ 0.8 and is suppressed at late times. We quantify its amplitude at z = 0 in Fig. 2, where we plot the ratio r_L=P_L,13/P⁽³⁾_L,13 which ranges between 0 and 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.

Figure 1

A comparison between P₁₄ (blue dashed line) and P₁₁ (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 P₁₄. For a fair comparison with P₁₁, which is not smoothed, P₁₄ has been redivided by the smoothing function.

Figure 2

The ratio between P_L,13 and P⁽³⁾_L,13: going from bottom to top, bin 1 (red), bin 2 (blue), bin 3 (green) and bin 4 (black).

5 ANALYTICAL TREATMENT OF NON-LINEARITIES

In this section, we deal with the non-linear evolution of the matter–halo system. We first compute one-loop corrections for the propagator and then perform the corresponding large-k resummation to all perturbative orders, thus extending the results presented by Crocce & Scoccimarro (2006a) 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

32

5.1 One-loop perturbation theory

The one-loop correction to the linear propagator (Crocce & Scoccimarro 2006a) is given by

33

(see Appendix A for its explicit expressions).

The one-loop correction to the two-point correlator (10) is given by the sum of two contributions

34

which are also known in the literature as ‘P₁₃’ and ‘P₂₂’, respectively.2 They are given by

35

36

with

37

The explicit expressions for ΔP^II_ab are given in Appendix A.

5.2 Large-k resummation for the propagator

In the large-k limit, the one-loop correction for the propagator, equation (33), grows as k², and eventually dominates over the (k-independent) linear propagator. Taking into account higher orders, the situation gets even worse. The two-loop correction grows as k⁴, the three-loop as k⁶ and so on. This is a manifestation of the perturbative expansion breakdown in cosmological PT, which appears not only in the computation of the propagator, but also in that 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-η limits can be resummed at all orders in PT, giving a well-behaved propagator.

The propagator connects the initial correlators with the cross-correlations between the final and the initial field configurations,

38

where the last term at the right-hand side comes from the initial non-Gaussianity of the matter and halo fields. At leading order, it is given by (see Appendix C)

39

where B_abc is the initial bispectrum at η = 0 (z=z_in) (see also Bernardeau, Crocce & Sefusatti 2010). In Section 6.1, we will use equation (38) to assess the validity of different approximation schemes for the propagator. In the large-k limit, G decays as

40

with

41

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

Although result (40) was obtained for the 2 × 2 propagator of the matter density–velocity system, it holds, taking into account proper modifications, also when haloes 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 one-loop result (equation 33) at low k and the Gaussian decay (equation 40) (with a modified prefactor) 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 2006a; 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 equation (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 Q_abcd = 0, so that the equations for P_ab and B_abc form a closed system:

42

which when 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. P_ab(k; η_i) and B_abc(k, −q, q−k; η_i).

From equation (32), we can observe that Ω₁₃ and Ω₂₃ 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 RESULTS

6.1 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 equation (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 do not have to deal with the shot-noise problem. We extract the cross-spectra from the simulation and compare them against those obtained using both the linear theory propagators P_PL≡g₃₁P⁰₁₁+g₃₂P⁰₂₁+g₃₁P⁰₃₁ and resummed propagators P_PR≡G₃₁P⁰₁₁+G₃₂P⁰₂₁+G₃₁P⁰₃₁+ΔP^NG₃₁; the result is shown in Fig. 3. We note that the linear model severely overpredicts the two-time cross-spectrum for k > 0.05 h Mpc⁻¹. 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 ΔP^NG₃₁ the effect of the initial non-Gaussianity of the halo field via its initial bispectra, computed as in equation (17). It turns out that the components giving a non-vanishing effect are of the B₁₁₃ type. Their contribution is suppressed by a D₊(z_in)/D₊(z = 0) factor with respect to that of a 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.

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 linear theory P_PL (blue dashed line) and resummed result P_PR (green solid line). The red dotted line shows the effect of neglecting the non-Gaussian term, i.e. P_PR−ΔP^NG₃₁. The vertical black dotted lines represent the limit up to which we expect our fluid approximation for the haloes to work.

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 Fig. 4 we show a comparison between the halo–matter cross-spectra extracted from the N-body simulation and 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 autospectrum (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 one-loop correction nor the TRG is 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 Fig. 5, where we plot the ratios between the spectra from the simulation and the theoretical results. The TRG gives a cross-spectrum within 5 per cent accuracy at least up to k = 0.1 h Mpc⁻¹ (barring bin 4), while linear theory does so up to k = 0.05 h Mpc⁻¹.

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 one-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.

Figure 5

The ratio between the spectrum from simulation and linear theory (blue squares), TRG without bispectra (red triangles) and 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.

As previously pointed out, however, the halo-density field is not Gaussian at z_in = 50. We use expressions (17) to account for initial bispectra in the TRG method (they are not present in the linear theory and at one-loop level). Fig. 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 P₁₃ is suppressed. To understand why, we can use PT; first, let us investigate the case of P₁₁. In the one-loop computation, the initial P₁₁ contributes to the final one with a ‘weight’ of [D₊(z = 0)/D₊(z_in)]² = e^2η, as it also happens in linear theory. The initial matter bispectrum B₁₁₁, instead, has a weight of e ^η, so it is suppressed by a factor of e ^−η.3 If we now consider the haloes, we have shown in equation (28) that there is a new decaying mode, responsible for the linear debiasing. This new mode, which involves only the halo field, carries an extra e ^−η suppression factor. We can now rank the contributions to P₁₃ according to their relevance:

• P₁₁;

• P₁₃, B₁₁₁;

• P₃₃, B₁₁₃;

• B₁₃₃;

• B₃₃₃.

Each item is suppressed by a factor of e ^−η with respect to the previous one. We can see that only B₁₁₁ has some relevance, while the other terms are highly suppressed. Even if the reasoning was based on PT, it is valid also for TRG, at a qualitative level. Incidentally, the fact that P₁₁ 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 Q₁₁₁₁ 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 further suppressed with respect to Q₁₁₁₁ 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 < z_in; 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 onwards, 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 Fig. 6 that the TRG is able to give a better prediction than full linear theory in equation (31), even neglecting the velocity bias that the linear theory accounts for.

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.

We can now look at the model predictions for the halo bias, defined as the ratio P₁₃/P₁₁. This quantity is plotted as a function of wavenumber in Fig. 7. While the linear-theory bias always increases with scale, irrespective 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.

Figure 7

The effective bias b≡P_mh/P_m 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 the TRG. The vertical black dotted lines represent the limit up to which we expect our fluid approximation for the haloes to work.

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 evolution of the bias as well. In Figs 8 and 9 we plot, respectively,

43

Again, our model is able to match accurately the trend of the simulation and improve upon linear theory, excluding bin 4. A key feature of the linear solution in equation (26) is the debiasing between halo and matter distributions with time. It is worth noting from Fig. 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.

Figure 8

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

Figure 9

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

7 CONCLUSIONS

We have presented a novel approach to modelling 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 protohaloes). Since the number of protohaloes is conserved over time, for sufficiently large scales (k < 0.2 h Mpc⁻¹), 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 one-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 the non-linear evolution using a semi-analytical procedure, namely an extension of the TRG.

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 as to reproduce the distribution of protohaloes 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.

• 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⁻¹. 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⁻¹ for a broad range of halo masses. This does not hold for very massive haloes (M > 10¹⁴h⁻¹ M_⊙), which 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 one-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.

1

Note that, adopting a local-bias model, δ_h(q) =b₁δ_m(q) +b₂δ²_m(q), provides a worse fit to simulation results. In this case, the model P₁₃ lacks power for all but the smallest wavenumbers.

2

They should not be confused with the P₁₃ and P₂₂ of our notation!

3

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

AE was supported through the SFB-Transregio 33 ‘The Dark Universe’ by the Deutsche Forschungsgemeinschaft (DFG).

AE and SK acknowledge the 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 2010 August.

Appendices

APPENDIX A: COMPLETE ANALYTIC EXPRESSIONS OF THE ONE-LOOP PROPAGATOR AND POWER SPECTRUM

Linear power spectra:

(A1)

A1 One-loop 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

(A2)

where

(A3)

while

(A4)

Looking at equation (33), we can consistently split the one-loop correction to the propagator as Δg_ab=Δ⁽¹⁾g_ab+Δ⁽²⁾g_ab, and we denote them as δ⁽ⁱ⁾g_ab, so that Δ⁽ⁱ⁾g_ab=∫dq δ⁽ⁱ⁾g_ab, for i = 1, 2:

(A5)

A2 One-loop power spectrum

From equation (37), one can see that the one-loop corrections require integrals over ds₁, ds₂ and d³q=q²sin θ dq dθ dϕ; the integration over q and over μ≡ cos θ cannot be performed analytically. The following expressions are therefore integral kernels, denoted by δ_qμP (the factor of q² is already taken into account). If cos θ=μ=k·q/(kq) and ⁠,

(A6)

(A7)

(A8)

APPENDIX B: THE RESUMMED PROPAGATOR

As it was discussed in Crocce & Scoccimarro (2006b), the leading contribution in k²exp(2η) to the propagator at a generic n-loop order contains a chain of propagators and vertices of the form

(B1)

The c_i indices have to be contracted in all possible pairings, by inserting n linear power spectra, each of the form

(B2)

The n-loop contribution is obtained by multiplying by and by the appropriate combinatoric factor, and then by integrating over

(B3)

Recall equations (A2)–(A4): if all the insertions are of the type, then the resummation goes exactly as in the standard case. Indeed,

(B4)

in the k≫q limit. Besides the explicit form for the vertices (see Section 1), we have used the composition property of the propagators,

(B5)

and have defined u = (1, 1, 1). Since each vertex, contracted by a vector, becomes proportional to a delta function in its first and third indices, the chain of propagators composes up to a single propagator g_ab(η) and the time integral can be easily performed. The momentum integrals factorize into n integrals of the type

(B6)

where σ² has been defined in equation (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 type is

(B7)

which resums to ⁠.

As for the insertions including the contribution to the linear power spectrum, the important point to realize is that, due to equation (A4) and the structure of the linear propagator and the vertices, a chain-like equation (B1) can be contracted by at most one power spectrum, the remaining ones being of the type ⁠. Moreover, these insertions only contribute to G₃₁ and G₃₂. Therefore, at nth order, we have only two types of contributions: those with all insertions, giving equation (B7), and those with one insertions and n− 1 ones, which can also be resummed in the large-k limit.

As a consequence, the complete resummed propagator in the large momentum limit is

(B8)

with f_b = (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 one-loop result at low k. This can be done, for instance for G₁₁, starting from its one-loop expression

(B9)

where the above approximation is exact in the large-η limit. Since

(B10)

the required interpolation is given by

(B11)

Proceeding analogously for the other components, we have

(B12)

where we have used the identities

(B13)

APPENDIX C: TAKING INTO ACCOUNT THE INITIAL HALO NON-GAUSSIANITY

While at z_in 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 Section 6.2. In order to obtain the effect of an initial non-Gaussianity on the cross-correlator at different times considered in equation (38), it suffices to compute the O(γ) correction to the linear evolution of the field, by inserting equation (22) in (5), to get

(C1)

The cross-correlator

(C2)

then includes the non-Gaussian expectation value

(C3)

and therefore gives equation (39).