© ESO, 2011

1. Introduction

In the standard cosmological scenario, the large-scale structures we observe in the recent Universe (galaxies, clusters, filaments, voids, etc.) have formed through the amplification by gravitational instability of small primordial perturbations generated in the primordial Universe by quantum fluctuations (Peebles 1980). Moreover, at the beginning of matter domination (z ~ 3000), the power increases on small scales for cold dark matter (CDM) models (Peebles 1982), which leads to a “hierarchical scenario” where smaller scales become nonlinear first. Then, increasingly large and massive structures form in the course of time, as small structures merge and larger scales turn nonlinear. Therefore, on large scales or at early times, it is possible to use linear theory or, more generally, perturbation theory, while on small scales or at late times one must use numerical simulations or phenomenological models. One very interesting quantity for characterizing the density fields built by these processes is the density power spectrum P(k), which is the Fourier transform of the density two-point correlation function. On very large linear scales this is actually sufficient for fully determining the statistical properties of the matter distribution if the initial conditions are (almost) Gaussian, as in standard scenarios of single-field inflation. On smaller nonlinear scales, higher order correlations are generated by the nonlinearities of the dynamics, but P(k) remains a standard tool for comparing observations with theoretical predictions (since higher order statistics are increasingly noisy or difficult to predict).

When one focuses on large scales, several observational probes, such as weak lensing surveys (Massey et al. 2007; Munshi et al. 2008), or galaxy redshift surveys provide rich cosmological information. For example, the signature of baryon acoustic oscillations (BAO) (Eisenstein et al. 1998, 2005), redshift-space distortions (Hamilton 1992; Percival & White 2009), non-Gaussianities in the primordial perturbations (Liguori et al. 2010; Desjacques & Seljak 2010), and the mass of neutrinos (Swanson et al. 2010; Saito et al. 2010) are sensitive to the clustering on linear to weakly nonlinear scales, where the departures from the linear regime are already noticeable but still modest. To match the accuracy of future observations a theoretical accuracy on the order of 1% is required, which is beyond the reach of simple phenomenological models or fits to numerical simulations (if one wishes to obtain robust predictions for a wide range of initial conditions and cosmological parameters). However, this regime should be within the range of validity of perturbation theories, which offer the advantage of systematic and reliable predictions, without fitting parameters. This has led to a renewed interest in perturbative approaches in recent years, as it may be possible to improve over the standard perturbation theory (Goroff et al. 1986; Bernardeau et al. 2002) by using resummation schemes that allow partial resummations of higher order terms (Crocce & Scoccimarro 2006b,a; Valageas 2007a; Matarrese & Pietroni 2007; Taruya & Hiramatsu 2008; Valageas 2008; Pietroni 2008; Matsubara 2008). In particular, this somewhat improves the accuracy of the predicted matter power spectrum on the large scales associated with BAO, as compared with standard perturbation theory truncated at the same order (Crocce & Scoccimarro 2008; Carlson et al. 2009; Taruya et al. 2009).

On the other hand, many observational probes, such as weak lensing surveys on smaller angular scales and galaxy surveys, are sensitive to highly nonlinear scales. There, systematic analytical approaches that have been proposed so far do not apply, in particular because of the importance of shell crossing effects that are not adequately taken into account (an exception is the formalism developed in Valageas 2004, which does not make the single-stream approximation since it deals with the phase-space distribution f(x,v;t), but leads to heavy computations). Then, one must resort to numerical simulations or to phenomenological models. One such model is the Lagrangian mapping presented in Hamilton et al. (1991; see Peacock & Dodds 1996, in Fourier space), that writes the nonlinear two-point correlation or power spectrum in terms of its linear counterpart on a different scale. A second model is the “halo model” (Scherrer & Bertschinger 1991; Cooray & Sheth 2002), where the matter distribution is described as a collection of halos. This yields two contributions for the power spectrum or two-point correlation function, a “1-halo” term associated with particle pairs that belong to the same halo, and a “2-halo” term associated with particles that belong to two different halos. Then, the 1-halo term, which dominates on small scales, is sensitive to the density profile and mass distribution of the halos, whereas the 2-halo term, which dominates on large scales, is sensitive to the correlation between halos. In practice, one often replaces the 2-halo term by the linear two-point correlation, or power spectrum, to reproduce linear theory on large scales. Previous works (Cooray & Sheth 2002; Smith et al. 2003, 2007; Giocoli et al. 2010) have shown that halo models provide a good match to numerical simulations, especially at high k, whereas Lagrangian mappings based on the stable clustering ansatz (Peebles 1980) do not fare as well.

It is clear that it is of great interest to develop unified models, which combine the benefits of systematic perturbation theories on large scales with the reasonable match to simulations on small scales provided by halo models. The goal of the present paper is precisely to investigate the building of such a model, that describes all scales, from the linear to the highly nonlinear regime. We can already note here two points that must be addressed to reach this goal: i) the increasingly large growth at high k of higher order terms obtained within standard perturbation theory, and ii) the nonzero limit for k → 0 of the 1-halo term (as written in previous works), which eventually becomes non-negligible (and even dominant) on large scales, since the linear power spectrum typically scales as P_L(k) ~ k at low k for CDM scenarios. Thus, the standard implementations of both frameworks are badly behaved beyond their regime of validity, where they should actually become negligible.

In this work we address these issues, and we perform a detailed comparison between different possible implementations of such a unified model, as well as with numerical simulations, as follows. We first consider the computation of the density power spectrum from a Lagrangian point of view in Sect. 2. This provides a decomposition over 2-halo and 1-halo terms that are slightly different from the versions found in previous works, and we emphasize their relationship with “perturbative” and “non-perturbative” terms. In particular, we explain why the 1-halo contribution contains a counterterm, missed in previous studies, that ensures that it becomes truly negligible on large scales, solving the second point ii) above. Next, we present a simple implementation of this model in Sect. 3, and we explain how the point i) above is solved by using resummation schemes instead of the standard perturbation theory. We describe in Sect. 4 the N-body simulations that we use to evaluate the accuracy of our models, and we perform detailed comparisons for the nonlinear density power spectrum in Sect. 5, from k ~ 0.02   h   Mpc^-1 up to k ~ 100   h   Mpc^-1, and for redshifts z = 0.35 up to z = 3. Then, we investigate in Sect. 6 the dependence of our results on various ingredients of the model, such as the properties of virialized halos or the choice of the perturbative resummation scheme. We also evaluate the importance of the 1-halo counterterm on large scales. Next, we consider the real-space two-point correlation function in Sect. 7. Finally, we estimate in Sect. 8 the accuracy of such unified models, for both the power spectrum and the two-point correlation, and we conclude in Sect. 9.

In this paper we neglect the impact of baryonic physics on the matter power spectrum, or more precisely we do not explicitly address this issue. Thus, we assume that it can be incorporated in an effective manner through the choice of the halo density profile that enters the model.

2. Decomposition of the density power spectrum from a Lagrangian point of view

We show in this section how the density power spectrum can be split into “perturbative” and “non-perturbative” terms, and how they are related to the 2-halo and 1-halo terms, from a Lagrangian point of view.

2.1. “Perturbative” and “non-perturbative” terms

