An improved fitting formula for the dark matter bispectrum

In this paper we present an improved fitting formula for the dark matter bispectrum motivated by the previous phenomenological approach of Scoccimarro&Couchman (2001). We use a set of LCDM simulations to calibrate the fitting parameters in the k-range of 0.03 h/Mpc<k<0.4 h/Mpc and in the redshift range of 0<z<1.5. This new proposed fit describes well the BAO-features although it was not designed to. The deviation between the simulations output and our analytic prediction is typically less than 5% and in the worst case is never above 10%. We envision that this new analytic fitting formula will be very useful in providing reliable predictions for the non-linear dark matter bispectrum for LCDM models.


Introduction
The dark matter and galaxy power spectrum have been widely used to study the growth of structure, to constrain cosmological parameters and galaxy bias models. These tools have proved very successful and have contributed to crystallize the current LCDM model e.g., [2] and refs therein. With ongoing and forthcoming galaxy surveys, like BOSS 1 and EUCLID 2 , the signal-to-noise of the data will increase and the uncertainties around this model will be reduced. Higher precision data will allow the use of not only the two-point correlation function, but also of higher-order statistics, in order to constrain and improve our theories and models. The bispectrum (the three-point correlation function in Fourier space) is naturally the next statistic to consider [3,4]. Using both the power spectrum and bispectrum we can improve our knowledge of the growth of structure and galaxy biasing [5][6][7][8][9][10][11][12][13], constrain possible departures from Gaussianity in the initial conditions of the matter density field [14][15][16][17][18] as well as constrain departures from GR e.g., [19,20].
From a theoretical point of view, perturbation theory and subsequent improvements such as renormalized perturbation theory [21], resummed perturbation theory or time-RG flow [22] is a physically well-motivated approach to study these statistical moments. Tree-level perturbation theory has demonstrated to describe well the behavior of the power spectrum and bispectrum at large scales. However, non-trivial computations are needed to obtain predictions at non-linear scales: the one-loop correction and beyond, for the power spectrum and bispectrum. For the power spectrum, however, other phenomenological approaches have been demonstrated to work better for a wide range of redshifts and different cosmologies e.g., [23][24][25]. For the bispectrum there are also simple phenomenological models that predict its behavior at non-linear scales [1], but they fail to accurately reproduce the BAO-features [26] and are only precise at the 20%-30% level. Therefore better analytical models are needed to describe the bispectrum at these non-linear scales.
In this paper we improve the phenomenological description presented by [1] (hereafter SC) more than 10 years ago. Using a set of modern simulations we fit the free parameters of our proposed analytic formula. Thus, we obtain an improved description for the bispectrum in the LCDM model scenario (including baryonic acoustic oscillations) in a range of 0.03 h/Mpc ≤ k ≤ 0.4 h/Mpc and for different redshifts, 0 ≤ z ≤ 1.5. This paper is organized as follows: in §2 we begin with a description of the density field statistics and different analytic approaches to the dark matter bispectrum. In §3 we describe the simulations we use to fit the parameters. In §4 we present our results, compare them with previous fitting formulae and with 1-loop corrections and discuss the differences. We finally conclude in §5. In Appendix A, we give details of how the bispectrum and its errors are computed from simulations. In Appendix B we test how our formula works for other non-standard LCDM models. In Appendix C we present a short description of 1-loop correction in Eulerian perturbation theory.

Power spectrum & bispectrum
The power spectrum P (k), the Fourier transform of the two-point correlation function, is one of the simplest statistics of interest one can extract from the dark matter overdensity field δ(k), The Dirac delta function ensures that the bispectrum is defined only for k-vector configurations that form closed triangles: i k i = 0. Note that δ(k 1 )δ(k 2 )δ(k 3 ) is in general a complex number, however once the average is taken, the imaginary part goes to zero. It is convenient to define the reduced bispectrum Q 123 ≡ Q(k 1 , k 2 , k 3 ) as, which takes away part of the dependence on scale and cosmology 3 . The reduced bispectrum is useful when comparing different models, since it only has a weak dependence on cosmology and one can thus break degeneracies between cosmological parameters in order to isolate the effects of gravity. The bispectrum for Gaussian initial conditions is zero and remains zero in linear theory, i.e. as long as the k-modes evolve independently 4 . However, when non-linearities start to play an important role, mode coupling is no longer negligible and the bispectrum becomes non-zero. Thus, by measuring the bispectrum one can extract information about how non-linear processes influence the evolution of dark matter clustering.

