Nonlinear effects of dark energy clustering beyond the acoustic scales
Abstract
We extend the resummation method of Anselmi & Pietroni (2012) to compute the total density power spectrum in models of quintessence characterized by a vanishing speed of sound. For standard CDM cosmologies, this resummation scheme allows predictions with an accuracy at the few percent level beyond the range of scales where acoustic oscillations are present, therefore comparable to other, common numerical tools. In addition, our theoretical approach indicates an approximate but valuable and simple relation between the power spectra for standard quintessence models and models where scalar field perturbations appear at all scales. This, in turn, provides an educated guess for the prediction of nonlinear growth in models with generic speed of sound, particularly valuable since no numerical results are yet available.
a,b]Stefano Anselmi, b,c]Diana López Nacir, d,b]and Emiliano Sefusatti Prepared for submission to JCAP
Nonlinear effects of dark energy clustering beyond the acoustic scales

Department of Physics/CERCA/ISO, Case Western Reserve University, Cleveland, OH 441067079 – USA

The Abdus Salam International Center for Theoretical Physics, strada costiera 11, I34151 Trieste – Italy

Departamento de Física and IFIBA, FCEyN UBA, Facultad de Ciencias Exactas y Naturales, Ciudad Universitaria, Pabellón I, 1428 Buenos Aires – Argentina

INAF  Osservatorio Astronomico di Brera, via E. Bianchi 46, I23807 Merate (LC) – Italy
Keywords: Cosmology, Large Scale Structure of the Universe
Contents
 1 Introduction
 2 Equations of motion and linear solutions
 3 Nonlinear solutions: an analytical insight
 4 The power spectrum in the mildly nonlinear regime
 5 Mapping smooth to clustering quintessence power spectra
 6 Conclusions
 A Nonlinear propagator and power spectrum in the large or eikonal limit
 B Evolution equations: factorization and interpolation at intermediate scales
 C Testing the AP resummation
