© ESO, 2011

1. Introduction

According to standard cosmological scenarios, the large-scale structures of the present Universe have formed through the amplification by gravitational instability of small almost-Gaussian primordial fluctuations (Peebles 1980). Then, from observations of the recent Universe, through galaxy surveys (Tegmark et al. 2006; Cole et al. 2005), weak-lensing studies (Massey et al. 2007; Munshi et al. 2008), measures of baryon acoustic oscillations (Eisenstein et al. 1998, 2005), one can derive constraints on the cosmological parameters (such as the mean matter and dark energy contents) and on the properties of the initial perturbations. The main statistical quantity used in this context is the density two-point correlation function, or power spectrum in Fourier space. This allows measuring the scale dependence and the amplitude of the initial conditions, and for Gaussian fields this fully determines all statistical properties of the matter distribution. However, to constrain the possible non-Gaussianities of the initial perturbations, to check the gravitational clustering scenario (through the measure of the non-Gaussianities generated by the nonlinear dynamics), and to break parameter degeneracies, it is useful to study higher order statistics. The three-point correlation function or bispectrum in Fourier space, which is the lowest order statistics beyond the Gaussian, is the main quantity studied in this context (Peebles 1980; Frieman & Gaztanaga 1994, 1999; Scoccimarro & Couchman 2001; Bernardeau et al. 2002b; Verde et al. 2002; Takada & Jain 2003; Kayo et al. 2004; Sefusatti et al. 2007; Sefusatti & Komatsu 2007; Nishimichi et al. 2007; Sefusatti et al. 2010; Nishimichi et al. 2010).

On large scales, or at early times, where the density fluctuations are small within cold dark matter (CDM) scenarios (Peebles 1982), one can use perturbation theory (Goroff et al. 1986; Scoccimarro 1997; Scoccimarro et al. 1998; Bernardeau et al. 2002a). In order to extend the validity of perturbative predictions to somewhat smaller scales and to increase their theoretical accuracy, many resummation schemes have been proposed recently (Crocce & Scoccimarro 2006b,a; Valageas 2007a; Matarrese & Pietroni 2007; Pietroni 2008; Taruya & Hiramatsu 2008; Matsubara 2008b; Taruya et al. 2009). Although these methods follow different routes and involve different expansion schemes, they are consistent with each other and with the standard perturbation theory up to the order of truncation, and only differ by higher order contributions. Thus, they complete the standard perturbation theory contributions (which consist of a finite number of Feynman diagrams) by (usually infinite) partial resummations of higher order diagrams. Most of these studies have focused on the matter power spectrum, except for Valageas (2008), who also considers three-point (and higher) density correlations, Bernardeau et al. (2008), who study higher order response functions (propagators) and the bispectrum, and Bartolo et al. (2010a), who consider both the power spectrum and bispectrum for Gaussian and non-Gaussian initial conditions.

On small scales, where the density fluctuations are large, one must use numerical simulations or phenomenological models, such as the Lagrangian mapping introduced in Hamilton et al. (1991) or the halo model (McClelland & Silk 1977; Scherrer & Bertschinger 1991; Cooray & Sheth 2002). Both approaches have already been used for the bispectrum. However, while the first method has only had some success on intermediate mildly nonlinear scales (Pan et al. 2007), the second method is also able to describe the highly nonlinear scales (Ma & Fry 2000b; Scoccimarro et al. 2001; Wang et al. 2004; Fosalba et al. 2005).

In order to compare theoretical predictions with observations, which go from large linear scales to small highly nonlinear scales, it is useful to build unified models, which combine for instance perturbation theories with halo models. This combines the high accuracy on large scales (~1%) provided by systematic perturbative expansions with the reasonably good accuracy on small scales (~10%) provided by the halo model. This is particularly useful on the weakly nonlinear scales associated with the baryon acoustic oscillations, where the matter power spectrum and bispectrum show small “wiggles” that would be difficult to reproduce with a high accuracy by a simple halo model or fits to simulations (unless one runs a simulation with the required cosmological parameters). As a first step in this direction, we have recently built such a unified model in Valageas & Nishimichi (2011), focusing on the matter power spectrum. In this paper, we extend this study to the matter bispectrum, using the same halo model and perturbative schemes. As for the power spectrum, we explain that the “1-halo” and “2-halo” contributions contain new counterterms that ensure a physical behavior on large scales. Then, we show that we obtain a good agreement with numerical simulations without introducing new parameters. We also note that it is possible to improve the predictions for the power spectrum on mildly nonlinear scales using the shape of the predicted bispectrum.

This paper is organized as follows. We first describe in Sect. 2 the computation of the density bispectrum from a Lagrangian point of view, making the link with halo models and perturbative expansions through the decomposition into 1-halo, 2-halo, and 3-halo contributions. Then, we present in Sect. 3 a simple implementation, based on the halo model and the perturbative approaches used in Valageas & Nishimichi (2011) for the power spectrum. We compare our results for the bispectrum with numerical simulations in Sect. 4, for equilateral and isosceles configurations, from linear to highly nonlinear scales. Next, we discuss the differences between our work and some previous studies in Sect. 5, and we investigate in Sect. 6 the dependence of the predictions on the choice of the perturbative scheme and on the halo parameters. We explain in Sect. 7 how to use the shape of the predicted reduced bispectrum to improve the predictions for the power spectrum and the bispectrum, and we estimate the accuracy of our unified model in Sect. 8 before concluding in Sect. 9.

2. Decomposition of the density bispectrum from a Lagrangian point of view

We show in this section how the density bispectrum can be split into 3-halo, 2-halo, and 1-halo terms from a Lagrangian point of view, and we make the connection with perturbation theory.

2.1. Lagrangian-space formulation

In a Lagrangian framework, one considers the trajectories x(q,t) of the particles, of initial Lagrangian coordinates q, and Eulerian coordinates x at time t. At any 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)We define the power spectrum, P(k), and the bispectrum, B(k₁,k₂,k₃), which are the Fourier transforms of the two-point and three-point density correlation functions, by $\begin{array}{ccc}& & \mathrm{⟨}\begin{array}{}\mathrm{˜}\\ \mathit{\delta }\end{array}\mathrm{\left(}{{k}}_{\mathrm{1}}\mathrm{\right)}\begin{array}{}\mathrm{˜}\\ \mathit{\delta }\end{array}\mathrm{\left(}{{k}}_{\mathrm{2}}\mathrm{\right)}\mathrm{⟩}\mathrm{=}{\mathit{\delta }}_{\mathit{D}}\mathrm{(}{{k}}_{\mathrm{1}}\mathrm{+}{{k}}_{\mathrm{2}}\mathrm{)}\hspace{0.17em}\mathit{P}\mathrm{\left(}{\mathit{k}}_{\mathrm{1}}\mathrm{\right)}\mathit{,}\\ & & \mathrm{⟨}\begin{array}{}\mathrm{˜}\\ \mathit{\delta }\end{array}\mathrm{\left(}{{k}}_{\mathrm{1}}\mathrm{\right)}\begin{array}{}\mathrm{˜}\\ \mathit{\delta }\end{array}\mathrm{\left(}{{k}}_{\mathrm{2}}\mathrm{\right)}\begin{array}{}\mathrm{˜}\\ \mathit{\delta }\end{array}\mathrm{\left(}{{k}}_{\mathrm{3}}\mathrm{\right)}\mathrm{⟩}\mathrm{=}{\mathit{\delta }}_{\mathit{D}}\mathrm{(}{{k}}_{\mathrm{1}}\mathrm{+}{{k}}_{\mathrm{2}}\mathrm{+}{{k}}_{\mathrm{3}}\mathrm{)}\hspace{0.17em}\mathit{B}\mathrm{\left(}{\mathit{k}}_{\mathrm{1}}\mathit{,}{\mathit{k}}_{\mathrm{2}}\mathit{,}{\mathit{k}}_{\mathrm{3}}\mathrm{\right)}\mathit{.}\end{array}$Here we used the fact that the system is statistically homogeneous, which gives rise to the Dirac prefactors associated with statistical invariance through translations, and statistically isotropic, which implies that P(k₁) and B(k₁,k₂,k₃) only depend on the lengths |k₁| and  {|k₁|,|k₂|,|k₃|} .

We have already studied the power spectrum in a previous work (Valageas & Nishimichi 2011), so that we here focus on the bispectrum. Our goal is to follow the same approach to combine the perturbation theory and the halo model and to build a simple unified model. Within this framework, we wish to decompose the bispectrum into “perturbative” and “non-perturbative” contributions. The former is associated with configurations of particle triplets (or p-uplets for the p-point correlation function) where the particles belong to different halos, whereas the latter is associated with configurations where several particles belong to the same halo. As in Valageas & Nishimichi (2011), our strategy is to evaluate the perturbative contribution through perturbation theory (either with the standard expansion or with resummation schemes) and the non-perturbative contribution through a halo model. Then, the Lagrangian framework allows us to recover the counterterms that arise in the halo-model contributions and that are necessary to ensure a good behavior on large scales. In addition, it is of interest to see how the standard expression for the bispectrum obtained from an Eulerian point of view (Cooray & Sheth 2002) can be recovered (with the addition of these new counterterms) from a Lagrangian point of view.

In Valageas & Nishimichi (2011), where we studied the power spectrum, we could factorize out the Dirac prefactor of Eq. (5) and write for the power spectrum the well-known expression (Schneider & Bartelmann 1995; Taylor & Hamilton 1996) $\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{-}\mathrm{i}{k}\mathrm{·}\mathrm{\Delta }{x}}\mathrm{-}{\mathrm{e}}^{\mathrm{i}{k}\mathrm{·}{q}}\mathrm{⟩}\mathit{,}$(7)where we introduced the Eulerian-space separation Δx = x(q) − x(0). Thus, Eq. (7) only depends on the relative displacements of the particles, as it must because uniform translations do not affect the structural properties of the density field. Then, Eq. (7) provides a convenient starting point because the Dirac prefactor of Eq. (5) has already been taken into account and it will not be put in danger by approximations used in later steps.