Analytic approaches in the literature
In order to understand the observational data, we need accurate theoretical predictions for B(k 1 , k 2 , k 3 ). A physically well-motivated analytic theory for doing this, is perturbation theory (PT hereafter) (see [27] for a review) or subsequent improvement such as renormalized PT, resummed PT etc.
In an Einstein de-Sitter Universe (hereafter EdS Universe) and at second order (tree-level) in Eulerian perturbation theory, the bispectrum is given by [28], where B 123 = B(k 1 , k 2 , k 3 ), P L i = P L (k i ) is the linear power spectrum, and the symmetrized twopoint kernel F s 2 is given by where θ ij is the angle between the vectors k i and k j . This formula is the second order perturbation theory contribution to the bispectrum which is the leading order contribution. On quasi-linear scales, this expression is a very good prediction but fails in the moderate non-linear regime. The dependence on cosmology of the two-point kernel F s 2 is very weak and hence the cosmology dependence of the bispectrum is almost completely contained in P L i . Because of this, in this work, we use the kernel of Eq. 2.5 even though we are dealing with the LCDM model.
One can improve the tree-level PT prediction by going one step further and including one-loop corrections. However, at this point the computation of the bispectrum becomes cumbersome. For an initially Gaussian δ-field this yields four additional terms to the tree-level contribution (see Appendix C for details).
An alternative way of reaching these non-linear scales, without using the one-loop correction, and to even push beyond the one-loop regime of validity, is with phenomenologically motivated models. Phenomenological formulae can give simpler expressions in the non-linear regime and accurate predictions for the bispectrum. However, their physical motivation is limited and they usually have free parameters that need to be calibrated using N-body simulations.
SC proposed a fitting formula based on the structure of the formula of Eq. 2.4. It consists in replacing the linear power spectrum by the non-linear one in Eq. 2.4 and the EdS two-point symmetrized kernel by where the functions a(n, k), b(n, k) and c(n, k) are chosen to interpolate between the tree-level results and the hyper-extended perturbation theory regime (HEPT) [29], a(n, k) = 1 + σ a6 c(n, k) = 1 + 4.5a 4 /[1.5 + (n + 3) 4 ](qa 5 ) n+3 1 + (qa 5 ) n+3.5 .
Here n is the slope of the linear power spectrum at k, and q ≡ k/k nl , where k nl is the scale where non-linearities start to be important and is defined as, a i are free parameters that must be fitted using data from simulations. In particular, SC propose the values, a 1 = 0.25, a 2 = 3.5, a 3 = 2, a 4 = 1, a 5 = 2, a 6 = −0.2 .
The function Q 3 (n) is given by (2.10) With all these changes, the SC approach reads, B 123 = 2F eff 2 (k 1 , k 2 )P 1 P 2 + cyc. perm., (2.11) where P i is the non-linear power spectrum at k i . On large scales, where the functions a, b and c → 1 we recover the tree-level PT formula for the bispectrum. On the other hand, on small scales a 2 → (7/10)Q 3 and b and c → 0 and we obtain Q 123 → Q 3 (n), which is the prediction of HEPT.
Another approach based on phenomenological formulae, is the one presented by [26]. The main idea is to rescale the linear formula of the bispectrum, by using some scale transformation in k. This way, the tree-level formula can easily be extended up to non-linear scales using the ansatz . This approach has by definition the drawback that it does not preserve the BAO-features of the bispectrum. In particular, the rescaling of k produces a spurious rescaling of the peaks and troughs of the BAO wiggles that do not match with the data, producing higher deviations than the SC approach. Because of that, we do not consider this approach in this paper.