1 Introduction
Ever since the discovery of the cosmic acceleration found by the observations of type Ia supernovae in 1998 [1, 2], a tremendous amount of observational and theoretical effort has been put in place in order to understand the cause of this phenomenon.
A plethora of alternative models attribute the observed acceleration to either the presence of a mysterious energy component dubbed dark energy (DE) which dominates at the current epoch, or to a modification of the gravitational laws on cosmological distances (see, e.g. [3] for an overview of proposed theoretical models). While the standard CDM model, where the acceleration is driven by a cosmological constant , provides a good fit to observations, several models cannot, so far, be ruledout. The next generation of largescale structure (LSS) surveys such as EUCLID [4] and DESI^{1}^{1}1http://desi.lbl.gov/ will provide an unprecedented determination of observables such as the weak lensing shear and galaxy power spectra over a very large range of scales, and will be therefore sensitive to the details of structure formation, representing a unique opportunity to discriminate between competing models. However, the improvement in the quality and quantity of observational data represents a challenge to the accuracy of theoretical predictions, particularly for nonstandard models.
A fundamental complication, in this respect, is given by the nonlinear nature of structure formation. Indeed many cosmological probes of matter density correlations will present small statistical errors on scales smaller than , where gravitational instability is responsible for relevant nonlinear corrections to linear theory predictions. The solution to the problem entails two goals: to understand and predict these corrections on one hand and, on the other, to achive a fast way (seconds) to compute numerically such predictions, in order to efficiently sample the theory parameter space in the data analysis. Our current understanding of the nonlinear regime, rather limited even for the standard paradigm, is very poor for non–standard cosmological models. In the following we will present an important step in this direction for the case of models with perturbations in the quintessence field.
Quintessence is one of the most popular models of DE, where the acceleration of the Universe is attributed to a scalar field with negative pressure [5]. Its standard formulation involves a minimally coupled canonical scalar field, whose fluctuations propagate at the speed of light . In this case, sound waves maintain the quintessence homogeneous on scales smaller than the Hubble radius [6]. The quintessence contribution to the total energy density affects the expansion of the Universe and the growth of dark matter perturbations, allowing observational constraints on the DE equation of state , where and are the pressure and energy density of the scalar field. In this setup, quintessence clustering effects take place only on scales of order , and their detection is therefore severely limited by cosmic variance. However, this may not be the case for more general quintessence models where the scalar field kinetic term is noncanonical^{2}^{2}2 The name “quintessence” is usually reserved for theories where the dark energy is given by a scalar field with a canonical kinetic term, however, following previous works on the subject, we will refer as quintessence to a generic minimally coupled scalar field, with either canonical or noncanonical kinetic term.. In fact, in [7, 8], following the approach of [9], the authors constructed the most general effective field theory describing quintessence perturbations around a given FriedmanRobertsonWalker background spacetime, showing the existence of theoretically consistent models where the fluctuations propagate at a speed . In particular, by a careful analysis of the theoretical consistency of the models, namely, the classical stability of the perturbations and the absence of ghosts, ref. [8] concluded that the region of parameter space with (the so called “phantom regime”) is in general excluded, unless the fluctuations are characterized by a practically vanishing, imaginary speed of sound, e.g. (note that observational limits do not exclude the case [10, 11]). In the latter case, the stability of the perturbations would be guaranteed by the presence of higher derivative operators with negligible effects on cosmological scales. As was emphasized in [7, 8], such a tiny speed of sound should not be considered as a fine tuning, since in the limit one recovers the ghost condensate model of ref. [12], which is invariant under a shift symmetry of the field, and hence a value of can be interpreted as a deformation of the limit where such symmetry is recovered.
What makes these models with of particular interest is that quintessence perturbations can cluster on all observable scales, opening up to a rich phenomenology. In fact, for such models the quintessence field is expected to affect not only the growth history of dark matter through the background evolution but also by actively participating to structure formation, both at the linear and nonlinear level. In the linear regime, the observational consequences of clustering quintessence have been investigated by several authors: see for instance [13, 14, 15, 16, 17, 18, 19] for large scale CMB anisotropies, [20] for galaxy redshift surveys, [21] for neutral hydrogen observations or [22, 23] for the crosscorrelation of the Integrated SachsWolfe effect in the CMB and the LSS. However, the possibility to detect quintessence perturbations effects focusing on observations in the linear regime alone is limited by large degeneracies with cosmological parameters, particularly when quintessence clusters on all relevant scales. In this case, in fact, such effects consist essentially in a change of the fluctuations amplitude. On the other hand the relation between the linear and nonlinear growth of density perturbations can be affected by clustering quintessence in a peculiar way, leading to specific signatures that could in principle allow to distinguish these models from CDM and homogeneous quintessence cosmologies [24, 25, 26].
Robust predictions of the nonlinear dynamics of the structure formation process ultimately require numerical simulations. The case of homogeneous quintessence (i.e., ), as it affects only the background evolution has been investigated by means of Nbody simulations in several articles (see for instance [27, 28] and references therein). The case of clustering quintessence, however, represents a much more difficult problem and numerical results are still lacking. The extension of common fitting formulas such as halofit [29, 30, 31], based in turn on numerical simulations, is also, a priori not possible. Precisely for this reason, the study of analytical approximations in cosmological perturbation theory (PT) to investigate the mildly nonlinear regime is particularly relevant. On one hand it can give insights in the physical processes at play. One of the main results of this work, in fact, consists in showing how fitting functions valid for tested, standard cosmologies can be used in clustering quintessence models. Such extensions have been already considered in the literature without a proper theoretical justification [19, 32, 33]. On the other hand, PT predictions can as well provide useful checks for future numerical results, which are likely to comprise significant approximations in their description of the dynamics of quintessence perturbations.
The nonlinear regime of structure formation in quintessence models with has been studied in previous literature. Several works focussed on the spherical collapse of structures in the presence of quintessence perturbations [34, 35, 36]. In addition, ref. [37] makes use of the spherical collapse model [38] (see [39, 40] for the case of arbitrary sound speed) along with the PressSchechter formalism [41], to provide predictions for the halo mass function. We will consider here, instead, analytical predictions in cosmological perturbation theory for the mildly nonlinear regime of the density power spectrum (see [42] for a review of standard PT).
In the context of standard, CDM cosmologies, perturbative techniques have witnessed a significant progress, motivated by the need for accurate predictions for Baryon Acoustic Oscillations measurements, following the seminal papers of Crocce & Scoccimarro [43, 44, 45]. Several different ÒresummedÓ perturbative schemes and approaches have been proposed in the last few years [46, 47, 48, 49, 50, 51, 52], along with public codes implementing some of the methods [53, 54]. The key improvement over standard PT consisted in reorganizing the perturbative series in a more efficient way by means of a new building block, the nonlinear propagator, a quantity measuring the response of the density and velocity perturbations to their initial conditions. The possibility to analytically compute this new object revolutionized indeed the cosmological perturbation theory approach.
A first extension of standard PT to the specific case of clustering quintessence can be found in [24], where the pressureless perfect fluid equations have been worked out, while in [25, 26] it was shown how to compute the nonlinear power spectrum (PS) by means of the TimeRenormalization Group (TRG) method of [25, 26]. However, the computation of the nonlinear propagator for these cosmologies was still lacking. This is the first achievement of the present work, allowing the extension of a large set of powerful resummed techniques to clustering Dark Energy models.
We compute the nonlinear PS exploiting the resummation scheme proposed by Anselmi & Pietroni [55] (hereafter AP). Such scheme takes into account the relevant diagrams in the perturbative expansion and proposes a new interpolation procedure between small and large scales. What makes AP so appealing is the possibility to extend the predictions, for the first time, to scales beyond the Baryonic Acoustic Oscillation (BAO) range, together with a fast and easy code implementation. It is important to stress that the AP approach, as other perturbative techniques, is intrinsic limited by the emergence of multistreaming effects at small scales, i.e. the generation of velocity dispersion due to orbit crossing. Clearly, the breakingdown of the pressureless perfect fluid (PPF) approximation employed as a starting point for the fluid equations has nothing to do with the perturbative scheme used. The AP prediction agrees with Nbody simulations at the 2% level accuracy on BAO scales while at smaller scales shows discrepancies only of a few percents up to at . As we expect, the level of the agreement is consistent with the findings of [56] where the authors estimated the departures from the PPF assumption by means of high resolution Nbody simulations. Clearly, it would be desirable to extend our results beyond the PPF description. However, so far, even in the standard scenario this is still subject of intense investigations (see [56, 57, 58, 59, 60, 61] for an incomplete list of contributions).
Taking advantage of the AP approach, we will show that additional nonlinear corrections due DE clustering become larger than for scales smaller than BAO for 10% variation from . Thus it allows us to finally test the approximation employed in the TRG predictions of [26] and the assumptions made in several works forecasting the detectability of dark energy models in presence of quintessence perturbations [19, 32, 33]. Furthermore, as a byproduct, we provide a theoretically motivated mapping from the nonlinear matter power spectrum in smooth quintessence models to the nonlinear, total density power spectrum in clustering quintessence models that could possibly serve as a consistency check for numerical results once these will become available. This is not of small significance, given the challenge posed by accurate simulations of structure formation in this kind of models.
This paper is organized as follows. In Section 2 we review the equations of motion for the quintessence–dark matter fluid. In Section 3 we recast the starting fluid equations in a compact form and make explicit the nonlinear behavior of perturbations. In Section 4, in order to compute the nonlinear density power spectrum, we apply the selected resummation scheme to the clustering quintessence fluid. In Section 5 we present our results and a well motivated mapping from the smooth to the clustering quintessence PS. Section 6 is devoted to our conclusions and possible future developments. Some details about the resummation scheme and a comparison with recent numerical simulations for a CDM cosmology are included as Appendices.
2 Equations of motion and linear solutions
The equations of motion of matter and dark energy perturbations for a quintessence field characterized by a vanishing speed of sound, along with their linear solutions, have been studied in several works [24, 25, 26]. In order to keep this work as much selfcontained as possible and to set the notation, in this section we will summarize these results and the relevant aspects of the solutions at the linear level.
We consider a spatially flat FRW spacetime with metric (where is the conformal time given by and the cosmic time), containing only cold dark matter (CDM) and dark energy. For sake of simplicity, we assume a DE equation of state with constant , so that the Hubble rate is given by
(2.1) 
where , with the subscript “0” denoting quantities evaluated today, and ( being the mean value of the matter (quintessence) energy density ().
At the level of perturbations, we consider the perfect fluid approximation in the Newtonian regime to describe both the CDM and DE components, assuming the system to interact only gravitationally. In addition, we restrict ourself to models where the rest frame sound speed vanishes, (for an analysis where is allowed to take values between and see [26]). In that case the two fluids are comoving and, following the notation of [24], the fluid equations will be given by the Euler equation for the common CDM and quintessence velocity field ,
(2.2) 
and by the continuity equations for the density contrast of dark matter and quintessence ,
(2.3)  
(2.4) 
with the conformal Hubble rate and the gravitational potential satisfying the Poisson equation,
(2.5) 
where . In the last equality we have used the convenient definition for the “total” density contrast,
(2.6) 
corresponding to the weighted sum of matter and quintessence perturbations, normalized such that for .
In Fourier space, by combining the previous equations and introducing the function
(2.7) 
it is straightforward to see that the evolution equations for the total density contrast and velocity divergence can be recast as
(2.8)  
(2.9) 
with
(2.10)  
(2.11) 
where we adopt the notation and .
Note that for and (i.e., ) these equations reduce to the usual equations that describe the matter density contrast in the smooth quintessence scenario (with ). If, in addition, we set , we reduce to the case of a cosmology. In other terms, the function captures all the corrections to the equations of motion induced by the clustering of quintessence. It is therefore crucial, for the interpretations of our results, to keep in mind the behavior of this function for different values of . This is shown in Fig. 1 (reproducing Fig. 1 of [24, 25]), where is plotted as a function of redshift for and different values of .
Following [24], the linear solutions for the total density contrast and for the velocity divergence can be written in the separable form
(2.12)  
(2.13) 
where is the linear growth function and the linear growth rate,
(2.14) 
For the smooth case () and constant equation of state parameter , the solution for the linear growing mode is known in terms of hypergeometric functions as [62, 63]
(2.15) 
with
(2.16) 
and where is normalized such that as , recovering the behavior at early time, during matter domination.
In the case of clustering quintessence with , for constant , the equations always admit an integral solution for the growing mode of the total density fluctuations, given by [24]
(2.17) 
It can also be expressed in terms of Hypergeometric functions as
(2.18) 
with defined as in eq. (2.16) above.
The linear growth factor for the matter perturbations alone can be obtained from the solution for the total density, given the relation [24]:
(2.19) 
These different solutions are shown in Fig. 2 as a function of redshift, for and . Top panels plot the ratio between the growth factor for the clustering case to the one for the CDM model, while bottom panels plot the ratio of the growth function for the clustering case to the corresponding solution for smooth quintessence^{3}^{3}3Note that while there are theoretical motivations for considering for , this is not the case for , due to the presence of ghost instabilities. However, we consider also the latter case just for comparison.. Left and right panels show, respectively, the results for the total and matter perturbations.
The upper left panel, in particular, illustrates how different and competing effects are at play in the clustering quintessence models. Since the fraction increases as when decreases, if ( leads to the opposite results) the redshift at which the acceleration becomes relevant (i.e. when ) is larger than for CDM, suppressing more the growth of density perturbations. However, in the case of clustering quintessence with ( ), the presence of quintessence fluctuations leads to a larger (smaller) growth for the total perturbations. This second effect dominates at low redshift for the total density perturbations, while it is quite subdominant for the matter perturbations. In fact, while at better that at for or .
3 Nonlinear solutions: an analytical insight
In this section we show how it is possible to obtain some insights on the nonlinear evolution of the total density field, simply by inspection of the formal solution to the equations of motion. Following the notation of [64], we rewrite the fluid equations using^{4}^{4}4From now on, unless explicitly stated, we will refer to and as the linear growth factors for total and matter perturbations respectively in the clustering DE scenario, i.e., we drop the “” subscript. as the time variable, and introduce the doublet (), defined by
(3.1) 
along with the vertex matrix, () whose only independent, nonvanishing, elements are
(3.2) 
and . Then, the two equations (2.8) and (2.9) can be recast as
(3.3) 
with
(3.4) 
and where repeated indexes are summed over and integration over momenta and is implied. Again, for and eq. (3.3), with eq. (3.4), reduces to the corresponding one for the case of smooth DE, after replacing the functions , and in eq. (3.1) by the appropriate ones.
As for CDM cosmologies, what complicates the implementation of a resummation technique to these equations is the fact that the matrix is time dependent. However, as for the CDM case [44], most of the time we have , and one can see that for the approximation works better than for CDM, while it worsen for . In other words, this means that at linear level almost all of the information about the cosmological parameters is encoded in the linear growth factor which has been rescaled out in eq. (3.1). Therefore, as commonly done for CDM cosmologies, in what follows we consider the approximation
(3.5) 
which is expected to be accurate at better than for , but quickly improving at larger redshift (see, in particular, Appendix A in [53] for a relevant test).
The linear solution can be expressed as
(3.6) 
where is the retarded linear propagator which obeys the equation
(3.7) 
with causal boundary conditions. Explicitly, it is given by [44]:
(3.8) 
with the stepfunction, and
(3.9) 
The initial conditions corresponding respectively to the growing and decaying modes can be selected by considering initial fields proportional to
(3.10) 
At this stage, before introducing statistics and computing the observables, we can extract the first insights about the non–linear behavior of the clustering quintessence fluid. The formal solution of eq. (3.3) reads
(3.11) 
where, assuming growing mode initial conditions and being the initial value of the field, we have . Following [43] we can expand in this fashion
(3.12) 
with
(3.13) 
where . Replacing (3.12) and (3.13) in (3.11) we find the kernel recursion relations
(3.14) 
where the rhs has to be symmetrized under interchange of any two wave vectors. For , . For a CDM cosmology, in the limit where the initial conditions are imposed in the infinite past we recover the well known SPT recursion relation [43]. On the other hand in the clustering quintessence case we can not analytically perform the time integrals. However, when we take into account just the propagator growing mode in eq. (3.14), the times integrals will contribute a factor ^{5}^{5}5The other terms contributing to , given by the propagator decaying mode, will be of the same order.
(3.15) 
Recalling that from eq. (2.19) we have
(3.16) 
where is the growth factor at some initial time . It follows that
(3.17) 
Therefore the order contribution of the total density contrast follows as was already pointed out at second order in [25]. Basically this means that the nonlinear interactions are driven by dark matter, while the effect of DE enters as a multiplicative factor (i.e. ), acting in the same way at both linear and non–linear level. The resummation we will perform confirms this behavior.
4 The power spectrum in the mildly nonlinear regime
We focus now on the nonlinear behavior of a fundamental quantity: the total density power spectrum. In this section we present the nonlinear evolution equations which we then solve numerically to compute the power spectrum. To this end we use of the resummation scheme proposed in [55, 65]. Although a detailed description of the resummation scheme can be found in these references, for the sake of completeness, in Appendices A and B we include a summary of the main relevant aspects of it applied to the clustering quintessence case. The generalization of this approach to the clustering quintessence case is quite straightforward because, in the absence of a sound horizon, there is no additional preferred scale in the system, and moreover, as described in the previous section, in the perfect fluid approximation the velocity field is the same for both DM and DE, allowing us to reduce the number of equations to those of a single fluid. Indeed, given the linear propagator approximation introduced in the previous section, the only difference in the implementation of the diagrammatic technique developed for the CDM case in [44], and of the functional method of [46], is that now each vertex comes multiplied by the time dependent function , as . In what follows, whenever we make use of the diagrammatic representations for specific perturbative contributions, we refer to the definition of the Feynman rules in Fig. 3, where is the linear power spectrum,
(4.1) 
which we define at the initial time in terms of the linear growing mode as , with the initial density power spectrum.
At nonlinear level, it is wellknown that either by exploring the diagrammatic structure of the contributions at arbitrary high order [43] or by applying functional methods [46], it is possible to express the exact evolution equations for the propagator and for the PS, as
(4.2) 
and
(4.3)  
Here the quantity , known as the “self energy”, is a causal kernel () given by the sum of oneparticleirreducible (1PI) diagrams ^{6}^{6}6 A 1PI diagram is one that cannot be separated into two disjoint parts by cutting one of its lines. connecting a dashed line at time with a continuous one at . The kernel is given by the the sum of 1PI diagrams that connect two dashed lines at and with no time ordering. The corresponding 1loop diagrams are shown in Fig. 4.
In particular, the starting point of [55, 65] is the equation for the propagator obtained by deriving eq. (4.2) with respect to
(4.4) 
where
(4.5) 
and the equation obtained deriving eq. (4.3) with respect to at equal times, i.e. after setting ,
(4.6) 
These are nonlocal integrodifferential equations which are very difficult to solve, unless some approximation is considered. The approximation considered in [65] relies on the observation that the equations simplify considerably in the limit of large and small values of , presenting a common structure in both limits. Indeed, it has been shown that in the large (or eikonal) limit, the kernel factorize as (see details in appendices A and B.1)
(4.7) 
where the bold indexes are not summed over, and
(4.8) 
with being the 1loop approximation to the full , given by (see also Fig. 4)
(4.9) 
Ref. [65] also shows that in the opposite limit , where linear theory is recovered, the same factorization applies to a very good approximation^{7}^{7}7More specifically, Ref. [65] shows that the factorization occurs for small values of but after contracting with one , . However, it argues that without contracting with the factorization works for the individual components of the propagators in the limit , and that for finite it is still a very good approximation. . Therefore, a natural way of interpolating between the small and large regimes consists in assuming eq. (4.7) to hold also for intermediate values of , leading to a local equation for the propagator
(4.10) 
In fact, the solution of this equation (which we denote as ) is exact in the large limit, and reduces to the 1loop propagator for . Computing this solution in our case we obtain, in the large limit (see Appendix A)
(4.11) 
with the superscript “eik” standing for the eikonal limit,
(4.12) 
where is defined by eq. (2.19), and
(4.13) 
Therefore, we recover the wellknown Gaussian decay found in [44] for the matter propagator in the standard scenario, which takes into account nonperturbatively the main infrared effects of the longwavelength modes. In the clustering quintessence case, a crucial difference is given by the fact that the derived propagator for the total density field, depends on the matter linear growth factor . This is exactly what we expected from the analysis of the previous section.
In [55], a similar procedure was followed to obtain an approximation for the evolution equation of the power spectrum. Of course, in the physical situation we are considering, the large limit is not a good approximation. In the CDM scenario, multi streaming effects, which are neglected in this framework, are certainly relevant for [56]. Moreover the separation of scales is not so large for . However, as for the propagator, a sensible choice is to use an interpolation between small and large regimes. The equation derived in [55] (see also Appendix B.2) takes the form
(4.14) 
where is given by the weighted sum of the diagrams in Fig. 11, and the explicit expression generalized to our case of clustering DE is given by (see Appendix B.2),
(4.15)  
Here is the 1loop approximation to given by (see the right panel of Fig. 4)
(4.16) 
The second term of eq. (4.15) contributes only for high values and it needs to be switched off for small values of the momenta, we do that by means of a filter function defined as (see Appendix B.2)
(4.17) 
where is taken to be the scale at which the two contribution in eq. (4.15) are equal at (which turns out to be of order ). Note that with this the 1loop result is trivially recovered at small values of up to corrections of .
In the rest of the paper we will consider eq.s (4.10) and (4.14), along with (4.8) and (4.15), and we solve them numerically ^{8}^{8}8 In eq. (4.14) one should in principle use the propagator , obtained after solving eq. (4.10). However, for the CDM case it has been checked that the approximation of replacing by in that contribution gives only sub percent differences [55] and the same happens in our case. Therefore, for the sake of simplicity, we also consider this approximation here. .
Notice that, as a by product of the resummation scheme and of the fact that at for viable models of quintessence, we expect the results for the power spectrum of the rescaled field in clustering quintessence models to be very close to the corresponding one in smooth, CDM models. This is one of the main results of this paper and it will be tested and exploited in the next section.
5 Mapping smooth to clustering quintessence power spectra
We now study the nonlinear solutions for the total density power spectrum both for smooth and clustering quintessence. We compare different models sharing the same values for the following cosmological parameters: , , and . In addition, the amplitude of the initial power spectrum is also the same and corresponds to a for the CDM case . On the other hand, we consider different values of the DE equation of state parameter , both for the extreme cases and . In other terms, while the initial linear PS [13] is the same for all cosmologies, varying results in a different linear and nonlinear growth. In order to solve the PS evolution equations we set the initial conditions at redshift . This is done evolving back the CDM linear PS from by means of linear theory in the Newtonian approximation. We verified that this prevents transient issues in the evolution equations.
In Appendix C we assess the range of validity of the approximation, comparing the results with measurements in Nbody simulations for the CDM model, since no simulations including DE perturbations are available at the moment. Our technique reproduces at level the Nbody power spectrum in the BAO regime, confirming previous results in [55]. A accuracy extends to at and to at higher redshifts. At , the agreement with simulations is at the level only for , and it worsens to the level up to . This is, at some extent, expected given that we are working in the singlestream approximation. In [56], the authors estimate the magnitude of corrections due to multistreaming effects, which appears to be comparable to the discrepancy between the AP and Nbody results^{9}^{9}9Notice that in [56] the authors assume : this high value enhances multistreaming corrections..