For the bispectrum, it is still possible to factorize out the Dirac prefactor of Eq. (6) by introducing the relative displacements. This yields $\begin{array}{ccc}\mathit{B}\mathrm{\left(}{\mathit{k}}_{\mathrm{1}}\mathit{,}{\mathit{k}}_{\mathrm{2}}\mathit{,}{\mathit{k}}_{\mathrm{3}}\mathrm{\right)}& \mathrm{=}& \int \frac{\mathrm{d}{{q}}_{\mathrm{2}}{{q}}_{\mathrm{3}}}{\mathrm{(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathrm{6}}}\hspace{0.17em}\mathrm{⟨}{\mathrm{e}}^{\mathrm{-}\mathrm{i}{{k}}_{\mathrm{2}}\mathrm{·}\mathrm{\Delta }{{x}}_{\mathrm{2}}\mathrm{-}\mathrm{i}{{k}}_{\mathrm{3}}\mathrm{·}\mathrm{\Delta }{{x}}_{\mathrm{3}}}\mathrm{⟩}\\ & & \mathrm{-}\hspace{0.17em}{\mathit{\delta }}_{\mathit{D}}\mathrm{\left(}{{k}}_{\mathrm{1}}\mathrm{\right)}\mathit{P}\mathrm{\left(}{\mathit{k}}_{\mathrm{2}}\mathrm{\right)}\mathrm{-}{\mathit{\delta }}_{\mathit{D}}\mathrm{\left(}{{k}}_{\mathrm{2}}\mathrm{\right)}\mathit{P}\mathrm{\left(}{\mathit{k}}_{\mathrm{3}}\mathrm{\right)}\mathrm{-}{\mathit{\delta }}_{\mathit{D}}\mathrm{\left(}{{k}}_{\mathrm{3}}\mathrm{\right)}\mathit{P}\mathrm{\left(}{\mathit{k}}_{\mathrm{1}}\mathrm{\right)}\\ & & \mathrm{-}\hspace{0.17em}{\mathit{\delta }}_{\mathit{D}}\mathrm{\left(}{{k}}_{\mathrm{1}}\mathrm{\right)}{\mathit{\delta }}_{\mathit{D}}\mathrm{\left(}{{k}}_{\mathrm{2}}\mathrm{\right)}\mathit{,}\end{array}$(8)where we introduced Δx₂ = x(q₂) − x(0) and Δx₃ = x(q₃) − x(0). However, the symmetry over  {k₁,k₂,k₃}  is no longer transparent in Eq. (8) and it is not very convenient to count the triplets that are within one, two, or three halos. Therefore, we prefer to work directly with the third-order moment (6), although this will require an additional approximation to those needed to compute the power spectrum in Valageas & Nishimichi (2011). Thus, using Eq. (4) we write $\begin{array}{ccc}\mathrm{⟨}\begin{array}{}\mathrm{˜}\\ \mathit{\delta }\end{array}\mathrm{\left(}{{k}}_{\mathrm{1}}\mathrm{\right)}\begin{array}{}\mathrm{˜}\\ \mathit{\delta }\end{array}\mathrm{\left(}{{k}}_{\mathrm{2}}\mathrm{\right)}\begin{array}{}\mathrm{˜}\\ \mathit{\delta }\end{array}\mathrm{\left(}{{k}}_{\mathrm{3}}\mathrm{\right)}\mathrm{⟩}& \mathrm{=}& {⟨}^{\mathrm{\int }}\frac{\mathrm{d}{{q}}_{\mathrm{1}}\mathrm{d}{{q}}_{\mathrm{2}}{{q}}_{\mathrm{3}}}{\mathrm{(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathrm{9}}}\hspace{0.17em}\mathrm{(}{\mathrm{e}}^{\mathrm{-}\mathrm{i}{{k}}_{\mathrm{1}}\mathrm{·}{{x}}_{\mathrm{1}}}\mathrm{-}{\mathrm{e}}^{\mathrm{-}\mathrm{i}{{k}}_{\mathrm{1}}\mathrm{·}{{q}}_{\mathrm{1}}}\mathrm{)}\\ & & \mathrm{\times }\hspace{0.17em}\mathrm{(}{\mathrm{e}}^{\mathrm{-}\mathrm{i}{{k}}_{\mathrm{2}}\mathrm{·}{{x}}_{\mathrm{2}}}\mathrm{-}{\mathrm{e}}^{\mathrm{-}\mathrm{i}{{k}}_{\mathrm{2}}\mathrm{·}{{q}}_{\mathrm{2}}}\mathrm{)}\mathrm{(}{\mathrm{e}}^{\mathrm{-}\mathrm{i}{{k}}_{\mathrm{3}}\mathrm{·}{{x}}_{\mathrm{3}}}\mathrm{-}{\mathrm{e}}^{\mathrm{-}\mathrm{i}{{k}}_{\mathrm{3}}\mathrm{·}{{q}}_{\mathrm{3}}}\mathrm{)}⟩\mathit{,}\end{array}$(9)where x_j = x(q_j). Next, we split the average (9) into three contributions, associated with the cases where the triplet  {q₁,q₂,q₃}  belongs to either one, two, or three halos: $\begin{array}{ccc}\mathrm{⟨}\begin{array}{}\mathrm{˜}\\ \mathit{\delta }\end{array}\mathrm{\left(}{{k}}_{\mathrm{1}}\mathrm{\right)}\begin{array}{}\mathrm{˜}\\ \mathit{\delta }\end{array}\mathrm{\left(}{{k}}_{\mathrm{2}}\mathrm{\right)}\begin{array}{}\mathrm{˜}\\ \mathit{\delta }\end{array}\mathrm{\left(}{{k}}_{\mathrm{3}}\mathrm{\right)}\mathrm{⟩}& \mathrm{=}& \mathrm{⟨}\begin{array}{}\mathrm{˜}\\ \mathit{\delta }\end{array}\mathrm{\left(}{{k}}_{\mathrm{1}}\mathrm{\right)}\begin{array}{}\mathrm{˜}\\ \mathit{\delta }\end{array}\mathrm{\left(}{{k}}_{\mathrm{2}}\mathrm{\right)}\begin{array}{}\mathrm{˜}\\ \mathit{\delta }\end{array}\mathrm{\left(}{{k}}_{\mathrm{3}}\mathrm{\right)}{\mathrm{⟩}}_{\mathrm{1}\mathrm{H}}\\ & & \mathrm{+}\mathrm{⟨}\begin{array}{}\mathrm{˜}\\ \mathit{\delta }\end{array}\mathrm{\left(}{{k}}_{\mathrm{1}}\mathrm{\right)}\begin{array}{}\mathrm{˜}\\ \mathit{\delta }\end{array}\mathrm{\left(}{{k}}_{\mathrm{2}}\mathrm{\right)}\begin{array}{}\mathrm{˜}\\ \mathit{\delta }\end{array}\mathrm{\left(}{{k}}_{\mathrm{3}}\mathrm{\right)}{\mathrm{⟩}}_{\mathrm{2}\mathrm{H}}\mathrm{+}\mathrm{⟨}\begin{array}{}\mathrm{˜}\\ \mathit{\delta }\end{array}\mathrm{\left(}{{k}}_{\mathrm{1}}\mathrm{\right)}\begin{array}{}\mathrm{˜}\\ \mathit{\delta }\end{array}\mathrm{\left(}{{k}}_{\mathrm{2}}\mathrm{\right)}\begin{array}{}\mathrm{˜}\\ \mathit{\delta }\end{array}\mathrm{\left(}{{k}}_{\mathrm{3}}\mathrm{\right)}{\mathrm{⟩}}_{\mathrm{3}\mathrm{H}}\mathit{.}\end{array}$(10)This assumes that all the matter is contained within halos, so that these three contributions sum up to Eq. (9), counting all particles.

2.2. “1-halo” contribution

Let us first consider the 1-halo contribution. It can be written as $\begin{array}{ccc}\mathrm{⟨}\begin{array}{}\mathrm{˜}\\ \mathit{\delta }\end{array}\mathrm{\left(}{{k}}_{\mathrm{1}}\mathrm{\right)}\begin{array}{}\mathrm{˜}\\ \mathit{\delta }\end{array}\mathrm{\left(}{{k}}_{\mathrm{2}}\mathrm{\right)}\begin{array}{}\mathrm{˜}\\ \mathit{\delta }\end{array}\mathrm{\left(}{{k}}_{\mathrm{3}}\mathrm{\right)}{\mathrm{⟩}}_{\mathrm{1}\mathrm{H}}& \mathrm{=}& {⟨}^{\mathrm{\int }}\underset{\mathit{j}\mathrm{=}\mathrm{1}}{\overset{\mathrm{3}}{\mathrm{􏽙}}}\frac{\mathrm{d}{{q}}_{\mathit{j}}}{\mathrm{(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathrm{3}}}\hspace{0.17em}\mathrm{(}{\mathrm{e}}^{\mathrm{-}\mathrm{i}{{k}}_{\mathit{j}}\mathrm{·}{{x}}_{\mathit{j}}}\mathrm{-}{\mathrm{e}}^{\mathrm{-}\mathrm{i}{{k}}_{\mathit{j}}\mathrm{·}{{q}}_{\mathit{j}}}\mathrm{)}\\ & & \mathrm{\times }\sum _{\mathit{\alpha }}\underset{\mathit{j}\mathrm{=}\mathrm{1}}{\overset{\mathrm{3}}{\mathrm{􏽙}}}\mathit{\theta }\mathrm{(}{{q}}_{\mathit{j}}\mathrm{\in }{\mathit{M}}_{\mathit{\alpha }}\mathrm{)}⟩\mathit{,}\end{array}$(11)where the sum runs over all halos, labeled by the index α, and the top-hat factors θ(q_j ∈ M_α) are unity if the particle q_j belongs to the halo α and zero otherwise. This clearly gives the contribution to Eq. (9) of the configurations where the triplet  {q₁,q₂,q₃}  belongs to a single halo. Equation (11) also reads as $\begin{array}{ccc}\mathrm{⟨}\begin{array}{}\mathrm{˜}\\ \mathit{\delta }\end{array}\mathrm{\left(}{{k}}_{\mathrm{1}}\mathrm{\right)}\begin{array}{}\mathrm{˜}\\ \mathit{\delta }\end{array}\mathrm{\left(}{{k}}_{\mathrm{2}}\mathrm{\right)}\begin{array}{}\mathrm{˜}\\ \mathit{\delta }\end{array}\mathrm{\left(}{{k}}_{\mathrm{3}}\mathrm{\right)}{\mathrm{⟩}}_{\mathrm{1}\mathrm{H}}& \mathrm{=}& {⟨}^{\mathrm{\int }}\underset{\mathit{j}\mathrm{=}\mathrm{1}}{\overset{\mathrm{3}}{\mathrm{􏽙}}}\frac{\mathrm{d}{{q}}_{\mathit{j}}}{\mathrm{(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathrm{3}}}\hspace{0.17em}\mathrm{(}{\mathrm{e}}^{\mathrm{-}\mathrm{i}{{k}}_{\mathit{j}}\mathrm{·}{{x}}_{\mathit{j}}}\mathrm{-}{\mathrm{e}}^{\mathrm{-}\mathrm{i}{{k}}_{\mathit{j}}\mathrm{·}{{q}}_{\mathit{j}}}\mathrm{)}\\ & & \mathrm{\times }\mathrm{\int }\mathrm{d}{{q}}^{\mathrm{c}}\mathrm{d}\mathit{M}\hspace{0.17em}\mathit{n̂}\mathrm{\left(}{{q}}^{\mathrm{c}}\mathit{,M}\mathrm{\right)}\hspace{0.17em}\underset{\mathit{j}\mathrm{=}\mathrm{1}}{\overset{\mathrm{3}}{\mathrm{􏽙}}}\mathit{\theta }\mathrm{(}{{q}}_{\mathit{j}}\mathrm{\in }\mathit{M}\mathrm{)}⟩\mathit{,}\end{array}$(12)where is the halo number density in a given realization of the initial conditions (as denoted by the hat, in contrast with its mean), $\mathit{n̂}\mathrm{\left(}{{q}}^{\mathrm{c}}\mathit{,M}\mathrm{\right)}\mathrm{=}\sum _{\mathit{\alpha }}{\mathit{\delta }}_{\mathit{D}}\mathrm{(}{{q}}^{\mathrm{c}}\mathrm{-}{{q}}_{\mathit{\alpha }}^{\mathrm{c}}\mathrm{)}\hspace{0.17em}{\mathit{\delta }}_{\mathit{D}}\mathrm{(}\mathit{M}\mathrm{-}{\mathit{M}}_{\mathit{\alpha }}\mathrm{)}\mathit{.}$(13)Here ${{q}}_{\mathit{\alpha }}^{\mathrm{c}}$ is the “center” of the halo α, which may be defined as the halo center of mass, in Lagrangian space (i.e., in the initial or linear density field). Note that in Eq. (12) we have chosen to count halos in Lagrangian space, by their Lagrangian center ${{q}}_{\mathit{\alpha }}^{\mathrm{c}}$, but we could as well count them in Eulerian space by using the Eulerian-space center of mass ${{x}}_{\mathit{\alpha }}^{\mathrm{c}}$. The mean of the halo number density does not depend on location, thanks to statistical homogeneity, and it is given by the halo mass function, $\mathrm{⟨}\mathit{n̂}\mathrm{\left(}{{q}}^{\mathrm{c}}\mathit{,M}\mathrm{\right)}\mathrm{⟩}\mathrm{=}\mathit{n}\mathrm{\left(}\mathit{M}\mathrm{\right)}\mathit{,}$(14)which we write 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)}}\mathrm{·}$(15)As usual, we have introduced in Eq. (15) the scaling function f(ν) and the reduced variable ν, where σ(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{,}$(16)and ${\mathit{\sigma }}^{\mathrm{2}}\mathrm{\left(}{\mathit{q}}_{\mathit{M}}\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{k}{\mathit{q}}_{\mathit{M}}{\mathrm{)}}^{\mathrm{2}}\mathit{,}$(17)where $\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{W}\end{array}\mathrm{\left(}\mathit{k}{\mathit{q}}_{\mathit{M}}\mathrm{\right)}$ is the Fourier transform of the top-hat of radius q_M, defined as $\begin{array}{}\mathrm{˜}\\ \mathit{W}\end{array}\mathrm{\left(}\mathit{k}{\mathit{q}}_{\mathit{M}}\mathrm{\right)}\mathrm{=}{\mathrm{\int }}_{{\mathit{V}}_{\mathit{M}}}\frac{\mathrm{d}{q}}{{\mathit{V}}_{\mathit{M}}}\hspace{0.17em}{\mathrm{e}}^{\mathrm{i}{k}\mathrm{·}{q}}\mathrm{=}\mathrm{3}\hspace{0.17em}\frac{\sin \mathrm{\left(}\mathit{k}{\mathit{q}}_{\mathit{M}}\mathrm{\right)}\mathrm{-}\mathit{k}{\mathit{q}}_{\mathit{M}}\cos \mathrm{\left(}\mathit{k}{\mathit{q}}_{\mathit{M}}\mathrm{\right)}}{\mathrm{(}\mathit{k}{\mathit{q}}_{\mathit{M}}{\mathrm{)}}^{\mathrm{3}}}\mathrm{·}$(18)In the second Eq. (15), the linear density contrast δ_L is related to the nonlinear density threshold δ that defines the halos through the spherical collapse dynamics, as δ = ℱ(δ_L), see Valageas (2009). Thus, Eq. (12) also writes as $\begin{array}{ccc}\mathrm{⟨}\begin{array}{}\mathrm{˜}\\ \mathit{\delta }\end{array}\mathrm{\left(}{{k}}_{\mathrm{1}}\mathrm{\right)}\begin{array}{}\mathrm{˜}\\ \mathit{\delta }\end{array}\mathrm{\left(}{{k}}_{\mathrm{2}}\mathrm{\right)}\begin{array}{}\mathrm{˜}\\ \mathit{\delta }\end{array}\mathrm{\left(}{{k}}_{\mathrm{3}}\mathrm{\right)}{\mathrm{⟩}}_{\mathrm{1}\mathrm{H}}& \mathrm{=}& \int \mathrm{d}{{q}}^{\mathrm{c}}\frac{\mathrm{d}\mathit{\nu }}{\mathit{\nu }}\hspace{0.17em}\frac{\overline{\mathit{\rho }}}{\mathit{M}}\mathit{f}\mathrm{\left(}\mathit{\nu }\mathrm{\right)}\\ & & \mathrm{\times }{⟨{\mathrm{\int }}_{{\mathit{V}}_{\mathit{M}}}\underset{\mathit{j}\mathrm{=}\mathrm{1}}{\overset{\mathrm{3}}{\mathrm{􏽙}}}\frac{\mathrm{d}{{q}}_{\mathit{j}}}{\mathrm{(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathrm{3}}}\hspace{0.17em}\mathrm{(}{\mathrm{e}}^{\mathrm{-}\mathrm{i}{{k}}_{\mathit{j}}\mathrm{·}{{x}}_{\mathit{j}}}\mathrm{-}{\mathrm{e}}^{\mathrm{-}\mathrm{i}{{k}}_{\mathit{j}}\mathrm{·}{{q}}_{\mathit{j}}}\mathrm{)}⟩}_{{{q}}^{\mathrm{c}}\mathit{,M}}\mathit{,}\end{array}$(19)where the average  ⟨ ... ⟩ _q^c,M is the conditional average under the constraint that the three particles q_j belong to the same halo of center q^c, mass M, and Lagrangian volume V_M. As in Valageas & Nishimichi (2011), we make the approximation of fully virialized halos, that is, the particle q has lost all memory of its initial location and velocity and it is located at random within the halo. This means that in Eulerian space the average is defined by the density profile of the halo, ρ_x^c,M(x), as $\mathrm{⟨}{\mathrm{e}}^{\mathrm{-}\mathrm{i}{k}\mathrm{·}{{x}}_{\mathit{j}}}{\mathrm{⟩}}_{{{q}}^{\mathrm{c}}\mathit{,M}}\mathrm{=}\frac{\mathrm{1}}{\mathit{M}}\mathrm{\int }\mathrm{d}{x}\hspace{0.17em}{\mathrm{e}}^{\mathrm{-}\mathrm{i}{k}\mathrm{·}{x}}\hspace{0.17em}{\mathit{\rho }}_{\mathit{M}}\mathrm{\left(}\mathrm{\right|}{x}\mathrm{-}{{x}}^{\mathrm{c}}\mathrm{\left|}\mathrm{\right)}\mathit{,}$(20)because q_j can be identified with a uniform probability to all the particles that make up the halo (approximation of complete relaxation). Here we also made the approximation of spherically symmetric halos, and introducing the normalized Fourier transform of the 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{,}$(21)we obtain $\mathrm{⟨}{\mathrm{e}}^{\mathrm{-}\mathrm{i}{k}\mathrm{·}{{x}}_{\mathit{j}}}{\mathrm{⟩}}_{{{q}}^{\mathrm{c}}\mathit{,M}}\mathrm{=}{\mathrm{e}}^{\mathrm{-}\mathrm{i}{k}\mathrm{·}{{x}}^{\mathrm{c}}}\hspace{0.17em}{\begin{array}{}\mathrm{˜}\\ \mathit{u}\end{array}}_{\mathit{M}}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathit{.}$(22)Making the approximation that these halos are also spherically symmetric in Lagrangian space, we can use Eq. (18) for the factors e^−ik_j·q_j and Eq. (19) reads as $\begin{array}{ccc}\mathrm{⟨}\begin{array}{}\mathrm{˜}\\ \mathit{\delta }\end{array}\mathrm{\left(}{{k}}_{\mathrm{1}}\mathrm{\right)}\begin{array}{}\mathrm{˜}\\ \mathit{\delta }\end{array}\mathrm{\left(}{{k}}_{\mathrm{2}}\mathrm{\right)}\begin{array}{}\mathrm{˜}\\ \mathit{\delta }\end{array}\mathrm{\left(}{{k}}_{\mathrm{3}}\mathrm{\right)}{\mathrm{⟩}}_{\mathrm{1}\mathrm{H}}& \mathrm{=}& \int \mathrm{d}{{q}}^{\mathrm{c}}\frac{\mathrm{d}\mathit{\nu }}{\mathit{\nu }}\hspace{0.17em}\frac{\overline{\mathit{\rho }}}{\mathit{M}}\mathit{f}\mathrm{\left(}\mathit{\nu }\mathrm{\right)}\hspace{0.17em}{\mathrm{e}}^{\mathrm{-}\mathrm{i}\mathrm{(}{{k}}_{\mathrm{1}}\mathrm{+}{{k}}_{\mathrm{2}}\mathrm{+}{{k}}_{\mathrm{3}}\mathrm{)}\mathrm{·}{{q}}^{\mathrm{c}}}\\ & & \mathrm{\times }{\left(\frac{\mathit{M}}{\overline{\mathit{\rho }}\mathrm{(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathrm{3}}}\right)}^{\mathrm{3}}\underset{\mathit{j}\mathrm{=}\mathrm{1}}{\overset{\mathrm{3}}{\mathrm{􏽙}}}\mathrm{(}{\mathrm{e}}^{\mathrm{-}\mathrm{i}{{k}}_{\mathit{j}}\mathrm{·}{\mathbf{\Psi }}^{\mathrm{c}}}{\begin{array}{}\mathrm{˜}\\ \mathit{u}\end{array}}_{\mathit{M}}\mathrm{\left(}{\mathit{k}}_{\mathit{j}}\mathrm{\right)}\mathrm{-}\begin{array}{}\mathrm{˜}\\ \mathit{W}\end{array}\mathrm{(}{\mathit{k}}_{\mathit{j}}{\mathit{q}}_{\mathit{M}}{\mathrm{)}}^{\mathrm{)}}\mathit{,}\end{array}$(23)where we introduced the displacement, Ψ^c = x^c − q^c, of the center of mass of the halo. We have written Eqs. (20) and (22), (23) as if the relation x^c(q^c) were deterministic, but a priori we should take the average over the displacement Ψ^c in Eq. (23). First, we note that  ⟨ e^−ik_j·Ψc ⟩  does not depend on q^c, thanks to statistical homogeneity, so that the integral over q^c in Eq. (23) gives the expected Dirac factor δ_D(k₁ + k₂ + k₃), in agreement with Eq. (6). Then, the displacement Ψ^c only appears in mixed terms, such as ${\mathrm{e}}^{\mathrm{-}\mathrm{i}\mathrm{(}{{k}}_{\mathrm{1}}\mathrm{+}{{k}}_{\mathrm{2}}\mathrm{)}\mathrm{·}{\mathbf{\Psi }}^{\mathrm{c}}}{\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{u}\end{array}}_{\mathit{M}}\mathrm{\left(}{\mathit{k}}_{\mathrm{1}}\mathrm{\right)}{\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{u}\end{array}}_{\mathit{M}}\mathrm{\left(}{\mathit{k}}_{\mathrm{2}}\mathrm{\right)}\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{W}\end{array}\mathrm{\left(}{\mathit{k}}_{\mathrm{3}}{\mathit{q}}_{\mathit{M}}\mathrm{\right)}$. Going back to Eq. (4), we can see that the term $\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{W}\end{array}\mathrm{\left(}{\mathit{k}}_{\mathrm{3}}{\mathit{q}}_{\mathit{M}}\mathrm{\right)}$ arises from a factor δ_D(k₃), so that in the regime where the 1-halo contribution is dominant, that is $\mathrm{⟨}\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{\delta }\end{array}\mathrm{\left(}{{k}}_{\mathrm{1}}\mathrm{\right)}\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{\delta }\end{array}\mathrm{\left(}{{k}}_{\mathrm{2}}\mathrm{\right)}\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{\delta }\end{array}\mathrm{\left(}{{k}}_{\mathrm{3}}\mathrm{\right)}\mathrm{⟩}\mathrm{\simeq }\mathrm{⟨}\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{\delta }\end{array}\mathrm{\left(}{{k}}_{\mathrm{1}}\mathrm{\right)}\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{\delta }\end{array}\mathrm{\left(}{{k}}_{\mathrm{2}}\mathrm{\right)}\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{\delta }\end{array}\mathrm{\left(}{{k}}_{\mathrm{3}}\mathrm{\right)}{\mathrm{⟩}}_{\mathrm{1}\mathrm{H}}$, it should reduce to δ_D(k₃). In Eq. (23) this corresponds to the fact that $\mathrm{\left|}\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{W}\end{array}\mathrm{\right(}{\mathit{k}}_{\mathrm{3}}{\mathit{q}}_{\mathit{M}}\mathrm{\left)}\mathrm{\right|}\mathrm{\ll }\mathrm{1}$ for k₃ ≫ 1/q_M. Then, thanks to the prefactor δ_D(k₁ + k₂ + k₃) we can see that in this regime |k₁ + k₂| → 0 and e^{−i(k₁ + k₂)·Ψc} → 1. In order to satisfy these properties, we simply make the approximation Ψ^c = 0, that is, we neglect the displacements of halos. Then, performing the integral over q^c, which allows the factorization of the Dirac factor δ_D(k₁ + k₂ + k₃), we obtain the 1-halo contribution to the bispectrum as $\begin{array}{ccc}{\mathit{B}}_{\mathrm{1}\mathrm{H}}\mathrm{\left(}{\mathit{k}}_{\mathrm{1}}\mathit{,}{\mathit{k}}_{\mathrm{2}}\mathit{,}{\mathit{k}}_{\mathrm{3}}\mathrm{\right)}& \mathrm{=}& \int \frac{\mathrm{d}\mathit{\nu }}{\mathit{\nu }}\hspace{0.17em}\mathit{f}\mathrm{\left(}\mathit{\nu }\mathrm{\right)}\hspace{0.17em}{\left(\frac{\mathit{M}}{\overline{\mathit{\rho }}\mathrm{(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathrm{3}}}\right)}^{\mathrm{2}}\\ & & \mathrm{\times }\hspace{0.17em}\underset{\mathit{j}\mathrm{=}\mathrm{1}}{\overset{\mathrm{3}}{\mathrm{􏽙}}}\mathrm{(}{\begin{array}{}\mathrm{˜}\\ \mathit{u}\end{array}}_{\mathit{M}}\mathrm{\left(}{\mathit{k}}_{\mathit{j}}\mathrm{\right)}\mathrm{-}\begin{array}{}\mathrm{˜}\\ \mathit{W}\end{array}\mathrm{(}{\mathit{k}}_{\mathit{j}}{\mathit{q}}_{\mathit{M}}{\mathrm{)}}^{\mathrm{)}}\mathit{.}\end{array}$(24)Thus, as for the power spectrum studied in Valageas & Nishimichi (2011), we recover the standard expression of the 1-halo contribution to the bispectrum derived within an Eulerian framework, see Cooray & Sheth (2002); Scoccimarro et al. (2001), with the addition of the new counterterms associated with the factors $\begin{array}{}{\mathrm{˜}}_{}\\ \mathit{W}\end{array}$. Note that for the power spectrum we obtained in Valageas & Nishimichi (2011) a factor $\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}\mathrm{\left(}\mathit{k}{\mathit{q}}_{\mathit{M}}{\mathrm{)}}^{\mathrm{2}}\mathrm{\right)}$, instead of the factor $\mathrm{\left[}{\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{u}\end{array}}_{\mathit{M}}\mathrm{\right(}\mathit{k}\mathrm{)}\mathrm{-}\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{W}\end{array}\mathrm{(}\mathit{k}{\mathit{q}}_{\mathit{M}}\mathrm{)}{\mathrm{]}}^{\mathrm{2}}$ that we would obtain using the approach described above. The difference arises because for the power spectrum we used Eq. (7), where the Dirac prefactor δ_D(k₁ + k₂) has already been factorized, whereas for the bispectrum we used Eq. (9) instead of Eq. (8). The advantage of Eqs. (7) and (8) is that they only involve relative positions, Δx, so that Eulerian position x^c of the halo never appears (for particles located in the same halo) and we do not need to introduce any approximation for the halo displacement Ψ^c.

However, for higher-order functions, such as the bispectrum or the trispectrum, the symmetric Eq. (9) (and its N-point generalization) provides a simpler starting point. In particular, the counting of the volumes over  {q₁,...,q_N}  associated with a single halo simply factorizes as ${\mathit{V}}_{\mathit{M}}^{\mathit{N}}$, whereas introducing relative positions makes the geometrical countings slightly more intricate.

The factors $\mathrm{\left[}{\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{u}\end{array}}_{\mathit{M}}\mathrm{\right(}{\mathit{k}}_{\mathit{j}}\mathrm{)}\mathrm{-}\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{W}\end{array}\mathrm{(}{\mathit{k}}_{\mathit{j}}{\mathit{q}}_{\mathit{M}}\mathrm{\left)}\mathrm{\right]}$ cannot be interpreted as modifications of the halo profiles ρ_M(x) under the form of compensated profiles (the result ${\mathit{\rho }}_{\mathit{M}}\mathrm{\left(}\mathit{x}\mathrm{\right)}\mathrm{-}\overline{\mathit{\rho }}$ would actually be negative at large x). Indeed, both quantities do not have the same radius in real space, and actually apply to two different spaces, namely the Eulerian and Lagrangian spaces.

A second important feature is that the expression (11) automatically ensures that the 1-halo contribution decays at least as ${\mathit{B}}_{\mathrm{1}\mathrm{H}}\mathrm{\left(}{\mathit{k}}_{\mathrm{1}}\mathit{,}{\mathit{k}}_{\mathrm{2}}\mathit{,}{\mathit{k}}_{\mathrm{3}}\mathrm{\right)}\mathrm{\propto }{\mathit{k}}_{\mathit{j}}^{\mathrm{2}}$ for any k_j → 0, since no linear dependence on k_j can remain because of statistical homogeneity and isotropy. This can be checked on the final expression (24), since ${\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{u}\end{array}}_{\mathit{M}}\mathrm{\left(}\mathrm{0}\mathrm{\right)}\mathrm{=}\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{W}\end{array}\mathrm{\left(}\mathrm{0}\mathrm{\right)}\mathrm{=}\mathrm{1}$ and the functions ${\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{u}\end{array}}_{\mathit{M}}\mathrm{\left(}\mathit{k}\mathrm{\right)}$ and $\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{W}\end{array}\mathrm{\left(}\mathit{k}{\mathit{q}}_{\mathit{M}}\mathrm{\right)}$ are regular functions of k² at the origin, ${\mathit{k}}_{\mathit{j}}\mathrm{\to }\mathrm{0}\mathrm{:}{\mathit{B}}_{\mathrm{1}\mathrm{H}}\mathrm{\left(}{\mathit{k}}_{\mathrm{1}}\mathit{,}{\mathit{k}}_{\mathrm{2}}\mathit{,}{\mathit{k}}_{\mathrm{3}}\mathrm{\right)}\mathrm{\propto }{\mathit{k}}_{\mathit{j}}^{\mathrm{2}}\mathit{.}$(25)This agrees with the requirements associated with a small-scale redistribution of matter. Here, there is a simplification as compared with the matter power spectrum P(k). Indeed, as discussed in Peebles (1974), small-scale redistributions of matter generically give a low-k tail $\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{\delta }\end{array}\mathrm{\left(}{k}\mathrm{\right)}\mathrm{\propto }\mathit{k}$, whence P(k) ∝ k², while taking into account momentum conservation implies the steeper decay $\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{\delta }\end{array}\mathrm{\left(}{k}\mathrm{\right)}\mathrm{\propto }{\mathit{k}}^{\mathrm{2}}$, whence P(k) ∝ k⁴. Thus, in Valageas & Nishimichi (2011) we recovered the expected k² tail, through the factor $\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}\mathrm{\left(}\mathit{k}{\mathit{q}}_{\mathit{M}}{\mathrm{)}}^{\mathrm{2}}\mathrm{\right)}$, since our analysis did not enforce momentum conservation. Using the approach described above, we would actually obtain a k⁴ tail, through the factor $\mathrm{\left[}{\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{u}\end{array}}_{\mathit{M}}\mathrm{\right(}\mathit{k}\mathrm{)}\mathrm{-}\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{W}\end{array}\mathrm{(}\mathit{k}{\mathit{q}}_{\mathit{M}}\mathrm{)}{\mathrm{]}}^{\mathrm{2}}$, once we make the approximation Ψ^c = 0. This clearly satisfies momentum conservation (the peculiar momentum of each halo is zero) but comes at the expense of an additional approximation and one should not give too much weight to this property. By contrast, for the bispectrum (and higher-order correlations) small-scale redistributions of matter lead to the same ${\mathit{k}}_{\mathit{j}}^{\mathrm{2}}$ tail, whether we take into account momentum conservation or not, because the linear term k_j vanishes by symmetry. Note that when all wavenumbers go to zero at the same rate, we have from Eq. (24) the decrease $\mathit{\lambda }\mathrm{\to }\mathrm{0}\mathrm{:}{\mathit{B}}_{\mathrm{1}\mathrm{H}}\mathrm{\left(}\mathit{\lambda }{\mathit{k}}_{\mathrm{1}}\mathit{,\lambda }{\mathit{k}}_{\mathrm{2}}\mathit{,\lambda }{\mathit{k}}_{\mathrm{3}}\mathrm{\right)}\mathrm{\propto }{\mathit{\lambda }}^{\mathrm{6}}\mathit{.}$(26)As noticed in previous works (Cooray & Sheth 2002), in the usual version of the halo model the 1-halo contribution to the bispectrum does not contain the counterterms $\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{W}\end{array}\mathrm{\left(}{\mathit{k}}_{\mathit{j}}{\mathit{q}}_{\mathit{M}}\mathrm{\right)}$, so that it goes to a strictly positive constant at large scales and dominates over the lowest-order prediction of perturbation theory. This led to rather inaccurate predictions on large scales (an overestimate of 20% at k ~ 0.01 h Mpc^-1), because for CDM cosmologies the perturbative prediction scales as B_pert ~ P_L(k)² ∝ k^2n_s with n_s ≃ 1 (for k_j ~ k for all j). As seen from Eqs. (25), (26), our approach solves this problem (we will return to this point in Sect. 5).

2.3. “2-halo” contribution

We now turn to the 2-halo contribution to the bispectrum, following the procedure described in the previous section for the 1-halo term. Thus, as in Eq. (12) we can write the 2-halo contribution as $\begin{array}{ccc}\mathrm{⟨}\begin{array}{}\mathrm{˜}\\ \mathit{\delta }\end{array}\mathrm{\left(}{{k}}_{\mathrm{1}}\mathrm{\right)}\begin{array}{}\mathrm{˜}\\ \mathit{\delta }\end{array}\mathrm{\left(}{{k}}_{\mathrm{2}}\mathrm{\right)}\begin{array}{}\mathrm{˜}\\ \mathit{\delta }\end{array}\mathrm{\left(}{{k}}_{\mathrm{3}}\mathrm{\right)}{\mathrm{⟩}}_{\mathrm{2}\mathrm{H}}& \mathrm{=}& {⟨}^{\mathrm{\int }}\underset{\mathit{j}\mathrm{=}\mathrm{1}}{\overset{\mathrm{3}}{\mathrm{􏽙}}}\frac{\mathrm{d}{{q}}_{\mathit{j}}}{\mathrm{(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathrm{3}}}\hspace{0.17em}\mathrm{(}{\mathrm{e}}^{\mathrm{-}\mathrm{i}{{k}}_{\mathit{j}}\mathrm{·}{{x}}_{\mathit{j}}}\mathrm{-}{\mathrm{e}}^{\mathrm{-}\mathrm{i}{{k}}_{\mathit{j}}\mathrm{·}{{q}}_{\mathit{j}}}\mathrm{)}\\ & & \mathrm{\times }\hspace{0.17em}{\mathrm{\int }}_{{{q}}_{\mathrm{1}}^{\mathrm{c}}\ne {{q}}_{\mathrm{2}}^{\mathrm{c}}}\mathrm{d}{{q}}_{\mathrm{1}}^{\mathrm{c}}\mathrm{d}{\mathit{M}}_{\mathrm{1}}\hspace{0.17em}\mathrm{d}{{q}}_{\mathrm{2}}^{\mathrm{c}}\mathrm{d}{\mathit{M}}_{\mathrm{2}}\hspace{0.17em}\mathit{n̂}\mathrm{\left(}{{q}}_{\mathrm{1}}^{\mathrm{c}}\mathit{,}{\mathit{M}}_{\mathrm{1}}\mathrm{\right)}\hspace{0.17em}\mathit{n̂}\mathrm{\left(}{{q}}_{\mathrm{2}}^{\mathrm{c}}\mathit{,}{\mathit{M}}_{\mathrm{2}}\mathrm{\right)}\\ & & \mathrm{\times }\hspace{0.17em}\mathit{\theta }\mathrm{(}{{q}}_{\mathrm{1}}\mathrm{\in }{\mathit{M}}_{\mathrm{1}}\mathrm{)}\hspace{0.17em}\mathit{\theta }\mathrm{(}{{q}}_{\mathrm{2}}\mathrm{\in }{\mathit{M}}_{\mathrm{2}}\mathrm{)}\hspace{0.17em}\mathit{\theta }\mathrm{(}{{q}}_{\mathrm{3}}\mathrm{\in }{\mathit{M}}_{\mathrm{2}}\mathrm{)}⟩\mathrm{+}\mathrm{2}\hspace{0.17em}\mathrm{cyc}\mathit{.}\mathit{,}\end{array}$(27)where we sum over all pairs of distinct halos, of mass M₁ and M₂, which contain one and two of the particles q_j. Here we have written the contribution where the halo M₁ contains the particle q₁, and the two additional contributions, noted “2 cyc.”, correspond to the cases where M₁ contains the particle q₂ or q₃. Using again the approximation (22) of fully relaxed halos, the average of Eq. (27) writes as $\begin{array}{ccc}\mathrm{⟨}\begin{array}{}\mathrm{˜}\\ \mathit{\delta }\end{array}\mathrm{\left(}{{k}}_{\mathrm{1}}\mathrm{\right)}\begin{array}{}\mathrm{˜}\\ \mathit{\delta }\end{array}\mathrm{\left(}{{k}}_{\mathrm{2}}\mathrm{\right)}\begin{array}{}\mathrm{˜}\\ \mathit{\delta }\end{array}\mathrm{\left(}{{k}}_{\mathrm{3}}\mathrm{\right)}{\mathrm{⟩}}_{\mathrm{2}\mathrm{H}}& \mathrm{=}& \mathrm{\int }\mathrm{d}{{q}}_{\mathrm{1}}^{\mathrm{c}}\frac{\mathrm{d}{\mathit{\nu }}_{\mathrm{1}}}{{\mathit{\nu }}_{\mathrm{1}}}\hspace{0.17em}\mathrm{d}{{q}}_{\mathrm{2}}^{\mathrm{c}}\frac{\mathrm{d}{\mathit{\nu }}_{\mathrm{2}}}{{\mathit{\nu }}_{\mathrm{2}}}\hspace{0.17em}\frac{{\overline{\mathit{\rho }}}^{\mathrm{2}}}{{\mathit{M}}_{\mathrm{1}}{\mathit{M}}_{\mathrm{2}}}\mathit{f}\mathrm{\left(}{\mathit{\nu }}_{\mathrm{1}}\mathrm{\right)}\mathit{f}\mathrm{\left(}{\mathit{\nu }}_{\mathrm{2}}\mathrm{\right)}\\ & & \mathrm{\times }\hspace{0.17em}\left(\mathrm{1}\mathrm{+}{\mathit{\xi }}_{\mathrm{12}}\right)\frac{{\mathit{M}}_{\mathrm{1}}{\mathit{M}}_{\mathrm{2}}^{\mathrm{2}}}{{\overline{\mathit{\rho }}}^{\mathrm{3}}\mathrm{(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathrm{9}}}\hspace{0.17em}{\mathrm{e}}^{\mathrm{-}\mathrm{i}{{k}}_{\mathrm{1}}\mathrm{·}{{q}}_{\mathrm{1}}^{\mathrm{c}}}\hspace{0.17em}{\mathrm{e}}^{\mathrm{-}\mathrm{i}\mathrm{(}{{k}}_{\mathrm{2}}\mathrm{+}{{k}}_{\mathrm{3}}\mathrm{)}\mathrm{·}{{q}}_{\mathrm{2}}^{\mathrm{c}}}\\ & & \mathrm{\times }\hspace{0.17em}\mathrm{(}{\mathrm{e}}^{\mathrm{-}\mathrm{i}{{k}}_{\mathrm{1}}\mathrm{·}{\mathbf{\Psi }}_{\mathrm{1}}^{\mathrm{c}}}\hspace{0.17em}{\begin{array}{}\mathrm{˜}\\ \mathit{u}\end{array}}_{{\mathit{M}}_{\mathrm{1}}}\mathrm{\left(}{\mathit{k}}_{\mathrm{1}}\mathrm{\right)}\mathrm{-}\begin{array}{}\mathrm{˜}\\ \mathit{W}\end{array}\mathrm{(}{\mathit{k}}_{\mathrm{1}}{\mathit{q}}_{\mathit{M}\mathrm{1}}{\mathrm{)}}^{\mathrm{)}}\\ & & \mathrm{\times }\hspace{0.17em}\underset{\mathit{j}\mathrm{=}\mathrm{2}}{\overset{\mathrm{3}}{\mathrm{􏽙}}}\mathrm{(}{\mathrm{e}}^{\mathrm{-}\mathrm{i}{{k}}_{\mathit{j}}\mathrm{·}{\mathbf{\Psi }}_{\mathrm{2}}^{\mathrm{c}}}\hspace{0.17em}{\begin{array}{}\mathrm{˜}\\ \mathit{u}\end{array}}_{{\mathit{M}}_{\mathrm{2}}}\mathrm{\left(}{\mathit{k}}_{\mathit{j}}\mathrm{\right)}\mathrm{-}\begin{array}{}\mathrm{˜}\\ \mathit{W}\end{array}\mathrm{(}{\mathit{k}}_{\mathit{j}}{\mathit{q}}_{\mathit{M}\mathrm{2}}{\mathrm{)}}^{\mathrm{)}}\mathrm{+}\mathrm{2}\hspace{0.17em}\mathrm{cyc}\mathit{.}\mathit{,}\end{array}$(28)where ${\mathit{\xi }}_{\mathrm{12}}\mathrm{\left(}\mathrm{\right|}{{q}}_{\mathrm{2}}^{\mathrm{c}}\mathrm{-}{{q}}_{\mathrm{1}}^{\mathrm{c}}\mathrm{\left|}\mathrm{\right)}$ is the two-point correlation function of halos of mass M₁ and M₂, in Lagrangian space. As a first step, let us neglect halo motions, Ψ^c = 0. Then, as expected, the contribution associated with the factor 1 in (1 + ξ₁₂) vanishes (because the integral over ${{q}}_{\mathrm{1}}^{\mathrm{c}}$ yields δ_D(k₁) and $\mathrm{\left[}{\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{u}\end{array}}_{{\mathit{M}}_{\mathrm{1}}}\mathrm{\right(}{\mathit{k}}_{\mathrm{1}}\mathrm{)}\mathrm{-}\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{W}\end{array}\mathrm{(}{\mathit{k}}_{\mathrm{1}}{\mathit{q}}_{\mathit{M}\mathrm{1}}\mathrm{\left)}\mathrm{\right]}\mathrm{=}\mathrm{0}$ for k₁ = 0), so that only the term associated with the large-scale correlation of halos remains. Writing this two-point correlation in terms of the halo power spectrum, ${\mathit{\xi }}_{\mathrm{12}}\mathrm{\left(}\mathrm{\right|}{{q}}_{\mathrm{2}}^{\mathrm{c}}\mathrm{-}{{q}}_{\mathrm{1}}^{\mathrm{c}}\mathrm{\left|}\mathrm{\right)}\mathrm{=}\mathrm{\int }\mathrm{d}{k}\hspace{0.17em}{\mathrm{e}}^{\mathrm{i}{k}\mathrm{·}\mathrm{(}{{q}}_{\mathrm{2}}^{\mathrm{c}}\mathrm{-}{{q}}_{\mathrm{1}}^{\mathrm{c}}\mathrm{)}}\hspace{0.17em}{\mathit{P}}_{\mathrm{12}}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathit{,}$(29)the integrals over ${{q}}_{\mathrm{1}}^{\mathrm{c}}$ and ${{q}}_{\mathrm{2}}^{\mathrm{c}}$ give the Dirac factors δ_D(k + k₁)δ_D(k − k₂ − k₃). Therefore, we recover the Dirac prefactor δ_D(k₁ + k₂ + k₃), which we can factorize out to write the 2-halo contribution to the bispectrum as $\begin{array}{ccc}{\mathit{B}}_{\mathrm{2}\mathrm{H}}\mathrm{\left(}{\mathit{k}}_{\mathrm{1}}\mathit{,}{\mathit{k}}_{\mathrm{2}}\mathit{,}{\mathit{k}}_{\mathrm{3}}\mathrm{\right)}& \mathrm{=}& \int \frac{\mathrm{d}{\mathit{\nu }}_{\mathrm{1}}}{{\mathit{\nu }}_{\mathrm{1}}}\frac{\mathrm{d}{\mathit{\nu }}_{\mathrm{2}}}{{\mathit{\nu }}_{\mathrm{2}}}\frac{{\mathit{M}}_{\mathrm{2}}}{\overline{\mathit{\rho }}\mathrm{(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathrm{3}}}\hspace{0.17em}\mathit{f}\mathrm{\left(}{\mathit{\nu }}_{\mathrm{1}}\mathrm{\right)}\mathit{f}\mathrm{\left(}{\mathit{\nu }}_{\mathrm{2}}\mathrm{\right)}\hspace{0.17em}{\mathit{P}}_{\mathrm{12}}\mathrm{\left(}{\mathit{k}}_{\mathrm{1}}\mathrm{\right)}\\ & & \mathrm{\times }\hspace{0.17em}\mathrm{(}{\begin{array}{}\mathrm{˜}\\ \mathit{u}\end{array}}_{{\mathit{M}}_{\mathrm{1}}}\mathrm{\left(}{\mathit{k}}_{\mathrm{1}}\mathrm{\right)}\mathrm{-}\begin{array}{}\mathrm{˜}\\ \mathit{W}\end{array}\mathrm{(}{\mathit{k}}_{\mathrm{1}}{\mathit{q}}_{\mathit{M}\mathrm{1}}{\mathrm{)}}^{\mathrm{)}}\underset{\mathit{j}\mathrm{=}\mathrm{2}}{\overset{\mathrm{3}}{\mathrm{􏽙}}}\mathrm{(}{\begin{array}{}\mathrm{˜}\\ \mathit{u}\end{array}}_{{\mathit{M}}_{\mathrm{2}}}\mathrm{\left(}{\mathit{k}}_{\mathit{j}}\mathrm{\right)}\mathrm{-}\begin{array}{}\mathrm{˜}\\ \mathit{W}\end{array}\mathrm{(}{\mathit{k}}_{\mathit{j}}{\mathit{q}}_{\mathit{M}\mathrm{2}}{\mathrm{)}}^{\mathrm{)}}\\ & & \mathrm{+}\mathrm{2}\hspace{0.17em}\mathrm{cyc}\mathit{.}\end{array}$(30)Then, as is customary, we could write the halo power spectrum as P₁₂(k) = b(M₁)b(M₂)P(k), where b(M) is a mass-dependent bias factor (of order unity for typical halos) that we approximate as scale-independent (thereby neglecting exclusion constraints between halos).

However, the expression (30) is not really satisfactory. Indeed, for k₁ → 0 it decreases as ${\mathit{k}}_{\mathrm{1}}^{\mathrm{2}}\mathit{P}\mathrm{\left(}{\mathit{k}}_{\mathrm{1}}\mathrm{\right)}$ because of the prefactor $\mathrm{\left[}{\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{u}\end{array}}_{{\mathit{M}}_{\mathrm{1}}}\mathrm{\right(}{\mathit{k}}_{\mathrm{1}}\mathrm{)}\mathrm{-}\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{W}\end{array}\mathrm{(}{\mathit{k}}_{\mathrm{1}}{\mathit{q}}_{\mathit{M}\mathrm{1}}\mathrm{\left)}\mathrm{\right]}$, while we would expect to recover a behavior proportional to P(k₁). Indeed, in this regime this 2-halo contribution should describe the power that arises from the large-scale correlation between density fluctuations at a position q₁ and at a position q₂, with |q₂ − q₁| ~ 2π/k₁, the fluctuations at q₂ being furthermore decomposed over smaller-scale fluctuations within single halos of mass M₂. In other words, the large-scale power P(k₁) should not be damped by the prefactor $\mathrm{\left[}{\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{u}\end{array}}_{{\mathit{M}}_{\mathrm{1}}}\mathrm{\right(}{\mathit{k}}_{\mathrm{1}}\mathrm{)}\mathrm{-}\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{W}\end{array}\mathrm{(}{\mathit{k}}_{\mathrm{1}}{\mathit{q}}_{\mathit{M}\mathrm{1}}\mathrm{\left)}\mathrm{\right]}$. Going back to Eq. (28), we can see that the product ${\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{u}\end{array}}_{{\mathit{M}}_{\mathrm{1}}}\mathrm{\left(}{\mathit{k}}_{\mathrm{1}}\mathrm{\right)}{\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{u}\end{array}}_{{\mathit{M}}_{\mathrm{2}}}\mathrm{\left(}{\mathit{k}}_{\mathrm{2}}\mathrm{\right)}{\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{u}\end{array}}_{{\mathit{M}}_{\mathrm{2}}}\mathrm{\left(}{\mathit{k}}_{\mathrm{3}}\mathrm{\right)}$ involves a prefactor ${\mathrm{e}}^{\mathrm{i}{{k}}_{\mathrm{1}}\mathrm{·}\mathrm{(}{\mathbf{\Psi }}_{\mathrm{2}}^{\mathrm{c}}\mathrm{-}{\mathbf{\Psi }}_{\mathrm{1}}^{\mathrm{c}}\mathrm{)}}$, which only depends on relative displacements and as such is physically meaningful and should be taken into account. This should be contrasted with the behavior obtained for the 1-halo contribution (23), where the product of three terms ${\begin{array}{}{\mathrm{˜}}_{}\\ \mathit{u}\end{array}}_{\mathit{M}}$ gave the prefactor e^{i(k₁ + k₂ + k₃)·Ψc}, which had to disappear because it does not depend on relative displacements, and was indeed irrelevant thanks to the Dirac prefactor δ_D(k₁ + k₂ + k₃). This means that the approximation Ψ^c = 0 is not as good for the 2-halo contribution as for the 1-halo contribution, because it neglects relative halo displacements, which did not appear in the latter case (i.e. its consequences are stronger in the former case).

In order to have a 2-halo contribution that behaves in a reasonable fashion, we choose to use instead of Eq. (30) the simple expression $\begin{array}{ccc}{\mathit{B}}_{\mathrm{2}\mathrm{H}}\mathrm{\left(}{\mathit{k}}_{\mathrm{1}}\mathit{,}{\mathit{k}}_{\mathrm{2}}\mathit{,}{\mathit{k}}_{\mathrm{3}}\mathrm{\right)}& \mathrm{=}& {\mathit{P}}_{\mathrm{L}}\mathrm{\left(}{\mathit{k}}_{\mathrm{1}}\mathrm{\right)}\mathrm{\int }\frac{\mathrm{d}\mathit{\nu }}{\mathit{\nu }}\frac{\mathit{M}}{\overline{\mathit{\rho }}\mathrm{(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathrm{3}}}\hspace{0.17em}\mathit{f}\mathrm{\left(}\mathit{\nu }\mathrm{\right)}\\ & & \mathrm{\times }\hspace{0.17em}\underset{\mathit{j}\mathrm{=}\mathrm{2}}{\overset{\mathrm{3}}{\mathrm{􏽙}}}\mathrm{(}{\begin{array}{}\mathrm{˜}\\ \mathit{u}\end{array}}_{\mathit{M}}\mathrm{\left(}{\mathit{k}}_{\mathit{j}}\mathrm{\right)}\mathrm{-}\begin{array}{}\mathrm{˜}\\ \mathit{W}\end{array}\mathrm{(}{\mathit{k}}_{\mathit{j}}{\mathit{q}}_{\mathit{M}}{\mathrm{)}}^{\mathrm{)}}\mathrm{+}\mathrm{2}\hspace{0.17em}\mathrm{cyc}\mathit{.}\end{array}$(31)Thus, we have made the approximation P₁₂(k) ≃ P_L(k), where P_L(k) is the linear matter density power spectrum, and we have set the spurious prefactor $\mathrm{\left[}{\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{u}\end{array}}_{{\mathit{M}}_{\mathrm{1}}}\mathrm{\right(}{\mathit{k}}_{\mathrm{1}}\mathrm{)}\mathrm{-}\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{W}\end{array}\mathrm{(}{\mathit{k}}_{\mathrm{1}}{\mathit{q}}_{\mathit{M}\mathrm{1}}\mathrm{\left)}\mathrm{\right]}$ to unity. In principle, instead of the linear power P_L(k) we could include higher orders of perturbation theory. However, in view of the approximate nature of Eq. (31) this is not really necessary, especially since the 2-halo contribution is subdominant on most scales of interest, as shown by the numerical results in Sect. 4.

In the large-scale limit we obtain from Eq. (31) the asymptotic behaviors ${\mathit{k}}_{\mathit{j}}\mathrm{\to }\mathrm{0}\mathrm{:}{\mathit{B}}_{\mathrm{2}\mathrm{H}}\mathrm{\left(}{\mathit{k}}_{\mathrm{1}}\mathit{,}{\mathit{k}}_{\mathrm{2}}\mathit{,}{\mathit{k}}_{\mathrm{3}}\mathrm{\right)}\mathrm{\propto }{\mathit{P}}_{\mathrm{L}}\mathrm{\left(}{\mathit{k}}_{\mathit{j}}\mathrm{\right)}\mathit{,}$(32)when only one wavenumber goes to zero, and $\mathit{\lambda }\mathrm{\to }\mathrm{0}\mathrm{:}{\mathit{B}}_{\mathrm{2}\mathrm{H}}\mathrm{\left(}\mathit{\lambda }{\mathit{k}}_{\mathrm{1}}\mathit{,\lambda }{\mathit{k}}_{\mathrm{2}}\mathit{,\lambda }{\mathit{k}}_{\mathrm{3}}\mathrm{\right)}\mathrm{\propto }{\mathit{\lambda }}^{\mathrm{4}}{\mathit{P}}_{\mathrm{L}}\mathrm{\left(}\mathit{\lambda }\mathrm{\right)}\mathit{,}$(33)when all wavenumbers decrease at the same rate. In Eq. (32) we used the fact that for CDM cosmologies we have P_L(k) ~ k^n_s with n_s ≃ 1 at low k. It is interesting to note that the behavior (32) agrees with perturbation theory, as shown in Appendix A and Eq. (A.5). There, we note that this scaling, which applies when only one of the wavenumbers k_j goes to zero, for instance k₁, is valid at all orders of perturbation theory. In this regime, the 2-halo contribution, which is non-perturbative (as discussed in Sect. 2.4 below), is therefore on the same order (i.e.,  ~P_L(k₁)) as the perturbative contribution. This can be understood from the arguments developed above to obtain the expression (31), as the non-perturbative effects associated with the formation of small-scale nonlinear structures around the nearby points x₂ and x₃ (or q₂ and q₃ in Lagrangian space) do not modify the larger-scale correlation with the very distant point x₁. That is, the non-perturbative effects associated with the formation of a halo correspond to a local redistribution of matter. Then, at leading order, non-perturbative terms simply lead to a “renormalization” of the coefficient that multiplies the factor P_L(k₁), by an amount that depends on how far in the nonlinear regime the wavenumbers k₂ and k₃ are.

When all wavenumbers go to zero, the property (33) again ensures that the 2-halo contribution becomes negligible as compared with the lowest order term of perturbation theory, which scales as P_L(λ)². Thus, as for the 1-halo contribution, the counterterms $\begin{array}{}{\mathrm{˜}}_{}\\ \mathit{W}\end{array}$ of Eq. (31) solve the problem encountered with the usual implementation of the halo model, which does not recover perturbation theory on very large scales.

2.4. “3-halo” contribution as a perturbative contribution

For the 3-halo contribution we follow the spirit of the approach described in Valageas & Nishimichi (2011), as we bypass the halo model to make contact with standard perturbation theory. Indeed, as we have seen above for the 2-halo term, it is not straightforward to take into account large-scale correlations within this Lagrangian implementation of the halo model, starting with the expression (9). This would require some modeling of relative displacements, which may not be worth the effort (if one wishes to improve in a significant manner over the approximation (34) used below). Thus, we take a much more direct and simpler route, noting that the 1-halo and 2-halo contributions vanish at all orders of perturbation theory. Here, we mean by perturbation theory any expansion over powers of the linear power spectrum P_L (for the Gaussian initial conditions that we consider in this article), such as the standard perturbation theory derived within the fluid approximation to the equations of motion (Goroff et al. 1986; Bernardeau et al. 2002a). Indeed, the halo mass function shows a large-mass tail of the form ${\mathrm{e}}^{\mathrm{-}{\mathit{\nu }}^{\mathrm{2}}\mathit{/}\mathrm{2}}\mathrm{=}{\mathrm{e}}^{\mathrm{-}{\mathit{\delta }}_{\mathrm{L}}^{\mathrm{2}}\mathit{/}\mathrm{\left(}\mathrm{2}{\mathit{\sigma }}^{\mathrm{2}}\mathrm{\right(}\mathit{M}\mathrm{\left)}\mathrm{\right)}}$, which vanishes exponentially for P_L → 0 and has a Taylor expansion over powers of P_L that is identically zero. Therefore, the remaining 3-halo contribution to Eq. (10) must be consistent with perturbation theory up to all orders, and we simply use the approximation ${\mathit{B}}_{\mathrm{3}\mathrm{H}}\mathrm{\left(}{\mathit{k}}_{\mathrm{1}}\mathit{,}{\mathit{k}}_{\mathrm{2}}\mathit{,}{\mathit{k}}_{\mathrm{3}}\mathrm{\right)}\mathrm{=}{\mathit{B}}_{\mathrm{pert}}\mathrm{\left(}{\mathit{k}}_{\mathrm{1}}\mathit{,}{\mathit{k}}_{\mathrm{2}}\mathit{,}{\mathit{k}}_{\mathrm{3}}\mathrm{\right)}\mathit{,}$(34)where B_pert is the bispectrum obtained from perturbation theory. To follow more closely the approach used in Valageas & Nishimichi (2011) for the power spectrum, we should multiply B_pert in Eq. (34) by a factor F_3H, with 0 < F_3H < 1, which counts the fraction of triplets that are located within three distinct halos. From the same arguments we have F_3H = 1 at all orders of perturbation theory, and it only differs from unity by non-perturbative terms such as ${\mathrm{e}}^{\mathrm{-}{\mathit{\delta }}_{\mathrm{L}}^{\mathrm{2}}\mathit{/}\mathrm{\left(}\mathrm{2}{\mathit{\sigma }}^{\mathrm{2}}\mathrm{\right(}\mathit{M}\mathrm{\left)}\mathrm{\right)}}$, where M is typically the mass scale associated with the maximum of  {2π/k_j} . For illustration purposes, we display in Fig. 1 the following estimate of F_3H, ${\mathit{F}}_{\mathrm{3}\mathrm{H}}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathrm{=}{\left({\mathrm{\int }}_{\mathrm{0}}^{{\mathit{\nu }}_{\mathit{k}}}\frac{\mathrm{d}\mathit{\nu }}{\mathit{\nu }}\hspace{0.17em}\mathit{f}\mathrm{\left(}\mathit{\nu }\mathrm{\right)}\right)}^{\mathrm{3}}\mathit{,}$(35)where ν_k = δ_L/σ(q_k) with q_k = 2π/k. From Eq. (15), this is the probability that three particles belong to halos of size smaller than q_k, neglecting halo correlations, and it should give an estimate of the range where F_3H ≃ 1 (for the equilateral bispectrum, which involves a single wavenumber k). As expected, the comparison with Fig. 5 below shows that the departure from unity of F_3H roughly corresponds to the scale where the 1-halo and 2-halo terms become of the same order as the 3-halo term. In view of the approximations involved in these 1-halo and 2-halo contributions, we do not try to include the deviations from unity of F_3H, which only play a role in the transition regime. Therefore, we make the approximation F_3H ≃ 1 and use the simple prescription (34).

In practice, the perturbative term B_pert is not exactly known (assuming the perturbative series is convergent), and one must truncate the perturbative expansion. There are many ways to do so, because a priori one can use the standard perturbation theory or any alternative expansion scheme, which corresponds to partial resummations of the standard diagrams. As shown in Valageas & Nishimichi (2011), for the power spectrum it happens that standard perturbation theory is not a good choice, because higher order terms grow increasingly fast on small scales and prevent a good description of the highly nonlinear regime. Then, one must use resummation schemes that agree with standard perturbation theory up to the required order while remaining well-behaved on small scales. As we will check in Sect. 4, for the bispectrum it turns out that at one-loop order the standard perturbation theory remains small at high k compared with the 1-halo and 2-halo contributions, which makes this an efficient choice; but we will also study the direct steepest-descent resummation scheme described in detail in Valageas & Nishimichi (2011) for the computation of the power spectrum.

In most studies that are based on the halo model, the 3-halo term is rather written as (Ma & Fry 2000b; Scoccimarro et al. 2001; Cooray & Sheth 2002) ${\mathit{B}}_{\mathrm{3}\mathrm{H}}^{\mathrm{h}\mathit{.}\mathrm{m}\mathit{.}}\mathrm{\left(}{\mathit{k}}_{\mathrm{1}}\mathit{,}{\mathit{k}}_{\mathrm{2}}\mathit{,}{\mathit{k}}_{\mathrm{3}}\mathrm{\right)}\mathrm{=}\left[\underset{\mathit{i}\mathrm{=}\mathrm{1}}{\overset{\mathrm{3}}{\mathrm{􏽙}}}\mathrm{⟨}\mathit{b}\mathrm{\left(}\mathit{M}\mathrm{\right)}{\begin{array}{}\mathrm{˜}\\ \mathit{u}\end{array}}_{\mathit{M}}\mathrm{\left(}{\mathit{k}}_{\mathit{i}}\mathrm{\right)}\mathrm{⟩}\right]{\mathit{B}}_{\mathrm{tree}}\mathrm{\left(}{\mathit{k}}_{\mathrm{1}}\mathit{,}{\mathit{k}}_{\mathrm{2}}\mathit{,}{\mathit{k}}_{\mathrm{3}}\mathrm{\right)}\mathit{,}$(36)where $\mathrm{⟨}\mathit{b}\mathrm{\left(}\mathit{M}\mathrm{\right)}{\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{u}\end{array}}_{\mathit{M}}\mathrm{⟩}$ is the average over mass, weighted by the halo mass function, of the bias parameter b and halo profile ${\begin{array}{}{\mathrm{˜}}_{}\\ \mathit{u}\end{array}}_{\mathit{M}}$; B_tree is the matter bispectrum obtained at lowest order through perturbation theory, given by Eq. (39) below. In our approach, we do not introduce this halo bias parameter because we simply write the 3-halo term as the perturbative matter bispectrum (34). This implies that we focus on the matter bispectrum because our approach does not contain the halo bispectrum ${\mathit{B}}_{\mathrm{3}\mathrm{H}}^{\mathrm{h}\mathit{.}\mathrm{m}\mathit{.}}\mathrm{\left(}{\mathit{k}}_{\mathrm{1}}\mathit{,}{\mathit{k}}_{\mathrm{2}}\mathit{,}{\mathit{k}}_{\mathrm{3}}\mathrm{;}{\mathit{M}}_{\mathrm{1}}\mathit{,}{\mathit{M}}_{\mathrm{2}}\mathit{,}{\mathit{M}}_{\mathrm{3}}\mathrm{\right)}$ as an intermediate tool. Therefore, to compute the bispectrum of halos of given masses we would need to add another prescription. On the other hand, this makes our model for the matter bispectrum simpler and more robust because it does not rely on bias parameters that are fairly difficult to compute in a systematic and well-controlled manner (because halos themselves do not appear in a natural fashion in the equations of motion). Moreover, Eq. (34) ensures that our results agree with perturbation theory up to the order of truncation. In this paper, we will go up to 1-loop order, while most previous studies that involved the halo model only used the tree-level contribution, as in Eq. (36).

2.5. Comparison with the Eulerian framework

As seen in the previous sections, the derivation of the bispectrum from a Lagrangian implementation of the halo model is not as straightforward nor as “clean” as for the power spectrum (Valageas & Nishimichi 2011). Indeed, we had to introduce some additional approximation regarding the halo displacements Ψ^c, and for the 2-halo term that mixes larger-scale and intra-halo wavenumbers we had to partially bypass the naive prediction of this halo model to recover the large-scale behavior (32).

More generally, it is not surprising that the Lagrangian implementation of the halo model requires more steps than the usual Eulerian version. Indeed, within an Eulerian framework we only need to describe the density field at the time of interest t. Within a Lagrangian framework, we need to add some information on the building of this density field, that is, the mapping between the initial coordinates q and the current positions x of the particles. This necessarily implies that a practical Lagrangian implementation requires more hypothesis or approximations (because there are more intermediate quantities, even though in principle the results would be the same if we made no approximation at all). However, we think it remains of interest to describe how N-point correlation functions or poly-spectra can be constructed within such a Lagrangian framework, based on the halo model, as we have done above. First, it is useful to see how identical or similar approximate models can be derived from different points of view, since this gives more weight to the results. Second, it provides additional insight on the underlying approximations and it could offer an alternative route to more precise modeling. Third, it provides a natural derivation of the counterterms $\begin{array}{}{\mathrm{˜}}_{}\\ \mathit{W}\end{array}$ in the 1-halo contribution (24), which ensure the large-scale behaviors (25), (26) that were missed in previous Eulerian implementations.

As explained above, for the 2-halo contribution we had to take some liberty with a naive implementation of the halo model to recover the asymptotic scaling (32). However, in view of the approximate nature of these Lagrangian halo models, we think it is best to be guided by physical arguments and to make sure physical constraints are satisfied, instead of strictly adhering to approximate models that can show a varying degree of accuracy, depending on the quantities under scrutiny.

However, as in Valageas & Nishimichi (2011), it is interesting to note that the counterterms of Eqs. (24) and (31) might also be guessed within an Eulerian framework, from the requirement that the bispectrum should vanish for a perfectly homogeneous universe. Indeed, with the usual halo-model approximations, we can still split a constant-density medium over the same set of “halos”, or cells, of radius q_M and center q^c (neglecting geometrical constraints associated with the impossibility of covering a 3D space with spherical cells). The only difference is that because these “halos” have not changed and have remained at the constant density $\overline{\mathit{\rho }}$, their Eulerian and Lagrangian properties are the same. In particular, the Eulerian radius x_M is equal to q_M, and the halo profiles are simply given by ${\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)}$ for this uniform system. Then, the counterterms in Eqs. (24) and (31) clearly satisfy the constraint B(k₁,k₂,k₃) = 0 for this homogeneous universe. However, by itself this argument is not sufficient to imply the precise form of Eq. (24) (nor of Eq. (31)), since the choice ${{}^{\mathrm{􏽑}}}_{\mathit{j}}{\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{u}\end{array}}_{\mathit{M}}\mathrm{\left(}{\mathit{k}}_{\mathit{j}}\mathrm{\right)}\mathrm{-}{{}^{\mathrm{􏽑}}}_{\mathit{j}}\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{W}\end{array}\mathrm{\left(}{\mathit{k}}_{\mathit{j}}{\mathit{q}}_{\mathit{M}}\mathrm{\right)}$ would also satisfy this constraint (but not the properties (25)).

Although the main reason for taking into account the counterterms $\begin{array}{}{\mathrm{˜}}_{}\\ \mathit{W}\end{array}$ in the 1-halo and 2-halo contributions is to satisfy physical requirements, this is also needed to reach high accuracy on weakly nonlinear scales. Thus, Guo & Jing (2009) find that on scales on the order of 0.1 h Mpc^-1 their 1-halo and 2-halo terms (which do not include these counterterms) are still too large and spoil the agreement of (tree-order) perturbation theory with numerical simulations. Giving a faster decay on large scales of these 1-halo and 2-halo contributions, the counterterms $\begin{array}{}{\mathrm{˜}}_{}\\ \mathit{W}\end{array}$ improve the agreement of the analytical predictions with the simulations, as we will check in Sect. 5.

3. A simple implementation

We briefly describe in this section the numerical implementation of our model for the matter density bispectrum. Further details can be found in Valageas & Nishimichi (2011).

3.1. Halo properties

For numerical computations we use the same halo model as the one used in Valageas & Nishimichi (2011) for the power spectrum. Thus, the halo mass function is given by (Valageas 2009) $\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{,}$(37)which is normalized to unity and provides a good match to numerical simulations. We choose the usual NFW profile for the halo density profile ρ_M(x) (Navarro et al. 1997), which also has an explicit form for its Fourier transform ${\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{u}\end{array}}_{\mathit{M}}\mathrm{\left(}\mathit{k}\mathrm{\right)}$ (Scoccimarro et al. 2001). In particular, we truncate the halo profiles at the density contrast δ(<r₂₀₀) = 200, which defines their radius r₂₀₀. For the concentration parameter we take $\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} {\mathit{M}}_{\mathrm{\odot }}}\right)}^{-0.1}\hspace{0.17em}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{z}{\mathrm{)}}^{-0.8}\mathit{,}$(38)which is similar to the behaviors measured in numerical simulations (Dolag et al. 2004; Duffy et al. 2008), and was found in Valageas & Nishimichi (2011) to provide a good tool for the density power spectrum (in this sense, it would also describe to some degree the effect of halo substructures).