Our analytic formula
Our approach in this paper is inspired by the SC approach. It consists of not only refitting the a i parameters from Eq. 2.7 but of also modifying their expression to make it more suitable for current precision N-body data and consider the redshift range of 0 ≤ z ≤ 1.5. In order to do that, we use simulations with more particles, larger box sizes, and more realizations (and thus higher precision, better statistics and better error-control) with respect to previous works; we also consider snapshots at different redshifts. In order to improve the fitting precision, we also add 3 more parameters to the original model. The modified functionsã(n, k),b(n, k),c(n, k) then read, Note that one recovers the original SC formulae in the limit of a 7 → 1 and a 8 , a 9 → 0. The original SC formula was not designed to reproduce the BAO features. Applying this formula to a power spectrum with BAOs produces unphysical oscillations. These oscillations are much larger than those observed in simulations (see black dashed line in the right panel of Fig. 1). These oscillations are caused by the oscillatory behavior of the slope parameter n. One solution to this problem is to "dewiggle" the linear power spectrum [30]. However here we want to preserve the BAO oscillations. We propose to smooth the oscillatory behavior of the parameter n by means of splines, as is shown in the blue solid line of the left panel of Fig. 1. This provides an improved fit to the BAO-features, as it is shown by the blue solid line in the right panel of Fig. 1. In order to smooth out n we calculate its spline by taking a number of points n(k), where the points are chosen to be in the middle of the amplitude of each wiggle, such that when the points are connected a smooth line would pass through them. These points are used in the spline routine, and their second order derivatives are calculated for each point k. This output is then fed into the spline routine, which returns a smoothed value of n for each value of k.

Simulations
The simulations in this paper consist of two different sets, namely A and B. Each simulation is characterized by the box size, L b , the number of particles, N p , and the number of independent   Table 1 for details on the simulations). Black dashed line is SC prediction without any spline in n(k) and blue solid line with the spline in n(k).  Table 1. Simulations details for simulations A and B. L b is the box size, Np is the number of particles, Nr is the number of independent realizations. A quarter of the Nyquist frequency, kN /4, is the upper threshold to which we trust the simulation results at the percent level. The force resolution is specified by the softening parameter and the Particle Mesh (PM) grid size. The short-range force accuracy is determined by α through the cell-opening criterion M l 2 > α|a old |r 4 , where M is the mass inside the cell, l its side length, a old the total acceleration of the particle in the previous time step, and r the distance between the particle and the cell. The remaining parameters set the time stepping in Gadget-2: the maximum global time step in the logarithm of the scale factor, max(∆ log a), and the parameter η in the individual time step criterion ∆a = 2η /|a|, where a is the acceleration of the individual particle.
runs, N r . Details about the two simulations are given in Table 1. As a rule of thumb, a maximum threshold in k for trusting the simulation data is set by a quarter of the Nyquist frequency, defined as . At this scale it has been observed that the power spectrum starts to deviate at the 1%-level with respect to higher resolution simulations [31]. We confirmed this result using our two sets of simulations. For all the plots and results shown in this paper this limit in k is never exceeded.
Both A and B simulations consist in a flat LCDM cosmology with cosmological parameters consistent with observational data. The cosmology used is Ω Λ = 0.73, Ω m = 0.27 h = 0.7, Ω b h 2 = 0.023, n s = 0.95 and σ 8 (z = 0) = 0.7913. The initial conditions were generated at z = 19 and z = 49 for simulations A and B respectively, by displacing the particles according to the second-order Lagrangian PT from their initial grid points. The initial power spectrum of the density fluctuations  Table 2. Best-fit parameters (according to Eq. 2.12) derived by combining data from simulations A and B, using different triangle configurations θ12/π = 0.1, 0.2, . . . , 0.9 and k2/k1 = 1.0, 1.5, 2.0, 2.5 and at different redshifts z = 0, 0.5, 1.0, 1.5. was computed by CAMB [32]. Taking only the gravitational interaction into account, the simulation was performed with GADGET-2 code [33].
As we estimate the error of the bispectrum from its dispersion among different realizations (see Eq. A.5 in Appendix A for details), and given that we only have 3 simulations of type B, we divide each of these 3 boxes into 8 sub-boxes. Each of these 24 sub-boxes is then treated as if it were an independent realization with smaller box size, L b = 937.5 h/Mpc, where each of these sub-boxes contains about 512 3 particles. The measurements of the bispectrum from sub-boxes suffer from two issues: a) the measurements are not completely independent and more importantly b) the sub-boxes are affected by modes larger than sub-box size. As a consequence of this, a new source of non-Gaussian errors arises for the power spectrum and bispectrum estimation, called beat-coupling effect [34][35][36]. However, by using the mean density measured in each sub-box instead of the global mean density for the normalization of the density contrast, δ ≡ ρ/ρ − 1, this effect gets strongly suppressed [37]. Hence, we expect that on overlapping scales the bispectrum errors estimated from simulation B to be slightly larger than those from A. This is shown in Appendix A.
In order to obtain the dark matter field from particles we discretize each box of simulation A and each sub-box of simulation B using 512 3 grid cells. Thus the size of the grid cells is 4.68 Mpc/h in A and 1.83 Mpc/h in B. We assign the particles to the cells using the count-in-cells prescription.
More details about the estimation of the bispectrum from simulations and the error bars computation are given respectively in Eq. A.4 and A.5 in Appendix A.