Let us recall that in a Lagrangian framework one considers the trajectories x(q,t) of all particles, of initial Lagrangian coordinates q and Eulerian coordinates x at time t. In particular, at any given time t, this defines a mapping, q → x, from Lagrangian to Eulerian space, which fully determines the Eulerian density field ρ(x) through the conservation of matter, $\mathit{\rho }\mathrm{\left(}{x}\mathrm{\right)}\hspace{0.17em}\mathrm{d}{x}\mathrm{=}\overline{\mathit{\rho }}\hspace{0.17em}\mathrm{d}{q}\mathit{,}$(1)where $\overline{\mathit{\rho }}$ is the mean comoving matter density of the Universe and we work in comoving coordinates. Then, defining the density contrast as $\mathit{\delta }\mathrm{\left(}{x}\mathit{,t}\mathrm{\right)}\mathrm{=}\frac{\mathit{\rho }\mathrm{\left(}{x}\mathit{,t}\mathrm{\right)}\mathrm{-}\overline{\mathit{\rho }}}{\overline{\mathit{\rho }}}\mathit{,}$(2)and its Fourier transform as $\begin{array}{}\mathrm{˜}\\ \mathit{\delta }\end{array}\mathrm{\left(}{k}\mathrm{\right)}\mathrm{=}\mathrm{\int }\frac{\mathrm{d}{x}}{\mathrm{(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathrm{3}}}\hspace{0.17em}{\mathrm{e}}^{\mathrm{-}\mathrm{i}{k}\mathrm{·}{x}}\hspace{0.17em}\mathit{\delta }\mathrm{\left(}{x}\mathrm{\right)}\mathit{,}$(3)we obtain from Eq. (1) $\begin{array}{}\mathrm{˜}\\ \mathit{\delta }\end{array}\mathrm{\left(}{k}\mathrm{\right)}\mathrm{=}\mathrm{\int }\frac{\mathrm{d}{q}}{\mathrm{(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathrm{3}}}\hspace{0.17em}\mathrm{(}{\mathrm{e}}^{\mathrm{-}\mathrm{i}{k}\mathrm{·}{x}\mathrm{\left(}{q}\mathrm{\right)}}\mathrm{-}{\mathrm{e}}^{\mathrm{-}\mathrm{i}{k}\mathrm{·}{q}}\mathrm{)}\mathit{.}$(4)Next, defining the density power spectrum as $\mathrm{⟨}\begin{array}{}\mathrm{˜}\\ \mathit{\delta }\end{array}\left({k}\mathrm{1}\right)\begin{array}{}\mathrm{˜}\\ \mathit{\delta }\end{array}{\left({k}\mathrm{2}\right)}^{\mathrm{⟩}}\mathrm{=}{\mathit{\delta }}_{\mathrm{D}}\left({k}\mathrm{1}\mathrm{+}{k}\mathrm{2}\right)\mathit{P}\mathrm{\left(}{\mathit{k}}_{\mathrm{1}}\mathrm{\right)}\mathit{,}$(5)we obtain from Eq. (4), using statistical homogeneity, $\mathit{P}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathrm{=}\mathrm{\int }\frac{\mathrm{d}{q}}{\mathrm{(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathrm{3}}}\hspace{0.17em}\mathrm{⟨}{\mathrm{e}}^{\mathrm{i}{k}\mathrm{·}\mathrm{\Delta }{x}}\mathrm{-}{\mathrm{e}}^{\mathrm{i}{k}\mathrm{·}{q}}\mathrm{⟩}\mathit{,}$(6)where we introduced the Eulerian-space separation Δx, $\mathrm{\Delta }{x}\mathrm{=}{x}\mathrm{\left(}{q}\mathrm{\right)}\mathrm{-}{x}\mathrm{\left(}\mathrm{0}\mathrm{\right)}\mathit{.}$(7)The expression (6) is fully general since it is a simple consequence of the matter conservation Eq. (1), and of statistical homogeneity. In particular, it holds for any dynamics, such as the one associated with the Zeldovich approximation (Zeldovich 1970). There, the mapping x(q) is given by the linear displacement field Ψ_L, as x(q) = q + Ψ_L(q), so that for Gaussian initial conditions (i.e. Ψ_L(q) is Gaussian) one can easily compute the average (6) and derive an explicit analytic expression for the density power spectrum (Schneider & Bartelmann 1995; Taylor & Hamilton 1996; Valageas 2010a).

The second term in Eq. (6), e^ik·q, merely gives a Dirac factor, δ_D(k), and can often be discarded in exact or systematic computations (e.g., perturbative expansions) if one consider k ≠ 0. Within perturbative expansions of P(k) over powers of the linear power spectrum, P_L(k), it simply cancels out the zeroth-order term so that P(k) = P_L(k) at lowest order. However, as we shall see below, it will prove important in the following analysis to keep track of this factor. Indeed, it will be split over two parts, along with the first term e^ik·Δx, and it is clear that one should either discard it or keep it in a consistent manner in both parts, contrary to what is implicitly done in usual versions of the halo model.

For the three-dimensional gravitational dynamics it is not possible to derive an explicit analytical expression for the average (6). Therefore, as in the Eulerian framework, one must resort to perturbative expansions (Bouchet et al. 1995; Matsubara 2008) or phenomenological models. In this paper, our goal is precisely to use Eq. (6) as a starting point to decompose the density power spectrum, P(k), over two parts, associated with “perturbative” and “non-perturbative” contributions.

Here and in the following, we call “perturbative” those contributions that can be obtained from perturbative expansions over powers of the linear power spectrum P_L(k) within the fluid approximation. In practice, these may be derived from the equations of motion, written in terms of the density and velocity fields in the usual hydrodynamical description of the system, by writing these nonlinear Eulerian fields as perturbative series over powers of the linear growing mode (Goroff et al. 1986; Bernardeau et al. 2002). Partial resummations of this standard perturbative series have been recently introduced (Crocce & Scoccimarro 2006b,a; Valageas 2007a; Matarrese & Pietroni 2007; Taruya & Hiramatsu 2008; Valageas 2008; Pietroni 2008) but remain within this hydrodynamical regime. Of course, it is possible to develop perturbative expansions, and associated resummations, in Lagrangian space (Bouchet et al. 1995; Matsubara 2008), by looking for an expansion of the displacement field Ψ(q) over powers of the linear displacement field. Taken at face value, these schemes imply some shell crossing. For instance, at lowest order one recovers the Zeldovich dynamics for the displacement field. However, this does not mean that they take into account in a systematic fashion the physics beyond shell crossing, which remains beyond the reach of current approaches of this kind. Technically, the gravitational potential is obtained from the Poisson equation, where the density contrast is obtained from the conservation of matter (1) as δ(x) =  |det(∂x/∂q)| ^-1 − 1, and one expands the Jacobian over powers of the displacement Ψ. In doing so, one discards the absolute value and misses the changes of sign that occur after shell crossing. On the other hand, these schemes give the same final expansion over powers of P_L for the nonlinear power spectrum as the Eulerian schemes, up to the order of truncation of the computation. Therefore, we consider both perturbative approaches as “perturbative” schemes, in the sense of expansions over powers of P_L that neglect shell crossing. By contrast, we call “non-perturbative” those contributions that cannot be obtained through these approaches, either because they are beyond the reach of perturbative expansions (such as factors of the form e^−1/P_L that cannot be uncovered by Taylor series over P_L) or because they arise from the behavior of the system beyond shell crossing (and require going beyond the single-stream approximation).

In Valageas (2010a) such a decomposition between perturbative and non-perturbative contributions to the density power spectrum was obtained within the framework of the Zeldovich dynamics. More precisely, a “sticky model” that only differs from the Zeldovich dynamics beyond shell crossing was introduced. Then, whereas the Zeldovich density power spectrum is exactly given by a resummation of perturbative terms (i.e., the non-perturbative contribution is zero), the “sticky model” power spectrum contains an additional non-perturbative contribution, while its perturbative part is identical to the one of the Zeldovich dynamics. Inspired by this simpler example, we perform a similar decomposition of Eq. (6) over perturbative and non-perturbative contributions. However, we no longer have a simple criterion to distinguish particle pairs that have undergone shell crossing, because for the gravitational dynamics even at the perturbative level (at third order) the displacement field develops rotational terms (Buchert 1994; Bernardeau & Valageas 2008). Therefore, we use a phenomenological halo model, in the sense that we assume that the matter distribution can be described as a collection of spherical halos defined by a given nonlinear density contrast δ. As usual, we can write the halo mass function as $\mathit{n}\mathrm{\left(}\mathit{M}\mathrm{\right)}\mathrm{d}\mathit{M}\mathrm{=}\frac{\overline{\mathit{\rho }}}{\mathit{M}}\mathit{f}\mathrm{\left(}\mathit{\nu }\mathrm{\right)}\frac{\mathrm{d}\mathit{\nu }}{\mathit{\nu }}\mathit{,}\mathrm{with}\mathit{\nu }\mathrm{=}\frac{{\mathit{\delta }}_{\mathrm{L}}}{\mathit{\sigma }\mathrm{\left(}\mathit{M}\mathrm{\right)}}\mathit{,}$(8)in terms of a scaling function f(ν) and of the reduced variable ν. Here, σ(M) is the rms linear density contrast at scale M, or Lagrangian radius q_M, with $\mathit{\sigma }\mathrm{\left(}\mathit{M}\mathrm{\right)}\mathrm{=}\mathit{\sigma }\mathrm{\left(}{\mathit{q}}_{\mathit{M}}\mathrm{\right)}\mathrm{with}\mathit{M}\mathrm{=}\overline{\mathit{\rho }}\frac{\mathrm{4}\mathit{\pi }}{\mathrm{3}}{\mathit{q}}_{\mathit{M}}^{\mathrm{3}}\mathit{,}$(9)and ${\mathit{\sigma }}^{\mathrm{2}}\mathrm{\left(}\mathit{q}\mathrm{\right)}\mathrm{=}\mathrm{4}\mathit{\pi }{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{\infty }}\mathrm{d}\mathit{k}\hspace{0.17em}{\mathit{k}}^{\mathrm{2}}{\mathit{P}}_{\mathrm{L}}\mathrm{\left(}\mathit{k}\mathrm{\right)}\begin{array}{}\mathrm{˜}\\ \mathit{W}\end{array}\mathrm{(}\mathit{kq}{\mathrm{)}}^{\mathrm{2}}\mathit{,}$(10)where $\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{W}\end{array}\mathrm{\left(}\mathit{kq}\mathrm{\right)}$ is the Fourier transform of the top-hat of radius q, defined as $\begin{array}{}\mathrm{˜}\\ \mathit{W}\end{array}\mathrm{\left(}\mathit{kq}\mathrm{\right)}\mathrm{=}{\mathrm{\int }}_{\mathit{V}}\frac{\mathrm{d}{q}}{\mathit{V}}\hspace{0.17em}{\mathrm{e}}^{\mathrm{i}{k}\mathrm{·}{q}}\mathrm{=}\mathrm{3}\hspace{0.17em}\frac{\sin \mathrm{\left(}\mathit{kq}\mathrm{\right)}\mathrm{-}\mathit{kq}\cos \mathrm{\left(}\mathit{kq}\mathrm{\right)}}{\mathrm{(}\mathit{kq}{\mathrm{)}}^{\mathrm{3}}}\mathrm{·}$(11)In the second Eq. (8) the linear density contrast δ_L is related to the nonlinear density threshold δ that defines the halos through the spherical collapse dynamics (Valageas 2009a), $\mathit{\delta }\mathrm{=}\mathrm{ℱ}\mathrm{\left(}{\mathit{\delta }}_{\mathrm{L}}\mathrm{\right)}\mathit{,}$(12)and for Gaussian initial conditions the large-mass tail shows the exponential falloff $\mathit{\nu }\mathrm{\to }\mathrm{\infty :}\ln \mathrm{\left[}\mathit{f}\mathrm{\right(}\mathit{\nu }\mathrm{\left)}\mathrm{\right]}\mathrm{~}\mathrm{-}\frac{{\mathit{\nu }}^{\mathrm{2}}}{\mathrm{2}}\mathit{,}$(13)whence $\mathit{M}\mathrm{\to }\mathrm{\infty :}\mathit{n}\mathrm{\left(}\mathit{M}\mathrm{\right)}\mathrm{~}{\mathrm{e}}^{\mathrm{-}{\mathit{\delta }}_{\mathrm{L}}^{\mathrm{2}}\mathit{/}\left(\mathrm{2}{\mathit{\sigma }}^{\mathrm{2}}\mathrm{\left(}\mathit{M}\mathrm{\right)}\right)}\mathit{.}$(14)This corresponds to the Press-Schechter falloff (Press & Schechter 1974), but with a somewhat lower threshold δ_L given by Eq. (12) (the usual Press-Schechter threshold actually corresponds to δ_L = ℱ^-1(∞), associated with full collapse to a point). However, the relation (12) only holds for δ < δ_vir, which typically gives δ < 200, as larger density contrasts are associated with inner shells where shell crossing plays a key role and modifies Eq. (12), associated with spherical dynamics at constant mass (Valageas 2009a). On the other hand, the halo mass function satisfies the normalization ${\mathrm{\int }}_{\mathrm{0}}^{\mathrm{\infty }}\frac{\mathrm{d}\mathit{\nu }}{\mathit{\nu }}\hspace{0.17em}\mathit{f}\mathrm{\left(}\mathit{\nu }\mathrm{\right)}\mathrm{=}\mathrm{1}\mathit{,}$(15)which ensures that all the mass is contained within such halos (for linear power spectra such that σ(q) grows to infinity on small scales). This also ensures that there is no overcounting (as would be the case if one used a mass function with a normalization greater than unity).

Then, in Lagrangian space, the probability dF for one particle q₁ to belong to a halo of mass in the range  [M,M + dM]  reads as $\mathrm{d}\mathit{F}\mathrm{=}\mathit{f}\mathrm{\left(}\mathit{\nu }\mathrm{\right)}\frac{\mathrm{d}\mathit{\nu }}{\mathit{\nu }}\mathrm{·}$(16)This is also the fraction of matter enclosed within such halos. Next, making the approximation that each halo comes from an initial spherical region in Lagrangian space, the probability for a second particle q₂, at distance q =  |q₂ − q₁| , to belong to the same halo reads as ${\mathit{F}}_{\mathit{M}}\mathrm{\left(}\mathit{q}\mathrm{\right)}\mathrm{=}\frac{{\mathrm{\int }}_{\mathit{V}}\mathrm{d}{q}\mathrm{1}\mathrm{\int }\mathrm{d}{\mathbf{\Omega }}_{\mathit{q}}\hspace{0.17em}\mathit{\theta }\left({q}\mathrm{2}\mathrm{\in }\mathit{V}\right)}{\frac{\mathrm{4}\mathit{\pi }}{\mathrm{3}}{\mathit{q}}_{\mathit{M}}^{\mathrm{3}}\hspace{0.17em}\mathrm{4}\mathit{\pi }}\mathrm{·}$(17)Here we integrate over all positions q₁ within the spherical volume V of radius q_M, and we integrate over all directions Ω_q of the Lagrangian vector q = q₂ − q₁. The top-hat factor θ(q₂ ∈ V) is unity if q₂ is located within the volume V, and zero otherwise. By isotropy the result only depends on the length q, and performing the integrations yields $\mathrm{0}\mathrm{\le }\mathit{q}\mathrm{\le }\mathrm{2}{\mathit{q}}_{\mathit{M}}\mathrm{:}{\mathit{F}}_{\mathit{M}}\mathrm{\left(}\mathit{q}\mathrm{\right)}\mathrm{=}\frac{{\left(\mathrm{2}{\mathit{q}}_{\mathit{M}}\mathrm{-}\mathit{q}\right)}^{\mathrm{2}}\left(\mathrm{4}{\mathit{q}}_{\mathit{M}}\mathrm{+}\mathit{q}\right)}{\mathrm{16}{\mathit{q}}_{\mathit{M}}^{\mathrm{3}}}\mathit{,}$(18)and F_M(q) = 0 for q > 2q_M. Therefore, combining Eqs. (16) and (18), we obtain the probability that a pair of particles of Lagrangian separation q belongs to the same halo of mass M as ${\mathit{q}}_{\mathit{M}}\mathit{>}\frac{\mathit{q}}{\mathrm{2}}\mathrm{:d}{\mathit{F}}_{\mathit{q}}\mathrm{\left(}\mathit{M}\mathrm{\right)}\mathrm{=}\frac{{\left(\mathrm{2}{\mathit{q}}_{\mathit{M}}\mathrm{-}\mathit{q}\right)}^{\mathrm{2}}\left(\mathrm{4}{\mathit{q}}_{\mathit{M}}\mathrm{+}\mathit{q}\right)}{\mathrm{16}{\mathit{q}}_{\mathit{M}}^{\mathrm{3}}}\hspace{0.17em}\mathit{f}\mathrm{\left(}\mathit{\nu }\mathrm{\right)}\frac{\mathrm{d}\mathit{\nu }}{\mathit{\nu }}\mathit{,}$(19)and dF_q(M) = 0 for q_M < q/2. In particular, the probability that the pair  { q₁,q₂ }  belongs to a single halo writes as ${\mathit{F}}_{\mathrm{1}\mathrm{H}}\mathrm{\left(}\mathit{q}\mathrm{\right)}\mathrm{=}{\mathrm{\int }}_{{\mathit{\nu }}_{\mathit{q}\mathit{/}\mathrm{2}}}^{\mathrm{\infty }}\frac{\mathrm{d}\mathit{\nu }}{\mathit{\nu }}\hspace{0.17em}\mathit{f}\mathrm{\left(}\mathit{\nu }\mathrm{\right)}\hspace{0.17em}\frac{{\left(\mathrm{2}{\mathit{q}}_{\mathit{M}}\mathrm{-}\mathit{q}\right)}^{\mathrm{2}}\left(\mathrm{4}{\mathit{q}}_{\mathit{M}}\mathrm{+}\mathit{q}\right)}{\mathrm{16}{\mathit{q}}_{\mathit{M}}^{\mathrm{3}}}\mathit{,}$(20)where ν_q/2 is defined as in Eq. (8), for the Lagrangian radius q_M = q/2. On the other hand, the probability F_2H(q) that the Lagrangian pair does not belong to a single halo (whence the two points belong to two different halos) reads as ${\mathit{F}}_{\mathrm{2}\mathrm{H}}\mathrm{\left(}\mathit{q}\mathrm{\right)}\mathrm{=}\mathrm{1}\mathrm{-}{\mathit{F}}_{\mathrm{1}\mathrm{H}}\mathrm{\left(}\mathit{q}\mathrm{\right)}\mathit{.}$(21)Then, we split the average in Eq. (6) over two terms, P_1H and P_2H, associated with pairs  { 0,q }  that belong to a single halo or to two different halos, $\mathit{P}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathrm{=}{\mathit{P}}_{\mathrm{1}\mathrm{H}}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathrm{+}{\mathit{P}}_{\mathrm{2}\mathrm{H}}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathit{,}$(22)with ${\mathit{P}}_{\mathrm{1}\mathrm{H}}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathrm{=}\mathrm{\int }\frac{\mathrm{d}{q}}{\mathrm{(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathrm{3}}}\hspace{0.17em}{\mathit{F}}_{\mathrm{1}\mathrm{H}}\mathrm{\left(}\mathit{q}\mathrm{\right)}\hspace{0.17em}{\mathrm{⟨}{\mathrm{e}}^{\mathrm{i}{k}\mathrm{·}\mathrm{\Delta }{x}}\mathrm{-}{\mathrm{e}}^{\mathrm{i}{k}\mathrm{·}{q}}\mathrm{⟩}}_{\mathrm{1}\mathrm{H}}$(23)and ${\mathit{P}}_{\mathrm{2}\mathrm{H}}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathrm{=}\mathrm{\int }\frac{\mathrm{d}{q}}{\mathrm{(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathrm{3}}}\hspace{0.17em}{\mathit{F}}_{\mathrm{2}\mathrm{H}}\mathrm{\left(}\mathit{q}\mathrm{\right)}\hspace{0.17em}{\mathrm{⟨}{\mathrm{e}}^{\mathrm{i}{k}\mathrm{·}\mathrm{\Delta }{x}}\mathrm{-}{\mathrm{e}}^{\mathrm{i}{k}\mathrm{·}{q}}\mathrm{⟩}}_{\mathrm{2}\mathrm{H}}\mathit{.}$(24)Here the averages  ⟨ ... ⟩ _1H and  ⟨ ... ⟩ _2H are the conditional averages, knowing that the pair of length q belongs to a single halo or to two halos. The decomposition (22) clearly corresponds to the 1-halo and 2-halo terms of the usual halo model (Cooray & Sheth 2002). Then, to make the connection with the distinction between perturbative and non-perturbative terms, we note that at a perturbative level F_1H is identically zero and F_2H unity, $\begin{array}{ccc}& & {\mathit{F}}_{\mathrm{1}\mathrm{H}}\mathrm{=}\mathrm{0} \mathrm{and} {\mathit{F}}_{\mathrm{2}\mathrm{H}}\mathrm{=}\mathrm{1}\mathit{,}\\ & & \mathrm{atallordersofperturbationtheory}\mathit{.}\end{array}$(25)Indeed, the large-mass tail (14) is actually a rare-event limit that holds both in the large-mass limit, at fixed linear density power spectrum (Valageas 2002b, 2009a), and in the quasi-linear limit at fixed mass, where the amplitude of the linear density power spectrum goes to zero (Valageas 2002a, 2009a). It is this second regime which corresponds to usual perturbation theories, where as recalled above we look for expansions over powers of the amplitude of the linear power spectrum. Then, because of the exponential decay of the form e ^{− 1/σ2(M)} we can see that the expansion over powers of P_L of F_1H(q) defined in Eq. (20), at fixed q, is identically zero. From Eq. (21) this also yields F_2H = 1 at all orders of perturbation theory. Therefore, we can see from Eqs. (23), (24) that the 1-halo contribution is a fully non-perturbative term, while the 2-halo contribution is (almost) the perturbative term multiplied by the factor F_2H(q).

The factor F_1H being non-perturbative is not mainly related to shell crossing, but simply to the fact that it cannot be recovered by a series expansion over powers of P_L. However, for halos that would be defined by a high density threshold, typically δ > 200, the exponential falloff (14) is modified by shell crossing (i.e. the factor δ_L is no longer given by Eq. (12), see Valageas 2009a) so that it would also be non-perturbative in this sense. On the other hand, even for lower threshold δ, the exact form of the factor F_1H(q), and in particular the low-mass tail of the halo mass function, depends on the behavior of the system beyond shell crossing.

The decomposition (22), with the Lagrangian-based interpretation (23) − (24), has the advantage to automatically satisfy the conservation of matter. Thus, thanks to Eq. (21) we count all particle pairs once. By contrast, in the Eulerian derivation of the halo model, where we first write the density field ρ(x) as a sum of halo profiles, we would need to pay attention to possible overlaps between halos, which arise when we use a spherical approximation. Thus, the splitting (22) is more easily expressed in this framework, and one can independently focus on the modelization of the averages  ⟨ e^ik·Δx − e^ik·q ⟩ . The latter also offers a closer link to the dynamics, through the mapping q → x. We shall not make much use of this relationship in the following, as we use simple approximations that allow us to recover the usual Eulerian expressions, with the addition of simple prefactors and counterterms, but this may provide a route to more accurate modeling.

2.2. “2-halo” contribution

We first consider the 2-halo contribution (24). Since the conditional average involves the constraint that the pair  { 0,q }  does not belong to the same halo, the mean of e^ik·Δx is “biased” as compared with a mean over all possible pairs. However, in order to simplify the computation of this term, we note that at a perturbative level F_2H = 1, as seen in (25), so that we can replace the average  ⟨ ... ⟩ _2H by the mean over all pairs, as given by perturbation theory, ${\mathit{P}}_{\mathrm{2}\mathrm{H}}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathrm{\simeq }\mathrm{\int }\frac{\mathrm{d}{q}}{\mathrm{(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathrm{3}}}\hspace{0.17em}{\mathit{F}}_{\mathrm{2}\mathrm{H}}\mathrm{\left(}\mathit{q}\mathrm{\right)}\hspace{0.17em}{\mathrm{⟨}{\mathrm{e}}^{\mathrm{i}{k}\mathrm{·}\mathrm{\Delta }{x}}\mathrm{-}{\mathrm{e}}^{\mathrm{i}{k}\mathrm{·}{q}}\mathrm{⟩}}_{\mathrm{pert}}\mathit{,}$(26)where the subscript “pert” denotes quantities obtained from standard perturbation theory (or equivalently its various resummation schemes). Thus Eq. (26) is still exact at all orders of perturbation theory. This expression is best suited for perturbation theories developed within the Lagrangian framework. Unfortunately, as we shall discuss in Sect. 6.5 below, Lagrangian perturbation theories built so far are not as efficient as their Eulerian counterparts, especially when we consider available resummation schemes. Then, in order to make contact with Eulerian perturbation theories we further approximate Eq. (26) as $\begin{array}{ccc}{\mathit{P}}_{\mathrm{2}\mathrm{H}}\mathrm{\left(}\mathit{k}\mathrm{\right)}& \mathrm{\simeq }& {\mathit{F}}_{\mathrm{2}\mathrm{H}}\mathrm{(}\mathrm{1}\mathit{/}\mathit{k}\mathrm{)}\mathrm{\int }\frac{\mathrm{d}{q}}{\mathrm{(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathrm{3}}}\hspace{0.17em}{\mathrm{⟨}{\mathrm{e}}^{\mathrm{i}{k}\mathrm{·}\mathrm{\Delta }{x}}\mathrm{-}{\mathrm{e}}^{\mathrm{i}{k}\mathrm{·}{q}}\mathrm{⟩}}_{\mathrm{pert}}\\ & \mathrm{=}& {\mathit{F}}_{\mathrm{2}\mathrm{H}}\mathrm{(}\mathrm{1}\mathit{/}\mathit{k}\mathrm{)}{\mathit{P}}_{\mathrm{pert}}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathit{,}\end{array}$(27)where we have replaced the q-dependent factor F_2H(q) by its value at a typical scale q ~ 1/k. Again, this is legitimate at a perturbative level, where F_2H = 1, so that Eq. (27) remains exact at all orders of perturbation theory. To obtain the second line we simply used the exact expression (6), which implies the same equality in terms of perturbative expansions.

Let us recall that from (25) the 1-halo contribution is zero at all orders of perturbation theory, so that the full power spectrum P(k) of Eq. (22) automatically agrees with perturbation theory at all orders, whether we use Eq. (26) or (27). Then, within these approximations the only effect of non-perturbative corrections to the 2-halo term is to multiply the perturbative power spectrum P_pert(k) by the prefactor F_2H(1/k).

Within the usual halo model the 2-halo term reads as (Cooray & Sheth 2002) $\begin{array}{ccc}{\mathit{P}}_{\mathrm{2}\mathrm{H}}^{\mathrm{h}\mathit{.}\mathrm{m}\mathit{.}}\mathrm{\left(}\mathit{k}\mathrm{\right)}& \mathrm{=}& \int \frac{\mathrm{d}{\mathit{\nu }}_{\mathrm{1}}\mathrm{d}{\mathit{\nu }}_{\mathrm{2}}}{{\mathit{\nu }}_{\mathrm{1}}{\mathit{\nu }}_{\mathrm{2}}}\mathit{f}\mathrm{(}{\mathit{\nu }}_{\mathrm{1}}\mathrm{)}\mathit{f}\mathrm{(}{\mathit{\nu }}_{\mathrm{2}}\mathrm{)}{\begin{array}{}\mathrm{˜}\\ \mathit{u}\end{array}}_{{\mathit{M}}_{\mathrm{1}}}\mathrm{\left(}\mathit{k}\mathrm{\right)}{\begin{array}{}\mathrm{˜}\\ \mathit{u}\end{array}}_{{\mathit{M}}_{\mathrm{2}}}\mathrm{\left(}\mathit{k}\mathrm{\right)}{\mathit{P}}_{{\mathit{M}}_{\mathrm{1}}{\mathit{M}}_{\mathrm{2}}}\mathrm{\left(}\mathit{k}\mathrm{\right)}\\ & \mathrm{\simeq }& {\mathit{P}}_{\mathrm{L}}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathit{,}\end{array}$(28)where ${\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{u}\end{array}}_{\mathit{M}}\mathrm{\left(}\mathit{k}\mathrm{\right)}$ is the normalized halo density profile, defined in Eq. (31) below, and P_M1M2(k) is the halo power spectrum. The second line (28) is obtained in the low-k limit, so that ${\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{u}\end{array}}_{\mathit{M}}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathrm{\to }\mathrm{1}$ and P_M1M2(k) ≃ b(M₁)b(M₂)P_L(k), with a halo bias b(M) that is normalized to unity. It is also possible to combine perturbation theory and nonlinear halo bias to make the expression above consistent with standard perturbation theory while building a model for the halo power spectrum itself, see Smith et al. (2007).

Of course, in order to describe the weakly nonlinear regime one can as well replace P_L(k) by P_pert(k) in Eq. (28), which gives an expression very similar to Eq. (27). Then, we can see that a first difference between the halo-model expression (28) and Eq. (27) is that we did not need to introduce any halo bias to derive Eq. (27). In fact, although the contribution (24) is associated with a 2-halo term, by making the simple approximation (26) we avoid any need to consider halo biasing, and as explained above this does not spoil the agreement with perturbation theory. The second difference is the prefactor F_2H(1/k) in Eq. (27). As seen from Eq. (6) and the splitting (22), this term is required by self-consistency, to ensure that the two averages F_1H(q)    ⟨ e^ik·Δx ⟩ _1H and F_2H(q)    ⟨ e^ik·Δx ⟩ _2H sum up to  ⟨ e^ik·Δx ⟩ , and in particular that they sum up to unity for k → 0. Within the usual halo model, this factor is implicitly set to unity by taking the large-scale limit in Eq. (28) and ignoring exclusion constraints on the halos. This is also valid within perturbation theory, as seen in (25). However, in order to describe the weakly nonlinear regime, where the 1-halo term is nonzero, it is best to keep the prefactor F_2H(1/k) in Eq. (27), to keep a consistent model and to avoid any overcounting.

2.3. “1-halo” contribution

From Eqs. (20) and (23) the 1-halo contribution to the density power spectrum reads as $\begin{array}{ccc}{\mathit{P}}_{\mathrm{1}\mathrm{H}}\mathrm{\left(}\mathit{k}\mathrm{\right)}& \mathrm{=}& \int \frac{\mathrm{d}{q}}{\mathrm{(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathrm{3}}}{\mathrm{\int }}_{{\mathit{\nu }}_{\mathit{q}\mathit{/}\mathrm{2}}}^{\mathrm{\infty }}\frac{\mathrm{d}\mathit{\nu }}{\mathit{\nu }}\hspace{0.17em}\mathit{f}\mathrm{\left(}\mathit{\nu }\mathrm{\right)}\hspace{0.17em}\frac{{\left(\mathrm{2}{\mathit{q}}_{\mathit{M}}\mathrm{-}\mathit{q}\right)}^{\mathrm{2}}\left(\mathrm{4}{\mathit{q}}_{\mathit{M}}\mathrm{+}\mathit{q}\right)}{\mathrm{16}{\mathit{q}}_{\mathit{M}}^{\mathrm{3}}}\\ & & \mathrm{\times }{\mathrm{⟨}{\mathrm{e}}^{\mathrm{i}{k}\mathrm{·}\mathrm{\Delta }{x}}\mathrm{-}{\mathrm{e}}^{\mathrm{i}{k}\mathrm{·}{q}}\mathrm{⟩}}_{\mathit{M}}\mathit{.}\end{array}$(29)In order to compute the average in Eq. (29), within a halo of mass M, we make the approximation of fully virialized spherical halos. Thus, we describe the halos as spherical objects, truncated at a radius r_M such that the mean density contrast within this radius is the nonlinear threshold δ used to define these objects, and the enclosed mass is equal to M, $\mathit{M}\mathrm{=}\overline{\mathit{\rho }}\frac{\mathrm{4}\mathit{\pi }}{\mathrm{3}}{\mathit{q}}_{\mathit{M}}^{\mathrm{3}}\mathrm{=}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{\delta }\mathrm{)}\overline{\mathit{\rho }}\frac{\mathrm{4}\mathit{\pi }}{\mathrm{3}}{\mathit{r}}_{\mathit{M}}^{\mathrm{3}}\mathit{.}$(30)We introduce as usual the normalized Fourier transform of this halo radial profile, ${\begin{array}{}\mathrm{˜}\\ \mathit{u}\end{array}}_{\mathit{M}}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathrm{=}\frac{\mathrm{\int }\mathrm{d}{x}\hspace{0.17em}{\mathrm{e}}^{\mathrm{-}\mathrm{i}{k}\mathrm{·}{x}}{\mathit{\rho }}_{\mathit{M}}\mathrm{\left(}\mathit{x}\mathrm{\right)}}{\mathrm{\int }\mathrm{d}{x}\hspace{0.17em}{\mathit{\rho }}_{\mathit{M}}\mathrm{\left(}\mathit{x}\mathrm{\right)}}\mathrm{=}\frac{\mathrm{1}}{\mathit{M}}\mathrm{\int }\mathrm{d}{x}\hspace{0.17em}{\mathrm{e}}^{\mathrm{-}\mathrm{i}{k}\mathrm{·}{x}}{\mathit{\rho }}_{\mathit{M}}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathit{,}$(31)where ρ_M(x) is the halo density profile. Next, using the approximation of fully virialized halos, that is, that the two Lagrangian particles “0” and “q” have lost all memory of their initial locations and are independently located at random within the halo, we write $\begin{array}{ccc}{\mathrm{⟨}{\mathrm{e}}^{\mathrm{i}{k}\mathrm{·}\mathrm{\Delta }{x}}\mathrm{⟩}}_{\mathit{M}}& \mathrm{=}& \frac{\mathrm{1}}{{\mathit{M}}^{\mathrm{2}}}\mathrm{\int }\mathrm{d}{x}\mathrm{1}\mathrm{d}{x}\mathrm{2}\hspace{0.17em}{\mathit{\rho }}_{\mathit{M}}\mathrm{(}{x}{\mathrm{1}}^{\mathrm{)}}{\mathit{\rho }}_{\mathit{M}}\mathrm{(}{x}{\mathrm{2}}^{\mathrm{)}}\hspace{0.17em}{\mathrm{e}}^{\mathrm{i}{k}\mathrm{·}\left({x}\mathrm{2}\mathrm{-}{x}\mathrm{1}\right)}\\ & \mathrm{=}& {\begin{array}{}\mathrm{˜}\\ \mathit{u}\end{array}}_{\mathit{M}}\mathrm{(}\mathit{k}{\mathrm{)}}^{\mathrm{2}}\mathit{.}\end{array}$Substituting into Eq. (29) and exchanging the order of integration gives $\begin{array}{ccc}{\mathit{P}}_{\mathrm{1}\mathrm{H}}\mathrm{\left(}\mathit{k}\mathrm{\right)}& \mathrm{=}& \int \begin{array}{}\mathrm{\infty }\\ \mathrm{0}\end{array}\frac{\mathrm{d}\mathit{\nu }}{\mathit{\nu }}\mathit{f}\mathrm{\left(}\mathit{\nu }\mathrm{\right)}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{2}{\mathit{q}}_{\mathit{M}}}\frac{\mathrm{d}{q}}{\mathrm{(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathrm{3}}}\hspace{0.17em}\frac{{\left(\mathrm{2}{\mathit{q}}_{\mathit{M}}\mathrm{-}\mathit{q}\right)}^{\mathrm{2}}\left(\mathrm{4}{\mathit{q}}_{\mathit{M}}\mathrm{+}\mathit{q}\right)}{\mathrm{16}{\mathit{q}}_{\mathit{M}}^{\mathrm{3}}}\\ & & \mathrm{\times }\mathrm{(}{\begin{array}{}\mathrm{˜}\\ \mathit{u}\end{array}}_{\mathit{M}}\mathrm{(}\mathit{k}{\mathrm{)}}^{\mathrm{2}}\mathrm{-}{\mathrm{e}}^{\mathrm{i}{k}\mathrm{·}{q}}\mathrm{)}\mathit{,}\end{array}$(34)and the integration over q yields ${\mathit{P}}_{\mathrm{1}\mathrm{H}}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathrm{=}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{\infty }}\frac{\mathrm{d}\mathit{\nu }}{\mathit{\nu }}\mathit{f}\mathrm{\left(}\mathit{\nu }\mathrm{\right)}\frac{\mathit{M}}{\overline{\mathit{\rho }}\mathrm{(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathrm{3}}}\mathrm{\left(}{\begin{array}{}\mathrm{˜}\\ \mathit{u}\end{array}}_{\mathit{M}}\mathrm{\right(}\mathit{k}{\mathrm{)}}^{\mathrm{2}}\mathrm{-}\begin{array}{}\mathrm{˜}\\ \mathit{W}\end{array}{\left(\mathit{k}{\mathit{q}}_{\mathit{M}}\right)}^{\mathrm{2}}\mathrm{)}\mathit{.}$(35)Therefore, we recover the usual 1-halo term of the halo model (Cooray & Sheth 2002), with the addition of the new counterterm $\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{W}\end{array}\mathrm{(}\mathit{k}{\mathit{q}}_{\mathit{M}}{\mathrm{)}}^{\mathrm{2}}$, where $\begin{array}{}{\mathrm{˜}}_{}\\ \mathit{W}\end{array}$ was defined in Eq. (11).

This counterterm originates from the second term in Eq. (6), which actually sums up to a Dirac factor δ_D(k). However, since we eventually split the nonlinear power spectrum as in Eq. (22) and we use different approximations for the 2-halo and 1-halo contributions, it is best to keep track of this factor in both contributions. In particular, as pointed out in Sect. 2.1, one should avoid keeping it in one contribution and disregarding it in the other one. Thus, it has been explicitly used in Eq. (27) to cancel out the zeroth-order term F_2Hδ_D(k), so that one should keep it in the 1-halo term (23). Within the usual halo model this term is usually missed because the splitting between the 2-halo and 1-halo terms is not treated so carefully. In particular, as noticed in Sect. 2.2, the usual halo model implicitly takes F_2H = 1 (which is valid in perturbation theory) while taking F_1H ≠ 0 by adding the 1-halo term, which is not fully consistent.

In fact, the counterterm of Eq. (35) can be recovered without going through the previous steps, associated with a Lagrangian point of view, by the following simple argument. Let us consider a uniform universe, where there are no density fluctuations and the density is everywhere equal to $\overline{\mathit{\rho }}$. Then, within the spirit of the usual (Eulerian) halo model, and making the same approximations (e.g., neglecting geometrical constraints, associated with exclusion constraints, and departures from spherical symmetry), we can consider that all the matter is within halos of constant density $\overline{\mathit{\rho }}$ with an arbitrary distribution of halo radii. That is, we can split this uniform system over an arbitrary set of cells, which we can call halos. Then, within such a halo model the usual 1-halo term for the full density correlation, ${\mathit{P}}_{\mathrm{1}\mathrm{H}}^{\mathit{\rho \rho }}\mathrm{\left(}\mathit{k}\mathrm{\right)}$ (defined as $\mathrm{⟨}\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{\rho }\end{array}\mathrm{\left(}{k}{}_{\mathrm{1}}\mathrm{\right)}\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{\rho }\end{array}\mathrm{\left(}{k}{}_{\mathrm{2}}\mathrm{\right)}\mathrm{⟩}\mathrm{=}{\mathit{\delta }}_{\mathrm{D}}\mathrm{(}{k}{}_{\mathrm{1}}\mathrm{+}{k}{}_{\mathrm{2}}\mathrm{)}{\mathit{P}}^{\mathit{\rho \rho }}\mathrm{\left(}{\mathit{k}}_{\mathrm{1}}\mathrm{\right)}$, as in Eq. (5)), reads as the counterterm of Eq. (35), with a positive sign and multiplied by ${\overline{\mathit{\rho }}}^{\mathrm{2}}$. Indeed, the density profile of such halos is simply the top-hat of radius q_M, since δ = 0 whence r_M = q_M, so that the normalized Fourier transform of the halo radial profile defined in Eq. (31) is equal to the Fourier transform (11), ${\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{u}\end{array}}_{\mathit{M}}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathrm{=}\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{W}\end{array}\mathrm{\left(}\mathit{k}{\mathit{q}}_{\mathit{M}}\mathrm{\right)}$. This means that within a halo model, for any halo mass function, the 1-halo contribution to the power spectrum associated with an uniform-density medium reads as the counterterm of Eq. (35). Since the density-contrast power spectrum can also be defined as ${\overline{\mathit{\rho }}}^{\mathrm{2}}\mathit{P}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathrm{=}{\mathit{P}}^{\mathit{\rho \rho }}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathrm{-}{\mathit{P}}^{\overline{\mathit{\rho }}\hspace{0.17em}\overline{\mathit{\rho }}}\mathrm{\left(}\mathit{k}\mathrm{\right)}$ (i.e. we consider the connected two-point correlation), we must subtract this uniform-density term as in Eq. (35), as the first term ${\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{u}\end{array}}_{\mathit{M}}\mathrm{(}\mathit{k}{\mathrm{)}}^{\mathrm{2}}$ corresponds to the total density power spectrum, P^ρρ(k).

In agreement with this result, we can note by an integration over k of the second term of Eq. (34) that the second term of Eq. (35) is indeed a well-normalized approximation to δ_D(k), where the Dirac peak is broadened over a size  ~1/q_M associated with the typical halo scale (thus, as could be expected, this becomes an increasingly good approximation to δ_D(k) as the halo size goes to infinity and the associated power is repelled to k → 0). Clearly, the fact that this does not yield exactly a Dirac factor, whence that it does not vanish for nonzero k, is due to the fact that it is only a partial contribution to the full power spectrum, since by construction within such halo models there is always a 2-halo term, and it is only the sum of both contributions that contains such a Dirac factor. (This also agrees with the fact that we recover a Dirac factor in the limit q_M → ∞, since at fixed scale 1/k the 2-halo term clearly disappears in this limit.)

In practice, the counterterm of Eq. (35) is small on most scales of interest, especially in the nonlinear regime. However, it plays an important role on very large scales. Indeed, since ${\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{u}\end{array}}_{\mathit{M}}\mathrm{\left(}\mathrm{0}\mathrm{\right)}\mathrm{=}\mathrm{1}$ we can see that the first term of Eq. (35), associated with the usual halo model, goes to a nonzero constant as k → 0. Since for CDM-like initial power spectra P_L(k) ∝ k^n_s with n_s ≃ 1 as k → 0, this implies that the usual 1-halo term dominates on very large scale, which is not correct. In fact, from the conservation of matter a small-scale redistribution of matter generically yields a P(k) ∝ k² tail at low k, while taking into account momentum conservation one obtains a k⁴ tail (Peebles 1974). To solve this problem the use of compensated halo profiles (i.e. with ${\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{u}\end{array}}_{\mathit{M}}\mathrm{\left(}\mathrm{0}\mathrm{\right)}\mathrm{=}\mathrm{0}$) was investigated in Cooray & Sheth (2002). However, they noticed that this recipe fails because it also spoils the 2-halo term, as seen from the first line in Eq. (28). Here we can see that this problem is automatically solved within our approach by the counterterm $\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{W}\end{array}\mathrm{(}\mathit{k}{\mathit{q}}_{\mathit{M}}{\mathrm{)}}^{\mathrm{2}}$ in Eq. (35). In particular we recover at low k the expected slope P_1H(k) ∝ k², associated with small-scale rearrangements (indeed, the 1-halo term only describes the redistribution within small clumps, and is not sensitive to the long-range correlations between clumps, that is described by the 2-halo term).

Looking at the steps of the derivation of Eq. (35) we can also see why we fail to recover the k⁴ tail due to momentum conservation. Indeed, in the approximation (33), where we have assumed full virialization, we have erased all memory of the initial positions and velocities of the two Lagrangian particles  { 0,q } , which clearly implies that we have disregarded the constraints associated with momentum conservation. Therefore, in order to obtain a k⁴ tail one should improve this approximation, by taking care at some level of the particle velocities. Of course, a simpler remedy would be to modify the counterterm $\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{W}\end{array}\mathrm{(}\mathit{k}{\mathit{q}}_{\mathit{M}}{\mathrm{)}}^{\mathrm{2}}$ in Eq. (35), by using a window ${\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{W}\end{array}}_{\mathrm{2}}\mathrm{(}\mathit{k}{\mathit{q}}_{\mathit{M}}{\mathrm{)}}^{\mathrm{2}}$ that cancels both the terms k⁰ and k² of ${\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{u}\end{array}}_{\mathit{M}}\mathrm{(}\mathit{k}{\mathrm{)}}^{\mathrm{2}}$. However, since the CDM linear power spectrum roughly behaves as P_L(k) ∝ k at low k the k² tail is sufficient to make the 1-halo term subdominant on large scales, as we shall check in Sect. 5 below. Therefore, we shall keep the simple approximation (33) in the following, and Eq. (35) for the 1-halo term, as it appears to be sufficient to reach a good agreement with numerical simulations.

2.4. Extension to the galaxy power spectrum

Although in this article we focus on the dark matter power spectrum, we briefly discuss in this section how the previous results can be used for the galaxy power spectrum, which is also a quantity of great practical interest. This necessarily involves additional ingredients, as we must define the relationship between the dark matter density field and the galaxies. As usual (Seljak 2000; Cooray & Sheth 2002), we keep the splitting (22) over 1-halo and 2-halo terms. Then, we simply write the galaxy 2-halo contribution as ${\mathit{P}}_{\mathrm{2}\mathrm{H}}^{\mathrm{gal}}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathrm{=}\mathrm{⟨}\mathit{b}{\mathrm{⟩}}^{\mathrm{2}}{\mathit{P}}_{\mathrm{2}\mathrm{H}}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathit{,}$(36)where  ⟨ b ⟩  is the mean bias of the galaxy population that we consider. This can be computed using one of the bias models that have been proposed in previous works (Mo & White 1996; Sheth & Tormen 1999; Sheth et al. 2001; Valageas 2009a, 2010b) or fits to numerical simulations (Tinker et al. 2010). Next, following the Eulerian interpretation of the counterterm of Eq. (35) discussed in the previous section, as arising from the difference  ⟨ n^galn^gal ⟩  −  ⟨ n^gal ⟩  ⟨ n^gal ⟩ , we write the 1-halo term as $\begin{array}{ccc}{\mathit{P}}_{\mathrm{1}\mathrm{H}}^{\mathrm{gal}}\mathrm{\left(}\mathit{k}\mathrm{\right)}& \mathrm{=}& \int \begin{array}{}\mathrm{\infty }\\ \mathrm{0}\end{array}\frac{\mathrm{d}\mathit{\nu }}{\mathit{\nu }}\mathit{f}\mathrm{\left(}\mathit{\nu }\mathrm{\right)}\hspace{0.17em}\frac{\mathit{M}}{\overline{\mathit{\rho }}\mathrm{(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathrm{3}}}\hspace{0.17em}\frac{{\overline{\mathit{\rho }}}^{\mathrm{2}}}{{\overline{\mathit{n}}}^{\mathrm{2}}}\hspace{0.17em}\frac{\mathrm{⟨}\mathit{N}\mathrm{(}\mathit{N}\mathrm{-}\mathrm{1}\mathrm{)}\mathrm{⟩}}{{\mathit{M}}^{\mathrm{2}}}\\ & & \mathrm{\times }\mathrm{(}{\begin{array}{}\mathrm{˜}\\ \mathit{u}\end{array}}_{\mathit{M}}\mathrm{(}\mathit{k}{\mathrm{)}}^{\mathit{p}}\mathrm{-}\begin{array}{}\mathrm{˜}\\ \mathit{W}\end{array}{\left(\mathit{k}{\mathit{q}}_{\mathit{M}}\right)}^{\mathit{p}}\mathrm{)}\mathit{,}\end{array}$(37)where $\overline{\mathit{n}}$ is the mean galaxy number density,  ⟨ N(N − 1) ⟩  is the mean squared number of galaxies (minus a factor N) in halos of mass M, and p = 1 or 2 depending on whether the former quantity is smaller or greater than unity (Seljak 2000; Cooray & Sheth 2002). Thus, we have simply subtracted to the usual expression the counterterm associated with galaxies set at random within the Lagrangian radius q_M (i.e., corresponding to a uniform density universe). Again, the quantity  ⟨ N(N − 1) ⟩  must be computed from a model defined for the population of galaxies that we consider.

The derivation of Eq. (37) is more phenomenological than the derivation of Eq. (35) from the Lagrangian point of view (29). However, since one always needs to introduce some phenomenological model to relate galaxies to the matter distribution (here through the two quantities  ⟨ b ⟩ (M) and  ⟨ N(N − 1) ⟩ (M)), this should be sufficient.

The contribution (37) to the galaxy power spectrum vanishes as k² in the large-scale limit, which again solves the unphysical nonzero limit obtained in previous implementations. This behavior remains valid even though galaxies are discrete objects, since this property only implies a constant shot-noise asymptote for P^gal(k) in the high-k limit. This shot-noise has actually been subtracted in Eq. (37), using the quantity  ⟨ N(N − 1) ⟩  instead of  ⟨ N² ⟩ , to follow the usual convention (Peebles 1980; Seljak 2000; Cooray & Sheth 2002). At low k, we are dominated by the 2-halo contribution and in case of constant large-scale bias we recover the slope P^gal(k) ∝ k^n_s of the primordial matter power spectrum¹.

3. A simple implementation

We describe in more detail in this section how we implement Eqs. (27) and (35) to compute the nonlinear density power spectrum.

3.1. “2-halo” contribution

We first consider the 2-halo contribution (27), and more precisely the perturbative part P_pert(k), as we shall discuss the non-perturbative prefactor F_2H(1/k) in Sect. 3.2 below.

Since we cannot compute and resum all terms of the perturbative expansion of P(k), contrary to the simpler Zeldovich dynamics, we must use for P_pert(k) an approximation associated with a truncation at a finite order. This leads to several possible choices, since a priori we can use the standard perturbation theory (Goroff et al. 1986; Bernardeau et al. 2002) as well as any of the various resummation schemes that have been introduced in the recent years (Crocce & Scoccimarro 2006b,a; Valageas 2007a; Matarrese & Pietroni 2007; Valageas 2008; Taruya & Hiramatsu 2008; Pietroni 2008). Indeed, all such methods give expressions for P_pert(k) that agree up to the order of the computation, and only differ by higher order terms.

However, it turns out that standard perturbation theory cannot suit our purposes because higher order terms grow increasingly fast at high k (Bernardeau et al. 2002). As we shall check in Sect. 6.4 below, at one-loop order it already gives a steep contribution that remains non-negligible at high k, because the prefactor F_2H(1/k) does not go to zero very fast (because for CDM power spectra σ(q) only grows logarithmically at small q if n_s = 1, and even remains finite if n_s < 1). This implies that to get a well-behaved 2-halo term we need to use other schemes, that provide a perturbative part P_pert(k) that remains small at high k as compared with the 1-halo term. Therefore, we must use one of the resummation schemes that have been recently developed and show a well-behaved high-k behavior.

It is interesting to point out that, although such approaches were devised in order to improve the accuracy of perturbation theory on large scale and to increase its range of validity (by improving its convergence through partial resummations of an infinite number of higher order terms), our study shows that a second important benefit of these schemes is to provide well-behaved expressions for P_pert(k). Here the point is not that the high-k behavior is accurately reproduced (since anyway perturbation theory does not apply in this range) but that it remains within reasonable bounds (typically close to or smaller than the linear power) so as to be subdominant. Thus, by giving systematic expansions of P_pert(k), which agree with perturbation theory up to the order of the computation while remaining well-behaved at high k, these schemes allow the construction of unified models such as those studied in this article, based on Eq. (22), that show a correct behavior at both low and high k. (In the same spirit, we have checked in Sect. 2.3 that the 1-halo term is well-behaved at low k.) This provides a second important motivation for such resummation schemes, as compared with standard perturbation theory.

In this article we focus on the “direct steepest-descent method” introduced in Valageas (2007a). In terms of diagrams it means that the two-point correlation C is given at one-loop order by the resummations shown in Fig. 1, where single lines are the linear correlation and response functions C_L and R_L, while the double lines are the nonlinear correlation and response functions C and R obtained at that order, see also Valageas (2008). At this one-loop order the standard perturbation theory consists in keeping only the three diagrams with zero or one bubble among the infinite series shown in the second equality in Fig. 1. We recall the details of the computation of the (perturbative) density power spectrum P_pert(k) within this approach in Appendix A.

This direct steepest-descent method is not necessarily the most accurate resummation scheme. In particular, it yields a response function that does not decay at high k or late times, but shows increasingly fast oscillations with an amplitude that follows the linear response function. This is not realistic, since one expects a Gaussian-like decay for Eulerian response functions, as can be seen from theoretical arguments and numerical simulations (Crocce & Scoccimarro 2006b,a; Valageas 2007b; Bernardeau & Valageas 2010a). However, the fast oscillations still provide an effective damping in a weak sense (that is when the response function is integrated over). The reason why we consider this direct steepest-descent method here is that it provides a simple and efficient method, which satisfies our requirements and proves to be reasonably accurate, as we shall check in Sect. 5 below. Indeed, while by construction it agrees with standard perturbation theory at one-loop order, it is well-behaved at high k as it remains close to the linear power spectrum on all scales (when we truncate the computation at one-loop order, as in this article), see Valageas (2007a).

It is also particularly efficient because the integrals involved in the computation of the power spectrum factorize (because the bubble in the upper right diagram of Fig. 1 involves a product of linear two-point correlation functions, whose time-dependence can be factorized). This factorization can be clearly seen in Eq. (A.29), as the two-time correlation C(k;η₁,η₂) can be written as the sum of three products, each of them being of the form ϕ(η₁) × ϕ(η₂). Therefore, there is no need to compute two-dimensional time integrals, over both η₁ and η₂. A second property is that the integro-differential equations satisfied by the nonlinear response R(k;η₁,η₂) (which is needed at intermediate stages to obtain the power spectrum) can be reduced to (third-order) ordinary differential equations, as seen in Eqs. (A.21), (A.22). This avoids the need to compute integrals over past times at each time step over η₁. These two properties allow a fast numerical computation, which can be a useful feature for practical purposes.

3.2. “1-halo” contribution

We now turn to the 1-halo contribution (35), which requires a model for the halo mass function and for the halo density profiles. In this article we consider the scaling function f(ν) given by $\mathit{f}\mathrm{\left(}\mathit{\nu }\mathrm{\right)}\mathrm{=}\mathrm{0.502}\mathrm{\left[}\mathrm{\right(}\mathrm{0.6}\hspace{0.17em}\mathit{\nu }{\mathrm{)}}^{\mathrm{2.5}}\mathrm{+}\mathrm{(}\mathrm{0.62}\hspace{0.17em}\mathit{\nu }{\mathrm{)}}^{\mathrm{0.5}}\mathrm{]}\hspace{0.17em}{\mathrm{e}}^{\mathrm{-}{\mathit{\nu }}^{\mathrm{2}}\mathit{/}\mathrm{2}}\mathit{.}$(38)It satisfies the normalization constraint (15) and it has already been shown to provide a good match to numerical simulations for halos defined by a density contrast δ = 200 (Valageas 2009a). Thus, substituting Eq. (38) into (8) we obtain the mass function of these halos, with a linear threshold δ_L = ℱ^-1(200) in the second Eq. (8). At redshift z = 0 this gives δ_L ≃ 1.59 for a ΛCDM cosmology, and it increases slightly to δ_L ≃ 1.6 at high z, see Valageas (2009a) and the fit of Eq. (12) therein. Next, the halo mass function also gives the probabilities F_1H(q) and F_2H(q) defined in Eqs. (20), (21). This yields in turn the prefactor F_2H(1/k) that enters the 2-halo term (27).

For the halo density profile we use the usual NFW profile (Navarro et al. 1997), ${\mathit{\rho }}_{\mathit{M}}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{=}\frac{{\mathit{\rho }}_{\mathit{s}}}{\mathit{x}\mathit{/}{\mathit{r}}_{\mathit{s}}{\mathrm{(}\mathrm{1}\mathrm{+}\mathit{x}\mathit{/}{\mathit{r}}_{\mathit{s}}\mathrm{)}}^{\mathrm{2}}}\mathit{,}$(39)which we truncate at the radius r₂₀₀, associated with the overall density contrast δ = 200. The scale radius r_s is defined in terms of the outer radius r₂₀₀ through the concentration c(M₂₀₀), $\mathit{c}\mathrm{\left(}{\mathit{M}}_{\mathrm{200}}\mathrm{\right)}\mathrm{=}\frac{{\mathit{r}}_{\mathrm{200}}}{{\mathit{r}}_{\mathit{s}}}\mathit{,}$(40)while the characteristic density is given by ${\mathit{\rho }}_{\mathit{s}}\mathrm{=}\overline{\mathit{\rho }}\hspace{0.17em}\frac{\mathrm{201}}{\mathrm{3}}\hspace{0.17em}\frac{{\mathit{c}}^{\mathrm{3}}}{\ln \mathrm{(}\mathrm{1}\mathrm{+}\mathit{c}\mathrm{)}\mathrm{-}\mathit{c}\mathit{/}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{c}\mathrm{)}}\mathrm{·}$(41)Therefore, it remains to specify the concentration parameter as a function of the halo mass, which we have labeled M₂₀₀ above to remind that it corresponds to the mass defined by the nonlinear density contrast δ = 200. Although we shall also consider fits to numerical simulations that have already been proposed in the literature, we shall find in Sect. 5 below that a better match to numerical simulations for the density power spectrum (especially in the high-k tail) is obtained with the following prescription, $\mathit{c}\mathrm{\left(}{\mathit{M}}_{\mathrm{200}}\mathrm{\right)}\mathrm{=}\mathrm{10.04}\hspace{0.17em}{\left(\frac{{\mathit{M}}_{\mathrm{200}}}{\mathrm{2}\mathrm{\times }{\mathrm{10}}^{\mathrm{12}} {\mathit{h}}^{-1}\hspace{0.17em}{\mathit{M}}_{\mathrm{\odot }}}\right)}^{-0.1}\hspace{0.17em}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{z}{\mathrm{)}}^{-0.8}\mathit{,}$(42)which is intermediate between the behaviors found in Dolag et al. (2004) and Duffy et al. (2008). However, this should not be trusted beyond the range studied in this paper (z ≤ 3).

Since we obtain Eq. (42) from a comparison with the density power spectrum measured in simulations, it can be understood as an “effective” concentration parameter that also includes the mean effect of halo substructures. Although it is possible to build more sophisticated halo models that explicitly describe substructures (Giocoli et al. 2010), we do not investigate such refinements here. In the same spirit, it could be possible to include the mean effect of baryons on the total matter profile of the halos through such an “effective” concentration parameter (or another choice for the shape of the halo profile than the NFW formula (39)). However, this is also beyond the scope of the present study, and we leave such investigations to future works. In any case, it is clear that any modeling of halo profiles (either at an “effective” level or through detailed explicit models) can be incorporated in the framework we describe in this paper. In particular, this could improve the accuracy at high k as simulations probe smaller scales and provide tighter (and more robust) constraints on halo properties.

For practical purposes it is convenient to use halo profiles that have an explicit expression for the Fourier transform ${\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{u}\end{array}}_{\mathit{M}}\mathrm{\left(}\mathit{k}\mathrm{\right)}$ defined in Eq. (31), especially if one is interested in high wavenumbers k where the integrand shows fast oscillations. Let us recall that for the NFW profile (39) one has (Scoccimarro et al. 2001) $\begin{array}{ccc}{\begin{array}{}\mathrm{˜}\\ \mathit{u}\end{array}}_{\mathit{M}}\mathrm{\left(}\mathit{k}\mathrm{\right)}& \mathrm{=}& {\left[\ln \mathrm{(}\mathrm{1}\mathrm{+}\mathit{c}\mathrm{)}\mathrm{-}\frac{\mathit{c}}{\mathrm{1}\mathrm{+}\mathit{c}}\right]}^{-1}\{\frac{\mathrm{-}\sin \mathrm{(}\mathit{ck}{\mathit{r}}_{\mathit{s}}\mathrm{)}}{\mathrm{(}\mathrm{1}\mathrm{+}\mathit{c}\mathrm{)}\mathit{k}{\mathit{r}}_{\mathit{s}}}\\ & & \mathrm{+}\cos \mathrm{(}\mathit{k}{\mathit{r}}_{\mathit{s}}\mathrm{)}[\mathrm{Ci}\left[\mathrm{(}\mathrm{1}\mathrm{+}\mathit{c}\mathrm{)}\mathit{k}{\mathit{r}}_{\mathit{s}}\right]\mathrm{-}\mathrm{Ci}\mathrm{(}\mathit{k}{\mathit{r}}_{\mathit{s}}\mathrm{)}]\\ & & \mathrm{+}\sin \mathrm{(}\mathit{k}{\mathit{r}}_{\mathit{s}}\mathrm{)}[\mathrm{Si}\left[\mathrm{(}\mathrm{1}\mathrm{+}\mathit{c}\mathrm{)}\mathit{k}{\mathit{r}}_{\mathit{s}}\right]\mathrm{-}\mathrm{Si}\mathrm{(}\mathit{k}{\mathit{r}}_{\mathit{s}}\mathrm{)}]\}\mathit{,}\end{array}$(43)where $\mathrm{Ci}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{=}\mathrm{-}{{}^{\mathrm{\int }}}_{\mathit{x}}^{\mathrm{\infty }}\cos \mathit{t}\hspace{0.17em}\mathrm{d}\mathit{t}\mathit{/}\mathit{t}$ and $\mathrm{Si}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{=}{{}^{\mathrm{\int }}}_{\mathrm{0}}^{\mathit{x}}\sin \mathit{t}\hspace{0.17em}\mathrm{d}\mathit{t}\mathit{/}\mathit{t}$ are the Cosine and Sine integrals (Abramowitz & Stegun 1970).

4. N-body simulations

In this section we describe the details of our N-body simulations, and we show the results of some tests of their possible systematic errors.

4.1. Set up

We adopt a flat ΛCDM model with cosmological parameters (Ω_m,Ω_b,h,σ₈,n_s) = (0.279,0.046035,0.701,0.817,0.96), which is consistent with 5 year observation of WMAP (Komatsu et al. 2009). We use a publicly available code, CAMB (Lewis et al. 2000), to compute the linear power spectrum. All the initial conditions for the simulations are created by adding displacements to particles put on a regular lattice by second-order Lagrangian perturbation theory (Scoccimarro 1998; Crocce et al. 2006) at z = 99, and evolved by a publicly available Tree Particle-Mesh code, GADGET2 (Springel 2005). We employ N = 2048³ particles in four periodic boxes with different box sizes (L_box = 512,1024,2048 and 4096   h^-1 Mpc). We call these runs as L9-N11, L10-N11, L11-N11 and L12-N11, respectively. The softening lengths for the tree force are set to be 5% of the mean inter-particle distance.

In addition to the four main simulations, we run five smaller simulations (L9-N8, L9-N9, L9-N10, L10-N9 and L11-N10) for convergence tests. The box sizes and number of particles of these simulations are summarized in Table 1.

4.2. Measuring power spectrum and two-point correlation function

We basically follow the ordinary FFT method to measure both the power spectrum and the two-point correlation function. Before assigning particles to the mesh, however, we fold the particle distribution into a smaller box to avoid the systematic error caused by the assignment on a finite number of grid points, which is important on small scales near the grid interval (see Colombi et al. 2009, for a similar discussion). Namely, we perform the following operation to the particle position, x: ${x}\mathrm{\to }{x}\hspace{0.17em}\mathrm{\%}\hspace{0.17em}\left({\mathit{L}}_{\mathrm{box}}\mathit{/}{\mathrm{2}}^{\mathit{n}}\right)\mathit{,}$(44)where the operation, a   %   b, gives the reminder of the division of a by b, and we vary the integer n from 0 to 6. After that, we assign particles on the 2048³ grid points of the folded small box with side length of L_box/2ⁿ using the CIC algorithm (Hockney & Eastwood 1981). We then Fourier transform the density contrast and take the average of $\mathrm{\left|}\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{\delta }\end{array}\mathrm{\right(}{k}\mathrm{)}{\mathrm{|}}^{\mathrm{2}}$ within k-bins to obtain the power spectrum. The binning is chosen to be logarithmic on small scale with 20 bins per decade (k > 0.3   h   Mpc^-1), but we adopt a linear binning with the interval Δk = 0.005   h   Mpc^-1 on large scale (k < 0.3   h   Mpc^-1) to see the feature of BAOs more clearly. While the folding procedure makes the systematic effect caused by assignment smaller, it reduces the number of available modes. We choose the value of n for the folding depending on the wavenumber so that the number of modes is large enough. We plot in Fig. 7 the resulting statistical error (1-σ level) estimated from $\frac{\mathrm{\Delta }\mathit{P}\mathrm{\left(}\mathit{k}\mathrm{\right)}}{\mathit{P}\mathrm{\left(}\mathit{k}\mathrm{\right)}}\mathrm{=}\frac{\mathrm{1}}{\sqrt{{\mathit{N}}_{\mathrm{mode}}\mathrm{\left(}\mathit{k}\mathrm{\right)}}}\mathit{,}$(45)where N_mode(k) denotes the number of independent mode in that bin. This is exact when $\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{\delta }\end{array}\mathrm{\left(}{k}\mathrm{\right)}$ is Gaussian (Feldman et al. 1994), and this can be shown to remain a good estimate in the current situation (see e.g., Takahashi et al. 2009). One may notice a dip at k ~ 0.3   h   Mpc^-1. This corresponds to the change of binning from linear to logarithmic. Another feature is a characteristic zig-zag pattern on smaller scale. This pattern corresponds to the change in the value of n used in the folding. Even after reducing the number of available modes by folding, however, we keep enough Fourier modes to achieve accuracy of  ~  sub-percent level on small scale $\mathrm{(}\mathit{k}{}^{\mathit{>}}\mathrm{~}_{\mathrm{~}}{}_{\mathrm{~}}\mathrm{1} \mathit{h}\hspace{0.17em}{\mathrm{Mpc}}^{-1}\mathrm{)}$.

The two-point correlation function is measured in a similar manner. We simply perform an inverse FFT to the three dimensional power spectrum calculated on grid points. We estimate the statistical error using the formula: ${\left|\mathrm{\Delta }\mathit{\xi }\mathrm{\left(}\mathit{x}\mathrm{\right)}\right|}^{\mathrm{2}}\mathrm{=}\mathrm{64}{\mathit{\pi }}^{\mathrm{4}}\mathrm{\int }\frac{\mathrm{d}\mathit{k}}{\mathit{V}}\hspace{0.17em}{\mathit{k}}^{\mathrm{2}}\hspace{0.17em}{\left\{{\mathit{j}}_{\mathrm{0}}\mathrm{\left(}\mathit{kx}\mathrm{\right)}\mathit{P}\mathrm{\left(}\mathit{k}\mathrm{\right)}\right\}}^{\mathrm{2}}\mathit{,}$(46)where the quantity j₀ is the spherical Bessel function of the first kind. We use the power spectrum measured from the same simulation in the integrand for consistency, see Taruya et al. (2009) for more details.

4.3. Convergence tests

One may naively think that one can model the power spectrum by combining simulations with different box sizes. However, it is not trivial to determine which simulation gives the best measurement on a given scale. The simulation in a larger volume lacks the power on small scale, while that in a smaller box does so on large scale.

A characteristic wavenumber for the finite box size is given by 2π/L_box, the fundamental mode. For the finiteness of mass/force resolution, one may define two important wavenumbers. First, the initial density fluctuation cannot be input on scales smaller than the mean inter-particle distance (2π/L_box × N^1/3 in wavenumber). Second, the force calculation is not accurate below the softening scale. This wavenumber is given by 2π/L_box × N^1/3 × 20 for our simulations (i.e., 5% of the mean inter-particle distance). One can basically trust the simulation results from the fundamental mode to the inter-particle distance, or at most to the softening scale. These wavenumbers are shown in the left panel of Fig. 2 for the four main simulations.

After nonlinear evolution, however, the initial small systematic error on one scale may propagate to different scales due to the coupling between different modes. This makes the situation complicated. To separate the two systematic effects caused by the finiteness of mass/force resolution and simulated volume, we perform two tests using additional simulations in this subsection. We will show that both finiteness effects lead to a lack of power on certain scales. We especially focus on the scales where these effects are more than 1%, and avoid to use the data points on these scales in later discussions.

4.3.1. Test 1: finite mass/force resolution

We first examine the effect of the finiteness of the mass/force resolution. We use L9-N8, L9-N9, L9-N10 and L9-N11 for this test (see Table 1). These simulations have the same volume but different number of particles, so that we can see the systematic effect purely from the finiteness of the resolution (see the middle panel of Fig. 2). We set the same random linear density field in creating the initial conditions of the four simulations to make the comparison clearer. Since the three additional simulations, L9-N8, L9-N9 and L9-N10, have the same resolution as L12-N11, L11-N11 and L10-N11, respectively, we can assess the impact of the systematic effect for the “main” simulations by this test.

The results are shown in Fig. 3. We plot the ratio of the power spectrum measured from the three lower resolution runs to that from the highest resolution one (i.e., L9-N11). One can clearly see that all curves blow up at large wavenumbers. This corresponds to the scale where the shot noise contribution becomes important. In addition to this feature, the curves show a lack of power on intermediate scales, which is more important for lower resolution runs at higher redshifts. We can interpret this as follows. At higher redshifts, the lack of power on around mean inter-particle scale survives. But the power generated by the pure nonlinear growth gradually dominates over the initial power and the difference between the four runs become smaller at lower redshifts.

We also test the effect of finite resolution on the two-point correlation function. Similar trends can be seen in Fig. 4 for the two-point correlation function.

4.3.2. Test 2: finite box size

We next discuss the effect of the finiteness of the simulated volume. We use L9-N8, L10-N9, L11-N10 and L12-N11 for this test (see Table 1). These runs have the same mass/force resolution but different box size (see the right panel of Fig. 2). Again, since the three smaller simulations have the same volume as the main runs, L9-N11, L10-N11 and L11-N11, we can test the possible systematics of these main simulations.

We show the results in Fig. 5. In contrast to the test for the finite resolution, the data points seem noisy. This is because we cannot set the same random realization of the linear density field for simulations with different box sizes. When one sees the result of L9-N8, the systematic effect is prominent at 0.3 ~ 2   h   Mpc^-1, and is growing with time. This is contrary to what we see in the finite resolution test, which as we showed was decaying. This growing nature of the systematics suggests that the effect comes from the mode coupling of large-scale modes. The systematic effect is not important (i.e., below 1% level) for the rest of the simulations.

We test the convergence of the two-point correlation function which is shown in Fig. 6. Unfortunately, however, we cannot determine a strict range of the trustable scales due to large statistical error. Since the two-point correlation function is a weighted sum of the power spectrum over all scales, a large error bar at a certain wavenumber may propagate to the two-point correlation function over a wide range of scales. Furthermore, the two-point correlation function is not positive definite unlike the power spectrum. Thus, the size of its fractional error can be very large on scales where the two-point correlation function is close to zero. Finally, we cannot control the size of error bars by changing the width of the bins, since it does not explicitly depend on this width [see Eq. (46)].

The ratio of the two-point correlation function on scales above  ≳5   h^-1 Mpc seems consistent with unity at all the three redshifts. As seen in Fig. 6, the finite resolution effect is important below  ≲5   h^-1 Mpc for L9-N8. This prevents us from a clear test for the finite box size since the four simulations used in this test have the same resolution as L9-N8. We conclude here that we do not find any clear evidence of systematic error on the two-point correlation function arising from the finiteness of volume, and we will determine which simulation to use following only the resolution test. We leave more complete tests for the finite volume to future works.

5. Comparison with numerical simulations

We first show in Fig. 8 the quantity F_2H(1/k), which is defined by Eq. (21) and enters the 2-halo term (27), at redshifts z = 0.35, 1, and 3. At fixed comoving scale, q = 1/k, it decreases at lower redshift, since the probability F_1H for a Lagrangian pair of distance q to belong to the same collapsed halo grows with time, as larger scales turn nonlinear, and F_2H = 1 − F_1H. We can note that it decreases very slowly at high k. In fact, since we consider an initial power spectrum with a primordial slope n_s = 0.96 that is smaller than unity, the rms linear density contrast σ(q) of Eq. (10) does not diverge to infinity as q → 0 but goes to a finite asymptote. This implies that the probability F_2H(1/k) does not go to zero at high k (i.e. on small scales), but reaches the nonzero asymptote ${\mathit{F}}_{\mathrm{2}\mathrm{H}}\mathrm{\left(}\mathrm{0}\mathrm{\right)}\mathrm{=}{{}^{\mathrm{\int }}}_{\mathrm{0}}^{{\mathit{\delta }}_{\mathrm{L}}\mathit{/}\mathit{\sigma }\mathrm{\left(}\mathrm{0}\mathrm{\right)}}\mathrm{d}\mathit{\nu f}\mathrm{\left(}\mathit{\nu }\mathrm{\right)}\mathit{/}\mathit{\nu }$. Although this value is unlikely to be exact, because of the approximations involved in the derivation of F_2H(q), it is possible that for such power spectra without significant power on small scales there always remains some fraction of matter which has never experienced shell crossing and remains described by the single-stream fluid approximation. Even though this is an interesting point from a theoretical perspective, regarding the dynamical properties of the system, this is beyond the subject of the present work and it plays no role for our purposes, since at high k the power spectrum is clearly dominated by the 1-halo term, as we shall check in Fig. 9 below. On the other hand, we can see that if we require a high accuracy, the deviation of F_2H from unity at k ~ 0.3   h   Mpc^-1 for z = 1 cannot be completely neglected and has some effect on the 2-halo term (27).

We now compare our model for the density power spectrum, that combines a systematic perturbative approach with a phenomenological halo model, with results obtained from numerical simulations. We show in Fig. 9 the power per logarithmic interval of k, defined as ${\mathrm{\Delta }}^{\mathrm{2}}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathrm{=}\mathrm{4}\mathit{\pi }{\mathit{k}}^{\mathrm{3}}\mathit{P}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathit{,}$(47)for the three redshifts z = 0.35, 1, and 3. The sudden rise of the N-body results at very high k (which is more apparent in the right panel at z = 3) is due to the finite resolution and shot noise, see Sect. 4. This marks the highest wavenumber where we can compare the simulations with our theoretical models.

As discussed in Sect. 3.1, we can check that the 2-halo term ${\mathrm{\Delta }}_{\mathrm{2}\mathrm{H}}^{\mathrm{2}}$ remains well-behaved at high k and is subdominant with respect to the 1-halo term ${\mathrm{\Delta }}_{\mathrm{1}\mathrm{H}}^{\mathrm{2}}$. It falls somewhat below the linear power ${\mathrm{\Delta }}_{\mathrm{L}}^{\mathrm{2}}$, contrary to the one-loop resummation of ${\mathrm{\Delta }}_{\mathrm{pert}}^{\mathrm{2}}$ given by the direct steepest-descent scheme, which becomes very close to ${\mathrm{\Delta }}_{\mathrm{L}}^{\mathrm{2}}$ as seen in Valageas (2007a), because of the prefactor F_2H(1/k) in Eq. (27).

In agreement with the discussion of Sect. 2.3, we can also check that the 1-halo term shows a fast decline at low k. At high k we can see that it is possible to reach a good agreement with the numerical simulations by using an appropriate prescription for the concentration parameter c(M), such as the one given in Eq. (42) and represented by the green solid line. The formula given in Eq. (42) is obtained by looking for values of the free parameters (the normalization and the two exponents) that provide a reasonable match with the power spectrum measured in the simulations at high k (looking among a few values close to the fits already proposed in the literature). Nevertheless, it is interesting to note that using fits for c(M) proposed in previous analysis of the halo profiles formed in numerical simulations, one obtains predictions that are either larger (Dolag et al. 2004) or smaller (Duffy et al. 2008) than the one associated with Eq. (42). This shows that the nonlinear power spectrum is fully consistent in this range with a simple halo model, such as (35), and with the properties of halos seen in numerical simulations. This also gives an estimate of the dependence of the power spectrum on the prescription used for c(M). In agreement with previous works (Huffenberger & Seljak 2003; Giocoli et al. 2010), we recover the fact that the nonlinear power Δ²(k) is larger and steeper at high k for concentration relations c(M) that have a larger normalization and a steeper dependence on M. A nice feature is that the nonlinear power is largely independent of the details of the halo model up to Δ² < 100, so that models such as the one studied in this article remain quite predictive. Even at higher k, we can see that up to Δ² < 10³, or k < 100   h   Mpc^-1, using any of these prescriptions for c(M) provides a reasonable estimate of the power spectrum and even fares better than the direct fit to P(k) that was obtained from older simulations (Smith et al. 2003). This suggests that models based on phenomenological ingredients such as the halo model may prove more robust than direct fits to numerical results. However, our model for P(k) should not be trusted beyond the domain where it has been checked, that is k ≤ 100   h   Mpc^-1 and z ≤ 3.

In order to obtain predictions at much higher redshifts and wavenumbers, one should use a prescription for c(M) that is based on some physical arguments rather than simple fits such as Eq. (42). This would probably lead to a loss of accuracy in the range tested in Fig. 9, as compared with the use of Eq. (42), but this is likely to be more robust as we extrapolate to other regimes. However, we shall not investigate the building of physical models for c(M) in this article, as this is a topic by itself.

We can see in Fig. 9 that the full nonlinear power spectrum (22) obtained by our approach, combining a perturbative expansion with a phenomenological halo model, is able to reproduce the results measured in numerical simulations, up to a reasonable accuracy. Let us point out that the perturbative part, P_pert(k), which dominates on large scales, contains no free parameter, since it is given by a systematic perturbation theory. The 1-halo contribution that dominates on small scales clearly contains some parameters, through the choice of the halo mass function and halo density profile. However, the mass function and the shape of the halo profile are already set by other measures from numerical simulations, so that the main free parameter is the concentration c(M), which is not as well constrained. However, as seen in Fig. 9, this uncertainty only affects the very high-k tail. Moreover, even in this region the resulting model is competitive with direct fits to the power spectrum measured in older simulations. Therefore, models such as the ones developed in this article should prove useful to obtain reliable predictions for P(k). Another advantage of such models, as compared with simpler fitting formulas, is that they should be quite accurate on large scales since they are consistent with perturbation theory, up to the order of truncation of the computation.

Thus, we plot in Fig. 10 the ratio of the power spectrum to a “no-wiggle” linear power spectrum P_Ls(k) which does not include baryonic acoustic oscillations, in order to see more clearly the low-k behavior. We can check that in this range our model, based on the direct steepest-descent resummation, shows a good match with the numerical simulations. As could be expected, it fares better than a direct fit from Smith et al. (2003) which was not designed to model this range with a high accuracy (but still remains surprisingly good). Furthermore, the use of the direct steepest-descent resummation proves to provide significantly more accurate results than the standard perturbation theory, truncated at one-loop order, given by the well-known expressions ${\mathit{P}}_{\mathrm{1}\mathrm{loop}}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathrm{=}{\mathit{P}}_{\mathrm{L}}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathrm{+}{\mathit{P}}_{\mathrm{22}}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathrm{+}{\mathit{P}}_{\mathrm{31}}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathit{,}$(48)where formulas for the terms P₂₂ and P₃₁ may be found for instance in Bernardeau et al. (2002). Let us recall here that the term P_pert(k) given by this direct steepest-descent scheme contains no free parameter, nor any interpolation procedure, and is consistent with standard perturbation theory up to one-loop order (i.e. the difference between P_1loop and P_pert is due to the partial resummation of higher order terms). We can see in the left panel, at z = 0.35, that around k ~ 0.14   h   Mpc^-1 the curve P_pert is slightly above the full nonlinear power spectrum P. This is due to the prefactor F_2H(1/k) in Eq. (27), which is slightly below unity. However, this is only a very small effect on these scales. On the other hand, at higher k (e.g., k > 0.26   h   Mpc^-1 at z = 0.35) the nonlinear power spectrum rises above P_pert and keeps growing, while P_pert remains close to P_L at high k as seen in Fig. 9. This is due to the 1-halo contribution, which starts being non-negligible. However, on these scales the dependence on the details of the halo model is extremely weak. Indeed, the three green curves plotted in Fig. 9, associated with the prescription (42) for the concentration c(M) and the two fits given by Dolag et al. (2004) and Duffy et al. (2008), are also plotted in Fig. 10. However, they almost exactly fall on top of each other and cannot be distinguished in the figure.

Thus, Figs. 9 and 10 show that by combining perturbation theories and halo models it is possible to obtain a good model for the nonlinear density power spectrum, both on quasi-linear and highly nonlinear scales. However, we can see in Fig. 9 that in the intermediate regime, where Δ² ~ 5, our predictions fall below the N-body results. This is also apparent in the high-k parts of Fig. 10, where the full nonlinear prediction P(k) starts to grow more slowly than the power measured in the numerical simulations. This regime corresponds to the transition between the 2-halo and 1-halo contributions (see Fig. 9) and as such it is at the limit of validity of the approximations used for both terms.

On the perturbative side, that is the 2-halo term, the discrepancy can be due to the truncation at one-loop order of the perturbative term. Indeed, we can expect that by going to higher orders we can extend the range of validity of the the perturbative term P_pert and push the downturn shown by the blue dotted curve in Fig. 10 to higher k. Within such resummation schemes this means that we include all diagrams up to n loops, and partial resummations for higher order terms. We can note that the discrepancy looks somewhat more severe at z = 3 in Fig. 9. More precisely, the range of k where there is a noticeable mismatch before the 1-halo term becomes dominant (larger than P_L) is somewhat more extended than at z = 0.35. This agrees with the results of Valageas (2010a), where it was found (within the Zeldovich framework) that the scope of perturbation theory is somewhat greater at higher redshift for CDM power spectra, in the sense that the range where higher order perturbative terms are important (i.e. larger than the non-perturbative correction associated with shell crossing effects) is wider and that the perturbative expansion makes sense up to higher orders. However, we shall not try to go beyond one-loop order in this paper and we leave such a task to future works.

On the non-perturbative side, it is clear that by definition halo models are only phenomenological models and do not provide a systematic tool. Even if we admit that it makes sense to “recognize” halos in the density field (even though this is a descriptive tool rather than a direct ingredient of the equations of motion), it is clear that realistic halos are not spherical nor fully virialized. In particular, apart from the inaccuracies introduced by the spherical approximations, outer halo radii are not fully virialized and particles do not instantly loose all memory of their initial conditions, so that the average in Eq. (33) is only approximate. As we have discussed in Sect. 2.3, such approximations also prevent us from exactly satisfying momentum conservation. Since these effects mostly apply to the outer parts of the halos they can be expected to have an impact on the moderately low-k part of the power spectrum, that is the transition region. Another effect that is likely to play a role is the truncation of the halo profiles at the radius r₂₀₀. This may be investigated in a rather straightforward manner by studying the change in P(k) as we modify this density threshold, and we shall consider in Sect. 6.1 below the results obtained for δ = 50.

On a more fundamental level, since the splitting between 2-halo and 1-halo contributions is necessarily approximate (it cannot be derived in a systematic fashion from the equations of motion and always involves some approximations, at least within analytical frameworks) it is not surprising that the match is not perfect, especially in the transition range.

Another approximation used in this work is the assumption of smooth halos, defined by the regular profile (39). As noticed in Sect. 3.2, it is possible to extend the halo model by including substructures (Giocoli et al. 2010), however this only has an effect on the very high-k tail of the power spectrum (as could be expected since these are small-scale modifications). We do not investigate such modifications in this work, as the agreement we obtain with numerical simulations is already reasonably good on these scales and we prefer not to introduce additional parameters. These effects are also degenerate with the mean concentration relation c(M), as seen in Fig. 9 by the comparison between the dashed and dotted green curves, associated with different prescriptions from Dolag et al. (2004) and Duffy et al. (2008). In particular, as pointed out in Sect. 3.2, the prescription (42) that we have obtained by requiring a good match at high k can also be seen as an effective model, that describes the combined effects of both the concentration of the mean radial profile and the presence of small substructures. At some stage, if one wants to build a model where the halo properties are independently obtained from the measure of individual halos in N-body simulations, so that c(M) can no longer be tuned, it may nevertheless prove necessary to explicitly include such effects.

6. Dependence on various ingredients of the model

We now investigate the impact of various ingredients of the model on the predictions obtained for the density power spectrum.

6.1. Halos defined by δ = 50

As we have discussed in the previous section, the lack of power on intermediate scales seen in Fig. 9, around the transition between the 2-halo and 1-halo terms, may be due in part to the truncation of the halo profiles at the radius r₂₀₀. Indeed, realistic halos extend to larger radii, which should generate some extra power on large scales. In order to take into account such outer radii, while keeping the mass M that enters the definition of the halo mass function as the halo mass within the truncation radius (so as to avoid overcounting and violations of mass normalizations), we consider defining halos by a smaller nonlinear density contrast, such as δ = 50. This lower threshold means that we include outer shells (as compared with the case δ = 200), while halos still remain well separated from the background.

For the halo mass function we still use the scaling function (38), but the linear density contrast is now δ_L = ℱ^-1(50), which gives for instance δ_L ≃ 1.50 at z = 1. This automatically satisfies the normalization (15) as well as the large-mass tail (14). It was shown in Valageas (2009a) that this recipe provides a reasonable match to numerical simulations (but not as good as for δ = 200) for halos defined by a density threshold δ = 100, and we shall assume that it still provides a reasonable model for δ = 50. This should be sufficient for our purposes, as the 1-halo term needed for the density power spectrum only depends on integrals over the halo mass function and we have seen that both the normalization and the large-mass tail are correctly reproduced.

We still use the NFW profile of Eq. (39), but we now need to obtain the scale radius r_s in terms of r₅₀. To do so, we now define the concentration parameter as $\mathit{c}\mathrm{\left(}{\mathit{M}}_{\mathrm{50}}\mathrm{\right)}\mathrm{=}\frac{{\mathit{r}}_{\mathrm{50}}}{{\mathit{r}}_{\mathit{s}}}\mathit{,}$(49)and the halo characteristic density now reads as ${\mathit{\rho }}_{\mathit{s}}\mathrm{=}\overline{\mathit{\rho }}\hspace{0.17em}\frac{\mathrm{51}}{\mathrm{3}}\hspace{0.17em}\frac{{\mathit{c}}^{\mathrm{3}}}{\ln \mathrm{(}\mathrm{1}\mathrm{+}\mathit{c}\mathrm{)}\mathrm{-}\mathit{c}\mathit{/}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{c}\mathrm{)}}\mathrm{·}$(50)We again look for a prescription for c(M₅₀) of the form of Eq. (42), and we find in Fig. 11 below that a reasonable agreement with the numerical simulations can be achieved by using $\mathit{c}\mathrm{\left(}{\mathit{M}}_{\mathrm{50}}\mathrm{\right)}\mathrm{=}\mathrm{13}\hspace{0.17em}{\left(\frac{{\mathit{M}}_{\mathrm{50}}}{\mathrm{2}\mathrm{\times }{\mathrm{10}}^{\mathrm{12}} {\mathit{h}}^{-1}\hspace{0.17em}{\mathit{M}}_{\mathrm{\odot }}}\right)}^{-0.1}\hspace{0.17em}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{z}{\mathrm{)}}^{-0.7}\mathit{.}$(51)This is mostly set by the behavior of the power spectrum at high k, where Δ² > 100, as seen in Fig. 9 from the comparison between various models for c(M) for halos defined by δ = 200. Thus, the properties of the density power spectrum on large scales, where Δ² < 10, are almost independent of the choice of parameters in Eq. (51) (restricted to a reasonable range).

We compare in Fig. 11 the power spectra obtained on nonlinear scales by using halos defined either by δ = 200 or δ = 50 (while the perturbative term P_pert(k) is not modified). As expected, truncating halo profiles at a larger radius yields some extra power at large k and improves somewhat the agreement with numerical simulations in the intermediate range, where Δ² ~ 5. However, the change is quite small and the match to the N-body results is not perfect yet. In fact, at redshift z = 0.35, we obtain a smaller power spectrum around Δ² ~ 100 by using halos defined by δ = 50. However, one may modify this by using a different concentration relation c(M₅₀) at that redshift. In any case, these results suggest that the truncation at r₂₀₀ is not the unique reason for the discrepancy between the model and N-body simulations in this transitory range, and as explained in Sect. 5 another possible source of inaccuracy is the partial resummation of high-order terms in the perturbative term P_pert.

Next, we compare in Fig. 12 the power spectra obtained on quasi-linear scales. We can see that on these very large scales, k < 0.3   h   Mpc^-1, the modification of the 1-halo term has a very weak effect. However, it improves somewhat the agreement with the simulations at low redshift, z = 1 (at k ~ 0.3   h   Mpc^-1, as the gap below the N-body results in the intermediate range is repelled to slightly higher k), although it leads to a slight overestimate for P(k) at z = 0.35. On the other hand, it appears to make almost no change over these scales at higher redshift, z = 3. In agreement with the discussion in Sect. 5 and with Valageas (2010a), this expresses the fact that for CDM power spectra the scope of perturbative expansions is somewhat broader at higher z, in the sense that the range of k where higher order terms play a role is wider. This corresponds to the interval  [ k_1loop,k_s.c. ]  where terms beyond linear order are already significant, as compared with the linear power P_L, but still larger than the non-perturbative correction associated for instance with shell crossing effects.

That our results for the density power spectrum only show a very weak dependence on the details of the halo definition (by δ = 200 or 50) can actually be considered as a reassuring sign, rather than a problem in the attempt to improve the agreement with simulations. Indeed, if the splitting (22) and all subsequent computations were exact, halos would only appear as an intermediate tool and would disappear in the final result, since the matter density field is independent of how we split it over halos at later stages. Therefore, the final result for P(k) should not depend on the details of the definition of the halos. In practice, this cannot be the case because the splitting (22) itself is only approximate as long as we consider spherical halos, and the subsequent treatment of both 2-halo and 1-halo contributions involves some approximations, as in Eqs. (27) and (32). This leads to some artificial dependence on the details of the halo model, but as shown in Figs. 11 and 12 this effect is rather small and argues for the robustness of the model.

Therefore, Figs. 11 and 12 suggest that there is still room for systematic improvement, by including higher order terms in the perturbative contribution P_pert. However, this may come at the price of heavier and slower computations, which we do not investigate here. On the other hand, they show that predictions on large scales are largely independent of the details of the underlying halo model.

6.2. Impact of the “1-halo” counterterm

We have seen in Sects. 5 and 6.1 that quantitative details of the 1-halo term can have some effect on the high-k tail (the concentration c(M) of the halo profiles) or on the intermediate range where Δ²(k) ~ 5 (the truncation radius of the halos). We investigate here the role of a more fundamental ingredient of the 1-halo contribution, namely the counterterm $\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{W}\end{array}\mathrm{(}\mathit{k}{\mathit{q}}_{\mathit{M}}{\mathrm{)}}^{\mathrm{2}}$ in Eq. (35). Thus, focusing on the case of halos defined by the nonlinear density contrast δ = 200, we compare in Fig. 13 the predictions of our fiducial model with those that are obtained by setting to zero the 1-halo counterterm $\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{W}\end{array}\mathrm{(}\mathit{k}{\mathit{q}}_{\mathit{M}}{\mathrm{)}}^{\mathrm{2}}$ of Eq. (35). Thus, we can see that, even though the 1-halo counterterm plays a negligible role at high k, it cannot be discarded at low k if we require a high accuracy. As explained in Sect. 2.3, this is also related to the fact that this counterterm is required to obtain the k² tail at low k expected for the 1-halo contribution (if we do not insist on momentum conservation). Discarding this counterterm leads to a nonzero asymptote at low k for P_1H(k), which eventually dominates over the linear power spectrum, which roughly decreases as P_L(k) ~ k for CDM cosmologies. On the scales probed in Fig. 13, k ~ 0.2   h   Mpc^-1, we have not entered yet the very low k regime where this would give rise to a spurious divergence. However, we can clearly see the spurious extra power that would be generated by the lack of this counterterm. Note that on these scales there are no fitting parameters that could be tuned to compensate for this extra power. Indeed, the details of the halo model, such as the concentration c(M) (see Fig. 10), or the truncation radius of the profiles (see Fig. 12), only have a weak impact in this domain and could not balance this extra power. On the other hand, the 2-halo term is also highly constrained, as the blue dashed and dotted curves in Fig. 10, associated with standard 1-loop perturbation theory and steepest-descent resummation, give an estimate of the possible scatter between the predictions that can be obtained from various perturbative schemes. This is most clearly seen around the peak at k ~ 0.13   h   Mpc^-1 for z = 0.35 (as could be expected the extra power associated with such an incorrect modelization of the 1-halo term grows at lower z as larger halos form). By contrast, the reasonably good agreement with simulations obtained by our complete model (solid lines) suggests that it does not miss important physical ingredients (as discussed in Sect. 5 difficulties mostly arise at higher k, associated with the transition between the 2-halo and 1-halo contributions).

6.3. Benefit of higher order perturbative terms

Before we investigate other choices for the perturbative term P_pert in the next sections, we first compare in Figs. 14 and 15 our model with the more standard approach, where one simply uses the linear power P_L for the 2-halo contribution. Thus, we replace Eq. (27) by the approximation P_2H(k) = P_L(k), while we keep the same 1-halo contribution (35), with its new counterterm that ensures a satisfactory behavior at low k. In agreement with previous works, we can see in Fig. 14 that on the large scales associated with the baryon acoustic oscillations the approximation P_2H(k) = P_L(k) is not sufficient to obtain a good match to the numerical simulations, contrary to the use of the one-loop perturbative term P_pert(k) obtained by the “steepest-descent” scheme recalled in Sect. 3.1. This means that on these scales the departure from the linear power is not yet due to the non-perturbative contribution associated with the 1-halo term, but to higher order perturbative terms. This agrees with the results of Valageas (2010a), based on the Zeldovich dynamics, which show that many orders of perturbation theory are relevant before the power spectrum is dominated by the non-perturbative contribution associated with shell crossing or halo formation (typically one can go up to order ${\mathit{P}}_{\mathrm{L}}^{\mathrm{9}}$ at z = 0 and ${\mathit{P}}_{\mathrm{L}}^{\mathrm{66}}$ at z = 3 for a ΛCDM cosmology, at least for this simpler case). This motivates the use of perturbative approaches that include higher order contributions.

On smaller scales, shown in Fig. 15, the difference between our model for P_2H(k) and the popular approximation P_2H(k) = P_L(k) is negligible. This agrees with Fig. 9, where we could see that at high k the one-loop resummed contribution to P_2H(k) is near the linear power P_L(k) and subdominant as compared with the 1-halo contribution.

6.4. Eulerian perturbation theories

We consider here two other choices for the perturbative term P_pert of the 2-halo contribution (27). A first choice is simply to use the standard perturbation theory and to write at one-loop order P_pert = P_1loop(k), where P_1loop(k) was given in Eq. (48). A second choice is to keep the diagrammatic expression given by the first line in Fig. 1, but to replace the external double lines by the Gaussian-decay response function (propagator) obtained by Crocce & Scoccimarro (2006b,a) instead of using the response function obtained by the steepest descent resummation scheme. This is not identical to the “RPT” resummation used in Crocce & Scoccimarro (2008), where the “bubble” in the upper right diagram of Fig. 1 is also replaced by resummed two-point correlations, obtained from this same nonlinear propagator (which corresponds to the partial resummation of further diagrams). Here we keep linear two-point correlation functions in this “bubble” to keep factorizable integrals, as in Eq. (A.29), since our goal is simply to estimate the dependence on the choice of resummation scheme. Moreover, it is interesting to evaluate various factorizable schemes of this kind as they allow faster computations. Two other differences with the prescription used in Crocce & Scoccimarro (2008) are the presence of the prefactor F_2H(1/k) in Eq. (27), and the fact that we do not introduce an extra tuning in the decay of the response function by replacing the linear velocity dispersion σ_v, which governs its high-k cutoff, by a nonlinear velocity dispersion obtained from fits to the nonlinear power spectrum from numerical simulations. Indeed, one of the main points of our approach is precisely that the term P_pert(k) in the 2-halo contribution (27) should be obtained from systematic perturbation theories, to ensure a good accuracy and robustness for a variety of cosmologies, so that we avoid introducing fitting parameters in this factor. Then, these two alternatives, “standard perturbation theory” (SPT) and “Gaussian-like decaying response function” (DR), and the complete steepest-descent resummation used in the previous sections, agree with each other up to one-loop order.

We first show in Fig. 16 the power spectra obtained on large scales for redshifts z = 0.35,1, and 3. We plot the results obtained by using either the 1-loop “standard perturbation theory”, or the “decaying response” function from Crocce & Scoccimarro (2006a) in the external double lines of the two upper diagrams of Fig. 1, for the term P_pert in the 2-halo contribution (27). In both cases the 1-halo term is the same as the one used in Sect. 5, so that these curves only differ from those shown in Fig. 10 through the term P_pert. As was already noticed in Fig. 10, and in agreement with previous works (Crocce & Scoccimarro 2008; Taruya & Hiramatsu 2008; Carlson et al. 2009), the 1-loop standard perturbation theory overestimates the power on these scales, especially at low redshift. The dashed curves shown in Fig. 16 are not identical to those shown in Fig. 10 and in other works, because the perturbative term P_pert is multiplied by the prefactor F_2H(1/k) in Eq. (27). Moreover, the results also include the (small) 1-halo contribution. However, on these scales the prefactor F_2H(1/k) remains very close to unity (see Fig. 8) and it is not sufficient to significantly improve the agreement with the N-body results. The use of the “decaying response” function leads to results that are very close to those obtained with standard perturbation theory and also overestimates the power spectrum on these large scales. Using the full “renormalized perturbation theory” (RPT) described in Crocce & Scoccimarro (2008) leads to a smaller prediction for P(k). This can be understood from the fact that inserting the “decaying response” function into the bubble of the upper right diagram of Fig. 1, as in “RPT”, clearly suppresses the nonlinear power and explicit computations show that this yields a good match to results from simulations (Crocce & Scoccimarro 2008; Carlson et al. 2009).

Next, we show in Fig. 17 the power per logarithmic interval of k, as in Fig. 9. We plot both the full nonlinear power Δ², as in Fig. 16, and the 2-halo contribution, for the “SPT” and “DR” cases. As already pointed out in Sect. 3.1, we can see that standard perturbation theory fails at a qualitative level as the 2-halo contribution keeps growing at high k and becomes very large, whereas the term P_pert being associated with particle pairs that have not collapsed within a single object it should remain near P_L at most. Moreover, this spurious growth is not suppressed by the prefactor F_2H(1/k), which does not decrease sufficiently fast at high k, as seen in Fig. 8. Even though this 2-halo contribution remains smaller than the 1-halo contribution on these scales, it leads to a significant extra power that worsens the agreement with the numerical simulations. As can be seen from Figs. 9 and 11, this extra power cannot be compensated by modifications to the underlying halo model (such as the concentration c(M) or the halo truncation radius), especially in the range 5 < Δ² < 200 where the dependence on the details of the halo model is rather weak. Of course, this problem would become increasingly severe as one includes higher order terms that grow increasingly fast at high k (Bernardeau et al. 2002; Crocce & Scoccimarro 2006b). This means that standard perturbation theory cannot be used to build systematic models for the density power spectrum, such as those studied in this paper. As pointed out in Sect. 3.1, this provides a further motivation to investigate resummation schemes that would be better behaved.

We can see that using the “decaying response” function in the external double lines of Fig. 1 leads to a 2-halo contribution that remains small and close to the linear power spectrum at high k, as for the direct steepest-descent scheme used in Fig. 9. This shows again the benefit of using resummation schemes instead of standard perturbation theory to obtain a well-behaved perturbative term P_pert(k) that also agrees with standard perturbation theory up to the order of truncation. In fact, it happens that using this “decaying response” function yields a somewhat larger perturbative term than the one obtained in Sect. 5 on large scales, as seen in Fig. 16, and this extra power actually improves the agreement with the N-body results in the range where 1 < Δ² < 100. Indeed, as discussed in Sect. 5, the implementation used in Fig. 9 underestimates the nonlinear power spectrum in this range, and the additional power due to the use of this “decaying response” function helps to bridge the gap. It actually provides a very good quantitative match, although at z = 3 we again underestimate somewhat the power spectrum in the transition range. Thus, for practical quantitative purposes this may provide a good variant of the model developed in this paper to compute the power spectrum, if one is interested in these transition scales. However, because on large scales this prescription actually gives a worse match to numerical simulations, as seen in Fig. 16, the better agreement on the transition scales is probably not more than a lucky coincidence. Thus, it is likely to compensate a meaningful lack of power (due to the various reasons discussed in Sect. 5, such as the neglect of some higher order perturbative terms or the limitations of the description in terms of spherical virialized halos) by an unrelated artificial overestimate.

6.5. Lagrangian perturbation theories

The Lagrangian re-interpretation of the halo model developed in this paper, leading to the expression (26) for the 2-halo contribution, would most naturally lead us to consider Lagrangian perturbation theories. It is possible to substitute such schemes into Eq. (26), which allows us to integrate in a self-consistent fashion over the q-dependent prefactor F_2H(q). Unfortunately, this does not yield very successful results while leading to somewhat more complex expressions than those encountered in Eulerian frameworks. Therefore, in this section we illustrate the problems faced by current Lagrangian schemes by using the simpler expression (27), where the term F_2H(q) has been factorized as the prefactor F_2H(1/k). Thus, as for the Eulerian perturbation theories used in the previous sections, we only need the prediction for the resummed power P_pert(k) to compute the 2-halo contribution. This simplifies matters and only makes a negligible change on large scales where F_2H ≃ 1.

We consider two schemes that are inspired by a Lagrangian framework. A first approach is the one presented in Matsubara (2008), which for our purposes amounts to expand P(k) after factorization of a Gaussian damping term ${\mathrm{e}}^{\mathrm{-}{\mathit{k}}^{\mathrm{2}}{\mathit{\sigma }}_{\mathit{v}}^{\mathrm{2}}}$ (see also Crocce & Scoccimarro 2006b; Valageas 2010a), where σ_v is the rms one-dimensional linear displacement (or the rms one-dimensional linear velocity, up to a time-dependent factor) given by ${\mathit{\sigma }}_{\mathit{v}}^{\mathrm{2}}\mathrm{=}\frac{\mathrm{4}\mathit{\pi }}{\mathrm{3}}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{\infty }}\mathrm{d}\mathit{k}\hspace{0.17em}{\mathit{P}}_{\mathrm{L}}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathit{.}$(52)Thus, at 1-loop order the density power spectrum reads as (Matsubara 2008) ${\mathit{P}}_{\mathrm{Mat}}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathrm{=}{\mathrm{e}}^{\mathrm{-}{\mathit{k}}^{\mathrm{2}}{\mathit{\sigma }}_{\mathit{v}}^{\mathrm{2}}}\mathrm{(}{\mathit{P}}_{\mathrm{L}}\mathrm{+}{\mathit{P}}_{\mathrm{22}}\mathrm{+}{\mathit{P}}_{\mathrm{31}}\mathrm{+}{\mathit{k}}^{\mathrm{2}}{\mathit{\sigma }}_{\mathit{v}}^{\mathrm{2}}{\mathit{P}}_{\mathrm{L}}\mathrm{)}\mathit{,}$(53)where P₂₂ and P₃₁ are the 1-loop terms obtained in the standard Eulerian perturbation theory, as in Eq. (48). Of course, Eq. (53) could be obtained at once by requiring consistency with Eq. (48) at order ${\mathit{P}}_{\mathrm{L}}^{\mathrm{2}}$. As discussed in Valageas (2007b, 2008, 2010a), and as can be seen from the computation of Lagrangian response functions (Bernardeau & Valageas 2008, 2010a), the prefactor ${\mathrm{e}}^{\mathrm{-}{\mathit{k}}^{\mathrm{2}}{\mathit{\sigma }}_{\mathit{v}}^{\mathrm{2}}}$ is associated with the coherent displacement of density structures by the long wavelengths of the velocity field and should only appear in different-time quantities, under the form ${\mathrm{e}}^{\mathrm{-}\mathrm{(}{\mathit{D}}_{\mathrm{1}}\mathrm{-}{\mathit{D}}_{\mathrm{2}}{\mathrm{)}}^{\mathrm{2}}{\mathit{k}}^{\mathrm{2}}{\mathit{\sigma }}_{\mathit{v}\mathrm{0}}^{\mathrm{2}}\mathit{/}\mathrm{2}}$, where D_i is the linear growth factor at redshift z_i and σ_v0 is the linear dispersion today, where D(z = 0) = 1. Thus, looking for an expansion under the form (53) leads to an artificially strong Gaussian damping at high k, which has no physical meaning. In order to remedy to this problem, which can also be seen as a consequence of the truncation to a finite number of terms between the parenthesis in Eq. (53) (once we insist on looking for an expansion of this form), we note that the nonlinear power spectrum generated by the Zeldovich dynamics (Zeldovich 1970) can also be written under this form, ${\mathit{P}}_{\mathrm{Z}}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathrm{=}{\mathrm{e}}^{\mathrm{-}{\mathit{k}}^{\mathrm{2}}{\mathit{\sigma }}_{\mathit{v}}^{\mathrm{2}}}\sum _{\mathit{n}\mathrm{=}\mathrm{1}}^{\mathrm{\infty }}{\mathit{P}}_{\mathrm{Z}}^{\mathrm{\left(}\mathit{n}\mathrm{\right)}}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathit{,}$(54)where each term ${\mathit{P}}_{\mathrm{Z}}^{\mathrm{\left(}\mathit{n}\mathrm{\right)}}\mathrm{\left(}\mathit{k}\mathrm{\right)}$ scales as (P_L)ⁿ. Moreover, both the full nonlinear expression of the Zeldovich power spectrum, P_Z(k), (i.e. the resummation of the series (54)), and the expression of the terms (P_L)ⁿ of all orders, are explicitly known (Crocce & Scoccimarro 2006b; Valageas 2007b, 2010a). Then, a possible recipe is to use for the nonlinear power spectrum (more precisely the perturbative term P_pert) the “modified Zeldovich” expression $\begin{array}{ccc}{\mathit{P}}_{\mathrm{mZ}}\mathrm{\left(}\mathit{k}\mathrm{\right)}& \mathrm{=}& {\mathrm{e}}^{\mathrm{-}{\mathit{k}}^{\mathrm{2}}{\mathit{\sigma }}_{\mathit{v}}^{\mathrm{2}}}\left({\mathit{P}}_{\mathrm{L}}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathrm{+}{\mathit{P}}^{\mathrm{\left(}\mathrm{2}\mathrm{\right)}}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathrm{+}\sum _{\mathit{n}\mathrm{=}\mathrm{3}}^{\mathrm{\infty }}{\mathit{P}}_{\mathrm{Z}}^{\mathrm{\left(}\mathit{n}\mathrm{\right)}}\mathrm{\left(}\mathit{k}\mathrm{\right)}\right)\\ & \mathrm{=}& {\mathit{P}}_{\mathrm{Z}}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathrm{+}{\mathrm{e}}^{\mathrm{-}{\mathit{k}}^{\mathrm{2}}{\mathit{\sigma }}_{\mathit{v}}^{\mathrm{2}}}\mathrm{(}{\mathit{P}}^{\mathrm{\left(}\mathrm{2}\mathrm{\right)}}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathrm{-}{\mathit{P}}_{\mathrm{Z}}^{\mathrm{\left(}\mathrm{2}\mathrm{\right)}}\mathrm{(}\mathit{k}{\mathrm{)}}^{\mathrm{)}}\mathit{,}\end{array}$where P⁽ⁿ⁾(k) is the term of order (P_L)ⁿ in the expansion of the form (54) for the gravitational power spectrum, whence ${\mathit{P}}^{\mathrm{\left(}\mathrm{2}\mathrm{\right)}}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathrm{=}{\mathit{P}}_{\mathrm{22}}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathrm{+}{\mathit{P}}_{\mathrm{31}}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathrm{+}{\mathit{k}}^{\mathrm{2}}{\mathit{\sigma }}_{\mathit{v}}^{\mathrm{2}}{\mathit{P}}_{\mathrm{L}}\mathrm{\left(}\mathit{k}\mathrm{\right)}$ from Eq. (53). In other words, we complete the expansion (53) by substituting the terms obtained for the Zeldovich power spectrum for all higher orders, which allows us to perform an explicit resummation as in the second line of Eq. (56). This is clearly a systematic procedure, that can be pushed up to any order, and agrees with standard perturbation theory up to the order of the substitution. Since the resummed Zeldovich power spectrum no longer shows the spurious Gaussian decay of Eq. (53) but a more physical power-law decay (Taylor & Hamilton 1996; Valageas 2007b, 2010a), one may hope that this could improve the agreement with results from numerical simulations.

We show in Fig. 18 the results we obtain using either Eq. (53) or (56), at redshifts z = 0.35, 1, and 3. As pointed out above, we can see that the Gaussian damping term of Eq. (53) leads to a fast decay for the nonlinear density power spectrum, which cannot be compensated by the 1-halo term. The “modified Zeldovich” expansion (56), which only implies a slower power-law decay at high k, fares better as the nonlinear power spectrum keeps growing on these scales. However, this is not sufficient to obtain a satisfactory match with the N-body results. Thus, the comparison with Fig. 10 shows that these Lagrangian-based expansions are not competitive with Eulerian schemes, and we do not investigate further these approaches in this paper. We note however that this is rather unfortunate, since in principles Lagrangian perturbation theories would seem better suited to compute the 2-halo term (26), and more importantly they should be better suited to computations in redshift space (Matsubara 2008). Therefore, it remains of interest to investigate whether other resummation schemes could be developed within the Lagrangian framework, but Fig. 18 suggests that some important new ingredients are required to make significant progress.

7. Real-space two-point correlation

Finally, we consider in this section the real-space two-point correlation function, ξ( | x₂ − x₁ | ) =  ⟨ δ(x₁)δ(x₂) ⟩ , that is predicted by our model. It is related to the Fourier-space power spectrum through $\begin{array}{ccc}\mathit{\xi }\mathrm{\left(}\mathit{x}\mathrm{\right)}& \mathrm{=}& \mathrm{4}\mathit{\pi }{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{\infty }}\mathrm{d}\mathit{k}\hspace{0.17em}{\mathit{k}}^{\mathrm{2}}\mathit{P}\mathrm{\left(}\mathit{k}\mathrm{\right)}\hspace{0.17em}\frac{\sin \mathrm{\left(}\mathit{kx}\mathrm{\right)}}{\mathit{kx}}\\ & \mathrm{=}& {\mathrm{\int }}_{\mathrm{0}}^{\mathrm{\infty }}\frac{\mathrm{d}\mathit{k}}{\mathit{k}}\hspace{0.17em}{\mathrm{\Delta }}^{\mathrm{2}}\mathrm{\left(}\mathit{k}\mathrm{\right)}\hspace{0.17em}\frac{\sin \mathrm{\left(}\mathit{kx}\mathrm{\right)}}{\mathit{kx}}\mathit{,}\end{array}$so that we obtain ξ(x) by integration over k of the nonlinear power spectrum (22). In particular, from the decomposition (22) we can write the nonlinear correlation as the sum of 1-halo and 2-halo contributions, $\mathit{\xi }\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{=}{\mathit{\xi }}_{\mathrm{1}\mathrm{H}}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{+}{\mathit{\xi }}_{\mathrm{2}\mathrm{H}}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathit{,}$(59)which are the Fourier transforms of the power spectra (23) and (24).

We show our results at redshifts z = 0.35, 1, and 3, in Figs. 19 and 20. Here we consider our fiducial model presented in Sects. 3 and 5, with the 2-halo perturbative term P_pert(k) given by the direct steepest-descent scheme detailed in Appendix A, and the 1-halo contribution obtained for halos defined by a density contrast δ = 200 and the mass-concentration relation (42).

In agreement with our results for the power spectrum shown in Fig. 9, we can see in Fig. 19 that our model provides a good match to numerical simulations on both small and large scales. As is well known, on large scales the two-point correlation is governed by the 2-halo contribution ξ_2H, whereas on small scales it is governed by the 1-halo contribution. As for the power spectrum, the transition scales where the 2-halo and 1-halo contributions are of the same order are more difficult to reproduce, and we again underestimate the power (here the real-space correlation) in this range. This is a direct consequence of the same feature observed in Fourier space in Fig. 9. In particular, we also find that the range and amplitude of this discrepancy is larger at z = 3 than at z = 0.35. Again, improving theoretical predictions on this domain would require pushing perturbation schemes to higher orders and/or improving the description of halo outer shells.

We focus in Fig. 20 on the large scales around the baryon acoustic peak. We can check that our model agrees very well with the N-body results and reproduces the well-known damping of the BAO peak by the nonlinearities of the dynamics (mode-coupling effect). This shows that for the purpose of BAO studies, the accuracy of our theoretical model is certainly sufficient. In particular, on these large scales, which are far from the transition region seen in Fig. 19, the 1-halo contribution is completely negligible and the two-point correlation is set by the perturbative term P_pert of the 2-halo contribution. Thus, it is given by a systematic perturbation theory, and we can check that the simple 1-loop direct steepest-descent scheme detailed in Appendix A already provides quite satisfactory results. For practical purposes, the main improvement that remains to be achieved is to extend these results to redshift space, but this is beyond the scope of this paper and we leave this point for future studies.

We compare in Fig. 21 the correlation functions obtained with halos defined by a density threshold δ = 200 or δ = 50. This is the real-space counterpart of Fig. 11. In agreement with Fig. 11, we can see that extending the truncation radius of the halo profiles increases somewhat the power in the transition range between the 2-halo and 1-halo contributions. However, this is not sufficient to obtain a perfect match to the numerical simulations, especially at z = 3. Again, this suggests that part of the discrepancy is due to the incomplete resummation of higher order terms of the perturbative contribution P_pert(k).

We do not plot the counterpart of Fig. 20, because on these large scales the 1-halo contribution is very small and we have checked that the curves obtained with halos defined by δ = 200 or δ = 50 cannot be distinguished in the figure: the two curves coincide within the thickness of the curves of Fig. 20. This also shows that the predictions on these large scales are quite robust, since they are independent of the underlying halo model and are given by systematic perturbation theory. The effect of the 1-halo counterterm is also very small so that we do not plot a specific figure to compare the curves obtained with and without this counterterm. Note that the difference between various curves in Figs. 12 and 13 was amplified by the fact that we divided by a “no-wiggle” linear power spectrum P_Ls(k).

8. Typical accuracy of combined models

In this last section before conclusion, as a summary of what can be achieved from simple analytical models that combine perturbation theories with halo models, we estimate the accuracy reached by our fiducial model. Thus, we plot in Fig. 22 the relative difference between our model, described in Sects. 3 and 5, and the results from the N-body simulations presented in Sect. 4, $\frac{\mathrm{\Delta }\mathit{P}}{\mathit{P}}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathrm{=}\frac{\left|{\mathit{P}}_{\mathrm{model}}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathrm{-}{\mathit{P}}_{\mathit{N}\mathrm{-}\mathrm{body}}\mathrm{\left(}\mathit{k}\mathrm{\right)}\right|}{{\mathit{P}}_{\mathit{N}\mathrm{-}\mathrm{body}}\mathrm{\left(}\mathit{k}\mathrm{\right)}}\mathrm{·}$(60)We also show the statistical error due to the finite number of modes in the simulation box, given by Eq. (45), and the shot-noise error given by $\mathrm{\Delta }{\mathit{P}}_{\mathrm{shot}\mathrm{-}\mathrm{noise}}\mathrm{=}{\mathit{L}}_{\mathrm{box}}^{\mathrm{3}}\mathit{/}\mathit{N}$. Since the number of modes within a bin of fixed size Δk scales as k²Δk, the relative statistical error decreases at higher k as 1/k, see Eq. (45). The sudden jumps are due to the folding procedure, see Sect. 4.2 and Fig. 7. On the other hand, the relative shot noise grows as 1/P(k) and dominates at high k. It seems that the simple approximation $\mathrm{\Delta }{\mathit{P}}_{\mathrm{shot}\mathrm{-}\mathrm{noise}}\mathrm{=}{\mathit{L}}_{\mathrm{box}}^{\mathrm{3}}\mathit{/}\mathit{N}$ overestimates somewhat the inaccuracy of the simulations, however we do not look for a better estimate here, as this is already sufficient to understand the high-k part of Fig. 22 and to mark the wavenumber where shot noise becomes dominant.

We can see that on the large scales associated with the BAO oscillations, 0.03 < k < 0.3   h   Mpc^-1, the analytical approach is competitive with the numerical simulations. It typically gives an accuracy of 1%, which however worsens somewhat at lower redshift to reach 2% at z = 0.35. However, a part of ΔP/P shown in Fig. 22 arises from the inaccuracies of the simulations themselves, so that the true value of relative error between the analytical model and the exact nonlinear power spectrum may be slightly smaller. Of course, on even larger scales, k < 0.03   h   Mpc^-1, the analytical model becomes exact as it agrees with linear theory and the relative difference ΔP/P is solely due to the inaccuracies of the simulation results (as can be checked by comparison with the relative statistical error of the simulations). In the 1-halo region, 3 < k < 100   h   Mpc^-1, we reach an accuracy on the order of 10%, and even better (4%) in the region where 10 < Δ²(k) < 100. Moreover, as we have seen in Sect. 6, the predictions of our model are rather robust in this range, since the dependence on the details of the underlying halo model is quite weak. At very high k the numerical simulations are strongly affected by shot noise, so that the curve ΔP/P mostly measures this source of error and the accuracy of the model may be better than what Fig. 22 suggests.

In agreement with the discussion in Sect. 5, a salient feature in Fig. 22 is the peak in the transition range between the 2-halo and 1-halo contributions, at k ~ 1   h   Mpc^-1. This is due to the underestimation noticed in Fig. 9, which can reach 30% around the peak. This feature shifts to higher k and somewhat broadens at higher redshift, in agreement with the discussion in Sect. 5.

Next, we show in Fig. 23 the relative accuracy of our model and our simulations for the real-space two-point correlation function. As in the Fourier-space Fig. 22, we plot the relative difference between our model, presented in Sect. 7 (for halos defined by δ = 200), and our simulations, $\frac{\mathrm{\Delta }\mathit{\xi }}{\mathit{\xi }}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{=}\frac{\left|{\mathit{\xi }}_{\mathrm{model}}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{-}{\mathit{\xi }}_{\mathit{N}\mathrm{-}\mathrm{body}}\mathrm{\left(}\mathit{x}\mathrm{\right)}\right|}{{\mathit{\xi }}_{\mathit{N}\mathrm{-}\mathrm{body}}\mathrm{\left(}\mathit{x}\mathrm{\right)}}\mathrm{·}$(61)We also show the statistical error, given by Eq. (46), and the shot-noise error given by $\mathrm{\Delta }{\mathit{\xi }}_{\mathrm{shot}\mathrm{-}\mathrm{noise}}\mathrm{=}\mathrm{[}{\mathit{L}}_{\mathrm{box}}^{\mathrm{3}}\mathit{/}\mathrm{(}\mathrm{4}\mathit{\pi }{\mathit{x}}^{\mathrm{3}}\mathit{/}\mathrm{3}\mathrm{\left)}\mathrm{\right]}\mathit{/}\mathit{N}$. As for the power spectrum, the relative statistical error grows on large scales, as  ∝ x, because of the smaller number of modes, and the sudden jumps are due to the folding procedure, see Sect. 4.2 and Fig. 7. On the other hand, the relative shot noise grows as 1/(x³ξ(x)) on small scales and dominates at low x.

We can check that Fig. 23 is consistent with Fig. 22. On large scales, x > 10   h^-1 Mpc, the analytical approach is again competitive with the numerical simulations and it reaches an accuracy of 1%. On very large scales the analytical model converges to the linear two-point correlation and becomes exact, and the relative difference Δξ/ξ is solely due to the inaccuracies of the simulations. On small scales dominated by the 1-halo contribution, x < 0.1   h^-1 Mpc, we reach an accuracy of 10%, and even better at z ≥ 1. On very small scales the difference Δξ/ξ mostly measures the level of shot noise in the simulations, so that the accuracy of the model could be better than shown in the figure.

Again, the most difficult region to reproduce is the transition between the 2-halo and 1-halo contributions, x ~ 1   h^-1 Mpc. The underestimation seen in Fig. 19 leads to a peak for the relative error of the analytical model that can reach about 30% and is somewhat greater at higher redshift. Of course, the boundaries of these various domains shift to smaller scales at higher redshifts.

Thus, it appears that analytical models such as the one studied in this paper can already provide reasonably good estimates for the density power spectrum and two-point correlation function over a large range of scales. We have not scanned all possible models in this paper (although we have investigated several variants in Sect. 6), but we expect similar approaches, based on a combination of perturbation theory and halo model, to yield similar levels of accuracy, for current resummation schemes and halo models. At low and high k, and at small and large x, the level of accuracy is already rather satisfactory, and may be improved in a straightforward and systematic manner by pushing perturbation theories to higher orders and measuring halo profiles down to smaller scales in simulations (and possibly taking into account additional effects such as substructures and baryon impact on gravitational collapse). An advantage of such models, especially on large scales, it to provide robust predictions that can be applied to a variety of initial conditions or cosmological parameters. This is obvious on the large scales described by systematic perturbation theories, but this remains true to some extent at high k and small x as the power spectrum is expressed in terms of more elementary quantities that satisfy physical constraints (such as the normalization or scaling of the mass function, and reasonable halo profiles). This ensures at the very least that the power spectra and two-point correlations obtained in this fashion make sense and scale in the appropriate fashion, even for (reasonable) cases where they have not been tested.

What may prove to be the most difficult regime for systematic improvement is the transition range. Although we have not investigated such an approach here, it is certainly possible to build fitting formulas that improve the agreement with simulations on this range. Thus, one may choose some interpolation formula between the low-k and high-k domains and obtain a good match by tuning some free parameters that govern the curvature of the interpolation curve in this regime. However, one must take some care not to spoil the good agreement on BAO scales and to keep a percent accuracy there. On the other hand, if we insist on avoiding such a direct interpolation procedure, so as to derive the power spectrum from more elementary quantities (such as perturbative expansions and halo profiles) without further modifications, it is not clear how much progress can be made on this range. Indeed, the description of the density field in terms of halos is by itself an approximate picture, especially for the moderate-density regions of space associated with outer shells or filaments that have already experienced shell crossing but have not fully virialized and are not enclosed within well-identified objects. This may prevent a perfect matching on the transition scales, unless it is somehow enforced a posteriori.

9. Conclusion

In this article we have explained how to build unified models to predict the matter density power spectrum, and the real-space two-point correlation function, by combining perturbation theories with halo models. First, we have shown how a Lagrangian point of view allows us to re-interpret the decomposition of the power spectrum into 2-halo and 1-halo contributions. This provides a convenient route to this splitting, which automatically ensures the conservation of mass and gives an explicit relationship with the dynamics of the system, through the direct Lagrangian map q → x. We have also emphasized the relationship between this decomposition and the separation into perturbative and non-perturbative terms. This explicitly shows how one can achieve consistency with perturbation theory through the 2-halo term, as the 1-halo term is fully non-perturbative and does not contribute to any order of perturbation theory.

Next, we have shown that the 1-halo contribution contains a counterterm, which was missed in previous studies, that gives the expected low-k tail P_1H(k) ∝ k² associated with the conservation of matter. This ensures that the 1-halo term becomes subdominant on large scales, as linear CDM power spectra typically behave as P_L(k) ∝ k at low k. On the other hand, we have explained why the k⁴-tail associated with momentum conservation is not recovered because of the approximations that we use (such as instantaneous virialization and loss of memory), but this is not a problem for practical purposes. If needed, in future works aiming at greater accuracy that include higher orders of perturbation theory, it may be possible to enforce the k⁴-tail by modifying this counterterm. Then, we have pointed out that the 2-halo contribution, and more precisely its term P_pert(k) that should be consistent with perturbation theory, cannot be computed by truncating standard perturbation theory at a finite order. Indeed, higher order terms of this perturbative expansion grow increasingly fast at high k (with large cancellations between various orders), so that an abrupt truncation at a given order leads to a 2-halo contribution that is non-negligible at high k (and would even be dominant on small scales if we go to two loops or higher orders), which is not physical and prevents a good agreement with numerical simulations. The remedy is to use resummation schemes that agree with standard perturbation theory, up to the required order, while remaining well-behaved at high k (typically close to or smaller than the linear power spectrum). In addition to the superior accuracy of such schemes on large scales, this is a second important motivation to develop such perturbative approaches. We have described in particular a simple implementation, based on the direct steepest-descent resummation at 1-loop order, which satisfies this property while being fast to compute, thanks to its factorized form. In this fashion, both the 1-halo and 2-halo contributions remain well-behaved beyond the domain where they dominate. This allows us to build a meaningful model from their combination, which can describe all scales and regimes.

We have compared a simple implementation of this model to N-body simulations. We can reach an accuracy of 1% on quasi-linear scales, associated with baryon acoustic oscillations, and 10% on highly nonlinear scales dominated by the 1-halo contribution. The agreement on large scale is actually better than the one obtained with standard 1-loop perturbation theory. Since it is based on systematic perturbation theory and contains no free parameter (except through the prefactor F_2H, which however is very close to unity on these scales), this gives a robust and reliable prediction that also improves over simple fitting formulas. On the other hand, the good agreement on highly nonlinear scales, Δ²(k) > 200, depends on the properties of the halos included in the underlying halo model, and more specifically on the mass-concentration relation c(M). We have given a simple model that provides a reasonable match to the N-body results, but this regime is clearly less reliable than the quasi-linear regime. In particular, to apply this framework to very high redshifts one should use a better constrained prescription c(M) or a physically motivated model. Nevertheless, in this high-k regime where 200 < Δ²(k) < 1000, we find that using various formulas for c(M) that have been proposed in the literature already improves over simple fitting formulas for P(k), in spite of the scatter of their predictions. Indeed, on these scales the dependence on the details of the model (the low-mass tail of c(M)) is not too large yet. Moreover, on somewhat less nonlinear scales where 10 < Δ² < 200, the predictions appear even more robust and also provide a good match to simulations.

The transition range where the 2-halo and 1-halo contributions are of the same order, roughly where 1 < Δ² < 10, proves to be the most difficult to reproduce. In agreement with previous works (Giocoli et al. 2010), we find that the model underestimates the power spectrum on these scales, and the inaccuracy is somewhat larger at higher z. Some of the discrepancy may arise from the truncation of halo profiles at a density contrast δ = 200, and we have shown that going to larger radii, defined by δ = 50, gives a slightly better match to the simulations, but this is not sufficient. Another route to improve the model would be to truncate the resummation scheme (involved in the 2-halo contribution) at a higher order, that is to include exactly all perturbative terms up to a higher loop order, while higher orders are only partially resummed. As seen in previous works (Valageas 2010a), the potential for improvement is greatest at higher z for CDM power spectra, since the range where perturbative expansions are relevant is somewhat broader and one can push to higher orders before non-perturbative terms become dominant. This would explain why the underestimate of the power spectrum appears to be more significant at z = 3 than at z = 0.35.

However, it is likely that even a complete resummation of the perturbative term would not provide a perfect match, because of the approximations involved by the splitting into 2-halo and 1-halo contributions. Indeed, even though the decomposition (22), with Eqs. (23), (24), can be made exact if we have a well-defined criterion to define halos, the computation of the 2-halo term involves some approximation as we replace the conditional average  ⟨ .. ⟩ _2H by the perturbative average  ⟨ .. ⟩ _pert. Going beyond this approximation requires going beyond perturbation theory (which does not know about halos). On the other hand, the 1-halo contribution also involves approximations that are difficult to improve, such as the instantaneous virialization, that is, the fact that particles are located at random in a halo, independently of their initial location. Nevertheless, this is a range of scales where there is still room for improvement.

Next, we have investigated the dependence of such predictions on the details of the model. We have shown that it is important to take into account the 1-halo counterterm in order to get a good accuracy on large scales, k ~ 0.1 to 0.3   h   Mpc^-1. Otherwise, the 1-halo contribution does not decrease fast enough at low k and spoils the good match with simulations (and would even become dominant on very large scales). On the perturbative side, we have shown that standard perturbation theory cannot be used to build unified model, because it yields higher order terms that grow increasingly fast at high k and prevent a good match to numerical simulations in the nonlinear regime. This cannot be compensated by details of the underlying model and would require at best the addition of some ad-hoc damping prefactor. Another perturbative scheme than the 1-loop direct steepest-descent resummation we have focused on in this paper, based on a more realistic “decaying response” function, gives similar results to the former, as it also provides a well-behaved 2-halo contribution at high k. It even improves somewhat the match with simulations on transition scales, but this is likely to be a coincidence as this comes at the price of an overestimate on the larger scales probed by BAO. Next, we have shown that current Lagrangian resummation schemes, which a priori could have been expected to be best suited to our purposes, actually fare much worse than Eulerian resummation schemes. Indeed, although they yield a well-behaved perturbative term at high k, the associated damping is so large that it leads to a significant underestimate for the power spectrum already on BAO scales, k ~ 0.1 to 0.3   h   Mpc^-1. However, if it is possible to devise new Lagrangian schemes that do not suffer from this effect they may provide a promising route to improve unified models for the power spectrum.

Finally, we have shown that our model also provides reliable predictions for the density two-point correlation function. As for the power spectrum, it describes rather well both small and large scales, but underestimates somewhat the correlation in the transition range between the 2-halo and 1-halo contributions. On very large scales, around the BAO peak, we obtain a very good agreement with numerical simulations, so that the simple 1-loop resummation scheme used in this work is probably sufficient to describe this feature.

Thus, we have shown that it is possible to build successful unified models for the density power spectrum and two-point correlation, that combine systematic and reliable perturbation theories on large scales with phenomenological halo models on small scales. Two issues where there is still room for improvement are the underestimate of the power spectrum and correlation function on transition scales and the building of accurate Lagrangian resummation schemes (but this is a more theoretical goal since for practical purposes Eulerian schemes are sufficient). A third issue is the modeling of the underlying virialized halos, which plays a key role at very high k. There, many points can be improved, such as taking into account substructures, deviations from spherical profiles, or the effect of baryons on the halo density profile. In particular, one should be cautious when using such models beyond the range where they have been tested (Δ²(k) < 1000 and z ≤ 3). There, the main uncertainty comes from the low-mass tail of the mass-concentration relation c(M), and it is probably safer to use physically motivated models for c(M) (even though they may not be very accurate), or at least to check that the formula used for c(M) is still reasonable.

The model presented in this paper could be extended to higher order statistics, such as the bispectrum or three-point correlation function. It should also be possible to handle the case of non-Gaussian initial conditions, at least as long as the halo density profiles are not strongly modified (or this would require a separate study to evaluate this effect). Finally, to allow a direct comparison with galaxy surveys, it would be desirable to extend the model to redshift space. This is likely to be a more involved task, that we leave for future works.

---

Acknowledgments

T.N. is supported by a Grant-in-Aid for Japan Society for the Promotion of Science (JSPS) Fellows and by World Premier International Research Center Initiative (WPI Initiative), MEXT, Japan. Numerical computations for the present work have been carried out in part on Cray XT4 at Center for Computational Astrophysics, CfCA, of National Astronomical Observatory of Japan, and in part under the Interdisciplinary Computational Science Program in Center for Computational Sciences, University of Tsukuba.

References

Appendix A: Computation of P_pert(k) with the “direct steepest-descent method”

We briefly recall here how the perturbative term P_pert(k) of the 2-halo contribution (27) is computed within the “direct steepest-descent method”, see also Valageas (2007a, 2008) for details. First, to simultaneously describe the density and velocity fields it is convenient to define the two-component vector ψ(k,η) as (Crocce & Scoccimarro 2006b), $\mathit{\psi }\mathrm{\left(}{k}\mathit{,\eta }\mathrm{\right)}\mathrm{=}\left(\begin{array}{c}\\ \mathit{\delta }\mathrm{\left(}{k}\mathit{,\eta }\mathrm{\right)}\\ \mathrm{-}\mathit{\theta }\mathrm{\left(}{k}\mathit{,\eta }\mathrm{\right)}\mathit{/}\left(\mathit{ȧ}\mathit{f}\right)\end{array}\right)\mathit{,}$(A.1)where θ is the divergence of the velocity field, θ = ∇·v, η = lnD(z) where D(z) is the linear growth factor (normalized by D(0) = 1 today), ȧ = da/dt where a(t) is the scale factor, and f = dlnD/dlna. Then, introducing the coordinate x = (k,η,i), where i = 1 or 2 is the discrete index of the two-component vectors, the equations of motion write as (Valageas 2007a) $\mathrm{𝒪}\left(\mathit{x,}{\mathit{x}}^{\mathrm{\prime }}\right)\mathrm{·}\mathit{\psi }\left({\mathit{x}}^{\mathrm{\prime }}\right)\mathrm{=}{\mathit{K}}_{\mathit{s}}\mathrm{\left(}\mathit{x}\mathrm{;}{\mathit{x}}_{\mathrm{1}}\mathit{,}{\mathit{x}}_{\mathrm{2}}\mathrm{\right)}\mathrm{·}\mathit{\psi }\mathrm{\left(}{\mathit{x}}_{\mathrm{1}}\mathrm{\right)}\mathit{\psi }\mathrm{\left(}{\mathit{x}}_{\mathrm{2}}\mathrm{\right)}\mathit{,}$(A.2)where the matrix reads as $\begin{array}{ccc}\mathrm{𝒪}\left(\mathit{x,}{\mathit{x}}^{\mathrm{\prime }}\right)& \mathrm{=}& {\mathit{\delta }}_{\mathrm{D}}\left({k}\mathrm{-}{k}\mathrm{\prime }\right)\hspace{0.17em}{\mathit{\delta }}_{\mathrm{D}}\left(\mathit{\eta }\mathrm{-}{\mathit{\eta }}^{\mathrm{\prime }}\right)\\ & & \mathrm{\times }\left(\begin{array}{c}\\ & \\ & \end{array}\right)\mathit{,}\end{array}$(A.3)whereas the symmetric vertex K_s(x;x₁,x₂) = K_s(x;x₂,x₁) writes as $\begin{array}{ccc}{\mathit{K}}_{\mathit{s}}\mathrm{(}\mathit{x}\mathrm{;}{\mathit{x}}_{\mathrm{1}}\mathit{,}{\mathit{x}}_{\mathrm{2}}\mathrm{)}& \mathrm{=}& {\mathit{\delta }}_{\mathrm{D}}\left({k}\mathrm{1}\mathrm{+}{k}\mathrm{2}\mathrm{-}{k}\right)\hspace{0.17em}{\mathit{\delta }}_{\mathrm{D}}\left({\mathit{\eta }}_{\mathrm{1}}\mathrm{-}\mathit{\eta }\right)\hspace{0.17em}{\mathit{\delta }}_{\mathrm{D}}\left({\mathit{\eta }}_{\mathrm{2}}\mathrm{-}\mathit{\eta }\right)\\ & & \mathrm{\times }{\mathit{\gamma }}_{\mathit{i}\mathrm{;}{\mathit{i}}_{\mathrm{1}}\mathit{,}{\mathit{i}}_{\mathrm{2}}}^{\mathit{s}}\left({k}\mathrm{1}\mathit{,}{k}\mathrm{2}\right)\mathit{,}\end{array}$(A.4)with $\begin{array}{ccc}{\mathit{\gamma }}_{\mathrm{1}\mathrm{;}\mathrm{1}\mathit{,}\mathrm{2}}^{\mathit{s}}\left({k}\mathrm{1}\mathit{,}{k}\mathrm{2}\right)& \mathrm{=}& \frac{\left({k}\mathrm{1}\mathrm{+}{k}\mathrm{2}\right)\mathrm{·}{k}\mathrm{2}}{\mathrm{2}{\mathit{k}}_{\mathrm{2}}^{\mathrm{2}}}\mathit{,}\\ {\mathit{\gamma }}_{\mathrm{1}\mathrm{;}\mathrm{2}\mathit{,}\mathrm{1}}^{\mathit{s}}\left({k}\mathrm{1}\mathit{,}{k}\mathrm{2}\right)& \mathrm{=}& \frac{\left({k}\mathrm{1}\mathrm{+}{k}\mathrm{2}\right)\mathrm{·}{k}\mathrm{1}}{\mathrm{2}{\mathit{k}}_{\mathrm{1}}^{\mathrm{2}}}\mathit{,}\\ {\mathit{\gamma }}_{\mathrm{2}\mathrm{;}\mathrm{2}\mathit{,}\mathrm{2}}^{\mathit{s}}\left({k}\mathrm{1}\mathit{,}{k}\mathrm{2}\right)& \mathrm{=}& \frac{{\left|{k}\mathrm{1}\mathrm{+}{k}\mathrm{2}\right|}^{\mathrm{2}}\left({k}\mathrm{1}\mathrm{·}{k}\mathrm{2}\right)}{\mathrm{2}{\mathit{k}}_{\mathrm{1}}^{\mathrm{2}}{\mathit{k}}_{\mathrm{2}}^{\mathrm{2}}}\mathit{,}\end{array}$and zero otherwise (Crocce & Scoccimarro 2006b). In addition to the usual two-point correlation function, $\mathit{C}\mathrm{\left(}{\mathit{x}}_{\mathrm{1}}\mathit{,}{\mathit{x}}_{\mathrm{2}}\mathrm{\right)}\mathrm{=}⟨\mathit{\psi }\mathrm{\left(}{\mathit{x}}_{\mathrm{1}}\mathrm{\right)}\mathit{\psi }\mathrm{\left(}{\mathit{x}}_{\mathrm{2}}\mathrm{\right)}⟩\mathit{,}$(A.8)resummation schemes often involve the response function (or propagator), defined as the functional derivative of the nonlinear field ψ with respect to an infinitesimal “noise” ζ that would be added to the right hand side of the equation of motion (A.2), $\mathit{R}\mathrm{\left(}{\mathit{x}}_{\mathrm{1}}\mathit{,}{\mathit{x}}_{\mathrm{2}}\mathrm{\right)}\mathrm{=}{⟨\frac{\mathrm{𝒟}\mathit{\psi }\mathrm{\left(}{\mathit{x}}_{\mathrm{1}}\mathrm{\right)}}{\mathrm{𝒟}\mathit{\zeta }\mathrm{\left(}{\mathit{x}}_{\mathrm{2}}\mathrm{\right)}}⟩}_{\mathit{\zeta }\mathrm{=}\mathrm{0}}\mathit{.}$(A.9)As described in Valageas (2007a, 2008), the nonlinear correlation and response functions obey the exact Schwinger-Dyson equations $\begin{array}{ccc}\mathrm{𝒪}\mathrm{\left(}\mathit{x,z}\mathrm{\right)}\mathrm{·}\mathit{C}\mathrm{\left(}\mathit{z,y}\mathrm{\right)}& \mathrm{=}& \mathrm{\Sigma }\mathrm{\left(}\mathit{x,z}\mathrm{\right)}\mathrm{·}\mathit{C}\mathrm{\left(}\mathit{z,y}\mathrm{\right)}\mathrm{+}\mathrm{\Pi }\mathrm{\left(}\mathit{x,z}\mathrm{\right)}\mathrm{·}{\mathit{R}}^{\mathit{T}}\mathrm{\left(}\mathit{z,y}\mathrm{\right)}\mathit{,}\\ \mathrm{𝒪}\mathrm{\left(}\mathit{x,z}\mathrm{\right)}\mathrm{·}\mathit{R}\mathrm{\left(}\mathit{z,y}\mathrm{\right)}& \mathrm{=}& {\mathit{\delta }}_{\mathrm{D}}\mathrm{(}\mathit{x}\mathrm{-}\mathit{y}\mathrm{)}\mathrm{+}\mathrm{\Sigma }\mathrm{\left(}\mathit{x,z}\mathrm{\right)}\mathrm{·}\mathit{R}\mathrm{\left(}\mathit{z,y}\mathrm{\right)}\mathit{,}\end{array}$where Σ(x,y) and Π(x,y) are “self-energy” terms, that can be written as diagrammatic expansions over the linear two-point functions C_L and R_L. Equation (A.10) can be integrated as $\mathit{C}\mathrm{=}\mathit{R}\mathrm{\times }{\mathit{C}}_{\mathrm{L}}\left({\mathit{\eta }}_{\mathrm{I}}\right)\mathrm{\times }{\mathit{R}}^{\mathit{T}}\mathrm{+}\mathit{R}\mathrm{·}\mathrm{\Pi }\mathrm{·}{\mathit{R}}^{\mathit{T}}\mathit{,}$(A.12)where C_L(η_I) is the linear correlation function at an early time η_I, with η_I →  − ∞, and the first product “ × C_L × ” does not contain any integration over time.

Then, at 1-loop order the direct steepest-descent scheme amounts to use for Σ and Π their lowest order terms, over powers of C_L and R_L. This means that Eq. (A.12) for the correlation C reads as the first line in Fig. 1, while a similar diagrammatic expression holds for the response R (see Fig. 8 in Valageas 2008). In particular, solving Eq. (A.11) as a perturbative series over powers of Σ gives a diagrammatic expansion for the response R in terms of “bubble” diagrams, and substituting into Eq. (A.12) this gives the series of diagrams shown in the second equality of Fig. 1. In practice, one does not compute the two-point function C(x,y), whence the density power spectrum, from these diagrammatic series, but from the integro-differential Eq. (A.11) and the explicit expression (A.12).

First, using in the following the approximation Ω_m/f² ≃ 1, the linear correlation and response functions read as ${\mathit{C}}_{\mathrm{L}}\mathrm{\left(}{\mathit{x}}_{\mathrm{1}}\mathit{,}{\mathit{x}}_{\mathrm{2}}\mathrm{\right)}\mathrm{=}{\mathrm{e}}^{{\mathit{\eta }}_{\mathrm{1}}\mathrm{+}{\mathit{\eta }}_{\mathrm{2}}}\hspace{0.17em}{\mathit{P}}_{\mathrm{L}\mathrm{0}}\mathrm{\left(}{\mathit{k}}_{\mathrm{1}}\mathrm{\right)}\hspace{0.17em}{\mathit{\delta }}_{\mathrm{D}}\left({k}\mathrm{1}\mathrm{+}{k}\mathrm{2}\right)\left(\begin{array}{c}\\ \mathrm{1}& \mathrm{1}\\ \mathrm{1}& \mathrm{1}\end{array}\right)\mathit{,}$(A.13)where P_L0(k) is the linear density power spectrum today (when η = 0), and $\begin{array}{ccc}{\mathit{R}}_{\mathrm{L}}\mathrm{(}{\mathit{x}}_{\mathrm{1}}\mathit{,}{\mathit{x}}_{\mathrm{2}}\mathrm{)}& \mathrm{=}& \mathit{\theta }\left({\mathit{\eta }}_{\mathrm{1}}\mathrm{-}{\mathit{\eta }}_{\mathrm{2}}\right)\hspace{0.17em}{\mathit{\delta }}_{\mathrm{D}}\left({k}\mathrm{1}\mathrm{-}{k}\mathrm{2}\right)\\ & & \mathrm{\times }\mathrm{[}{\mathrm{e}}^{{\mathit{\eta }}_{\mathrm{1}}\mathrm{-}{\mathit{\eta }}_{\mathrm{2}}}{\mathit{R}}_{\mathrm{0}}^{\mathrm{+}}\mathrm{+}{\mathrm{e}}^{\mathrm{-}\mathrm{3}\left({\mathit{\eta }}_{\mathrm{1}}\mathrm{-}{\mathit{\eta }}_{\mathrm{2}}\right)\mathit{/}\mathrm{2}}{{\mathit{R}}_{\mathrm{0}}^{\mathrm{-}}}^{\mathrm{]}}\mathit{,}\end{array}$(A.14)with ${\mathit{R}}_{\mathrm{0}}^{\mathrm{+}}\mathrm{=}\left(\begin{array}{c}\\ \mathrm{3}\mathit{/}\mathrm{5}& \mathrm{2}\mathit{/}\mathrm{5}\\ \mathrm{3}\mathit{/}\mathrm{5}& \mathrm{2}\mathit{/}\mathrm{5}\end{array}\right)\mathrm{and}{\mathit{R}}_{\mathrm{0}}^{\mathrm{-}}\mathrm{=}\left(\begin{array}{c}\\ \mathrm{2}\mathit{/}\mathrm{5}& \mathrm{-}\mathrm{2}\mathit{/}\mathrm{5}\\ \mathrm{-}\mathrm{3}\mathit{/}\mathrm{5}& \mathrm{3}\mathit{/}\mathrm{5}\end{array}\right)\mathit{.}$(A.15)In Eq. (A.14) the Heaviside factor θ(η₁ − η₂) expresses the causality of the response function. Then, the “self-energy” terms read at lowest order as (Valageas 2007a) $\begin{array}{ccc}{\mathrm{\Sigma }}_{\mathrm{0}}\mathrm{(}{\mathit{x}}_{\mathrm{1}}\mathit{,}{\mathit{x}}_{\mathrm{2}}\mathrm{)}& \mathrm{=}& \mathit{\theta }\left({\mathit{\eta }}_{\mathrm{1}}\mathrm{-}{\mathit{\eta }}_{\mathrm{2}}\right)\hspace{0.17em}{\mathit{\delta }}_{\mathrm{D}}\left({k}\mathrm{1}\mathrm{-}{k}\mathrm{2}\right)\\ & & \mathrm{\times }\mathrm{[}{\mathrm{e}}^{\mathrm{2}{\mathit{\eta }}_{\mathrm{1}}}{\mathrm{\Sigma }}_{\mathrm{0}}^{\mathrm{+}}\mathrm{(}{\mathit{k}}_{\mathrm{1}}\mathrm{)}\mathrm{+}{\mathrm{e}}^{\mathrm{-}{\mathit{\eta }}_{\mathrm{1}}\mathit{/}\mathrm{2}\mathrm{+}\mathrm{5}{\mathit{\eta }}_{\mathrm{2}}\mathit{/}\mathrm{2}}{\mathrm{\Sigma }}_{\mathrm{0}}^{\mathrm{-}}{\mathrm{(}{\mathit{k}}_{\mathrm{1}}\mathrm{)}}^{\mathrm{]}}\mathit{,}\end{array}$(A.16)with $\begin{array}{ccc}{\mathrm{\Sigma }}_{\mathrm{0}\mathrm{;}{\mathit{i}}_{\mathrm{1}}{\mathit{i}}_{\mathrm{2}}}^{\mathrm{\pm }}\mathrm{\left(}\mathit{k}\mathrm{\right)}& \mathrm{=}& \mathrm{4}\sum _{{\mathit{j}}_{\mathrm{1}}{\mathit{j}}_{\mathrm{2}}{\mathit{l}}_{\mathrm{1}}{\mathit{l}}_{\mathrm{2}}}\mathrm{\int }\mathrm{d}{q}{\mathit{\gamma }}_{{\mathit{i}}_{\mathrm{1}}\mathrm{;}{\mathit{j}}_{\mathrm{1}}\mathit{,}{\mathit{j}}_{\mathrm{2}}}^{\mathit{s}}\left({k}\mathrm{-}{q}\mathit{,}{q}\right)\\ & & \mathrm{\times }{\mathit{\gamma }}_{{\mathit{l}}_{\mathrm{1}}\mathrm{;}{\mathit{i}}_{\mathrm{2}}\mathit{,}{\mathit{l}}_{\mathrm{2}}}^{\mathit{s}}\mathrm{(}{k}\mathit{,}\mathrm{-}{q}\mathrm{)}\hspace{0.17em}{\mathit{P}}_{\mathrm{L}\mathrm{0}}\mathrm{\left(}\mathit{q}\mathrm{\right)}\hspace{0.17em}{\mathit{R}}_{\mathrm{0}\mathrm{;}{\mathit{j}}_{\mathrm{1}}{\mathit{l}}_{\mathrm{1}}}^{\mathrm{\pm }}\mathit{,}\end{array}$(A.17)and ${\mathrm{\Pi }}_{\mathrm{0}}\mathrm{\left(}{\mathit{x}}_{\mathrm{1}}\mathit{,}{\mathit{x}}_{\mathrm{2}}\mathrm{\right)}\mathrm{=}{\mathit{\delta }}_{\mathrm{D}}\left({k}\mathrm{1}\mathrm{+}{k}\mathrm{2}\right)\hspace{0.17em}{\mathrm{e}}^{\mathrm{2}{\mathit{\eta }}_{\mathrm{1}}\mathrm{+}\mathrm{2}{\mathit{\eta }}_{\mathrm{2}}}\hspace{0.17em}{\mathrm{\Pi }}_{\mathrm{0}}\mathrm{\left(}{\mathit{k}}_{\mathrm{1}}\mathrm{\right)}\mathit{,}$(A.18)with $\begin{array}{ccc}{\mathrm{\Pi }}_{\mathrm{0}\mathrm{;}{\mathit{i}}_{\mathrm{1}}{\mathit{i}}_{\mathrm{2}}}\mathrm{\left(}\mathit{k}\mathrm{\right)}& \mathrm{=}& \mathrm{2}\sum _{{\mathit{j}}_{\mathrm{1}}{\mathit{j}}_{\mathrm{2}}{\mathit{l}}_{\mathrm{1}}{\mathit{l}}_{\mathrm{2}}}\mathrm{\int }\mathrm{d}{q}{\mathit{\gamma }}_{{\mathit{i}}_{\mathrm{1}}\mathrm{;}{\mathit{j}}_{\mathrm{1}}\mathit{,}{\mathit{j}}_{\mathrm{2}}}^{\mathit{s}}\left({q}\mathit{,}{k}\mathrm{-}{q}\right)\\ & & \mathrm{\times }{\mathit{\gamma }}_{{\mathit{i}}_{\mathrm{2}}\mathrm{;}{\mathit{l}}_{\mathrm{1}}\mathit{,}{\mathit{l}}_{\mathrm{2}}}^{\mathit{s}}\left(\mathrm{-}{q}\mathit{,}\mathrm{-}{k}\mathrm{+}{q}\right){\mathit{P}}_{\mathrm{L}\mathrm{0}}\mathrm{\left(}\mathit{q}\mathrm{\right)}{\mathit{P}}_{\mathrm{L}\mathrm{0}}\mathrm{\left(}\mathrm{\right|}{k}\mathrm{-}{q}\mathrm{\left|}\mathrm{\right)}\mathit{.}\end{array}$(A.19)These “self-energy” terms also satisfy the symmetries ${\mathrm{\Sigma }}_{\mathrm{0}}^{\mathrm{-}}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathrm{=}\left(\begin{array}{c}\\ {\mathrm{\Sigma }}_{\mathrm{0}\mathrm{;}\mathrm{22}}^{\mathrm{+}}& \mathrm{-}{\mathrm{\Sigma }}_{\mathrm{0}\mathrm{;}\mathrm{12}}^{\mathrm{+}}\\ \mathrm{-}{\mathrm{\Sigma }}_{\mathrm{0}\mathrm{;}\mathrm{21}}^{\mathrm{+}}& {\mathrm{\Sigma }}_{\mathrm{0}\mathrm{;}\mathrm{11}}^{\mathrm{+}}\end{array}\right)\mathit{,}{\mathrm{\Pi }}_{\mathrm{0}\mathrm{;}\mathit{ij}}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathrm{=}{\mathrm{\Pi }}_{\mathrm{0}\mathrm{;}\mathit{ji}}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathit{.}$(A.20)Thanks to the simple time dependence of Σ₀(x₁,x₂), which arises from the fact that we have expanded Σ over the linear two-point functions C_L and R_L and only kept the lowest order term, it is possible to eliminate the integral terms in the Schwinger-Dyson Eq. (A.11) by taking successive derivatives with respect to η₁. This yields $\begin{array}{ccc}\frac{{\mathit{\partial }}^{\mathrm{3}}{\mathit{R}}_{\mathrm{1}}}{\mathit{\partial }{\mathit{\eta }}_{\mathrm{1}}^{\mathrm{3}}}\mathrm{-}\frac{\mathrm{3}}{\mathrm{2}}\frac{{\mathit{\partial }}^{\mathrm{2}}{\mathit{R}}_{\mathrm{1}}}{\mathit{\partial }{\mathit{\eta }}_{\mathrm{1}}^{\mathrm{2}}}\mathrm{-}\frac{{\mathit{\partial }}^{\mathrm{2}}{\mathit{R}}_{\mathrm{2}}}{\mathit{\partial }{\mathit{\eta }}_{\mathrm{1}}^{\mathrm{2}}}\mathrm{-}\frac{\mathit{\partial }{\mathit{R}}_{\mathrm{1}}}{\mathit{\partial }{\mathit{\eta }}_{\mathrm{1}}}\mathrm{+}\frac{\mathrm{3}}{\mathrm{2}}\frac{\mathit{\partial }{\mathit{R}}_{\mathrm{2}}}{\mathit{\partial }{\mathit{\eta }}_{\mathrm{1}}}\mathrm{+}{\mathit{R}}_{\mathrm{2}}\mathrm{=}& & \\ & & {\mathrm{e}}^{\mathrm{2}{\mathit{\eta }}_{\mathrm{1}}}\left[\mathrm{(}{\mathrm{\Sigma }}_{\mathrm{0}\mathrm{;}\mathrm{11}}^{\mathrm{+}}\mathrm{+}{{\mathrm{\Sigma }}_{\mathrm{0}\mathrm{;}\mathrm{22}}^{\mathrm{+}}}^{\mathrm{)}}\frac{\mathit{\partial }{\mathit{R}}_{\mathrm{1}}}{\mathit{\partial }{\mathit{\eta }}_{\mathrm{1}}}\mathrm{+}\frac{\mathrm{5}}{\mathrm{2}}{\mathrm{\Sigma }}_{\mathrm{0}\mathrm{;}\mathrm{11}}^{\mathrm{+}}{\mathit{R}}_{\mathrm{1}}\mathrm{+}\frac{\mathrm{5}}{\mathrm{2}}{\mathrm{\Sigma }}_{\mathrm{0}\mathrm{;}\mathrm{12}}^{\mathrm{+}}{\mathit{R}}_{\mathrm{2}}\right]\end{array}$(A.21)$\begin{array}{ccc}\frac{{\mathit{\partial }}^{\mathrm{3}}{\mathit{R}}_{\mathrm{2}}}{\mathit{\partial }{\mathit{\eta }}_{\mathrm{1}}^{\mathrm{3}}}\mathrm{-}\frac{\mathrm{3}}{\mathrm{2}}\frac{{\mathit{\partial }}^{\mathrm{2}}{\mathit{R}}_{\mathrm{1}}}{\mathit{\partial }{\mathit{\eta }}_{\mathrm{1}}^{\mathrm{2}}}\mathrm{-}\frac{{\mathit{\partial }}^{\mathrm{2}}{\mathit{R}}_{\mathrm{2}}}{\mathit{\partial }{\mathit{\eta }}_{\mathrm{1}}^{\mathrm{2}}}\mathrm{+}\frac{\mathrm{9}}{\mathrm{4}}\frac{\mathit{\partial }{\mathit{R}}_{\mathrm{1}}}{\mathit{\partial }{\mathit{\eta }}_{\mathrm{1}}}\mathrm{-}\frac{\mathrm{7}}{\mathrm{4}}\frac{\mathit{\partial }{\mathit{R}}_{\mathrm{2}}}{\mathit{\partial }{\mathit{\eta }}_{\mathrm{1}}}\mathrm{+}\frac{\mathrm{3}}{\mathrm{2}}{\mathit{R}}_{\mathrm{1}}\mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{R}}_{\mathrm{2}}\mathrm{=}& & \\ & & {\mathrm{e}}^{\mathrm{2}{\mathit{\eta }}_{\mathrm{1}}}\left[\mathrm{(}{\mathrm{\Sigma }}_{\mathrm{0}\mathrm{;}\mathrm{11}}^{\mathrm{+}}\mathrm{+}{{\mathrm{\Sigma }}_{\mathrm{0}\mathrm{;}\mathrm{22}}^{\mathrm{+}}}^{\mathrm{)}}\frac{\mathit{\partial }{\mathit{R}}_{\mathrm{2}}}{\mathit{\partial }{\mathit{\eta }}_{\mathrm{1}}}\mathrm{+}\frac{\mathrm{5}}{\mathrm{2}}{\mathrm{\Sigma }}_{\mathrm{0}\mathrm{;}\mathrm{21}}^{\mathrm{+}}{\mathit{R}}_{\mathrm{1}}\mathrm{+}\frac{\mathrm{5}}{\mathrm{2}}{\mathrm{\Sigma }}_{\mathrm{0}\mathrm{;}\mathrm{22}}^{\mathrm{+}}{\mathit{R}}_{\mathrm{2}}\right]\end{array}$(A.22)where the pair (R₁,R₂) corresponds either to (R₁₁,R₂₁) or (R₁₂,R₂₂). These equations apply to η₁ > η₂ and the initial conditions at η₁ = η₂ are $\begin{array}{ccc}{\mathit{X}}^{\mathit{T}}& \mathrm{=}& \left(\begin{array}{c}\\ & & & & & \end{array}\right)\\ & \mathrm{=}& \left(\begin{array}{c}\\ & & & & & \end{array}\right)\end{array}$for the pair (R₁₁,R₂₁), and ${\mathit{X}}^{\mathit{T}}\mathrm{=}\left(\begin{array}{c}\\ \mathrm{0}& \mathrm{1}& \mathrm{1}& \mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}& \mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}& \frac{\mathrm{7}}{\mathrm{4}}\mathrm{+}{\mathrm{e}}^{\mathrm{2}{\mathit{\eta }}_{\mathrm{1}}}\mathrm{(}{\mathrm{\Sigma }}_{\mathrm{0}\mathrm{;}\mathrm{11}}^{\mathrm{+}}\mathrm{+}{{\mathrm{\Sigma }}_{\mathrm{0}\mathrm{;}\mathrm{22}}^{\mathrm{+}}}^{\mathrm{)}}\end{array}\right)$(A.25)for the pair (R₁₂,R₂₂). Thus, we obtain the response function by solving the differential Eqs. (A.21), (A.22), with the initial conditions (A.24), (A.25), using an ordinary Runge-Kutta algorithm. A great simplification is that there are no longer integral terms over past times in the right hand side of Eqs. (A.21), (A.22), so that numerical computations are very fast.