3.2. Perturbative contribution

For the perturbative contribution (34) to the bispectrum, we investigate both the standard perturbation theory and the “direct steepest-descent” resummation scheme introduced (among other perturbative expansions) in Valageas (2007a). We only go up to 1-loop order for both methods, so that the first approach only contains all 1-loop diagrams, while the second approach also includes partial resummations of diagrams at all orders.

A detailed description of these methods for two-point and three-point functions is given in Valageas (2008). To facilitate the theoretical comparison between both approaches it is convenient to write them with the same diagrammatic language. Then, the diagrams associated with standard perturbation theory are shown in Fig. 2, where the solid lines are the linear two-point functions, either the linear correlation (lines without arrows) or the linear response (lines with an arrow that shows the flow of time). The three-leg vertex is the kernel K_s that appears in the equation of motion, which can be written as , where  is a linear operator and ψ is a two-component vector that describes both the density and velocity fields. The diagrams shown in Fig. 2 are not those associated with the usual description of standard perturbation theory, which is described in Appendix A, but of course they are equivalent. The expansion of Fig. 2 applies to the three-point correlation C₃ =  ⟨ ψψψ ⟩ , which contains the density and velocity three-point functions, as well as their cross correlations, but in this article we only consider the density three-point correlation, that is, the matter density bispectrum in Fourier space.