Results
In order to find the best-fit parameters from Eq. 2.12, namely a i , we minimize using a set of triangle configurations: k 2 /k 1 = 1.0, 1.5, 2.0, 2.5 and θ 12 /π = 0.1, 0.2, . . . , 0.9, at different redshifts: z = 0, 0.5, 1.0, 1.5 5 . The algorithm used for the minimization is amoeba [38]. In our analysis we neglect that the errors of the data points are correlated. However, since our errors are small (typically less than 5%) we expect that error correlations do not play an important role for this method: the dominating source of the error of our fitting formula given in Eq. 2.12 comes from the imperfection of the functional form of the fitting formula and not from the uncertainties in the simulation data. We have checked that the result converges from different starting points. The resulting best-fit values are shown in Table 2.
We have also checked that this is a very good fit not only for the reduced bispectrum Q but also for the bispectrum B.
In Fig. 2 and 3 we show the results of our fit and also the predictions of two other models for different triangles configurations for z = 0 (Fig. 2) and for z = 1 (Fig. 3). These two models are 1-loop Eulerian PT (see Appendix C) and the SC method (Eq. 2.11) + smoothed-n. In each plot we show the reduced bispectrum Q vs. k 1 for: N-body data (black circles for simulations A and black squares for simulations B), 1-loop correction (red solid line corresponding to data of simulations A and red dashed line to data of simulations B), SC formula (green solid line for A and green dashed line for B) and our model (blue solid and dashed line for A and B respectively). In the bottom part of each panel we show the deviation of these models with respect to the N-body data: red symbols depict the deviation of the 1-loop prediction with respect to the data, green symbols are the deviation of the SC formula and blue symbols are the deviation of our model. Circles are the deviation with respect to simulation A and squares with respect to simulation B. The error bars show the error of the bispectrum measured from the simulations (see Appendix A for details). Each panel shows a different triangle configuration: from left to right k 2 /k 1 = 1.0, 1.5, 2.0 and from top to bottom θ 12 /π = 0.2, 0.4, 0.6, 0.8. In order to avoid sample variance effects, we only use data points with k i > 0.03 h/Mpc for simulation A and k i > 0.09 h/Mpc for simulation B.
For some triangle configurations, we observe a mismatch between the data points in the overlapping region of simulations A and B . This is due to the fact that the data points do not exactly correspond to the same triangle configurations. The central value for k i is the same for both simulations in each panel. However, since A and B have different fundamental frequencies k f ≡ 2π/L b and we take the bin width to be ∆k = 3k f , the k-space over which we average the bispectrum for each bin is different (see Appendix A for further explanation). Hence, when one computes the effectivek i using Eq. A.8 one obtains that, especially for elongated triangle configurations, the simulations do not represent the same triangle configuration. Since the theoretical predictions (both 1-loop, SC and our model) are computed fromk i there is also a mismatching between solid and dashed lines for the same reason. This effect is also noted in [30].
We have checked that this mismatch is not due to the fact that the two simulations have different resolution but to the fact that different triangle configurations are sampled. When the triangle configurations and scales coincide, we do not observe any mismatch. Note also that in the cases where a mismatch appears -due to different configurations being sampled-, it is also present (and quantitatively similar) in the the 1-loop theory prediction (see for example top left panel of Fig. 3).
At z = 0, the deviation of our model from the data is typically less than 5% and always less than 10% for k < 0.4 h/Mpc, whereas for SC the deviation reaches values of up to 20% and 1-loop clearly breaks down for k 0.1 h/Mpc. We also observe that for triangles close to equilateral, both the SC approach and our work present maximum differences to the N-body data. This might have to do with the fact that for equilateral (or close to equilateral) triangles the 3 sides enter the nonlinear regime at the same time, and thus, non-linearities play a stronger role than for other triangle configurations, where each side enters the non-linear regime at different redshifts. Because of this, for other configurations the differences are smaller and remain within 5% deviation.
For z = 1 all models work better because non-linearities do not play such an important role. However, even in this regime, our model works better than the other two. Again, for triangle configurations close to equilateral, our model reaches its maximum deviation of about 10%. At z = 1, all other triangle configurations typically have errors within 5%.
As a cross-check, in Appendix B we compare our model with another set of simulations of nonstandard LCDM model. In particular for equilateral configurations, our formula reaches deviations up to 10%. However for the scalene configurations k 2 = 2k 1 the deviations are only of order 3%. In the cases studied here our model improves significantly the SC fitting formula.

