Issue 
A&A
Volume 532, August 2011



Article Number  A4  
Number of page(s)  23  
Section  Cosmology (including clusters of galaxies)  
DOI  https://doi.org/10.1051/00046361/201116638  
Published online  11 July 2011 
Combining perturbation theories with halo models for the matter bispectrum
^{1}
Institut de Physique Théorique, CEA Saclay,
91191
GifsurYvette,
France
email: patrick.valageas@cea.fr
^{2}
Institute for the Physics and Mathematics of the Universe,
University of Tokyo, Kashiwa, Chiba
2778568,
Japan
Received:
3
February
2011
Accepted:
30
May
2011
Aims. We investigate how unified models should be built to be able to predict the matterdensity bispectrum (and power spectrum) from very large to small scales and that are at the same time consistent with perturbation theory at low k and with halo models at high k.
Methods. We use a Lagrangian framework to decompose the bispectrum into “3halo”, “2halo”, and “1halo” contributions, related to “perturbative” and “nonperturbative” terms. We describe a simple implementation of this approach and present a detailed comparison with numerical simulations.
Results. We show that the 1halo and 2halo contributions contain counterterms that ensure their decay at low k, as required by physical constraints, and allow a better match to simulations. Contrary to the power spectrum, the standard 1loop perturbation theory can be used for the perturbative 3halo contribution because it does not grow too fast at high k. Moreover, it is much simpler and more accurate than two resummation schemes investigated in this paper. We obtain a good agreement with numerical simulations on both large and small scales, but the transition scales are poorly described by the simplest implementation. This cannot be amended by simple modifications to the halo parameters, but we show how it can be corrected for the power spectrum and the bispectrum through a simple interpolation scheme that is restricted to this intermediate regime. Then, we reach an accuracy on the order of 10% on mildly and highly nonlinear scales, while an accuracy on the order of 1% is obtained on larger weakly nonlinear scales. This also holds for the realspace twopoint correlation function.
Key words: largescale structure of Universe
© ESO, 2011
1. Introduction
According to standard cosmological scenarios, the largescale structures of the present Universe have formed through the amplification by gravitational instability of small almostGaussian primordial fluctuations (Peebles 1980). Then, from observations of the recent Universe, through galaxy surveys (Tegmark et al. 2006; Cole et al. 2005), weaklensing 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 twopoint 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 nonGaussianities of the initial perturbations, to check the gravitational clustering scenario (through the measure of the nonGaussianities generated by the nonlinear dynamics), and to break parameter degeneracies, it is useful to study higher order statistics. The threepoint 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 threepoint (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 nonGaussian 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 “1halo” and “2halo” 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 1halo, 2halo, and 3halo 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 3halo, 2halo, and 1halo terms from a Lagrangian point of view, and we make the connection with perturbation theory.
2.1. Lagrangianspace 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, (1)where is the mean comoving matter density of the Universe and we work in comoving coordinates. Then, defining the density contrast as (2)and its Fourier transform as (3)we obtain from Eq. (1) (4)We define the power spectrum, P(k), and the bispectrum, B(k_{1},k_{2},k_{3}), which are the Fourier transforms of the twopoint and threepoint density correlation functions, by 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_{1}) and B(k_{1},k_{2},k_{3}) only depend on the lengths k_{1} and {k_{1},k_{2},k_{3}} .
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 “nonperturbative” contributions. The former is associated with configurations of particle triplets (or puplets for the ppoint 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 nonperturbative contribution through a halo model. Then, the Lagrangian framework allows us to recover the counterterms that arise in the halomodel 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 wellknown expression (Schneider & Bartelmann 1995; Taylor & Hamilton 1996) (7)where we introduced the Eulerianspace 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 (8)where we introduced Δx_{2} = x(q_{2}) − x(0) and Δx_{3} = x(q_{3}) − x(0). However, the symmetry over {k_{1},k_{2},k_{3}} 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 thirdorder 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 (9)where x_{j} = x(q_{j}). Next, we split the average (9) into three contributions, associated with the cases where the triplet {q_{1},q_{2},q_{3}} belongs to either one, two, or three halos: (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. “1halo” contribution
Let us first consider the 1halo contribution. It can be written as (11)where the sum runs over all halos, labeled by the index α, and the tophat 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_{1},q_{2},q_{3}} belongs to a single halo. Equation (11) also reads as (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), (13)Here 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 , but we could as well count them in Eulerian space by using the Eulerianspace center of mass . 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, (14)which we write as (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 (16)and (17)where is the Fourier transform of the tophat of radius q_{M}, defined as (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 (19)where the average ⟨ ... ⟩ _{qc,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, ρ_{xc,M}(x), as (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, (21)we obtain (22)Making the approximation that these halos are also spherically symmetric in Lagrangian space, we can use Eq. (18) for the factors e^{−ikj·qj} and Eq. (19) reads as (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^{−ikj·Ψ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_{1} + k_{2} + k_{3}), in agreement with Eq. (6). Then, the displacement Ψ^{c} only appears in mixed terms, such as . Going back to Eq. (4), we can see that the term arises from a factor δ_{D}(k_{3}), so that in the regime where the 1halo contribution is dominant, that is , it should reduce to δ_{D}(k_{3}). In Eq. (23) this corresponds to the fact that for k_{3} ≫ 1/q_{M}. Then, thanks to the prefactor δ_{D}(k_{1} + k_{2} + k_{3}) we can see that in this regime k_{1} + k_{2} → 0 and e^{−i(k1 + k2)·Ψ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_{1} + k_{2} + k_{3}), we obtain the 1halo contribution to the bispectrum as (24)Thus, as for the power spectrum studied in Valageas & Nishimichi (2011), we recover the standard expression of the 1halo 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 . Note that for the power spectrum we obtained in Valageas & Nishimichi (2011) a factor , instead of the factor 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_{1} + k_{2}) 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 higherorder functions, such as the bispectrum or the trispectrum, the symmetric Eq. (9) (and its Npoint generalization) provides a simpler starting point. In particular, the counting of the volumes over {q_{1},...,q_{N}} associated with a single halo simply factorizes as , whereas introducing relative positions makes the geometrical countings slightly more intricate.
The factors cannot be interpreted as modifications of the halo profiles ρ_{M}(x) under the form of compensated profiles (the result 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 1halo contribution decays at least as 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 and the functions and are regular functions of k^{2} at the origin, (25)This agrees with the requirements associated with a smallscale redistribution of matter. Here, there is a simplification as compared with the matter power spectrum P(k). Indeed, as discussed in Peebles (1974), smallscale redistributions of matter generically give a lowk tail , whence P(k) ∝ k^{2}, while taking into account momentum conservation implies the steeper decay , whence P(k) ∝ k^{4}. Thus, in Valageas & Nishimichi (2011) we recovered the expected k^{2} tail, through the factor , since our analysis did not enforce momentum conservation. Using the approach described above, we would actually obtain a k^{4} tail, through the factor , 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 higherorder correlations) smallscale redistributions of matter lead to the same 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 (26)As noticed in previous works (Cooray & Sheth 2002), in the usual version of the halo model the 1halo contribution to the bispectrum does not contain the counterterms , so that it goes to a strictly positive constant at large scales and dominates over the lowestorder 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)^{2} ∝ k^{2ns} 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. “2halo” contribution
We now turn to the 2halo contribution to the bispectrum, following the procedure described in the previous section for the 1halo term. Thus, as in Eq. (12) we can write the 2halo contribution as (27)where we sum over all pairs of distinct halos, of mass M_{1} and M_{2}, which contain one and two of the particles q_{j}. Here we have written the contribution where the halo M_{1} contains the particle q_{1}, and the two additional contributions, noted “2 cyc.”, correspond to the cases where M_{1} contains the particle q_{2} or q_{3}. Using again the approximation (22) of fully relaxed halos, the average of Eq. (27) writes as (28)where is the twopoint correlation function of halos of mass M_{1} and M_{2}, 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 + ξ_{12}) vanishes (because the integral over yields δ_{D}(k_{1}) and for k_{1} = 0), so that only the term associated with the largescale correlation of halos remains. Writing this twopoint correlation in terms of the halo power spectrum, (29)the integrals over and give the Dirac factors δ_{D}(k + k_{1})δ_{D}(k − k_{2} − k_{3}). Therefore, we recover the Dirac prefactor δ_{D}(k_{1} + k_{2} + k_{3}), which we can factorize out to write the 2halo contribution to the bispectrum as (30)Then, as is customary, we could write the halo power spectrum as P_{12}(k) = b(M_{1})b(M_{2})P(k), where b(M) is a massdependent bias factor (of order unity for typical halos) that we approximate as scaleindependent (thereby neglecting exclusion constraints between halos).
However, the expression (30) is not really satisfactory. Indeed, for k_{1} → 0 it decreases as because of the prefactor , while we would expect to recover a behavior proportional to P(k_{1}). Indeed, in this regime this 2halo contribution should describe the power that arises from the largescale correlation between density fluctuations at a position q_{1} and at a position q_{2}, with q_{2} − q_{1} ~ 2π/k_{1}, the fluctuations at q_{2} being furthermore decomposed over smallerscale fluctuations within single halos of mass M_{2}. In other words, the largescale power P(k_{1}) should not be damped by the prefactor . Going back to Eq. (28), we can see that the product involves a prefactor , 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 1halo contribution (23), where the product of three terms gave the prefactor e^{i(k1 + k2 + k3)·Ψc}, which had to disappear because it does not depend on relative displacements, and was indeed irrelevant thanks to the Dirac prefactor δ_{D}(k_{1} + k_{2} + k_{3}). This means that the approximation Ψ^{c} = 0 is not as good for the 2halo contribution as for the 1halo 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 2halo contribution that behaves in a reasonable fashion, we choose to use instead of Eq. (30) the simple expression (31)Thus, we have made the approximation P_{12}(k) ≃ P_{L}(k), where P_{L}(k) is the linear matter density power spectrum, and we have set the spurious prefactor 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 2halo contribution is subdominant on most scales of interest, as shown by the numerical results in Sect. 4.
In the largescale limit we obtain from Eq. (31) the asymptotic behaviors (32)when only one wavenumber goes to zero, and (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^{ns} 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_{1}, is valid at all orders of perturbation theory. In this regime, the 2halo contribution, which is nonperturbative (as discussed in Sect. 2.4 below), is therefore on the same order (i.e., ~P_{L}(k_{1})) as the perturbative contribution. This can be understood from the arguments developed above to obtain the expression (31), as the nonperturbative effects associated with the formation of smallscale nonlinear structures around the nearby points x_{2} and x_{3} (or q_{2} and q_{3} in Lagrangian space) do not modify the largerscale correlation with the very distant point x_{1}. That is, the nonperturbative effects associated with the formation of a halo correspond to a local redistribution of matter. Then, at leading order, nonperturbative terms simply lead to a “renormalization” of the coefficient that multiplies the factor P_{L}(k_{1}), by an amount that depends on how far in the nonlinear regime the wavenumbers k_{2} and k_{3} are.
When all wavenumbers go to zero, the property (33) again ensures that the 2halo contribution becomes negligible as compared with the lowest order term of perturbation theory, which scales as P_{L}(λ)^{2}. Thus, as for the 1halo contribution, the counterterms 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. “3halo” contribution as a perturbative contribution
For the 3halo 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 2halo term, it is not straightforward to take into account largescale 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 1halo and 2halo 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 largemass tail of the form , which vanishes exponentially for P_{L} → 0 and has a Taylor expansion over powers of P_{L} that is identically zero. Therefore, the remaining 3halo contribution to Eq. (10) must be consistent with perturbation theory up to all orders, and we simply use the approximation (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 nonperturbative terms such as , 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}, (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 1halo and 2halo terms become of the same order as the 3halo term. In view of the approximations involved in these 1halo and 2halo 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).
Fig. 1 Probability F_{3H} that the triplets associated with the equilateral bispectrum B_{eq}(k) belong to three different halos, from Eq. (35), at redshifts z = 0.35,1, and 3. 
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 wellbehaved on small scales. As we will check in Sect. 4, for the bispectrum it turns out that at oneloop order the standard perturbation theory remains small at high k compared with the 1halo and 2halo contributions, which makes this an efficient choice; but we will also study the direct steepestdescent 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 3halo term is rather written as (Ma & Fry 2000b; Scoccimarro et al. 2001; Cooray & Sheth 2002) (36)where is the average over mass, weighted by the halo mass function, of the bias parameter b and halo profile ; 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 3halo 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 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 wellcontrolled 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 1loop order, while most previous studies that involved the halo model only used the treelevel 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 2halo term that mixes largerscale and intrahalo wavenumbers we had to partially bypass the naive prediction of this halo model to recover the largescale 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 Npoint correlation functions or polyspectra 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 in the 1halo contribution (24), which ensure the largescale behaviors (25), (26) that were missed in previous Eulerian implementations.
As explained above, for the 2halo 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 halomodel approximations, we can still split a constantdensity 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 , 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 for this uniform system. Then, the counterterms in Eqs. (24) and (31) clearly satisfy the constraint B(k_{1},k_{2},k_{3}) = 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 would also satisfy this constraint (but not the properties (25)).
Although the main reason for taking into account the counterterms in the 1halo and 2halo 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 1halo and 2halo terms (which do not include these counterterms) are still too large and spoil the agreement of (treeorder) perturbation theory with numerical simulations. Giving a faster decay on large scales of these 1halo and 2halo contributions, the counterterms 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) (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 (Scoccimarro et al. 2001). In particular, we truncate the halo profiles at the density contrast δ(<r_{200}) = 200, which defines their radius r_{200}. For the concentration parameter we take (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
Fig. 2 Diagrams associated with the standard perturbation theory as obtained from a pathintegral formalism up to 1loop order for the threepoint correlation C_{3}. Although they are different from those obtained by the standard approach, their sum at each order is identical to the sum of the standard diagrams of the same order over P_{L}(k). The lines are the linear twopoint correlation C_{L} (blue solid line) and the linear response R_{L} (red solid line with an arrow). The black dots are the threeleg vertex K_{s} that enters the quadratic term of the equation of motion. The numbers are the multiplicity factors of each diagram. The treeorder diagram a) gives Eq. (39), for the density bispectrum, while the 1loop diagrams b), ..., i), give the contribution of order . 
For the perturbative contribution (34) to the bispectrum, we investigate both the standard perturbation theory and the “direct steepestdescent” resummation scheme introduced (among other perturbative expansions) in Valageas (2007a). We only go up to 1loop order for both methods, so that the first approach only contains all 1loop diagrams, while the second approach also includes partial resummations of diagrams at all orders.
A detailed description of these methods for twopoint and threepoint 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 twopoint functions, either the linear correlation (lines without arrows) or the linear response (lines with an arrow that shows the flow of time). The threeleg 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 twocomponent 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 threepoint correlation C_{3} = ⟨ ψψψ ⟩ , which contains the density and velocity threepoint functions, as well as their cross correlations, but in this article we only consider the density threepoint 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 Npoint 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 pathintegral 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 threeleg 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 largek tail P_{L}(k) ∝ k^{n} 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 lowestorder contribution, or “treeorder” term, is given by diagram (a), which yields the wellknown result (39)The explicit expressions of the 1loop order diagrams (b), ..., (i), which are on the order of (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 steepestdescent” 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 steepestdescent” resummation scheme, the diagrams for the bispectrum up to 1loop order are those shown in Fig. 3. The difference with Fig. 2 is that we now have doubleline twopoint functions. They correspond to the nonlinear twopoint correlation and response functions, as given by the same method up to the same 1loop order. The single lines are the linear twopoint functions, as in Fig. 2.
Fig. 3 Diagrams obtained at 1loop order for the threepoint correlation C_{3} by the “direct steepestdescent” resummation scheme. The double lines are the nonlinear twopoint functions C and R, while the single lines are the linear twopoint functions C_{L} and R_{L} as in Fig. 2. The black dots are again the threeleg vertex K_{s}. 
Within this resummation scheme, nonlinear twopoint functions at “1loop” order actually contain an infinite number of “bubble” diagrams, as shown in Figs. 8 and 9 of Valageas (2008). There, the label “1loop” does not mean that we only include 1loop diagrams, as in the standard perturbation theory of Fig. 2, but that the diagrammatic expansion is only complete up to 1loop (i.e., although we include some diagrams at all orders, we miss some contributions at 2loop and higher orders). Then, substituting the expression of the nonlinear twopoint 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 twopoint functions are replaced by their linear counterparts.
At the order , which corresponds to the 1loop 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 , and only differ by the number of higher order terms that are included.
Fig. 4 Diagrams obtained at 1loop order for the threepoint correlation C_{3} within a “mixed” case. We keep the resummed diagram (a) of Fig. 3 but we replace diagrams (f), (g), (h), and (i) of Fig. 3 by their leading contributions, given by diagrams (f), (g), (h), and (i) of Fig. 2. 
Thus, for the perturbative (3halo) 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_{1},k_{2},k_{3}) contains no free parameter and is exact, up to 1loop 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 treeorder bispectrum (39) for the 3halo 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 Nbody simulations in a ΛCDM universe, consistent with the fiveyear observation of the WMAP satellite (Komatsu et al. 2009), which are described in Valageas & Nishimichi (2011). We analyze the four main simulations (N = 2048^{3}), together with supplementary smaller simulations (N = 1024^{3},512^{3} and 256^{3}) for convergence tests.
We measure the bispectrum from a fast Fourier transformation of the density field obtained by the cloudincells interpolation of the Nbody 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_{1}, k_{2}, and k_{3}, 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), (40)where the factor s_{123} 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_{1},k_{2},k_{3}).
We do not correct for the shotnoise 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 interparticle 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), for the power spectrum and the bispectrum, respectively. Note that the shot noise is included in P(k_{1}), P(k_{2}), and P(k_{3}) of Eq. (42). Instead of subtracting these contributions from the measured spectra, we assess the shotnoise 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 (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 shotnoise 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_{1}, k_{2}, k_{3})
4.2.1. Equilateral triangles
Fig. 5 Bispectrum B_{eq}(k) = B(k,k,k) for equilateral configurations, at redshifts z = 0.35,1, and 3. The symbols are the results from the numerical simulations. We show the bispectrum obtained at tree order (black solid line), standard 1loop order (blue dotdashed line), using the full model (red solid line), as well as the 2halo (cyan dotted line) and 1halo (cyan dashed line) contributions. The magenta dotdashed line labeled “Pan” is the phenomenological model of Pan et al. (2007). 
We show in Fig. 5 our results for the bispectrum, from linear to highly nonlinear scales, for equilateral configurations. Here we take the perturbative 3halo contribution equal to the standard 1loop 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 2halo and 1halo 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 twopoint 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.
Fig. 6 Bispectrum B(k_{1},k_{2},k_{3}) for isosceles configurations, k_{1} = k_{2} = k and k_{3} = k/α, at redshifts z = 0.35,1, and 3. We consider the cases α = 2 (upper row) and α = 4 (lower row). The symbols are the same as in Fig. 5. 
We clearly see in Fig. 5 the decay on large scales of the 1halo and 2halo contributions, in agreement with Eqs. (26) and (33). As explained in Sect. 2, this is due to the new counterterms 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 1loop perturbative contribution significantly extends the domain of validity of the 3halo perturbative term, as compared with the treelevel contribution.
4.2.2. Isosceles triangles
Fig. 7 Bispectrum B(k_{1},k_{2},k_{3}) for isosceles configurations, k_{1} = k_{2}, at redshifts z = 0.35,1, and 3 (from top to bottom in each panel). We show our results as a function of k_{1} = k_{2} at fixed k_{3} (upper row), and as a function of k_{3} at fixed k_{1} = k_{2} (lower row). The symbols are the same as in Fig. 5. 
We show in Fig. 6 our results for isosceles triangles as a function of scale k and for the fixed shapes defined by k_{1} = k_{2} = k and k_{3} = k/α, with constant ratios α. We consider the cases α = 2 and α = 4. As expected, as α increases (i.e., the triangles get “squeezed”), the 2halo contribution becomes more important compared with the 1halo contribution, because the scale 1/k_{3} 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 2halo contribution always remains subdominant because the perturbative (i.e., “3halo”) contribution is still dominant when the 1halo term becomes larger than the 2halo term. For more “squeezed” shapes, an intermediate regime would appear, dominated by the 2halo 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_{1},k_{2},k_{3}} , for isosceles configurations k_{1} = k_{2}.
First, we consider in the upper row the dependence on k_{1} = k_{2} = k at fixed k_{3}, so that the two equallength sides run from k_{3}/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_{1} = k_{2}. The value of the fixed side k_{3} 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_{3}/2 < k_{1} = k_{2} < 10k_{3}). In the left panel, we are dominated by the 3halo perturbative contribution, and since k_{3} remains small, for high values of k_{1} = k_{2} the 2halo term is larger than the 1halo term. In the middle panel, we are already dominated by the 1halo term (except at low k at z = 3 where the 3halo perturbative contribution is larger), and the the 2halo term is subdominant. In the right panel the bispectrum is completely dominated by the 1halo contribution.
Fig. 8 Reduced bispectrum, Q_{eq}(k) = B_{eq}(k)/ [3P(k)^{2}] , for equilateral configurations, at redshifts z = 0.35,1, and 3. The symbols are the same as in Fig. 5. The vertical arrow in the upper right part shows the wavenumber beyond which the simulation shot noise is greater than 10%. 
Second, we consider in the lower row the dependence on k_{3} at fixed value of the common sides k_{1} = k_{2} = k, so that the third side runs from 0 (squeezed limit) to 2k_{1} (flat triangle that reduces to a line). Again, the bispectrum decreases as the dependent wavenumber, here k_{3}, grows. However, the dependence is much weaker than the one found in the upper row. In the left panel, we are dominated by the linearorder perturbative contribution, so that all curves agree with each other. The too high values of the numerical points at low k_{3} 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 3halo and 1halo regimes. The 2halo term is never dominant on these scales. In the right panel we are mostly dominated by the 1halo 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 (44)Indeed, this ratio goes to a constant on large scales, as seen from the treeorder 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 steepestdescent resummation for the perturbative 2halo 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 shotnoise 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 highk 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 largescale plateau associated with the treelevel 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 realspace coefficients , reach a constant or keep growing on small scales, as the density variance 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 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^{1} 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.
Fig. 9 Ratio B(k_{1},k_{2},k_{3})/ [P_{S}(k_{1})P_{S}(k_{2}) + 2cyc.] , where P_{S}(k) is the fit to the power spectrum from Smith et al. (2003). The red solid line is our model, as in Figs. 5 to 8, the black dotted line is the approximation B_{3H} = B_{tree}, the green dotdashed line corresponds to no counterterms in the 1halo and 2halo contributions (), and the blue dashed line makes both approximations. Upper row: equilateral configurations, at redshifts z = 0.35,1, and 3, as in Fig. 5. Lower row: various isosceles configurations at redshift z = 0.35, as in left panels of Figs. 6 and 7. 
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 2halo and 3halo terms, or introducing a highmass 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 selfconsistent 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 1loop perturbation theory instead of the treelevel contribution (39) for the 3halo term; and ii) the counterterms in Eqs. (24) and (31) in the 1halo and 2halo 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 treelevel contribution (39) or the 1halo contribution with a negligible counterterm.
First, we can see in Fig. 9 that discarding the contribution of 1loop 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 1loop 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 1halo term is dominant. At z = 3, the 1loop contribution also extends up to k ~ 0.4 h Mpc^{1} the good agreement with simulations, while using the treelevel 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 1halo 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 2halo and 1halo contributions (while keeping the 1loop contribution in the 3halo 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 2halo and 1halo contributions spoil the good agreement with simulations.
Making simultaneously both approximations, discarding 1loop diagrams and halo counterterms, as in usual implementations of the halo model, also yields significantly worse results. Although both effects partly compensate, neglecting 1loop 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 2halo and 1halo contributions), our approach gives a better accuracy on weakly nonlinear scales. This shows the importance of including both 1loop contributions and the 2halo and 1halo 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
Fig. 10 Ratio B(k_{1},k_{2},k_{3})/ [P_{S}(k_{1})P_{S}(k_{2}) + 2cyc.] , as in Fig. 9. We show the results obtained by using for the 3halo term each of the three 1loop perturbative expansions given in Figs. 2−4. In each case, we show both the 3halo contribution alone (which decays at high k, upper line in front of the labels) and the full result (i.e., the sum of the 3halo, 2halo, and 1halo terms, lower line in front of the labels). We also plot the treelevel perturbative result (black solid line). Upper row: equilateral configurations, at redshifts z = 0.35,1, and 3, as in Fig. 5. Lower row: various isosceles configurations at redshift z = 0.35, as in left panels of Figs. 6 and 7. 
We investigate in this section the dependence of our results on the perturbative scheme used for the 3halo 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 3halo contribution either the standard 1loop result of Fig. 2, the full steepestdescent resummation of Fig. 3, or the mixed case of Fig. 4. For each choice we show both the 3halo term itself (blue curves that decrease at high k) and the full prediction, where we sum the 3halo, 2halo, and 1halo contributions. (The 2halo and 1halo 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 3halo contribution, yields a slightly better agreement with numerical simulations around 1 h Mpc^{1}. However, this is not conclusive, especially since the 1halo term is already nonnegligible on these scales and it is not sufficient to significantly improve the match to simulations. The complete steepestdescent resummation of Fig. 3 yields a 3halo 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 1halo 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 1loop 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 1loop perturbation theory for the perturbative 3halo 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 1loop steepestdescent 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 1loop perturbation theory yields a 2halo 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 1halo term).
We do not have this smallscale problem for the bispectrum, as seen in Figs. 5, 6, and 8. Indeed, the 1loop 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 1halo term in this regime. This is best seen in Fig. 8, which shows that the 1loop 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)^{2}. Therefore, contrary to the power spectrum, it is now possible to use the 1loop 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 3halo 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 loworder lowk behavior and the highk 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 wellunderstood “sweeping effect” associated with the random transport of density structures by largescale velocity flows (Valageas 2007b; Bernardeau & Valageas 2008, 2010). Another advantage of this approach is that its extensions to highorder quantities (Bernardeau et al. 2008), to nonGaussian 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 “largeN 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 selfenergy 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 manybody 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 lowerorder correlations. This method has already been used to study the effect of massive neutrinos (Lesgourgues et al. 2009) and of primordial nonGaussianities (Bartolo et al. 2010b). An advantage of this approach is that it does not involve propagators, but only the usual manybody correlations also encountered in standard perturbation theory, and it is written in terms of singletime 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 wellsuited to the approach described in this paper and presents the advantage that it provides direct extensions to redshiftspace 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
Fig. 11 Reduced bispectrum (upper panel, as in Fig. 8) and the bispectrum scaled by 3P_{S}(k)^{2} (lower panel, as in the upper row of Fig. 10) obtained for modified halo properties on equilateral triangles at z = 0.35. We plot our fiducial model (red solid line, as in previous figures), the impact of a high mass cutoff (M < 10^{15} h^{1} M_{⊙}, blue dashed line; M < 5 × 10^{14} h^{1} M_{⊙}, green dotdashed line), and the effect of defining halos by the lower nonlinear density threshold δ = 50 (black dotted line). 
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.
Fig. 12 Reduced bispectrum (upper panel) and the bispectrum scaled by 3P_{S}(k)^{2} (lower panel) obtained for modified halo properties on equilateral triangles, as in Fig. 11 but at z = 3. 
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^{15} or M < 5 × 10^{14} 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^{14} 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 3halo term includes 1loop perturbative contributions, which are still important on these scales as shown in Figs. 8 and 10, and that the 1halo and 2halo contributions contain the counterterms , so that the relative importance of the 2halo and 1halo 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 highmass 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^{15} h^{1} M_{⊙}, so that the cutoff at M < 5 × 10^{14} 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.
Fig. 13 Reduced bispectrum, Q_{eq}(k) = B_{eq}(k)/ [3P(k)^{2}] , for equilateral configurations, at redshifts z = 0.35,1, and 3. The points are the results from the numerical simulations and the solid line is our “direct” model, as in Fig. 8. The black dashed line is the lower tangent to the “direct” curve Q_{eq}(k), on the transition region. It defines the two contact points of abscissa k_{−} and k_{+}, shown by the two small vertical lines, between which we introduce the two modifications “tang.” (green dotdashed line) and “flat” (blue dashed line), explained in the main text. 
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^{13} and M < 10^{13} 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 highmass 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 highorder perturbative terms that have not been included in the perturbative expansions used here, which are only complete up to oneloop 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 halomodel 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 welldefined 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 (45)The modified power “tang.” is obtained from the “direct” prediction described in Valageas & Nishimichi (2011) by drawing in the (log k,log Δ^{2}) plane the upper tangent line on the interval [k_{−},k_{+}] that runs through the left point (log k_{−},log Δ^{2}(k_{−})). For the cases shown in Fig. 14 the right contact point has an abscissa (by construction we have ), and the power spectrum is only modified on the interval . As seen in Fig. 14, in this fashion we correct most of the artificial “dip” shown by our model without modifying the largescale and smallscale regimes where the “direct” predictions are satisfactory^{2}.
Fig. 14 Power per logarithmic interval of k, as defined in Eq. (45), at redshifts z = 0.35, 1, and 3 (from top to bottom). We show the “direct” model (red solid line), studied in Valageas & Nishimichi (2011), and the geometrical modification “tang.” (green dotdashed line), given by the upper tangent that runs through the point of abscissa k_{−}. The small dashed vertical line shows the location of the contact point of this tangent line with the “direct” curve. 
Fig. 15 Equilateral bispectrum at redshifts z = 0.35, 1, and 3 (from top to bottom). We show the “direct” model (red solid line), already shown in Fig. 5, and the geometrical modification “tang.” (green dotdashed line) given by the upper tangent that runs through the point of abscissa k_{−}. The small dashed vertical line shows the location of the contact point of this tangent line with the “direct” curve. 
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 , 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 , with fixed values that would mark the transition regime (such as and ). 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 redshiftindependent 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 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 . 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 curve over the interval [k_{−},k_{+}] , starting from the right boundary k_{+}, and imposing a monotonic decrease as k decreases. Thus, and are identical over most of [k_{−},k_{+}] (and on all larger and smaller scales), except under the remaining “bump” of , where is constant and equal to the local minimum of 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 welldefined 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 realspace twopoint 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
Fig. 16 Accuracy of our models and our numerical simulations at redshifts z = 0.35, 1, and 3 for the power spectrum. The red solid line “direct” is the relative difference (46) between the simulations and our model as described in Valageas & Nishimichi (2011). The green dashed line “tang.” corresponds to the geometrical modification shown in Fig. 14. The black dashed line is the error level of the simulations, including both the relative statistical error (which grows at low k) and the relative shotnoise (which grows at high k). 
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, (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 Δ^{2}(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.
Fig. 17 Accuracy of our model and our numerical simulations at redshifts z = 0.35, 1, and 3 for the bispectrum on equilateral configurations. We plot our “direct” model (red solid line) and its geometrical modification “tang.” (green dotdashed line), shown in Fig. 15. The black dashed line shows the statistical and shotnoise errors of the simulations. 
Next, we show in Fig. 17 the relative difference between our model and the numerical simulations for the equilateral bispectrum, (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 1loop 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 , 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 3halo term and the nonperturbative 2halo and 1halo terms is the main source of error because of the intrinsic limitations of a description in terms of relaxed spherical halos.
Fig. 18 Accuracy of our model and our numerical simulations at redshifts z = 0.35, 1, and 3 for the reduced bispectrum on equilateral configurations. We plot our “direct” model (red solid line) and its geometrical modifications “tang.” (green dotdashed line) and “flat” (blue dashed line), shown in Fig. 13. The black dashed line shows the statistical and shotnoise errors of the simulations. The vertical arrow in the upper right part shows the wavenumber beyond which the simulation shot noise is greater than 10%. 
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 and 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 3halo, 2halo, and 1halo contributions, which we relate to perturbative and nonperturbative terms. This explains how one can build a model that agrees with perturbation theory at all orders because the 1halo and 2halo terms are nonperturbative corrections that vanish at all orders over P_{L}. Moreover, we explained how new counterterms appear in the 1halo and 2halo contributions. This ensures that these contributions vanish at low k, as required by physical constraints on the power generated by smallscale 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 1loop order, which can be used for the perturbative 3halo 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 1loop 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 steepestdescent 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 1loop 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 1loop 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 highk cutoff. This problem does not occur for the bispectrum because the standard 1loop contribution is now negligible on small scales compared with the 1halo 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 highk 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 Δ^{2}(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 Δ^{2}(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 realspace twopoint 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 1loop 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 nonperturbative 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 nonGaussian 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 redshiftspace distortions (Smith et al. 2008). However, we leave these tasks to future studies.
The halo model can be made to recover the stableclustering ansatz predictions if the mass function scales at low mass as n(M)dM ∝ dM/M^{2} (Valageas 1999; Ma & Fry 2000b). This unrealistic formal limit, where the apparent amount of matter per unit volume is infinite, corresponds to a multiple counting of “halos”, which contain an infinite hierarchy of substructures that are also counted in the mass function, in agreement with a fractal model.
Although we show in Fig. 14 this geometrical construction for the logarithmic power Δ^{2}(k), applying the same construction to the power P(k), that is, in the (log k,log P) plane, yields the same results (since the abscissa of the contact point with the upper tangent line is the same).
Acknowledgments
We thank S. Colombi for useful discussions about the folding procedure for the bispectrum measurements. T.N. is supported by a GrantinAid 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
 Anselmi, S., Matarrese, S., & Pietroni, M. 2011, JCAP, 06, 015 [NASA ADS] [Google Scholar]
 Bartolo, N., Almeida, J. P. B., Matarrese, S., Pietroni, M., & Riotto, A. 2010a, JCAP, 3, 11 [NASA ADS] [CrossRef] [Google Scholar]
 Bartolo, N., Beltrán Almeida, J. P., Matarrese, S., Pietroni, M., & Riotto, A. 2010b, J. Cosmology Astropart. Phys., 3, 11 [Google Scholar]
 Bernardeau, F., & Valageas, P. 2008, Phys. Rev. D, 78, 083503 [NASA ADS] [CrossRef] [Google Scholar]
 Bernardeau, F., & Valageas, P. 2010, Phys. Rev. D, 81, 043516 [NASA ADS] [CrossRef] [Google Scholar]
 Bernardeau, F., Colombi, S., Gaztañaga, E., & Scoccimarro, R. 2002a, Phys. Rep., 367, 1 [NASA ADS] [CrossRef] [EDP Sciences] [MathSciNet] [Google Scholar]
 Bernardeau, F., Mellier, Y., & van Waerbeke, L. 2002b, A&A, 389, L28 [Google Scholar]
 Bernardeau, F., Crocce, M., & Scoccimarro, R. 2008, Phys. Rev. D, 78, 103521 [NASA ADS] [CrossRef] [Google Scholar]
 Bernardeau, F., Crocce, M., & Sefusatti, E. 2010, Phys. Rev. D, 82, 083507 [NASA ADS] [CrossRef] [Google Scholar]
 Cole, S., Percival, W. J., Peacock, J. A., et al. 2005, MNRAS, 362, 505 [NASA ADS] [CrossRef] [Google Scholar]
 Colombi, S., Jaffe, A., Novikov, D., & Pichon, C. 2009, MNRAS, 393, 511 [NASA ADS] [CrossRef] [Google Scholar]
 Cooray, A., & Sheth, R. 2002, Phys. Rep., 372, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Crocce, M., & Scoccimarro, R. 2006a, Phys. Rev. D, 73, 063520 [NASA ADS] [CrossRef] [Google Scholar]
 Crocce, M., & Scoccimarro, R. 2006b, Phys. Rev. D, 73, 063519 [NASA ADS] [CrossRef] [Google Scholar]
 Crocce, M., & Scoccimarro, R. 2008, Phys. Rev. D, 77, 023533 [NASA ADS] [CrossRef] [Google Scholar]
 Dolag, K., Bartelmann, M., Perrotta, F., et al. 2004, A&A, 416, 853 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Duffy, A. R., Schaye, J., Kay, S. T., & Vecchia, C. D. 2008, MNRAS, 390, L64 [NASA ADS] [CrossRef] [Google Scholar]
 Eisenstein, D. J., Hu, W., & Tegmark, M. 1998, ApJ, 504, L57 [NASA ADS] [CrossRef] [Google Scholar]
 Eisenstein, D. J., Zehavi, I., Hogg, D. W., et al. 2005, ApJ, 633, 560 [NASA ADS] [CrossRef] [Google Scholar]
 Fosalba, P., Pan, J., & Szapudi, I. 2005, ApJ, 632, 29 [NASA ADS] [CrossRef] [Google Scholar]
 Frieman, J. A., & Gaztanaga, E. 1994, ApJ, 425, 392 [NASA ADS] [CrossRef] [Google Scholar]
 Frieman, J. A., & Gaztanaga, E. 1999, ApJ, 521, L83 [NASA ADS] [CrossRef] [Google Scholar]
 Giocoli, C., Bartelmann, M., Sheth, R. K., & Cacciato, M. 2010, MNRAS, 408, 300 [NASA ADS] [CrossRef] [Google Scholar]
 Goroff, M. H., Grinstein, B., Rey, S.J., & Wise, M. B. 1986, ApJ, 311, 6 [NASA ADS] [CrossRef] [Google Scholar]
 Guillet, T., Teyssier, R., & Colombi, S. 2010, MNRAS, 405, 525 [NASA ADS] [Google Scholar]
 Guo, H., & Jing, Y. P. 2009, ApJ, 698, 479 [NASA ADS] [CrossRef] [Google Scholar]
 Hamilton, A. J. S., Kumar, P., Lu, E., & Matthews, A. 1991, ApJ, 374, L1 [NASA ADS] [CrossRef] [Google Scholar]
 Jain, B., & Bertschinger, E. 1996, ApJ, 456, 43 [NASA ADS] [CrossRef] [Google Scholar]
 Jing, Y. P., & Suto, Y. 2002, ApJ, 574, 538 [NASA ADS] [CrossRef] [Google Scholar]
 Kayo, I., Suto, Y., Nichol, R. C., et al. 2004, Publ. Astron. Soc. Japan, 56, 415 [Google Scholar]
 Komatsu, E., Dunkley, J., Nolta, M. R., et al. 2009, ApJS, 180, 330 [NASA ADS] [CrossRef] [Google Scholar]
 Lesgourgues, J., Matarrese, S., Pietroni, M., & Riotto, A. 2009, J. Cosmology Astropart. Phys., 6, 17 [Google Scholar]
 Ma, C.P., & Fry, J. N. 2000b, ApJ, 543, 503 [NASA ADS] [CrossRef] [Google Scholar]
 Massey, R., Rhodes, J., Leauthaud, A., et al. 2007, ApJS, 172, 239 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Matarrese, S., & Pietroni, M. 2007, JCAP, 6, 26 [NASA ADS] [Google Scholar]
 Matsubara, T. 2008a, Phys. Rev. D, 78, 083519 [NASA ADS] [CrossRef] [Google Scholar]
 Matsubara, T. 2008b, Phys. Rev. D, 77, 063530 [NASA ADS] [CrossRef] [Google Scholar]
 McClelland, J., & Silk, J. 1977, ApJ, 217, 331 [NASA ADS] [CrossRef] [Google Scholar]
 Munshi, D., Valageas, P., van Waerbeke, L., & Heavens, A. 2008, Phys. Rep., 462, 67 [NASA ADS] [CrossRef] [Google Scholar]
 Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 490, 493 [NASA ADS] [CrossRef] [Google Scholar]
 Nishimichi, T., Kayo, I., Hikage, C., et al. 2007, Publ. Astron. Soc. Japan, 59, 93 [Google Scholar]
 Nishimichi, T., Shirata, A., Taruya, A., et al. 2009, Publ. Astron. Soc. Japan, 61, 321 [Google Scholar]
 Nishimichi, T., Taruya, A., Koyama, K., & Sabiu, C. 2010, JCAP, 7, 2 [NASA ADS] [Google Scholar]
 Pan, J., Coles, P., & Szapudi, I. 2007, MNRAS, 382, 1460 [NASA ADS] [Google Scholar]
 Peebles, P. J. E. 1974, A&A, 32, 391 [NASA ADS] [Google Scholar]
 Peebles, P. J. E. 1980, The large scale structure of the universe (Princeton: Princeton university press) [Google Scholar]
 Peebles, P. J. E. 1982, ApJ, 263, L1 [NASA ADS] [CrossRef] [Google Scholar]
 Pietroni, M. 2008, JCAP, 10, 36 [Google Scholar]
 Scherrer, R. J., & Bertschinger, E. 1991, ApJ, 381, 349 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Schneider, P., & Bartelmann, M. 1995, MNRAS, 273, 475 [NASA ADS] [Google Scholar]
 Scoccimarro, R. 1997, ApJ, 487, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Scoccimarro, R., & Couchman, H. M. P. 2001, MNRAS, 325, 1312 [NASA ADS] [CrossRef] [Google Scholar]
 Scoccimarro, R., & Frieman, J. 1996, ApJS, 105, 37 [NASA ADS] [CrossRef] [Google Scholar]
 Scoccimarro, R., Colombi, S., Fry, J. N., et al. 1998, ApJ, 496, 586 [NASA ADS] [CrossRef] [Google Scholar]
 Scoccimarro, R., Sheth, R. K., Hui, L., & Jain, B. 2001, ApJ, 546, 20 [NASA ADS] [CrossRef] [Google Scholar]
 Sefusatti, E., & Komatsu, E. 2007, Phys. Rev. D, 76, 083004 [NASA ADS] [CrossRef] [Google Scholar]
 Sefusatti, E., Crocce, M., Pueblas, S., & Scoccimarro, R. 2007, Phys. Rev. D, 74, 023522 [Google Scholar]
 Sefusatti, E., Crocce, M., & Desjacques, V. 2010, MNRAS, 406, 1014 [NASA ADS] [Google Scholar]
 Sheth, R. K., & Jain, B. 2003, MNRAS, 345, 529 [NASA ADS] [CrossRef] [Google Scholar]
 Smith, R. E., Peacock, J. A., Jenkins, A., et al. 2003, MNRAS, 341, 1311 [NASA ADS] [CrossRef] [Google Scholar]
 Smith, R. E., Watts, P. I. R., & Sheth, R. K. 2006, MNRAS, 365, 214 [NASA ADS] [CrossRef] [Google Scholar]
 Smith, R. E., Sheth, R. K., & Scoccimarro, R. 2008, Phys. Rev. D, 78, 023523 [NASA ADS] [CrossRef] [Google Scholar]
 Takada, M., & Jain, B. 2003, MNRAS, 340, 580 [NASA ADS] [CrossRef] [Google Scholar]
 Taruya, A., & Hiramatsu, T. 2008, ApJ, 674, 617 [NASA ADS] [CrossRef] [Google Scholar]
 Taruya, A., Nishimichi, T., Saito, S., & Hiramatsu, T. 2009, Phys. Rev. D, 80, 123503 [NASA ADS] [CrossRef] [Google Scholar]
 Taylor, A. N., & Hamilton, A. J. S. 1996, MNRAS, 282, 767 [NASA ADS] [Google Scholar]
 Tegmark, M., Eisenstein, D. J., Strauss, M. A., et al. 2006, Phys. Rev. D, 74, 123507 [Google Scholar]
 Valageas, P. 1999, A&A, 347, 757 [NASA ADS] [Google Scholar]
 Valageas, P. 2002, A&A, 382, 477 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Valageas, P. 2004, A&A, 421, 23 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Valageas, P. 2007a, A&A, 465, 725 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Valageas, P. 2007b, A&A, 476, 31 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Valageas, P. 2008, A&A, 484, 79 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Valageas, P. 2009, A&A, 508, 93 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Valageas, P. 2011, A&A, 526, A67 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Valageas, P., & Nishimichi, T. 2011, A&A, 527, A87 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Verde, L., Heavens, A. F., Percival, W. J., et al. 2002, MNRAS, 335, 432 [Google Scholar]
 Wang, Y., Yang, X., Mo, H. J., van den Bosch, F. C., & Chu, Y. 2004, MNRAS, 353, 287 [NASA ADS] [CrossRef] [Google Scholar]
Appendix A: Behavior of the bispectrum for one low wavenumber in perturbation theory
We show in this appendix how the largescale behavior (32) also arises at all orders of perturbation theory. As is well known, within the standard perturbation theory the density contrast field is written as a perturbative expansion over powers of the linear density contrast , (A.1)with (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_{1},...,k_{n}} . Then, the bispectrum can be obtained within perturbation theory from the threepoint correlation (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 (A.4)where we used a shorthand notation for the diagram shown in Fig. A.1, and the Dirac factors contain sums over the wavenumbers . Of course, the three Dirac factors can be combined to factorize out a Dirac prefactor δ_{D}(k_{1} + k_{2} + k_{3}), 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_{12} = n_{13} = n_{22} = 0 is excluded).
Fig. A.1 A perturbative contribution to the threepoint connected correlation , as in Eq. (A.4). The big red circles corresponds to the three nonlinear density contrasts , , and . The small black dots represent the linear fields and the Gaussian average amounts to connect them by the linear correlation , shown by the blue solid lines. There are n_{11}, n_{22}, and n_{33} internal lines within the circles , , and . There are n_{12}, n_{13}, and n_{23} lines connecting the three circles. Here we have {n_{11},n_{22},n_{33}} = {2,2,3} and {n_{12},n_{13},n_{23}} = {3,2,4} . This is part of , which involves the kernels , , and . 
We are interested in the limit k_{1} → 0, at fixed k_{2} and k_{3}. A wellknown property of the kernels , which arises from momentum conservation (Peebles 1974; Goroff et al. 1986), is that as k = k_{1} + ... + 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 behaves as if 2n_{11} + n_{12} + n_{13} ≥ 2. Therefore, in the limit k_{1} → 0 the dominant contributions come from the diagrams with 2n_{11} + n_{12} + n_{13} = 1, since F_{1} ≡ 1 and for CDM cosmologies with n_{s} ≃ 1 at low k_{1}. This corresponds to either {n_{11} = 0,n_{12} = 1,n_{13} = 0} or {n_{11} = 0,n_{12} = 0,n_{13} = 1} because we only take the connected part of (A.4). Thus, for k_{1} → 0 the dominant contribution to the perturbative bispectrum reads as (A.5)where we note by “1 perm” the symmetric contribution with respect to {k_{2} ↔ k_{3}} . There is a subtlety in Eq. (A.5), because the limit at low k_{1} of the kernel contains divergent terms of the form (Goroff et al. 1986; Scoccimarro & Frieman 1996). However, taking the limit k_{1} → 0 elsewhere, we obtain k_{3} = −k_{2} so that the symmetric term with respect to {k_{2} ↔ k_{3}} actually corresponds to a change of sign of all wavenumbers. This cancels the divergences of the form so that the limit k_{1} → 0 is finite. The fact that such infrared divergences cancel out for equaltime 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 largescale limit (A.5) through the nonlinear corrections to and . In contrast, as we have seen above, only the linear term 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 highorder perturbative terms and nonperturbative corrections.
It is interesting to note that the “renormalization” of the prefactor to the P(k_{1}) 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^{2}P_{L}(k) at low k. This can also be understood from the fact that for the power spectrum there are no partial largescale limits. On the other hand, if all wavenumbers go to zero, as k_{1} ~ k_{2} ~ k_{3} ~ k with k → 0, higher order contributions to the bispectrum scale at least as k^{2}P_{L}(k)^{2} while the “treeorder” result scales as P_{L}(k)^{2}, so that in the simultaneous largescale limit we also recover the lowest order result.
Appendix B: Twopoint correlation function
Fig. B.1 Realspace twopoint correlation function ξ(x), at redshifts z = 0.35, 1, and 3 (from top to bottom). We show the “direct” model (red solid line) and its geometrical modification “tang” (green dotdashed line). They are the Fourier transforms of the corresponding power spectra shown in Fig. 14. 
Fig. B.2 Realspace twopoint correlation function ξ(x) at redshifts z = 0.35, 1, and 3 (from top to bottom) as in Fig. B.1 but on larger scales. In addition to the “direct” model (red solid line) and its geometrical modification “tang” (green dotdashed line), we also show the linear twopoint correlation (black solid line). 
Fig. B.3 Accuracy of our models and our numerical simulations at redshifts z = 0.35, 1, and 3 for the realspace twopoint correlation function. As in Figs. B.1 and B.2, the red solid line “direct” is the relative difference (46) between the simulations and our model as described in Valageas & Nishimichi (2011), while the green dashed line “tang.” corresponds to the geometrical modification shown in Fig. 14. The black dashed line shows the relative statistical error (which grows at large x) and shotnoise error (which grows at small x) of the simulations. 
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 realspace twopoint correlation function, which is given by Therefore, we compare in this appendix the twopoint 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 twopoint 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 Δ^{2}(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 shown in Fig. 14 was only continuous, with a discontinuous derivative at k_{−}, the correlation ξ_{tang}(x) has a welldefined and finite derivative over all x > 0.
Contrary to the power spectrum, the modification “tang” to the twopoint 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 1loop 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
Fig. 1 Probability F_{3H} that the triplets associated with the equilateral bispectrum B_{eq}(k) belong to three different halos, from Eq. (35), at redshifts z = 0.35,1, and 3. 

In the text 
Fig. 2 Diagrams associated with the standard perturbation theory as obtained from a pathintegral formalism up to 1loop order for the threepoint correlation C_{3}. Although they are different from those obtained by the standard approach, their sum at each order is identical to the sum of the standard diagrams of the same order over P_{L}(k). The lines are the linear twopoint correlation C_{L} (blue solid line) and the linear response R_{L} (red solid line with an arrow). The black dots are the threeleg vertex K_{s} that enters the quadratic term of the equation of motion. The numbers are the multiplicity factors of each diagram. The treeorder diagram a) gives Eq. (39), for the density bispectrum, while the 1loop diagrams b), ..., i), give the contribution of order . 

In the text 
Fig. 3 Diagrams obtained at 1loop order for the threepoint correlation C_{3} by the “direct steepestdescent” resummation scheme. The double lines are the nonlinear twopoint functions C and R, while the single lines are the linear twopoint functions C_{L} and R_{L} as in Fig. 2. The black dots are again the threeleg vertex K_{s}. 

In the text 
Fig. 4 Diagrams obtained at 1loop order for the threepoint correlation C_{3} within a “mixed” case. We keep the resummed diagram (a) of Fig. 3 but we replace diagrams (f), (g), (h), and (i) of Fig. 3 by their leading contributions, given by diagrams (f), (g), (h), and (i) of Fig. 2. 

In the text 
Fig. 5 Bispectrum B_{eq}(k) = B(k,k,k) for equilateral configurations, at redshifts z = 0.35,1, and 3. The symbols are the results from the numerical simulations. We show the bispectrum obtained at tree order (black solid line), standard 1loop order (blue dotdashed line), using the full model (red solid line), as well as the 2halo (cyan dotted line) and 1halo (cyan dashed line) contributions. The magenta dotdashed line labeled “Pan” is the phenomenological model of Pan et al. (2007). 

In the text 
Fig. 6 Bispectrum B(k_{1},k_{2},k_{3}) for isosceles configurations, k_{1} = k_{2} = k and k_{3} = k/α, at redshifts z = 0.35,1, and 3. We consider the cases α = 2 (upper row) and α = 4 (lower row). The symbols are the same as in Fig. 5. 

In the text 
Fig. 7 Bispectrum B(k_{1},k_{2},k_{3}) for isosceles configurations, k_{1} = k_{2}, at redshifts z = 0.35,1, and 3 (from top to bottom in each panel). We show our results as a function of k_{1} = k_{2} at fixed k_{3} (upper row), and as a function of k_{3} at fixed k_{1} = k_{2} (lower row). The symbols are the same as in Fig. 5. 

In the text 
Fig. 8 Reduced bispectrum, Q_{eq}(k) = B_{eq}(k)/ [3P(k)^{2}] , for equilateral configurations, at redshifts z = 0.35,1, and 3. The symbols are the same as in Fig. 5. The vertical arrow in the upper right part shows the wavenumber beyond which the simulation shot noise is greater than 10%. 

In the text 
Fig. 9 Ratio B(k_{1},k_{2},k_{3})/ [P_{S}(k_{1})P_{S}(k_{2}) + 2cyc.] , where P_{S}(k) is the fit to the power spectrum from Smith et al. (2003). The red solid line is our model, as in Figs. 5 to 8, the black dotted line is the approximation B_{3H} = B_{tree}, the green dotdashed line corresponds to no counterterms in the 1halo and 2halo contributions (), and the blue dashed line makes both approximations. Upper row: equilateral configurations, at redshifts z = 0.35,1, and 3, as in Fig. 5. Lower row: various isosceles configurations at redshift z = 0.35, as in left panels of Figs. 6 and 7. 

In the text 
Fig. 10 Ratio B(k_{1},k_{2},k_{3})/ [P_{S}(k_{1})P_{S}(k_{2}) + 2cyc.] , as in Fig. 9. We show the results obtained by using for the 3halo term each of the three 1loop perturbative expansions given in Figs. 2−4. In each case, we show both the 3halo contribution alone (which decays at high k, upper line in front of the labels) and the full result (i.e., the sum of the 3halo, 2halo, and 1halo terms, lower line in front of the labels). We also plot the treelevel perturbative result (black solid line). Upper row: equilateral configurations, at redshifts z = 0.35,1, and 3, as in Fig. 5. Lower row: various isosceles configurations at redshift z = 0.35, as in left panels of Figs. 6 and 7. 

In the text 
Fig. 11 Reduced bispectrum (upper panel, as in Fig. 8) and the bispectrum scaled by 3P_{S}(k)^{2} (lower panel, as in the upper row of Fig. 10) obtained for modified halo properties on equilateral triangles at z = 0.35. We plot our fiducial model (red solid line, as in previous figures), the impact of a high mass cutoff (M < 10^{15} h^{1} M_{⊙}, blue dashed line; M < 5 × 10^{14} h^{1} M_{⊙}, green dotdashed line), and the effect of defining halos by the lower nonlinear density threshold δ = 50 (black dotted line). 

In the text 
Fig. 12 Reduced bispectrum (upper panel) and the bispectrum scaled by 3P_{S}(k)^{2} (lower panel) obtained for modified halo properties on equilateral triangles, as in Fig. 11 but at z = 3. 

In the text 
Fig. 13 Reduced bispectrum, Q_{eq}(k) = B_{eq}(k)/ [3P(k)^{2}] , for equilateral configurations, at redshifts z = 0.35,1, and 3. The points are the results from the numerical simulations and the solid line is our “direct” model, as in Fig. 8. The black dashed line is the lower tangent to the “direct” curve Q_{eq}(k), on the transition region. It defines the two contact points of abscissa k_{−} and k_{+}, shown by the two small vertical lines, between which we introduce the two modifications “tang.” (green dotdashed line) and “flat” (blue dashed line), explained in the main text. 

In the text 
Fig. 14 Power per logarithmic interval of k, as defined in Eq. (45), at redshifts z = 0.35, 1, and 3 (from top to bottom). We show the “direct” model (red solid line), studied in Valageas & Nishimichi (2011), and the geometrical modification “tang.” (green dotdashed line), given by the upper tangent that runs through the point of abscissa k_{−}. The small dashed vertical line shows the location of the contact point of this tangent line with the “direct” curve. 

In the text 
Fig. 15 Equilateral bispectrum at redshifts z = 0.35, 1, and 3 (from top to bottom). We show the “direct” model (red solid line), already shown in Fig. 5, and the geometrical modification “tang.” (green dotdashed line) given by the upper tangent that runs through the point of abscissa k_{−}. The small dashed vertical line shows the location of the contact point of this tangent line with the “direct” curve. 

In the text 
Fig. 16 Accuracy of our models and our numerical simulations at redshifts z = 0.35, 1, and 3 for the power spectrum. The red solid line “direct” is the relative difference (46) between the simulations and our model as described in Valageas & Nishimichi (2011). The green dashed line “tang.” corresponds to the geometrical modification shown in Fig. 14. The black dashed line is the error level of the simulations, including both the relative statistical error (which grows at low k) and the relative shotnoise (which grows at high k). 

In the text 
Fig. 17 Accuracy of our model and our numerical simulations at redshifts z = 0.35, 1, and 3 for the bispectrum on equilateral configurations. We plot our “direct” model (red solid line) and its geometrical modification “tang.” (green dotdashed line), shown in Fig. 15. The black dashed line shows the statistical and shotnoise errors of the simulations. 

In the text 
Fig. 18 Accuracy of our model and our numerical simulations at redshifts z = 0.35, 1, and 3 for the reduced bispectrum on equilateral configurations. We plot our “direct” model (red solid line) and its geometrical modifications “tang.” (green dotdashed line) and “flat” (blue dashed line), shown in Fig. 13. The black dashed line shows the statistical and shotnoise errors of the simulations. The vertical arrow in the upper right part shows the wavenumber beyond which the simulation shot noise is greater than 10%. 

In the text 
Fig. A.1 A perturbative contribution to the threepoint connected correlation , as in Eq. (A.4). The big red circles corresponds to the three nonlinear density contrasts , , and . The small black dots represent the linear fields and the Gaussian average amounts to connect them by the linear correlation , shown by the blue solid lines. There are n_{11}, n_{22}, and n_{33} internal lines within the circles , , and . There are n_{12}, n_{13}, and n_{23} lines connecting the three circles. Here we have {n_{11},n_{22},n_{33}} = {2,2,3} and {n_{12},n_{13},n_{23}} = {3,2,4} . This is part of , which involves the kernels , , and . 

In the text 
Fig. B.1 Realspace twopoint correlation function ξ(x), at redshifts z = 0.35, 1, and 3 (from top to bottom). We show the “direct” model (red solid line) and its geometrical modification “tang” (green dotdashed line). They are the Fourier transforms of the corresponding power spectra shown in Fig. 14. 

In the text 
Fig. B.2 Realspace twopoint correlation function ξ(x) at redshifts z = 0.35, 1, and 3 (from top to bottom) as in Fig. B.1 but on larger scales. In addition to the “direct” model (red solid line) and its geometrical modification “tang” (green dotdashed line), we also show the linear twopoint correlation (black solid line). 

In the text 
Fig. B.3 Accuracy of our models and our numerical simulations at redshifts z = 0.35, 1, and 3 for the realspace twopoint correlation function. As in Figs. B.1 and B.2, the red solid line “direct” is the relative difference (46) between the simulations and our model as described in Valageas & Nishimichi (2011), while the green dashed line “tang.” corresponds to the geometrical modification shown in Fig. 14. The black dashed line shows the relative statistical error (which grows at large x) and shotnoise error (which grows at small x) of the simulations. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.