The usual description of standard perturbation theory arises from the expansion of the density and velocity fields over powers of the linear field δ_L, as in Eqs. (A.1), (A.2). Then, the N-point correlations are obtained by taking the Gaussian average of the product of these expansions, as in Eq. (A.3), using Wick’s theorem as in Eq. (A.4). This yields diagrams with vertices F_n that have an increasing number n of legs as one goes to higher orders, see also Fig. A.1. The diagrams of Fig. 2 are obtained from a path-integral formalism, where the averages over the initial conditions are already taken and one directly works with correlation and response functions. Then, expanding over powers of the nonlinear interaction vertex K_s gives the expansion of Fig. 2. It looks quite different from the standard diagrams, since only one vertex appears, i.e., the three-leg kernel K_s, whatever the order of the expansion. On the other hand, in addition to the linear power spectrum these diagrams also involve the linear response function. Then, the order of the expansion corresponds to the number of loops of the diagrams.

Since both expansions can also be written as expansions over powers of the linear power spectrum P_L(k), they are actually equivalent, at each order. However, at a given order, there can be a different number of diagrams between both expansions, and it is only the two sums over all diagrams of that order over P_L(k) that coincide. In particular, it can be seen that diagrams (h) and (i) in Fig. 2 show an ultraviolet divergence for linear power spectra with a large-k tail P_L(k) ∝ kⁿ with n ≥ −3. However, these two divergences cancel out and the sum is finite for n < −1, as in the case of the standard diagrams (see Valageas 2008, for details).