Conclusions
In this paper we propose a new simple formula to compute the dark matter bispectrum in the moderate non-linear regime (k < 0.4 h/Mpc) and for redshifts z ≤ 1.5. Our method is inspired by the approach presented in [1], but includes a modification of the original formulae, namely Eq. 2.12, and a prescription to better describe the BAO oscillations. The cosmology dependence of the reduced bispectrum is known to be very weak, and that of the bispectrum is almost completely contained in the power spectrum. Given that the cosmological model today is well constrained by observations we have considered a single cosmology here.
Using LCDM simulations we fit the free parameters of our model. We end up with a simple analytic formula that is able to predict accurately the bispectrum for a LCDM Universe including the effects of BAO. Our main results are summarized by Eq. 2.11 where the kernel is given by Eq.   The main conclusions of our work are listed below.
1. Our method is able to predict the dark matter bispectrum for a wide range of triangle configurations up to k = 0.4 h/Mpc and for a redshift range 0 ≤ z ≤ 1.5. In particular, for the reduced bispectrum, our fitting formula agrees within 5% with N-body data for most of the triangle configurations and always within 10% for the worst cases. This presents a considerable improvement over previous phenomenological approaches and over the prediction of Eulerian perturbation theory.
2. The equilateral and quasi-equilateral configurations are the ones for which our model deviates most strongly from N-body data. We interpret this as being due to the fact that when the 3 sides of the triangle are similar, non-linearities start to play a role at the same time, and thus, the effect on the bispectrum is stronger than when the non-linearities enter at different times, i.e. for elongated triangles. Other methods, like the one described by SC show the same behavior.
3. We have checked that our model also works well for non-standard LCDM cosmologies (see Appendix B). In particular, we have checked that for k 2 = 2k 1 the deviation between N-body data and our model is never higher than 3% and for equilateral triangles reaches 10%. Also in these non-standard LCDM cases studied here, our model works better than the SC fitting formula.
We envision that this new analytic fitting formula will be very useful in providing a reliable prediction for the non-linear dark matter bispectrum for LCDM models. In particular, simple analytic predictions with high accuracy will be needed for the data analysis in the forthcoming era of precision data.