Finally, it is possible to take advantage of the factorized time dependence of Π₀ in Eq. (A.18) to factorize Eq. (A.12). To compute the first product in Eq. (A.12) we define the two-component vector ${\mathit{\varphi }}_{\mathit{i}}^{\mathrm{I}}\mathrm{\left(}\mathit{k,\eta }\mathrm{\right)}\mathrm{=}\sum _{\mathit{j}\mathrm{=}\mathrm{1}}^{\mathrm{2}}{\mathit{R}}_{\mathit{i,j}}\left(\mathit{k}\mathrm{;}\mathit{\eta ,}{\mathit{\eta }}_{\mathrm{I}}\right)\hspace{0.17em}{\mathrm{e}}^{{\mathit{\eta }}_{\mathrm{I}}}\hspace{0.17em}\sqrt{{\mathit{P}}_{\mathrm{L}\mathrm{0}}\mathrm{\left(}\mathit{k}\mathrm{\right)}}\mathit{,}$(A.26)while for the second product we diagonalize the symmetric matrix Π₀ as ${\mathrm{\Pi }}_{\mathrm{0}}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathrm{=}{\mathit{p}}^{\mathrm{+}}{\mathit{\pi }}^{\mathrm{+}}\mathrm{·}{\mathit{\pi }}^{\mathrm{+}\mathit{T}}\mathrm{+}{\mathit{p}}^{\mathrm{-}}{\mathit{\pi }}^{\mathrm{-}}\mathrm{·}{\mathit{\pi }}^{\mathrm{-}\mathit{T}}\mathit{,}$(A.27)where p ^±  are the eigenvalues of Π₀ and π ^±  its normalized eigenvectors, and we define the two-component vectors, ${\mathit{\varphi }}^{\mathrm{\pm }}\mathrm{\left(}\mathit{k,\eta }\mathrm{\right)}\mathrm{=}{\mathrm{\int }}_{\mathrm{-}\mathrm{\infty }}^{\mathit{\eta }}\mathrm{d}{\mathit{\eta }}^{\mathrm{\prime }}\hspace{0.17em}{\mathrm{e}}^{\mathrm{2}{\mathit{\eta }}^{\mathrm{\prime }}}\hspace{0.17em}\mathit{R}\left(\mathit{k}\mathrm{;}\mathit{\eta ,}{\mathit{\eta }}^{\mathrm{\prime }}\right)\mathrm{·}{\mathit{\pi }}^{\mathrm{\pm }}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathit{.}$(A.28)Then, Eq. (A.12) writes as $\begin{array}{ccc}\mathit{C}\left(\mathit{k}\mathrm{;}{\mathit{\eta }}_{\mathrm{1}}\mathit{,}{\mathit{\eta }}_{\mathrm{2}}\right)& \mathrm{=}& {\mathit{\varphi }}^{\mathrm{I}}\left({\mathit{\eta }}_{\mathrm{1}}\right)\mathrm{·}{\mathit{\varphi }}^{\mathrm{I}}{\left({\mathit{\eta }}_{\mathrm{2}}\right)}^{\mathit{T}}\\ & & \mathrm{+}{\mathit{p}}^{\mathrm{+}}\hspace{0.17em}{\mathit{\varphi }}^{\mathrm{+}}\left({\mathit{\eta }}_{\mathrm{1}}\right)\mathrm{·}{\mathit{\varphi }}^{\mathrm{+}}{\left({\mathit{\eta }}_{\mathrm{2}}\right)}^{\mathit{T}}\mathrm{+}{\mathit{p}}^{\mathrm{-}}\hspace{0.17em}{\mathit{\varphi }}^{\mathrm{-}}\left({\mathit{\eta }}_{\mathrm{1}}\right)\mathrm{·}{\mathit{\varphi }}^{\mathrm{-}}{\left({\mathit{\eta }}_{\mathrm{2}}\right)}^{\mathit{T}}\mathit{.}\end{array}$(A.29)In particular, the nonlinear power spectrum is given by the equal-time upper-left component of the correlation matrix C, $\mathit{P}\mathrm{\left(}\mathit{k,\eta }\mathrm{\right)}\mathrm{=}{\mathit{C}}_{\mathrm{11}}\mathrm{\left(}\mathit{k}\mathrm{;}\mathit{\eta ,\eta }\mathrm{\right)}\mathit{.}$(A.30)The factorized form (A.29), which avoids to compute two-dimensional integrals over past times as in generic expressions of Eq. (A.12), thanks to the factorized time dependence of Π₀ in Eq. (A.18), again allows fast numerical computations.

All Tables

All Figures