The lowest-order contribution, or “tree-order” term, is given by diagram (a), which yields the well-known result $\begin{array}{ccc}{\mathit{B}}^{\mathrm{\left(}\mathit{a}\mathrm{\right)}}\mathrm{\left(}{\mathit{k}}_{\mathrm{1}}\mathit{,}{\mathit{k}}_{\mathrm{2}}\mathit{,}{\mathit{k}}_{\mathrm{3}}\mathrm{\right)}& \mathrm{=}& {\mathit{P}}_{\mathrm{L}}\mathrm{\left(}{\mathit{k}}_{\mathrm{2}}\mathrm{\right)}{\mathit{P}}_{\mathrm{L}}\mathrm{\left(}{\mathit{k}}_{\mathrm{3}}\mathrm{\right)}[\frac{\mathrm{10}}{\mathrm{7}}\mathrm{+}\left(\frac{{\mathit{k}}_{\mathrm{2}}}{{\mathit{k}}_{\mathrm{3}}}\mathrm{+}\frac{{\mathit{k}}_{\mathrm{3}}}{{\mathit{k}}_{\mathrm{2}}}\right)\frac{{{k}}_{\mathrm{2}}\mathrm{·}{{k}}_{\mathrm{3}}}{{\mathit{k}}_{\mathrm{2}}{\mathit{k}}_{\mathrm{3}}}\\ & & \mathrm{+}\hspace{0.17em}\frac{\mathrm{4}\mathrm{(}{{k}}_{\mathrm{2}}\mathrm{·}{{k}}_{\mathrm{3}}{\mathrm{)}}^{\mathrm{2}}}{\mathrm{7}{\mathit{k}}_{\mathrm{2}}^{\mathrm{2}}{\mathit{k}}_{\mathrm{3}}^{\mathrm{2}}}]\mathrm{+}\mathrm{2}\hspace{0.17em}\mathrm{cyc}\mathit{.}\end{array}$(39)The explicit expressions of the 1-loop order diagrams (b), ..., (i), which are on the order of ${\mathit{P}}_{\mathrm{L}}^{\mathrm{3}}$ (as shown by the three blue solid lines which they contain), are given in Appendix A of Valageas (2008). For numerical computations it is convenient to perform analytically integrations over angles instead of directly implementing these expressions into the codes. Although this requires some care (since angular integrations may yield trigonometric or hyperbolic functions, depending on the amplitude of the wavenumbers), there is no fundamental difficulty.

Of course, it is also possible (and more common) to use the standard diagrams for the computation of the bispectrum within standard perturbation theory (Scoccimarro 1997; Scoccimarro et al. 1998). As noticed above, for our purposes the interest of the description of Fig. 2 is that it clarifies the link with the “direct steepest-descent” resummation scheme that we also investigate in this article. We focus on this specific resummation scheme to be consistent with our previous study for the power spectrum in Valageas & Nishimichi (2011). In addition, as discussed in that article (in particular in its Appendix A), this method allows a fast numerical implementation. Moreover, its predictions for the bispectrum (and higher order correlations), and more precisely the structure of its diagrammatic expansion, have already been studied in Valageas (2008), so that no further theoretical work is required. Then, within this “direct steepest-descent” resummation scheme, the diagrams for the bispectrum up to 1-loop order are those shown in Fig. 3. The difference with Fig. 2 is that we now have double-line two-point functions. They correspond to the nonlinear two-point correlation and response functions, as given by the same method up to the same 1-loop order. The single lines are the linear two-point functions, as in Fig. 2.

Within this resummation scheme, nonlinear two-point functions at “1-loop” order actually contain an infinite number of “bubble” diagrams, as shown in Figs. 8 and 9 of Valageas (2008). There, the label “1-loop” does not mean that we only include 1-loop diagrams, as in the standard perturbation theory of Fig. 2, but that the diagrammatic expansion is only complete up to 1-loop (i.e., although we include some diagrams at all orders, we miss some contributions at 2-loop and higher orders). Then, substituting the expression of the nonlinear two-point functions in terms of the linear functions, we can check that the diagram (a) of Fig. 3 actually contains the five diagrams (a), ..., (e), of Fig. 2, as well as an infinite number of higher order diagrams. On the other hand, the diagrams (f), ..., (i), of Fig. 3 contain their counterparts of Fig. 2: each diagram among (f), ..., (i), of Fig. 2 is the lowest order contribution to the corresponding diagram in Fig. 3, where nonlinear two-point functions are replaced by their linear counterparts.

At the order ${\mathit{P}}_{\mathrm{L}}^{\mathrm{3}}$, which corresponds to the 1-loop diagrams of Fig. 2, we can also consider the “mixed” case shown in Fig. 4, where we keep the resummed diagram (a) of Fig. 3, but we use for diagrams (f), ..., (i), their lowest order terms shown in Fig. 2. Of course, these three choices, drawn in Figs. 2−4, agree up to order ${\mathit{P}}_{\mathrm{L}}^{\mathrm{3}}$, and only differ by the number of higher order terms that are included.

Thus, for the perturbative (3-halo) contribution (34) we will consider the three alternatives shown in Figs. 2−4. Contrary to the numerical computations used in Valageas (2008), we do not introduce any approximation for the computation of these diagrams. Therefore, the perturbative term B_pert(k₁,k₂,k₃) contains no free parameter and is exact, up to 1-loop order (or to the truncation order of the perturbative scheme). This is an improvement over most previous studies involving the halo model (Ma & Fry 2000b; Fosalba et al. 2005), which only used the tree-order bispectrum (39) for the 3-halo term, with the addition of bias factors (which we do not introduce in our approach).

4. Comparison with numerical simulations

4.1. Numerical procedure

We use a set of large N-body simulations in a ΛCDM universe, consistent with the five-year observation of the WMAP satellite (Komatsu et al. 2009), which are described in Valageas & Nishimichi (2011). We analyze the four main simulations (N = 2048³), together with supplementary smaller simulations (N = 1024³,512³ and 256³) for convergence tests.

We measure the bispectrum from a fast Fourier transformation of the density field obtained by the cloud-in-cells interpolation of the N-body particles, at z = 0.35, 1 and 3. In doing so, we use the folding scheme to speed up the measurements at large wavenumber bins without systematic error from the interpolation (see Valageas & Nishimichi 2011; and also Colombi et al. 2009). We set the wavenumber bins for k₁, k₂, and k₃, logarithmically with 2 bins per factor 2.

The statistical uncertainties (at 1-σ level) of the measurements are estimated assuming that they mainly come from their Gaussian part (Scoccimarro et al. 1998), $\begin{array}{ccc}{\left[\mathrm{\Delta }\mathit{B}\mathrm{\left(}{\mathit{k}}_{\mathrm{1}}\mathit{,}{\mathit{k}}_{\mathrm{2}}\mathit{,}{\mathit{k}}_{\mathrm{3}}\mathrm{\right)}\right]}^{\mathrm{2}}\mathrm{=}\frac{{\mathit{s}}_{\mathrm{123}}\mathit{V}}{\mathrm{(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathrm{3}}\mathrm{2}{\mathit{N}}_{\mathrm{tri}}}\mathit{P}\mathrm{\left(}{\mathit{k}}_{\mathrm{1}}\mathrm{\right)}\mathit{P}\mathrm{\left(}{\mathit{k}}_{\mathrm{2}}\mathrm{\right)}\mathit{P}\mathrm{\left(}{\mathit{k}}_{\mathrm{3}}\mathrm{\right)}\mathit{,}& & \end{array}$(40)where the factor s₁₂₃ is 6 for equilateral triangles, 2 for isosceles triangles and 1 otherwise. In the above, V denotes the volume of the simulations and N_tri is the number of Fourier space triangles that fall into the bin of (k₁,k₂,k₃).

We do not correct for the shot-noise contributions to the power spectrum and the bispectrum because the naive expectation for the shot noise obtained by assuming a Poisson process does not seem very accurate, especially on large scales. Because we start the simulation from a regular lattice, the discreteness noise is boxed in on scales smaller than the inter-particle distance at the beginning. After gravitational evolution, the density fluctuations neatly follow the linear growth rate at k ≲ 0.1 h Mpc^-1, hence we conclude that we do not have any sign of shot noise on those large scales. On small scales, however, we can see that our measurements from simulations approach the Poisson noise (e.g., Verde et al. 2002), $\begin{array}{ccc}& & {\mathit{P}}_{\mathrm{shot}}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathrm{(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathrm{3}}}\frac{{\mathit{L}}_{\mathrm{box}}^{\mathrm{3}}}{\mathit{N}}\mathit{,}\\ & & {\mathit{B}}_{\mathrm{shot}}\mathrm{\left(}{\mathit{k}}_{\mathrm{1}}\mathit{,}{\mathit{k}}_{\mathrm{2}}\mathit{,}{\mathit{k}}_{\mathrm{3}}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathrm{(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathrm{6}}}[\frac{{\mathit{L}}_{\mathrm{box}}^{\mathrm{3}}}{\mathit{N}}\left(\mathit{P}\mathrm{\left(}{\mathit{k}}_{\mathrm{1}}\mathrm{\right)}\mathrm{+}\mathit{P}\mathrm{\left(}{\mathit{k}}_{\mathrm{2}}\mathrm{\right)}\mathrm{+}\mathit{P}\mathrm{\left(}{\mathit{k}}_{\mathrm{3}}\mathrm{\right)}\right)\\ & & \mathrm{-}\hspace{0.17em}\mathrm{2}{\left(\frac{{\mathit{L}}_{\mathrm{box}}^{\mathrm{3}}}{\mathit{N}}\right)}^{\mathrm{2}}]\mathit{,}\end{array}$for the power spectrum and the bispectrum, respectively. Note that the shot noise is included in P(k₁), P(k₂), and P(k₃) of Eq. (42). Instead of subtracting these contributions from the measured spectra, we assess the shot-noise level from Eqs. (41) and (42), and then plot their relative importance in Figs. 16 and 17 below. For the reduced bispectrum, we estimate the shot noise level from $\begin{array}{ccc}{\mathit{Q}}_{\mathrm{eq}\mathit{,}\hspace{0.17em}\mathrm{shot}}\mathrm{\left(}\mathit{k}\mathrm{\right)}& \mathrm{=}& \end{array}$(43)which is plotted in Fig. 18 below, where again P(k) and B_eq(k) are the measured power spectrum and bispectrum (and thus include the shot-noise contribution and other numerical errors, as compared with the theoretical values).

We now compare our model, described in Sects. 2 and 3, with these numerical simulations.

4.2. Bispectrum B(k₁, k₂, k₃)

4.2.1. Equilateral triangles

We show in Fig. 5 our results for the bispectrum, from linear to highly nonlinear scales, for equilateral configurations. Here we take the perturbative 3-halo contribution equal to the standard 1-loop result, which also corresponds to Fig. 2. We can see that we obtain a reasonably good agreement with the numerical simulations over all these scales. This is remarkable since our model contains no free parameters. Indeed, the 2-halo and 1-halo terms are fully determined by the halo mass function and density profiles that were used in Valageas & Nishimichi (2011) for the power spectrum. This provides a significant improvement over the phenomenological model presented in Pan et al. (2007), based on a generalization of the scale transformation introduced for the two-point correlation function by Hamilton et al. (1991), which breaks down on highly nonlinear scales.

However, we can see that at high redshift (right panel at z = 3) we underestimate the bispectrum in the transition range between linear and nonlinear scales. The same behavior appears for the power spectrum, see Fig. 9 in Valageas & Nishimichi (2011). This is likely to be due in part to higher order perturbative terms, which play a greater role at z = 3 than at z = 0.35, in agreement with the detailed study performed in Valageas (2011) that compares perturbative and nonperturbative contributions. On the other hand, it is natural to expect the transition range to be the most difficult to reproduce by models of the kind studied in this paper. Indeed, this domain is already beyond perturbation theory but does not yet correspond to the inner relaxed cores of virialized halos. Therefore, it is at the limit of validity of the two ingredients (perturbation theory and halo model) used in our approach. We will discuss in more detail this transition range and possible improvements on these scales in Sects. 6.2 and 7.

We clearly see in Fig. 5 the decay on large scales of the 1-halo and 2-halo contributions, in agreement with Eqs. (26) and (33). As explained in Sect. 2, this is due to the new counterterms $\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{W}\end{array}\mathrm{\left(}{\mathit{k}}_{\mathit{j}}{\mathit{q}}_{\mathit{M}}\mathrm{\right)}$ of Eqs. (24) and (31), which ensure a physically meaningful behavior. On the other hand, in agreement with Sefusatti et al. (2010), we can see that taking into account the 1-loop perturbative contribution significantly extends the domain of validity of the 3-halo perturbative term, as compared with the tree-level contribution.

4.2.2. Isosceles triangles

We show in Fig. 6 our results for isosceles triangles as a function of scale k and for the fixed shapes defined by k₁ = k₂ = k and k₃ = k/α, with constant ratios α. We consider the cases α = 2 and α = 4. As expected, as α increases (i.e., the triangles get “squeezed”), the 2-halo contribution becomes more important compared with the 1-halo contribution, because the scale 1/k₃ grows and it is more likely to have two distinct halos separated by such a large distance rather than a single halo of this size. However, for these moderate values of α the 2-halo contribution always remains subdominant because the perturbative (i.e., “3-halo”) contribution is still dominant when the 1-halo term becomes larger than the 2-halo term. For more “squeezed” shapes, an intermediate regime would appear, dominated by the 2-halo term.

We can check that our model agrees well with the numerical simulations and is able to describe both large and small scales. Again, this extends the analytical predictions beyond the scales described by the phenomenological model of Pan et al. (2007). However, at z = 3 (right panels) we can see the underestimation on transition scales already noticed in Fig. 5.

Next, we show in Fig. 7 the evolution of the bispectrum with the shape of the triangle formed by  {k₁,k₂,k₃} , for isosceles configurations k₁ = k₂.

First, we consider in the upper row the dependence on k₁ = k₂ = k at fixed k₃, so that the two equal-length sides run from k₃/2 (flat triangle that reduces to a line) to infinity (squeezed limit). In agreement with the behavior observed in Figs. 5 and 6, the bispectrum decreases for higher values of k₁ = k₂. The value of the fixed side k₃ grows from the left to the right panel, so that we go from linear to highly nonlinear scales (in each case we consider the typical range k₃/2 < k₁ = k₂ < 10k₃). In the left panel, we are dominated by the 3-halo perturbative contribution, and since k₃ remains small, for high values of k₁ = k₂ the 2-halo term is larger than the 1-halo term. In the middle panel, we are already dominated by the 1-halo term (except at low k at z = 3 where the 3-halo perturbative contribution is larger), and the the 2-halo term is subdominant. In the right panel the bispectrum is completely dominated by the 1-halo contribution.

Second, we consider in the lower row the dependence on k₃ at fixed value of the common sides k₁ = k₂ = k, so that the third side runs from 0 (squeezed limit) to 2k₁ (flat triangle that reduces to a line). Again, the bispectrum decreases as the dependent wavenumber, here k₃, grows. However, the dependence is much weaker than the one found in the upper row. In the left panel, we are dominated by the linear-order perturbative contribution, so that all curves agree with each other. The too high values of the numerical points at low k₃ are a numerical artifact caused by the finite boxsize and the numerical binning used to compute the bispectrum, because linear theory must be increasingly accurate on larger scales. In the middle panel, we can see the crossover between the 3-halo and 1-halo regimes. The 2-halo term is never dominant on these scales. In the right panel we are mostly dominated by the 1-halo term.

In all cases we can see that our model agrees well with the numerical simulations. However, at z = 3, in the upper middle panel and in the lower right panel we can see a trace of the underestimation on the transition scales found in the right panels of Figs. 5 and 6.

4.3. Reduced bispectrum

As seen in the previous figures, the bispectrum varies by many orders of magnitude on the scales of interest. Therefore, it is customary to introduce the “reduced bispectrum” defined as $\mathit{Q}\mathrm{\left(}{\mathit{k}}_{\mathrm{1}}\mathit{,}{\mathit{k}}_{\mathrm{2}}\mathit{,}{\mathit{k}}_{\mathrm{3}}\mathrm{\right)}\mathrm{=}\frac{\mathit{B}\mathrm{\left(}{\mathit{k}}_{\mathrm{1}}\mathit{,}{\mathit{k}}_{\mathrm{2}}\mathit{,}{\mathit{k}}_{\mathrm{3}}\mathrm{\right)}}{\mathit{P}\mathrm{\left(}{\mathit{k}}_{\mathrm{1}}\mathrm{\right)}\mathit{P}\mathrm{\left(}{\mathit{k}}_{\mathrm{2}}\mathrm{\right)}\mathrm{+}\mathrm{2}\mathrm{cyc}\mathit{.}}\mathrm{·}$(44)Indeed, this ratio goes to a constant on large scales, as seen from the tree-order expression (39), while it shows a moderate growth on small scales (Bernardeau et al. 2002a).

We compare in Fig. 8 our results for the reduced equilateral bispectrum, Q_eq(k) = Q(k,k,k), with numerical simulations. Here, in the denominator of Eq. (44) we use for the theoretical predictions the power spectrum obtained with our model described in Valageas & Nishimichi (2011) (that is, using the same halo model as in this paper and the direct steepest-descent resummation for the perturbative 2-halo term), while for the numerical curves we use the power spectrum measured in the simulations.

The vertical arrow in the upper right part of each panel of Fig. 8 shows the wavenumber where the shot-noise level of the numerical simulations becomes larger than 10%. We can see that this also roughly corresponds to the scale where the measured reduced bispectrum is maximum. Therefore, the downturn and the decrease of the reduced bispectrum at higher k found in the simulations is only a numerical artifact caused by the limited resolution, and these points must be discarded. In particular, at lower redshift where the simulations are able to probe deeper into the nonlinear regime, we clearly see in the left panel at z = 0.35 the fast growth of Q_eq(k) in the highly nonlinear regime, after a small plateau on the transition scales. As shown by the two left panels, this high-k fast increase of Q_eq(k) is well reproduced by our model.

We can see that we obtain a reasonably good agreement with the simulations on both large and small scales. This is consistent with our previous results, shown in Fig. 5 and in Valageas & Nishimichi (2011), where we found that both the bispectrum and the power spectrum are well reproduced on quasilinear and highly nonlinear scales. As expected, we recover the large-scale plateau associated with the tree-level limit (39), and the early rise owing to higher order perturbative contributions. At very high k, in agreement with previous studies (Ma & Fry 2000b; Cooray & Sheth 2002; Fosalba et al. 2005), the halo model leads to a continuous growth of the reduced bispectrum, which is a marked difference with the prediction of the stable clustering ansatz Peebles (1980). This agrees with our numerical simulations up to the scales that are well resolved. As shown in Sect. 2 of Valageas (1999), the positivity of the matter density actually implies that the reduced bispectrum, and more generally the real-space coefficients ${\mathit{S}}_{\mathit{p}}\mathrm{=}\mathrm{⟨}{\mathit{\rho }}_{\mathit{R}}^{\mathit{p}}{\mathrm{⟩}}_{\mathrm{c}}\mathit{/}\mathrm{⟨}{\mathit{\rho }}_{\mathit{R}}^{\mathrm{2}}{\mathrm{⟩}}_{\mathrm{c}}^{\mathit{p}\mathrm{-}\mathrm{1}}$, reach a constant or keep growing on small scales, as the density variance $\mathrm{⟨}{\mathit{\rho }}_{\mathit{R}}^{\mathrm{2}}{\mathrm{⟩}}_{\mathrm{c}}$ increases. The limiting behavior of the constant reduced bispectrum and the constant coefficients S_p is achieved for the stable clustering ansatz, and more generally for multifractal models such that $\mathrm{⟨}{\mathit{\rho }}_{\mathit{R}}^{\mathit{p}}{\mathrm{⟩}}_{\mathrm{c}}$ are governed by a single fractal exponent α for the values of p that are considered. The density field described by a halo model clearly violates these conditions¹ because it does not display this scale invariance, with a characteristic nonlinear mass associated with the falloff of the halo mass function and reasonably smooth profiles that depend on the mass scale (through their concentration parameter). This implies that Q_eq(k) has to grow on small scales, as checked in Fig. 8.

On the transition scales, we underestimate both the bispectrum and the power spectrum, and it appears that the latter effect is dominant, which leads to the “bump” seen in Fig. 8. As seen from the numerical results, this feature is not physical and only caused by the shortcomings of the theoretical model on this range. Such a feature was also noticed in previous studies (Scoccimarro et al. 2001; Takada & Jain 2003; Fosalba et al. 2005), which found that by using a larger halo radius, taking into account exclusion constraints in the 2-halo and 3-halo terms, or introducing a high-mass cutoff in the halo mass function, one might remove this “bump” and reach a better agreement with numerical simulations. On the other hand, Guo & Jing (2009) find that no self-consistent halo model is able to reproduce well both the matter power spectrum and bispectrum on these transition scales. We will return to this problem and these modifications in Sect. 6.2 and will discuss in Sect. 7 an alternative procedure to improve the model on these scales, using this artificial “bump” as a signature of insufficient accuracy and introducing simple interpolations that resolve most of this problem, while keeping large and small scales unchanged.

5. Comparison with previous halo models

We show in Fig. 9 the impact of two features of our model that differ from previous implementations of the halo model (Scoccimarro et al. 2001; Fosalba et al. 2005; Guo & Jing 2009), i) the use of 1-loop perturbation theory instead of the tree-level contribution (39) for the 3-halo term; and ii) the counterterms $\begin{array}{}{\mathrm{˜}}_{}\\ \mathit{W}\end{array}$ in Eqs. (24) and (31) in the 1-halo and 2-halo terms.