A Appendix: bispectrum estimator & error bars
Here we present details on the computation of the bispectrum and its error bars from N-body simulations. Moreover, we compare our error estimates with the Gaussian analytic predictions and discuss the differences.
We start by defining the estimator for the bispectrum as, f is the volume of the fundamental cell, k f . The integration is defined over the bin k i − ∆k i /2 < q i < k i + ∆k i /2. In this paper we always take ∆k = 3k f . V B is the six-dimensional volume of triangles defined by the triangle sizes k 1 , k 2 and k 3 with uncertainty ∆k. Its value can be approximated by which is good enough for not too small values of k i . The variance associated to this estimator depends on higher-order correlation functions: up to the 6-point connected correlation function. However, the main contribution to the variance is given by the power spectrum. Assuming that the fields are Gaussian, the variance associated to estimator presented above is [39], where the symmetry factor is s B = 6, 2, 1 for equilateral, isosceles or scalene configurations. The factor (2π) 3 comes from our definition of the power spectrum and bispectrum in Eq. 2.1 and 2.2.
On the other hand, the discretized version of this estimator used in this paper is, B where N tri is the number of random triangle configurations used to compute the bispectrum; and j runs over these triangle configurations. For this work we use a number of random triangles that increases with k in the same way as the number of fundamental triangles: ∼Ṽ B /V 2 f . It reaches up to N tri ∼ 10 9 for scales k ∼ 0.4 h/Mpc. We have checked that increasing the number of random triangles beyond this value has no effect neither in the value of the bispectrum nor in its error. The index d in the δ field stands for a discrete and dimensionless quantity. Therefore the quantity Re δ d j (k 1 )δ d j (k 2 )δ d j (k 3 ) needs to be rescaled with the factor L 6 b to make B matching with the definition of the bispectrum in Eq. 2.2. We compute the variance of this estimator B by the sample variance derived from the N r realizations, where B i is the bispectrum derived from the realization i and B is the mean over all realizations, The error on the mean B is then simply given by When comparing the measured N-body bispectrum with theoretical models and also when comparing Eq. A.3 with Eq. A.5, it is important to take into account the effect on the finite size of the triangle bins: each configuration is defined in terms of the sides of the triangle k i ± ∆k/2. In this case, we are assuming ∆k = 3k f and a large number of fundamental triangles fit into this bin. For certain configurations, it turns out that we have more triangles with k larger than the central value, k i . Because of that, one must correct the sides of the triangles bỹ where i = 1, 2, 3 for each dimension and the sum is taken over all random triangle generated in the bin. This correction is extremely important at large scales and for very squeezed triangles, and less important for equilateral configuration. In Fig 4 we present a comparison of the error estimation from theoretical models (Eq. A.3) and simulations (Eq. A.5). In the left panel we show the error of the bispectrum associated to a volume of 1 single box for ∆k = 3k f : using the theoretical model and the simulations for the case of k 2 /k 1 = 1 and θ 12 = 0.6π triangle configurations. The blue line shows the theoretical model prediction for the error of the bispectrum of one single realization using the non-linear power spectra from simulations A and B (Eq. A.3); whereas the red squares and green circles show the dispersion among the runs of simulations A and B respectively (Eq. A.5). In the right panel the ratio between the errors according to the simulations and the Gaussian prediction is plotted for simulations A (red squares) and B (green circles). The error estimates of simulations A agrees well with the theoretical model. On the other hand, on small scales the error estimates of simulations B is larger than the theoretical model and further increase with decreasing scale. Similar results were found by [40]. These differences are due to the fact that Eq. A.3 neglects any higher-order contributions (because it assumes Gaussianity). However, at small scales this is no longer a good approximation as it is shown in [36]. Furthermore, the errors of simulations B have been estimated by dividing each of the 3 simulation boxes into 8 sub-boxes. This introduces extra non-Gaussian terms [36] that are not taken into account in Eq. A.3.  Here we test how our model works with different LCDM simulations to those we have used to fit the a i parameters. In particular, we test our model with LCDM simulations with a f (R)-like power spectrum. For a full description of the simulations and the f (R) gravity we refer the reader to [19]. These simulations were run with ENZO code and have slightly different cosmology than the ones used in the rest of this paper: Ω Λ = 0.76, Ω m = 0.24, Ω b = 0.04181, h = 0.73. They consist of 6 realizations that contain 256 3 particles in a box of 400 Mpc/h per side. The one-quarter Nyquist frequency is k N /4 = 0.5 h/Mpc. In Fig. 5 the reduced bispectrum Q is shown: in the right panel as a function of the angle between k 1 and k 2 , namely θ 12 , for k 2 = 2k 1 ; in the left panel as a function of k 1 for equilateral configuration, both for z = 0. All panels correspond to a LCDM model with no BAOs whose initial conditions make their power spectrum look like a f (R)-like one. In order to do that a running index has been adopted in the initial conditions (see Table 1 in [19] for details). Panels correspond to LCDM simulations that match with f (R) models whose |f R0 | parameter 6 is: 10 −4 (top panels), 10 −5 (middle panels) and 10 −6 (bottom panels ).
Black points are data from simulations, green line is the SC prediction and blue line the prediction of our model. In the right panel we only compare data points for 0.4 < θ 12 /π < 0.9 and in the left panel 0.1 h/Mpc < k < 0.5 h/Mpc in order to ensure that all the k i are smaller than a quarter of the Nyquist frequency.
Considering the right panels of Fig. 5 (k 2 = 2k 1 = 0.4 h/Mpc), we see that our model describes the data within about 3%. In particular for small scales (θ 12 < 0.7π), both SC and our model agree with the simulations data well; however at large scales (θ 12 > 0.7π), our model fits the data points better. In the left panel of Fig. 5 (equilateral configuration), we see that our model shows deviations up to 10%. As for the standard LCDM model, the equilateral configuration is the one with the largest deviations. However even in this case, our formula behaves better than the SC model, especially at small scales.
Therefore we conclude that our formula is general enough to be applied also to some non-standard LCDM models and in particular works better than SC at small scales.
C Appendix: one-loop correction terms for the power spectrum and bispectrum Here we present a short description of the equations used to compute the one-loop correction in Eulerian perturbation theory for the bispectrum shown in Fig. 2 and 3. For a detailed description of PT see [27,39]. 6 Figure 5. In the left panels: Q vs. k1 for equilateral configuration. In the right panels Q vs θ12/π for k2 = 2k1 = 0.4 h/Mpc. From top to bottom we show different non-standard LCDM models. All the panels correspond to LCDM models with f (R)-like power spectrum from [19]. Top panels correspond to matching to |fR0| = 10 −4 ; middle panels to |fR0| = 10 −5 ; bottom panels to |fR0| = 10 −6 . The sub-panel shows the corresponding ratio between simulations and different theoretical models: SC (green points) and our work (blue points). In the right sub-panels, dashed lines mark 2.5% deviation, whereas in the left panels mark the 5% deviation. The bispectrum and its error are estimated by Eq. A.6 and A.7 (see Appendix A).
For Gaussian initial conditions, the one-loop term consist of two terms 7 , P 22 = 2 (2π) 3 d 3 q F s 2 2 (q, k − q)P L (q)P L (|k − q|) (C.2) P 13 = 6 (2π) 3 P L (k) d 3 qF s 3 (k, q, −q)P L (q) (C.3) P 22 accounts for the mode coupling between waves with wave-vectors k − q and q, whereas P 13 can be interpreted as the one-loop correction to the linear propagator.