To distinguish more clearly between various models, that is, to reduce the size of the vertical axis, we plot in Fig. 9 the effective reduced bispectrum defined as in Eq. (44), but where we use the power spectrum P_S(k) from Smith et al. (2003) in the denominator, for the models as well as for the numerical simulations. Thus, by dividing the bispectra by the same denominator we avoid introducing additional errors through the denominator and can make a clear comparison between the various bispectra and the simulations. In other words, Fig. 9 is not meant to show the actual reduced bispectrum (44), which was displayed in Fig. 8 for equilateral triangles, but only to show the bispectrum B on a more convenient scale. We focus on quasilinear scales where the different models can be distinguished, since on very large or small scales they converge towards the tree-level contribution (39) or the 1-halo contribution with a negligible counterterm.

First, we can see in Fig. 9 that discarding the contribution of 1-loop diagrams (the diagrams (b) to (i) in the representation of Fig. 2) leads to a significant underestimate of the bispectrum on weakly nonlinear scales, k ~ 0.2 h Mpc^-1 at z ≤ 3, in agreement with previous works based on perturbation theory alone (Scoccimarro et al. 1998; Sefusatti et al. 2010). In particular, at z = 0.35 the 1-loop contribution is necessary to obtain a good match to the simulations and appears to be sufficient to make the bridge to the smaller scales where the 1-halo term is dominant. At z = 3, the 1-loop contribution also extends up to k ~ 0.4 h Mpc^-1 the good agreement with simulations, while using the tree-level contribution alone yields significant discrepancies as soon as k > 0.2 h Mpc^-1. However, it is no longer sufficient to make the bridge with the highly nonlinear scales dominated by the 1-halo term, because we now underestimate the bispectrum on the transition scales, k ~ 1 h Mpc^-1. This behavior was already noticed in Figs. 5 and 6.

Second, setting the counterterms to zero in the 2-halo and 1-halo contributions (while keeping the 1-loop contribution in the 3-halo term) obviously increases the bispectrum and worsens the agreement with simulations. This is most clearly seen on the quasilinear scales, which are already well described by standard perturbation theory so that the extra power associated with the unphysical constant asymptotes at low k of the uncorrected 2-halo and 1-halo contributions spoil the good agreement with simulations.

Making simultaneously both approximations, discarding 1-loop diagrams and halo counterterms, as in usual implementations of the halo model, also yields significantly worse results. Although both effects partly compensate, neglecting 1-loop contributions dominates so that the resulting bispectrum is significantly too low on weakly nonlinear scales. Therefore, in addition to a more satisfactory physical behavior (systematic agreement on large scales, up to higher order, and decay of 2-halo and 1-halo contributions), our approach gives a better accuracy on weakly nonlinear scales. This shows the importance of including both 1-loop contributions and the 2-halo and 1-halo counterterms.

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

6.1. Dependence on the perturbative scheme

We investigate in this section the dependence of our results on the perturbative scheme used for the 3-halo contribution. As in Fig. 9, the bispectra obtained from the analytical models and the numerical simulations are divided by the same power spectra, from Smith et al. (2003), to make a clear comparison.

We compare in Fig. 10 the results obtained using for the perturbative 3-halo contribution either the standard 1-loop result of Fig. 2, the full steepest-descent resummation of Fig. 3, or the mixed case of Fig. 4. For each choice we show both the 3-halo term itself (blue curves that decrease at high k) and the full prediction, where we sum the 3-halo, 2-halo, and 1-halo contributions. (The 2-halo and 1-halo contributions are the same for all three choices.)

Let us first consider the upper row, which corresponds to equilateral triangles. (At z = 3 we can see the underestimate of the bispectrum in the transition range that we had already noticed in Fig. 5.) The standard and mixed cases give very close results. At z = 3, the mixed case, giving a slightly greater 3-halo contribution, yields a slightly better agreement with numerical simulations around 1 h Mpc^-1. However, this is not conclusive, especially since the 1-halo term is already non-negligible on these scales and it is not sufficient to significantly improve the match to simulations. The complete steepest-descent resummation of Fig. 3 yields a 3-halo contribution that decays significantly faster at high k. At z = 0.35 and z = 1, two numerical points in the range 0.1 < k < 0.2 h Mpc^-1 suggest that this may be a true improvement on these quasilinear scales. This improvement might be important for the analysis of baryon acoustic oscillations; it has been shown that resummation schemes give better predictions for the power spectrum on this scale (e.g., Crocce & Scoccimarro 2008; Nishimichi et al. 2009; Valageas & Nishimichi 2011), and we could expect the same thing for the bispectrum. However, beyond 0.2 h Mpc^-1 the faster decay significantly worsens the agreement with the simulations.

The upper row of Fig. 10 does not allow us to clearly discriminate between the various schemes because they mainly differ in the transition regime, where the 1-halo contribution is already important and all models tend to underestimate the bispectrum. However, looking now at the lower row, associated with isosceles triangles, and especially at the two lower right panels that show the evolution of the bispectrum with the shape of the triangle, we can see that the standard 1-loop perturbation theory yields a significantly better agreement with simulations than the two other perturbative expansions. This is also the simplest scheme for practical purposes.

Thus, from the analysis of both equilateral and isosceles triangles, we can conclude that using the standard 1-loop perturbation theory for the perturbative 3-halo term is the most efficient and accurate scheme among the three approaches studied in this paper. This is quite different from the power spectrum, studied in Valageas & Nishimichi (2011), where the 1-loop steepest-descent resummation is clearly better than standard perturbation theory. There, it is more accurate on weakly nonlinear scales but it is also the only choice (with other resummation schemes) that allows a combined model for both large and small scales. Indeed, the standard 1-loop perturbation theory yields a 2-halo contribution that keeps growing on small scales (for the logarithmic power of Eq. (45) below) and that worsens the agreement with simulations (even though it is smaller than the 1-halo term).

We do not have this small-scale problem for the bispectrum, as seen in Figs. 5, 6, and 8. Indeed, the 1-loop standard perturbation theory gives a contribution to the bispectrum that quickly decreases at high k and does not spoil the good agreement with simulations shown by the 1-halo term in this regime. This is best seen in Fig. 8, which shows that the 1-loop standard perturbation theory contribution to the reduced equilateral bispectrum Q_eq(k) itself also decreases at high k, in spite of the denominator 3P(k)². Therefore, contrary to the power spectrum, it is now possible to use the 1-loop standard perturbation theory to build a combined model that covers all scales. Moreover, as shown by Fig. 10, it happens that this is also more accurate than the two other perturbative schemes investigated in this paper, which we presented in Figs. 3 and 4.

It is possible that other perturbative schemes that we did not study here provide more accurate results. On the other hand, higher orders of standard perturbation theory may yield contributions that are increasingly large on small scales, so that beyond a certain order they lead to a 3-halo term that is unphysically large and spoils the good agreement with simulations on small scales. Then, one would need to use other perturbative approaches, such as the resummation schemes investigated here, or to add some ad hoc cutoff on small scales.

For completeness, we mention a few different approaches. A first method, introduced in Crocce & Scoccimarro (2006b,a) from a diagrammatic approach, reorganizes the standard perturbation theory and performs partial resummations by introducing both correlation functions and propagators (or response functions), in a fashion somewhat similar to the approach described in Sect. 3.2. However, a key step in this method is the matching between the low-order low-k behavior and the high-k asymptote of the propagator. This ensures a good behavior of this quantity in the nonlinear regime, which has been checked against numerical simulations (Crocce & Scoccimarro 2006a; Bernardeau et al. 2008), and expresses a well-understood “sweeping effect” associated with the random transport of density structures by large-scale velocity flows (Valageas 2007b; Bernardeau & Valageas 2008, 2010). Another advantage of this approach is that its extensions to high-order quantities (Bernardeau et al. 2008), to non-Gaussian initial conditions (Bernardeau et al. 2010), and to higher perturbative orders (Anselmi et al. 2011), have already been studied.

A second approach, developed in Valageas (2007a,b) as a “large-N 2PI” expansion, and in Taruya & Hiramatsu (2008); Taruya et al. (2009) as a “closure approximation”, is quite similar to the one described in Sect. 3.2. Although in principle it may give more accurate results, its full implementation is more complex because self-energy terms depend on the nonlinear quantities that are looked for. This leads to coupled nonlinear equations, which depend on scales and times, and to heavier numerical computations (Valageas 2007a). Because for practical purposes one requires fast algorithms this presents a significant disadvantage, unless one introduces further simplifications (Taruya et al. 2009).

A third method, introduced in Pietroni (2008), directly obtains the hierarchy of equations obeyed by the many-body correlation functions from the equations of motion, in a fashion similar to the standard BBGKY hierarchy. Then, truncation at a given order (at the trispectrum in Pietroni (2008)) defines an approximation for lower-order correlations. This method has already been used to study the effect of massive neutrinos (Lesgourgues et al. 2009) and of primordial non-Gaussianities (Bartolo et al. 2010b). An advantage of this approach is that it does not involve propagators, but only the usual many-body correlations also encountered in standard perturbation theory, and it is written in terms of single-time quantities. This is a great simplification that should allow efficient numerical computations.

A fourth approach, developed in Matsubara (2008b,a), is based on a Lagrangian framework. This would be well-suited to the approach described in this paper and presents the advantage that it provides direct extensions to redshift-space statistics. Unfortunately, as noticed in Valageas & Nishimichi (2011), in its present form this Lagrangian perturbation theory does not fare as well as its Eulerian counterparts when tested against numerical simulations. However, this approach remains interesting, especially in view of applications to redshift space.

In any case, it would not be too difficult to combine any of these methods with halo models by using the framework described in this paper, which is more general than the use of the standard perturbation theory or the two resummation schemes described in Sect. 3.2.

6.2. Dependence on halo properties

We now investigate the dependence of our results on the details of the halo model. As we recalled in Sect. 4.3, the “bump” shown by the reduced bispectrum in Fig. 8 was also found in previous works, which noticed that this feature may be cured to some extent by tuning halo parameters (Scoccimarro et al. 2001; Takada & Jain 2003; Fosalba et al. 2005). Therefore, we show in Fig. 11 the reduced bispectrum (upper panel, as in Fig. 8) and the scaled bispectrum (lower panel, as in Fig. 10) for several modifications of halo parameters at redshift z = 0.35.

Following Scoccimarro et al. (2001); Wang et al. (2004); Fosalba et al. (2005), we investigate the impact of truncating the halo mass function below M < 10¹⁵ or M < 5  ×  10¹⁴ h^-1 M_⊙. In agreement with these studies, removing large halos decreases both the matter power spectrum and bispectrum on mildly nonlinear scales (since smaller scales are associated with small halos), and the net effect on the reduced bispectrum is also a small decrease of the artificial “bump”. However, the upper panel shows that even truncating at M < 5  ×  10¹⁴ h^-1 M_⊙ is not sufficient to erase the “bump” and to provide a good agreement with simulations for the reduced bispectrum. Two differences with previous studies are that the 3-halo term includes 1-loop perturbative contributions, which are still important on these scales as shown in Figs. 8 and 10, and that the 1-halo and 2-halo contributions contain the counterterms $\begin{array}{}{\mathrm{˜}}_{}\\ \mathit{W}\end{array}$, so that the relative importance of the 2-halo and 1-halo contributions is somewhat smaller than in previous models on the weakly nonlinear scales associated with the early rise of the “bump”. On the other hand, because both the power spectrum and the bispectrum simultaneously decrease when we add such a high-mass cutoff and these effects partly compensate in the reduced bispectrum, before we obtain a significant improvement for the latter quantity both the power spectrum and the bispectrum have already been decreased by a large amount that disagrees with numerical results, as shown in the lower panel for the bispectrum. Moreover, in our large simulations we find halos up to 3  ×  10¹⁵ h^-1 M_⊙, so that the cutoff at M < 5  ×  10¹⁴ h^-1 M_⊙ is already too low to be fully justified by the finite box size of the simulations.

Then, following Fosalba et al. (2005), we investigate the impact of larger halo radii. More precisely, we consider the impact of defining halos by a nonlinear density threshold δ = 50 instead of δ = 200, following the procedure described in Sect. 6.1 of Valageas & Nishimichi (2011). This also involves a slightly different halo mass – concentration parameter relation, so as to keep a satisfactory match to numerical results for the power spectrum. As seen in Fig. 11, this modification does not bring a significant improvement either.

We find similar results at higher redshifts, as seen in Fig. 12 at z = 3 (where we consider the smaller mass thresholds M < 3  ×  10¹³ and M < 10¹³ h^-1 M_⊙).

Therefore, we reach the same conclusion as Guo & Jing (2009), that tuning these halo parameters cannot simultaneously provide accurate results for the power spectrum and bispectrum, especially compared with large simulations where massive halos are present, which removes the justification for introducing severe high-mass cutoffs.

The discrepancies found on these transition scales are not so surprising because one expects the transition range to be the most difficult to describe in a systematic fashion, since both the perturbative expansion (that neglects shell crossing) and the halo model (which assumes spherical and relaxed objects in our implementation) may break down in this intermediate regime. Another possibility is that we are missing important high-order perturbative terms that have not been included in the perturbative expansions used here, which are only complete up to one-loop order. Indeed, for the power spectrum perturbative terms up to order  ~9 (at z = 0) or  ~66 (at z = 3) are likely to be relevant (Valageas 2011). On the halo-model side, because in reality the density field on these transition scales shows a crossover from relaxed inner halo shells to outer infalling shells and filaments that cannot any longer be described in terms of well-defined halos, one cannot expect a systematic convergence to the numerical results by adding such simple modifications to halo properties. Then, one is probably limited to some extent by the intrinsic limitations of the halo model itself.

7. Improving the predictions for P(k) and B_eq(k) using the reduced bispectrum Q_eq(k)

The artificial “bump” found in Fig. 8 for the reduced equilateral bispectrum Q_eq(k), which is mostly caused by the underestimation of the power P(k), suggests a simple trick to improve the predictions for P(k) and B_eq(k). The idea is to use the unphysical “bump” shown by the predicted reduced bispectrum Q_eq(k) to automatically detect the range  [k₋,k₊]  where the model is not sufficiently accurate. The procedure that we investigate in this paper is to draw in the (log k,Q_eq) plane the lower tangent line to the predicted curve that was plotted in Fig. 8. We show this construction in Fig. 13, where we again plot the prediction of our model, described in the previous sections and labeled “direct” in this figure, as well as the lower tangent on the region where the “bump” appears. The two contact points between the model curve and its tangent line define the two wavenumbers, k₋ and k₊, between which the artificial “bump” arises and the model needs to be improved.

First, we modify the density power spectrum as shown in Fig. 14, where we plot 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{.}$(45)The modified power “tang.” is obtained from the “direct” prediction described in Valageas & Nishimichi (2011) by drawing in the (log k,log Δ²) plane the upper tangent line on the interval  [k₋,k₊]  that runs through the left point (log k₋,log Δ²(k₋)). For the cases shown in Fig. 14 the right contact point has an abscissa ${\mathit{k}}_{\mathrm{+}}^{\mathrm{\prime }}\mathit{<}{\mathit{k}}_{\mathrm{+}}$ (by construction we have ${\mathit{k}}_{\mathrm{-}}\mathit{<}{\mathit{k}}_{\mathrm{+}}^{\mathrm{\prime }}\mathrm{\le }{\mathit{k}}_{\mathrm{+}}$), and the power spectrum is only modified on the interval $\mathrm{\left[}{\mathit{k}}_{\mathrm{-}}\mathit{,}{\mathit{k}}_{\mathrm{+}}^{\mathrm{\prime }}\mathrm{\right]}$. As seen in Fig. 14, in this fashion we correct most of the artificial “dip” shown by our model without modifying the large-scale and small-scale regimes where the “direct” predictions are satisfactory².

Second, for equilateral triangles we modify the bispectrum in the same manner. Thus, as shown in Fig. 15, we again draw in the (log k,log B_eq) plane the upper tangent line, on the interval  [k₋,k₊] , that runs through the left point (log k₋,log B_eq(k₋)). This yields another right contact point at ${\mathit{k}}_{\mathrm{+}}^{\mathrm{\prime \prime }}$, which in the cases shown in Fig. 15 is located to the left of k₊ as there is again a knee in the shape of the bispectrum. This ensures that we follow the “direct” prediction beyond the knee, where it shows a good match to the numerical simulations, while correcting most of the artificial “dip” shown by the model.

As explained above, the interest of this procedure, based on the reduced equilateral bispectrum Q_eq(k), is to automatically define the wavenumbers k₋ and k₊. One could imagine defining these boundaries in a simpler manner from the power spectrum itself, for instance by ${\mathrm{\Delta }}^{\mathrm{2}}\mathrm{\left(}{\mathit{k}}_{\mathrm{\pm }}\mathrm{\right)}\mathrm{=}{\mathrm{\Delta }}_{\mathrm{\pm }}^{\mathrm{2}}$, with fixed values ${\mathrm{\Delta }}_{\mathrm{\pm }}^{\mathrm{2}}$ that would mark the transition regime (such as ${\mathrm{\Delta }}_{\mathrm{-}}^{\mathrm{2}}\mathrm{=}\mathrm{1}$ and ${\mathrm{\Delta }}_{\mathrm{+}}^{\mathrm{2}}\mathrm{=}\mathrm{30}$). However, owing to the curvature of the linear CDM power spectrum, the interval  [k₋,k₊]  where the model departs from the numerical simulations is not given by these redshift-independent conditions, as can be seen in Fig. 14. In particular, at higher redshift (i.e., lower value of the effective slope n of the linear power spectrum on the relevant scales) the threshold ${\mathrm{\Delta }}_{\mathrm{-}}^{\mathrm{2}}$ is seen to decrease. This is partly captured by our procedure, described in Fig. 13, which is sensitive to the shape of the initial conditions, through the shape of the predicted reduced bispectrum Q_eq(k).

Finally, going back to Fig. 13, from the improved model “tang.” constructed for the power spectrum and the bispectrum, shown in Figs. 14 and 15, we can build a new prediction “tang.” for the reduced bispectrum from Eq. (44). For equilateral triangles this reads as ${\mathit{Q}}_{\mathrm{eq}}^{\mathrm{tang}\mathit{.}}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathrm{=}{\mathit{B}}_{\mathrm{eq}}^{\mathrm{tang}\mathit{.}}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathit{/}\mathrm{\left[}\mathrm{3}{\mathit{P}}^{\mathrm{tang}\mathit{.}}\mathrm{\right(}\mathit{k}{\mathrm{)}}^{\mathrm{2}}\mathrm{]}$. The result, plotted in Fig. 13, shows a significant improvement over the “direct” prediction, which was already plotted in Fig. 8. However, we can see that although it is much smaller, the artificial “bump” has not been fully removed by the modifications to the power spectrum and the bispectrum. If one is interested in the reduced bispectrum Q_eq(k), one can introduce a last improvement by replacing this small “bump” by a flat plateau. This “flat” model is obtained from the “tang.” curve by running down the ${\mathit{Q}}_{\mathrm{eq}}^{\mathrm{tang}\mathit{.}}\mathrm{\left(}\mathit{k}\mathrm{\right)}$ curve over the interval  [k₋,k₊] , starting from the right boundary k₊, and imposing a monotonic decrease as k decreases. Thus, ${\mathit{Q}}_{\mathrm{eq}}^{\mathrm{flat}}$ and ${\mathit{Q}}_{\mathrm{eq}}^{\mathrm{tang}\mathit{.}}$ are identical over most of  [k₋,k₊]  (and on all larger and smaller scales), except under the remaining “bump” of ${\mathit{Q}}_{\mathrm{eq}}^{\mathrm{tang}\mathit{.}}$, where ${\mathit{Q}}_{\mathrm{eq}}^{\mathrm{flat}}$ is constant and equal to the local minimum of ${\mathit{Q}}_{\mathrm{eq}}^{\mathrm{tang}\mathit{.}}$ to the right of this “bump”.

The geometrical improvements described in Figs. 13−15 are not as elegant as one would wish for. Indeed, by combining perturbative expansions with halo models, as in this article and in Valageas & Nishimichi (2011), one would hope to build a model that shows a good match to numerical simulations on all scales, without further modifications. Unfortunately, as noticed above, at the current stage there remain some discrepancies on the transition scales. As discussed in Sect. 6.2, this is not surprising because transition scales may be at the limit of validity of both perturbation theory and halo models. Indeed, shell crossing (which is beyond the reach of the perturbative approaches studied here, based on the fluid approximation) is already important, and these scales do not correspond to relaxed halo inner shells, but rather to outer infalling clumps and to filaments (which cannot be described as isolated objects with a well-defined profile). In particular, we have seen in Fig. 11 that tuning halo parameters does not easily provide a significantly better agreement with simulations. Moreover, this would require some ad hoc modifications and new free parameters, which make the model less predictive. Then, the geometrical modifications introduced above can be seen as a simple procedure to obtain a good matching between the weakly and highly nonlinear regimes, and are no less satisfactory than more “algebraic” matchings. They offer the advantage to bypass all these complicating matters and the free parameters they would involve, but they share with approaches based on modifications to the halo parameters the lack of a systematic method toward increasingly high accuracy.

In addition to the automatic definition of the interval  [k₋,k₊] , associated with the transition range, an important property of this procedure is that the power spectrum and the equilateral bispectrum are not modified outside of this interval. Therefore, we keep the good match obtained on larger and smaller scales. In particular, large scales are still obtained by systematic perturbative expansions and small scales by phenomenological halo models, so that the final result is still a combination of these two approaches and keeps their distinct benefits.

We describe in Appendix B the impact on the real-space two-point correlation function of this modification to the power spectrum. In particular, we check that we keep an accuracy on the order of 1% on weakly nonlinear scales, x > 10 h^-1 Mpc, while reaching an accuracy on the order of 10% on nonlinear scales.

8. Typical accuracy of combined models

We now consider the accuracy of the combined model described in the previous sections and in Valageas & Nishimichi (2011). We first plot in Fig. 16 the relative difference between our model and the numerical simulations for the power spectrum, $\frac{\mathrm{\Delta }\mathit{P}}{\mathit{P}}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathrm{=}\frac{\mathrm{\left|}{\mathit{P}}_{\mathrm{model}}\mathrm{\right(}\mathit{k}\mathrm{)}\mathrm{-}{\mathit{P}}_{\mathit{N}\mathrm{-}\mathrm{body}}\mathrm{(}\mathit{k}\mathrm{\left)}\mathrm{\right|}}{{\mathit{P}}_{\mathit{N}\mathrm{-}\mathrm{body}}\mathrm{\left(}\mathit{k}\mathrm{\right)}}\mathrm{·}$(46)The “direct” curve corresponds to the model described in Valageas & Nishimichi (2011), without the geometrical modifications introduced in the previous Sect. 7, and was already displayed in Fig. 22 of Valageas & Nishimichi (2011). We can see that the modification P^tang(k), shown in Fig. 14 in terms of the logarithmic power Δ²(k), provides a significant improvement on the transition scales, especially for z ≤ 1. Thus, we obtain an accuracy on the order of 10% on the transition scales, and 1% on weakly nonlinear scales. On very large and very small scales the curves in Fig. 16 are dominated by the error bars of the numerical simulations.

Next, we show in Fig. 17 the relative difference between our model and the numerical simulations for the equilateral bispectrum, $\frac{\mathrm{\Delta }{\mathit{B}}_{\mathrm{eq}}}{{\mathit{B}}_{\mathrm{eq}}}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathrm{=}\frac{\mathrm{\left|}{\mathit{B}}_{\mathrm{model}}\mathrm{\right(}\mathit{k,k,k}\mathrm{)}\mathrm{-}{\mathit{B}}_{\mathit{N}\mathrm{-}\mathrm{body}}\mathrm{(}\mathit{k,k,k}\mathrm{\left)}\mathrm{\right|}}{{\mathit{B}}_{\mathit{N}\mathrm{-}\mathrm{body}}\mathrm{\left(}\mathit{k,k,k}\mathrm{\right)}}\mathrm{·}$(47)We can see that we reach a typical accuracy of 10% on most scales. At low k we are dominated by the error introduced by the finite size of the simulation box, as shown by the rise of the statistical error, and the theoretical predictions are actually more accurate than appears in the figure (they actually become increasingly good on larger scales where 1-loop perturbation theory is increasingly accurate). At high k we are dominated by the error from the finite resolution of the simulations. As for the power spectrum, studied in Valageas & Nishimichi (2011) and shown in Fig. 16, the accuracy of the “direct” model is worst on transition scales, especially at z = 3. This agrees with the behavior found in Figs. 5 and 8. The modified bispectrum ${\mathit{B}}_{\mathrm{eq}}^{\mathrm{tang}\mathit{.}}\mathrm{\left(}\mathit{k}\mathrm{\right)}$, shown in Fig. 15, provides a modest improvement on the transition scales, except at z = 0.35 where in our case it happens to give a slightly larger error that remains on the order of 10%, which is typical in this range, so that the change caused by the geometrical modification to the bispectrum is not meaningful here.

Although this 10% accuracy (and better on very large scales) is already a satisfactory result, it is significantly below what can be achieved for the power spectrum, where an accuracy on the order of 1% is reached in the weakly nonlinear regime, as seen in Fig. 16. From the comparison between different resummation schemes discussed in Sect. 6.1, it is not clear whether going to higher orders of perturbation theory would provide a significant improvement. It may happen that on these scales, especially in the transition regime, the bridge between the perturbative 3-halo term and the non-perturbative 2-halo and 1-halo terms is the main source of error because of the intrinsic limitations of a description in terms of relaxed spherical halos.

On the other hand, the error bars of the simulations are much larger for the bispectrum than for the power spectrum, and at k < 0.1 h Mpc^-1 the level of 10% seen in Fig. 17 is mostly set by the simulations. The comparison with Fig. 16 suggests that on these large scales the accuracy of the theoretical model is actually much better, on the order of 1%, because it is determined by the systematic perturbation theory. Thus, our model appears to be competitive with numerical simulations because its fares as well or better on both large and small scales, but worse on the intermediate mildly nonlinear scales. Of course, for practical purposes, the main advantage of the analytical model is that it allows much faster computations, as well as a greater flexibility.

Finally, we plot in Fig. 18 the relative difference between our model and the numerical simulations for the equilateral reduced bispectrum, (ΔQ_eq/Q_eq)(k), defined in a fashion similar to Eqs. (46) and (47). As in Fig. 8, the vertical arrow in the upper right part shows the wavenumber where the shot noise of the simulations becomes large.

As for the bispectrum, we obtain a typical accuracy on the order of 10% (and better on larger scales), except on the transition scales where the prediction given by the “direct” model can be larger than the numerical results by up to a factor two. This corresponds to the artificial “bump” found in Fig. 8. The modified predictions ${\mathit{Q}}_{\mathrm{eq}}^{\mathrm{tang}\mathit{.}}$ and ${\mathit{Q}}_{\mathrm{eq}}^{\mathrm{flat}}$ allow us to recover an accuracy on the order of 10% in this transition regime. On smaller scales the error bars of the numerical simulations are too large to obtain a precise estimate of the accuracy of these models because the rise seen in Fig. 18 is due to the finite resolution.

9. Conclusion

Extending a previous work dedicated to the matter power spectrum, we have explained how to combine perturbation theories with halo models to build unified models that can describe all scales, from large linear scales to small highly nonlinear scales. Starting again from a Lagrangian point of view, instead of the usual Eulerian point of view, we have shown how to recover the decomposition into 3-halo, 2-halo, and 1-halo contributions, which we relate to perturbative and non-perturbative terms. This explains how one can build a model that agrees with perturbation theory at all orders because the 1-halo and 2-halo terms are non-perturbative corrections that vanish at all orders over P_L. Moreover, we explained how new counterterms appear in the 1-halo and 2-halo contributions. This ensures that these contributions vanish at low k, as required by physical constraints on the power generated by small-scale redistributions of matter. This improves previous models that displayed an unphysical constant asymptote at low k, and allows us to reach a higher accuracy. In addition to standard perturbation theory, we described two alternative perturbative schemes, also complete up to 1-loop order, which can be used for the perturbative 3-halo contribution. They contain infinite partial resummations of higher order diagrams.

Combining the halo model used in Valageas & Nishimichi (2011) for the matter power spectrum with the 1-loop standard perturbation theory, we obtain a good agreement with numerical simulations for the bispectrum without any new free parameter. We consider the bispectrum as well as the reduced bispectrum, using for the latter the power spectrum predicted by the same approach (but with the more accurate steepest-descent resummation instead of standard perturbation theory). We checked that this reasonably good match to simulations holds for equilateral as well as isosceles triangles, from large to small scales. However, as for the power spectrum, the intermediate mildly nonlinear scales are not as well reproduced by this direct implementation of our approach.

The comparison between the three perturbative schemes investigated in this article shows that the standard 1-loop perturbation theory is actually the most accurate one. Because it is also simpler and faster to compute, this is also the most efficient one. This is quite different from the power spectrum, where the standard 1-loop perturbation theory is not the most accurate on large scales and behaves badly on small scales (it grows too fast at high k), so that it cannot be used in unified models unless one adds at least an external high-k cutoff. This problem does not occur for the bispectrum because the standard 1-loop contribution is now negligible on small scales compared with the 1-halo contribution. However, if one pushes the perturbative contribution to higher orders, it may start being too large at high k so that one would need to resort to alternative, better behaved, perturbative approaches, or to add high-k cutoffs.

Next, we have shown how to improve our predictions on the transition scales for the matter power spectrum and bispectrum, using a simple interpolation scheme instead of modifications to halo parameters, which would involve new parameters and do not allow significant improvements. This method automatically detects this transition regime from the shape of the reduced equilateral bispectrum, and in this fashion adapts to the change of initial conditions. Thus, for CDM linear power spectra that are not pure power laws, the transition interval  [k₋,k₊]  shifts somewhat with redshift (i.e., with the local slope n of the linear power spectrum on the transition scale) with respect to the interval that would be defined by constant thresholds, such as Δ²(k_±) = 1 and 30. Moreover, the interpolation through tangent lines adapts to the characteristic bends of the power spectrum and bispectrum, seen for CDM initial conditions around Δ²(k) ~ 30. Since this only modifies our model on the transition interval  [k₋,k₊] , large scales are still determined by systematic perturbation theory and small scales by the halo model. Then, we obtain an accuracy on the order of 10% for the power spectrum and the bispectrum on nonlinear scales, and 1% on larger weakly nonlinear scales. The same levels of improvement and final accuracy are obtained for the real-space two-point correlation function.

Our model can still be improved in various manners. First, one may investigate other perturbative approaches because other resummation schemes may prove more accurate than the standard 1-loop perturbation theory. However, to be more efficient, they should not be much more difficult to compute than the standard perturbation theory. Alternatively, one may go to higher orders. For the power spectrum, higher orders are indeed relevant because various resummation schemes have already been shown to be more accurate (Crocce & Scoccimarro 2008; Taruya et al. 2009; Valageas & Nishimichi 2011) and non-perturbative contributions only dominate after many perturbative orders have become involved (Valageas 2011). For the bispectrum, the failure of the two resummation schemes investigated in this paper to improve over standard perturbation theory suggests that it may be more difficult to reach significant improvements.

Second, one may improve the halo model used in our unified approach. For instance, one could take into account substructures (Sheth & Jain 2003; Giocoli et al. 2010), deviations from spherical profiles (Jing & Suto 2002; Smith et al. 2006), or the effect of baryons (Guillet et al. 2010).

Our model could be extended to other initial conditions, especially non-Gaussian ones for which the bispectrum is a very useful and direct probe (Sefusatti et al. 2010). It would also be interesting to use this approach to describe velocity fields, and to take into account redshift-space distortions (Smith et al. 2008). However, we leave these tasks to future studies.

---

Acknowledgments

We thank S. Colombi for useful discussions about the folding procedure for the bispectrum measurements. 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: Behavior of the bispectrum for one low wavenumber in perturbation theory

We show in this appendix how the large-scale behavior (32) also arises at all orders of perturbation theory. As is well known, within the standard perturbation theory the density contrast field $\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{\delta }\end{array}\mathrm{\left(}{k}\mathrm{\right)}$ is written as a perturbative expansion over powers of the linear density contrast ${\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{\delta }\end{array}}_{\mathrm{L}}\mathrm{\left(}{k}\mathrm{\right)}$, $\begin{array}{}\mathrm{˜}\\ \mathit{\delta }\end{array}\mathrm{\left(}{k}\mathrm{\right)}\mathrm{=}\sum _{\mathit{n}\mathrm{=}\mathrm{1}}^{\mathrm{\infty }}{\begin{array}{}\mathrm{˜}\\ \mathit{\delta }\end{array}}^{\mathrm{\left(}\mathit{n}\mathrm{\right)}}\mathrm{\left(}{k}\mathrm{\right)}\mathit{,}$(A.1)with $\begin{array}{ccc}{\begin{array}{}\mathrm{˜}\\ \mathit{\delta }\end{array}}^{\mathrm{\left(}\mathit{n}\mathrm{\right)}}\mathrm{\left(}{k}\mathrm{\right)}& \mathrm{=}& \int \mathrm{d}{{k}}_{\mathrm{1}}\mathit{...}\mathrm{d}{{k}}_{\mathit{n}}\hspace{0.17em}{\mathit{\delta }}_{\mathit{D}}\mathrm{(}{{k}}_{\mathrm{1}}\mathrm{+}\mathit{...}\mathrm{+}{{k}}_{\mathit{n}}\mathrm{-}{k}\mathrm{)}\hspace{0.17em}{\begin{array}{}\mathrm{˜}\\ \mathit{F}\end{array}}_{\mathit{n}}\mathrm{\left(}{{k}}_{\mathrm{1}}\mathit{,}\mathit{...}\mathit{,}{{k}}_{\mathit{n}}\mathrm{\right)}\\ & & \mathrm{\times }\hspace{0.17em}{\begin{array}{}\mathrm{˜}\\ \mathit{\delta }\end{array}}_{\mathrm{L}}\mathrm{\left(}{{k}}_{\mathrm{1}}\mathrm{\right)}\mathit{...}{\begin{array}{}\mathrm{˜}\\ \mathit{\delta }\end{array}}_{\mathrm{L}}\mathrm{\left(}{{k}}_{\mathit{n}}\mathrm{\right)}\mathit{.}\end{array}$(A.2)The kernels F_n can be obtained from recursion relations, which are derived from the equations of motion (Goroff et al. 1986; Bernardeau et al. 2002a), and can be chosen as symmetric over the wavenumbers  {k₁,...,k_n} . Then, the bispectrum can be obtained within perturbation theory from the three-point correlation $\mathrm{⟨}\begin{array}{}\mathrm{˜}\\ \mathit{\delta }\end{array}\mathrm{\left(}{{k}}_{\mathrm{1}}\mathrm{\right)}\begin{array}{}\mathrm{˜}\\ \mathit{\delta }\end{array}\mathrm{\left(}{{k}}_{\mathrm{2}}\mathrm{\right)}\begin{array}{}\mathrm{˜}\\ \mathit{\delta }\end{array}\mathrm{\left(}{{k}}_{\mathrm{3}}\mathrm{\right)}{\mathrm{⟩}}_{\mathrm{c}}\mathrm{=}\sum _{{\mathit{n}}_{\mathrm{1}}\mathit{,}{\mathit{n}}_{\mathrm{2}}\mathit{,}{\mathit{n}}_{\mathrm{3}}}\mathrm{⟨}{\begin{array}{}\mathrm{˜}\\ \mathit{\delta }\end{array}}^{\mathrm{\left(}{\mathit{n}}_{\mathrm{1}}\mathrm{\right)}}\mathrm{\left(}{{k}}_{\mathrm{1}}\mathrm{\right)}{\begin{array}{}\mathrm{˜}\\ \mathit{\delta }\end{array}}^{\mathrm{\left(}{\mathit{n}}_{\mathrm{2}}\mathrm{\right)}}\mathrm{\left(}{{k}}_{\mathrm{2}}\mathrm{\right)}{\begin{array}{}\mathrm{˜}\\ \mathit{\delta }\end{array}}^{\mathrm{\left(}{\mathit{n}}_{\mathrm{3}}\mathrm{\right)}}\mathrm{\left(}{{k}}_{\mathrm{3}}\mathrm{\right)}{\mathrm{⟩}}_{\mathrm{c}}\mathit{,}$(A.3)where the subscript c recalls that we only take the connected part of the Gaussian average. Using Wick’s theorem and defining the linear power spectrum as in Eq. (5) gives $\begin{array}{ccc}\mathrm{⟨}\begin{array}{}\mathrm{˜}\\ \mathit{\delta }\end{array}\mathrm{\left(}{{k}}_{\mathrm{1}}\mathrm{\right)}\begin{array}{}\mathrm{˜}\\ \mathit{\delta }\end{array}\mathrm{\left(}{{k}}_{\mathrm{2}}\mathrm{\right)}\begin{array}{}\mathrm{˜}\\ \mathit{\delta }\end{array}\mathrm{\left(}{{k}}_{\mathrm{3}}\mathrm{\right)}{\mathrm{⟩}}_{\mathrm{c}}& \mathrm{=}& \sum _{{\mathit{n}}_{\mathrm{\ast }}}\mathrm{\int }\underset{\mathit{j}}{\mathrm{􏽙}}\mathrm{d}{{k}}_{\mathit{j}}^{\mathrm{\left(}{\mathit{n}}_{\mathrm{\ast }}\mathrm{\right)}}{\mathit{\delta }}_{\mathit{D}}\mathrm{(}{{k}}_{\mathit{j}}^{\mathrm{\left(}{\mathit{n}}_{\mathrm{12}}\mathrm{\right)}}\mathrm{+}{{k}}_{\mathit{j}}^{\mathrm{\left(}{\mathit{n}}_{\mathrm{13}}\mathrm{\right)}}\mathrm{-}{{k}}_{\mathrm{1}}\mathrm{)}\\ & & \mathrm{\times }\hspace{0.17em}{\mathit{\delta }}_{\mathit{D}}\mathrm{(}\mathrm{-}{{k}}_{\mathit{j}}^{\mathrm{\left(}{\mathit{n}}_{\mathrm{12}}\mathrm{\right)}}\mathrm{+}{{k}}_{\mathit{j}}^{\mathrm{\left(}{\mathit{n}}_{\mathrm{23}}\mathrm{\right)}}\mathrm{-}{{k}}_{\mathrm{2}}\mathrm{)}{\mathit{\delta }}_{\mathit{D}}\mathrm{(}\mathrm{-}{{k}}_{\mathit{j}}^{\mathrm{\left(}{\mathit{n}}_{\mathrm{13}}\mathrm{\right)}}\mathrm{-}{{k}}_{\mathit{j}}^{\mathrm{\left(}{\mathit{n}}_{\mathrm{23}}\mathrm{\right)}}\mathrm{-}{{k}}_{\mathrm{3}}\mathrm{)}\\ & & \mathrm{\times }\hspace{0.17em}{\begin{array}{}\mathrm{˜}\\ \mathit{F}\end{array}}_{\mathrm{2}{\mathit{n}}_{\mathrm{11}}\mathrm{+}{\mathit{n}}_{\mathrm{12}}\mathrm{+}{\mathit{n}}_{\mathrm{13}}}\mathrm{(}{{k}}_{\mathit{j}}^{\mathrm{\left(}{\mathit{n}}_{\mathrm{11}}\mathrm{\right)}}\mathit{,}\mathrm{-}{{k}}_{\mathit{j}}^{\mathrm{\left(}{\mathit{n}}_{\mathrm{11}}\mathrm{\right)}}\mathit{,}{{k}}_{\mathit{j}}^{\mathrm{\left(}{\mathit{n}}_{\mathrm{12}}\mathrm{\right)}}\mathit{,}{{k}}_{\mathit{j}}^{\mathrm{\left(}{\mathit{n}}_{\mathrm{13}}\mathrm{\right)}}\mathrm{)}\\ & & \mathrm{\times }\hspace{0.17em}{\begin{array}{}\mathrm{˜}\\ \mathit{F}\end{array}}_{\mathrm{2}{\mathit{n}}_{\mathrm{22}}\mathrm{+}{\mathit{n}}_{\mathrm{12}}\mathrm{+}{\mathit{n}}_{\mathrm{23}}}\mathrm{(}{{k}}_{\mathit{j}}^{\mathrm{\left(}{\mathit{n}}_{\mathrm{22}}\mathrm{\right)}}\mathit{,}\mathrm{-}{{k}}_{\mathit{j}}^{\mathrm{\left(}{\mathit{n}}_{\mathrm{22}}\mathrm{\right)}}\mathit{,}\mathrm{-}{{k}}_{\mathit{j}}^{\mathrm{\left(}{\mathit{n}}_{\mathrm{12}}\mathrm{\right)}}\mathit{,}{{k}}_{\mathit{j}}^{\mathrm{\left(}{\mathit{n}}_{\mathrm{23}}\mathrm{\right)}}\mathrm{)}\\ & & \mathrm{\times }\hspace{0.17em}{\begin{array}{}\mathrm{˜}\\ \mathit{F}\end{array}}_{\mathrm{2}{\mathit{n}}_{\mathrm{33}}\mathrm{+}{\mathit{n}}_{\mathrm{13}}\mathrm{+}{\mathit{n}}_{\mathrm{23}}}\mathrm{(}{{k}}_{\mathit{j}}^{\mathrm{\left(}{\mathit{n}}_{\mathrm{33}}\mathrm{\right)}}\mathit{,}\mathrm{-}{{k}}_{\mathit{j}}^{\mathrm{\left(}{\mathit{n}}_{\mathrm{33}}\mathrm{\right)}}\mathit{,}\mathrm{-}{{k}}_{\mathit{j}}^{\mathrm{\left(}{\mathit{n}}_{\mathrm{13}}\mathrm{\right)}}\mathit{,}\mathrm{-}{{k}}_{\mathit{j}}^{\mathrm{\left(}{\mathit{n}}_{\mathrm{23}}\mathrm{\right)}}\mathrm{)}\\ & & \mathrm{\times }\hspace{0.17em}\underset{\mathit{j}}{\mathrm{􏽙}}{\mathit{P}}_{\mathrm{L}}\mathrm{\left(}{\mathit{k}}_{\mathit{j}}^{\mathrm{\left(}{\mathit{n}}_{\mathrm{\ast }}\mathrm{\right)}}\mathrm{\right)}\mathit{,}\end{array}$(A.4)where we used a short-hand notation for the diagram shown in Fig. A.1, and the Dirac factors contain sums over the wavenumbers ${{k}}_{\mathit{j}}^{\mathrm{\left(}{\mathit{n}}_{\mathrm{\ast }}\mathrm{\right)}}$. Of course, the three Dirac factors can be combined to factorize out a Dirac prefactor δ_D(k₁ + k₂ + k₃), in agreement with Eq. (6). Taking the connected part means that each of the three circles in Fig. A.1 is connected to at least one other circle (for instance the case n₁₂ = n₁₃ = n₂₂ = 0 is excluded).

We are interested in the limit k₁ → 0, at fixed k₂ and k₃. A well-known property of the kernels ${\begin{array}{}{\mathrm{˜}}_{}\\ \mathit{F}\end{array}}_{\mathit{n}}$, which arises from momentum conservation (Peebles 1974; Goroff et al. 1986), is that ${\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{F}\end{array}}_{\mathit{n}}\mathrm{(}{{k}}_{\mathrm{1}}\mathrm{+}\mathit{...}\mathrm{+}{{k}}_{\mathit{n}}\mathrm{)}\mathrm{\propto }{\mathit{k}}^{\mathrm{2}}$ as k = k₁ + ... + k_n goes to zero while the individual k_j do not, for all n ≥ 2. Then, from the first Dirac factor in Eq. (A.4) we can see that the first kernel $\begin{array}{}{\mathrm{˜}}_{}\\ \mathit{F}\end{array}$ behaves as ${\mathit{k}}_{\mathrm{1}}^{\mathrm{2}}$ if 2n₁₁ + n₁₂ + n₁₃ ≥ 2. Therefore, in the limit k₁ → 0 the dominant contributions come from the diagrams with 2n₁₁ + n₁₂ + n₁₃ = 1, since F₁ ≡ 1 and for CDM cosmologies ${\mathit{P}}_{\mathrm{L}}\mathrm{\left(}{\mathit{k}}_{\mathrm{1}}\mathrm{\right)}\mathrm{\propto }{\mathit{k}}_{\mathrm{1}}^{{\mathit{n}}_{\mathit{s}}}$ with n_s ≃ 1 at low k₁. This corresponds to either  {n₁₁ = 0,n₁₂ = 1,n₁₃ = 0}  or  {n₁₁ = 0,n₁₂ = 0,n₁₃ = 1}  because we only take the connected part of (A.4). Thus, for k₁ → 0 the dominant contribution to the perturbative bispectrum reads as $\begin{array}{ccc}{\mathit{B}}_{\mathrm{pert}}\mathrm{\left(}{\mathit{k}}_{\mathrm{1}}\mathit{,}{\mathit{k}}_{\mathrm{2}}\mathit{,}{\mathit{k}}_{\mathrm{3}}\mathrm{\right)}& \mathrm{~}& {\mathit{P}}_{\mathrm{L}}\mathrm{\left(}{\mathit{k}}_{\mathrm{1}}\mathrm{\right)}\hspace{0.17em}\underset{{\mathit{k}}_{\mathrm{1}}\mathrm{\to }\mathrm{0}}{\lim }\{\sum _{{\mathit{n}}_{\mathrm{\ast }}}\mathrm{\int }\underset{\mathit{j}}{\mathrm{􏽙}}\mathrm{d}{{k}}_{\mathit{j}}^{\mathrm{\left(}{\mathit{n}}_{\mathrm{\ast }}\mathrm{\right)}}\\ & & \mathrm{\times }\hspace{0.17em}{\mathit{\delta }}_{\mathit{D}}\mathrm{(}{{k}}_{\mathit{j}}^{\mathrm{\left(}{\mathit{n}}_{\mathrm{23}}\mathrm{\right)}}\mathrm{-}{{k}}_{\mathrm{2}}\mathrm{)}\hspace{0.17em}{\begin{array}{}\mathrm{˜}\\ \mathit{F}\end{array}}_{\mathrm{2}{\mathit{n}}_{\mathrm{33}}\mathrm{+}{\mathit{n}}_{\mathrm{23}}}\mathrm{(}{{k}}_{\mathit{j}}^{\mathrm{\left(}{\mathit{n}}_{\mathrm{33}}\mathrm{\right)}}\mathit{,}\mathrm{-}{{k}}_{\mathit{j}}^{\mathrm{\left(}{\mathit{n}}_{\mathrm{33}}\mathrm{\right)}}\mathit{,}\mathrm{-}{{k}}_{\mathit{j}}^{\mathrm{\left(}{\mathit{n}}_{\mathrm{23}}\mathrm{\right)}}\mathrm{)}\\ & & \mathrm{\times }\hspace{0.17em}{\begin{array}{}\mathrm{˜}\\ \mathit{F}\end{array}}_{\mathrm{2}{\mathit{n}}_{\mathrm{22}}\mathrm{+}\mathrm{1}\mathrm{+}{\mathit{n}}_{\mathrm{23}}}\mathrm{(}{{k}}_{\mathit{j}}^{\mathrm{\left(}{\mathit{n}}_{\mathrm{22}}\mathrm{\right)}}\mathit{,}\mathrm{-}{{k}}_{\mathit{j}}^{\mathrm{\left(}{\mathit{n}}_{\mathrm{22}}\mathrm{\right)}}\mathit{,}\mathrm{-}{{k}}_{\mathrm{1}}\mathit{,}{{k}}_{\mathit{j}}^{\mathrm{\left(}{\mathit{n}}_{\mathrm{23}}\mathrm{\right)}}\mathrm{)}\\ & & \mathrm{\times }\hspace{0.17em}\underset{\mathit{j}}{\mathrm{􏽙}}{\mathit{P}}_{\mathrm{L}}\mathrm{\left(}{\mathit{k}}_{\mathit{j}}^{\mathrm{\left(}{\mathit{n}}_{\mathrm{\ast }}\mathrm{\right)}}\mathrm{\right)}\mathrm{+}\mathrm{1}\mathrm{perm}\mathit{.}\}\mathit{,}\end{array}$(A.5)where we note by “1 perm” the symmetric contribution with respect to  {k₂ ↔ k₃} . There is a subtlety in Eq. (A.5), because the limit at low k₁ of the kernel ${\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{F}\end{array}}_{\mathrm{2}{\mathit{n}}_{\mathrm{22}}\mathrm{+}\mathrm{1}\mathrm{+}{\mathit{n}}_{\mathrm{23}}}\mathrm{(}{{k}}_{\mathit{j}}^{\mathrm{\left(}{\mathit{n}}_{\mathrm{22}}\mathrm{\right)}}\mathit{,}\mathrm{-}{{k}}_{\mathit{j}}^{\mathrm{\left(}{\mathit{n}}_{\mathrm{22}}\mathrm{\right)}}\mathit{,}\mathrm{-}{{k}}_{\mathrm{1}}\mathit{,}{{k}}_{\mathit{j}}^{\mathrm{\left(}{\mathit{n}}_{\mathrm{23}}\mathrm{\right)}}\mathrm{)}$ contains divergent terms of the form $\mathrm{(}{{k}}_{\mathrm{1}}\mathrm{·}{{k}}_{\mathit{j}}\mathrm{)}\mathit{/}{\mathit{k}}_{\mathrm{1}}^{\mathrm{2}}$ (Goroff et al. 1986; Scoccimarro & Frieman 1996). However, taking the limit k₁ → 0 elsewhere, we obtain k₃ = −k₂ so that the symmetric term with respect to  {k₂ ↔ k₃}  actually corresponds to a change of sign of all wavenumbers. This cancels the divergences of the form $\mathrm{(}{{k}}_{\mathrm{1}}\mathrm{·}{{k}}_{\mathit{j}}\mathrm{)}\mathit{/}{\mathit{k}}_{\mathrm{1}}^{\mathrm{2}}$ so that the limit k₁ → 0 is finite. The fact that such infrared divergences cancel out for equal-time statistics can be traced to the Galilean invariance of the equations of motion (Jain & Bertschinger 1996; Valageas 2004).

Terms at all orders of perturbation theory contribute to the partial large-scale limit (A.5) through the nonlinear corrections to $\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{\delta }\end{array}\mathrm{\left(}{{k}}_{\mathrm{2}}\mathrm{\right)}$ and $\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{\delta }\end{array}\mathrm{\left(}{{k}}_{\mathrm{3}}\mathrm{\right)}$. In contrast, as we have seen above, only the linear term ${\begin{array}{}{}_{\mathrm{˜}}\\ \mathit{\delta }\end{array}}_{\mathrm{L}}\mathrm{\left(}{{k}}_{\mathrm{1}}\mathrm{\right)}$ contributes. The behavior (A.5) can also be understood from the physical arguments discussed in Sect. 2.3, and the explicit result (A.5) confirms that this asymptotic behavior should be preserved by both high-order perturbative terms and non-perturbative corrections.

It is interesting to note that the “renormalization” of the prefactor to the P(k₁) tail by the higher order terms in Eq. (A.5) never occurs for the power spectrum P(k), where we recover P(k) → P_L(k) at low k. There, as seen from the analysis described in Valageas (2002), perturbative terms beyond linear order scale at least as k²P_L(k) at low k. This can also be understood from the fact that for the power spectrum there are no partial large-scale limits. On the other hand, if all wavenumbers go to zero, as k₁ ~ k₂ ~ k₃ ~ k with k → 0, higher order contributions to the bispectrum scale at least as k²P_L(k)² while the “tree-order” result scales as P_L(k)², so that in the simultaneous large-scale limit we also recover the lowest order result.

Appendix B: Two-point correlation function

Although the main focus of this paper is the computation of the bispectrum, we have seen in Sects. 7 and 8 that the shape of the reduced bispectrum can be used to improve the model devised for the power spectrum. This leads in turn to an improved model for the real-space two-point correlation function, which is given by $\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}{\mathit{k}}^{\mathrm{2}}\hspace{0.17em}\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}}{\mathrm{\Delta }}^{\mathrm{2}}\mathrm{\left(}\mathit{k}\mathrm{\right)}\hspace{0.17em}\frac{\sin \mathrm{\left(}\mathit{kx}\mathrm{\right)}}{\mathit{kx}}\mathrm{·}\end{array}$Therefore, we compare in this appendix the two-point correlation functions obtained either with our “direct” model, already studied in Valageas & Nishimichi (2011), or with the geometrical modification “tang” to the power spectrum explained in Fig. 14. Thus, we plot in Figs. B.1 and B.2 the two-point correlation functions defined by the two power spectra of Fig. 14 through the Fourier transform (B.2). As expected, we can see in Fig. B.1 that, as for the power Δ²(k) in Fig. 14, the geometrical modification “tang” corrects most of the artificial “dip” that was shown by the “direct” model on transition scales (i.e., the mildly nonlinear regime). The curves even look closer to the simulation results and more regular because the integral (B.2) has regularized the geometrical modification. Indeed, while the power ${\mathrm{\Delta }}_{\mathrm{tang}}^{\mathrm{2}}\mathrm{\left(}\mathit{k}\mathrm{\right)}$ shown in Fig. 14 was only continuous, with a discontinuous derivative at k₋, the correlation ξ_tang(x) has a well-defined and finite derivative over all x > 0.

Contrary to the power spectrum, the modification “tang” to the two-point correlation is not restricted to a finite range  [x₋,x₊]  because the modification of the power on the wavenumber interval  [k₋,k₊]  contributes to all distances x through the integral (B.2). However, as expected, we can check in Fig. B.1 that on very small and very large scales the modified correlation ξ_tang(x) becomes increasingly close to the original model prediction ξ_direct(x), and as such shows a good agreement with numerical simulations.

Nevertheless, to check in more detail that the modification “tang” does not spread on large scales, where the initial power is much weaker, we show in Fig. B.2 the correlation functions of Fig. B.1 on larger scales, around the baryon acoustic oscillation. We can see that both correlations, ξ_tang(x) and ξ_direct(x), are very close and within the error bars of the numerical simulations. In particular, they capture very well the departure from the linear correlation, mostly due to the 1-loop perturbative contribution. Therefore, the modified correlation ξ_tang(x) provides a good model from small to very large scales.

As in Fig. 16 for the power spectrum, we show in Fig. B.3 the relative accuracy of our analytical models and numerical simulations. (The curve obtained for the “direct” model was already shown in Valageas & Nishimichi (2011).) In agreement with Figs. 16 and B.1, we can see that the modification “tang” provides a significant improvement on transition scales. Thus, it ensures an accuracy of 10% or better on nonlinear scales, while an accuracy on the order of 1% is again reached on quasilinear scales, associated for instance with the baryon acoustic oscillation of Fig. B.2.

All Figures