Perturbation theory approach for the power spectrum: from dark matter in real space to massive haloes in redshift space

We investigate the accuracy of Eulerian perturbation theory for describing the matter and galaxy power spectra in real and redshift space in light of future observational probes for precision cosmology. Comparing the analytical results with a large suite of N-body simulations (160 independent boxes of 13.8 (Gpc/h)^3 volume each, which are publicly available), we find that re-summing terms in the standard perturbative approach predicts the real-space matter power spectrum with an accuracy of<2% for k<0.20 h/Mpc at redshifts z<1.5. This is obtained following the widespread technique of writing the resummed propagator in terms of 1-loop contributions. We show that the accuracy of this scheme increases by considering higher-order terms in the resummed propagator. By combining resummed perturbation theories with several models for the mappings from real to redshift space discussed in the literature, the multipoles of the dark-matter power spectrum can be described with sub-percent deviations from N-body results for k<0.15h/Mpc at z<1. As a consequence, the logarithmic growth rate, f, can be recovered with sub-percent accuracy on these scales. Extending the models to massive dark-matter haloes in redshift space, our results describe the monopole term from N-body data within 2% accuracy for scales k<0.15 h/Mpc at z<0.5; here f can be recovered within<5% when the halo bias is known. We conclude that these techniques are suitable to extract cosmological information from future galaxy surveys.


Introduction
Galaxy clustering is a key observational probe to study the large-scale structure of the Universe. The shape of the galaxy power spectrum, bispectrum and higher-order moments, contain information about the matter content of the Universe, about gravity and also about possible non-Gaussian initial conditions. In galaxy surveys, the galaxy distribution is distorted along the line of sight due to peculiar velocities that cause Doppler shifts, namely redshift space distortions (RSD). As these distortions depend on the growth of structure, they offer a complementary technique (to studies of the cosmic expansion history) to measure the matter content or to test gravity e.g., [1][2][3]. On very large scales and at higher redshifts, the RSD can be described by linear theory. However, on smaller scales and at later epochs, non-linearities start to play an important role and must be taken into account to accurately estimate the cosmological parameters from observational data.
Current surveys like BOSS 1 [4], and future missions like EUCLID 2 [5] will soon provide unprecedented datasets about the distribution of galaxies on large scales. In order to extract useful information from data of this quality we need more accurate theoretical models of structure formation. Standard perturbation theory (SPT) is the straightforward way of proceeding and has been the workhorse in the field for decades. However, its practical applications to statistics of the density field provide limited accuracy in both real and redshift space. A number of studies have investigated the fidelity of different theoretical models to describe the redshift space distortions and the possibility of extracting cosmological information from them [6][7][8][9][10][11][12][13][14][15][16][17][18][19][20][21]. This work explores the potential of describing the redshift-space distortions combining different mappings between real and redshift space with resummed perturbation theories. This is the basis to model the redshift-space distorted power spectrum of dark matter and dark-matter haloes. We then focus on the systematic errors that these models produce when they recover the logarithmic growth rate. For future surveys with forecasted statistical error on this quantity at the % level, the accuracy of the analytic description of clustering must be such that residual systematic errors due to modeling is well below this level. A fast, analytic description of clustering, calibrated on N-body simulations is an approach fully complementary to one based entirely on N-body simulations e.g., [22,23].
In particular, we start by describing the real-space power spectrum using 1-and 2-loop standard perturbation theory and the so-called renormalized (or resummed) perturbation theory [24,25]. We combine these models with different methods to obtain the redshift-space power spectrum: i) the Kaiser model [26], ii) the Scoccimarro model [27] and iii) the Taruya et al. model [28] (TNS model, from now on). For dark matter, we additionally account for the Finger-of-God (FoG) effect produced by virial motions on small scales by introducing a phenomenological damping term. This term is not needed in the case of haloes, as we only consider isolated haloes which are not part of a larger host halo. We compare these different methods with the results of a large suite of N-body simulations, focusing on the multipole expansion of the redshift-space power spectrum. Our suite of N-body simulations sums up to a volume of 2212 (Gpc/h) 3 , which is much larger than the volume surveyed by any current or planned experiment, ensuring that statistical errors in the simulations are negligible.
This paper is organized as follows. In §2 we present the basic theory of redshift-space distortions. In particular, in §2.1 we review standard perturbation theory and renormalized perturbation theory, extending the latter to higher-order propagators than previously considered. In §2.2 we present different RSD models while Finger-of-God effects are discussed in §2. 3. In §3 we provide details of the N-body simulations and describe the halo catalogues used in this paper. In §4 we compare the results of our models for the real-space power spectrum and compare them with N-body simulations for the dark-matter case. We also consider the RSD models mentioned above, focusing on their capacity to recover the logarithmic growth rate f , both for dark matter and massive haloes. Finally in §5, we discuss and summarize the obtained results. Appendices A and B contain details about standard and resummed perturbation theory providing a justification of the formulae presented in §2.1. Note that our discussion of resummed theories does not assume a field-theory background and is supposed to be accessible to all readers.

Theory
The matter-matter real-space power spectrum P δδ (k), the Fourier transform of the two-point correlation function, is the simplest statistic of interest one can extract from the dark matter overdensity field in Fourier space, δ(k), where δ D denotes the Dirac delta function and . . . the ensemble average. Under the assumption of an isotropic Universe, the power spectrum in real space does not depend on the direction of the k-vector. Since we only have one observable Universe, under the hypothesis of ergodicity the average . . . can be taken over all different directions for each k-vector.
The mapping between the radial coordinate in real space and the radial coordinate in redshift space is given by the Hubble flow and the Doppler effect due to the peculiar velocities, v. Since only the radial distance is computed from the measured redshift, the two angular coordinates remain the same in both real and redshift space. In this paper, we adopt the distant observer approximation, i.e., we assume that all line-of-sights are virtually parallel to each other. If we identify this direction with the third axis of our coordinate system, the mapping from real-space coordinates x to redshift-space coordinates s reads where H(a) is the Hubble parameter at the scale factor a andx 3 denotes the unit vector of the third axis. Using the scaled velocity field u ≡ −v/[H(a)af (a)] where f is the logarithmic derivative of the linear growth factor with respect to the scale factor, f ≡ d ln D(a)/d ln a, we can write this mapping as 3) The 2-point correlation function in redshift space is then defined by ξ s 2 (r) = δ(s + r)δ(s) (2.4) and its Fourier transform, P s δδ (k, µ) = d 3 r ξ s 2 (r) exp(−ik · r), (2.5) is the power spectrum in redshift space. Here, µ ≡k ·k 3 is the cosine of the angle between the vector k and the line-of-sight. The redshift space power spectrum, can also be computed from Eq. 2.2. In this mapping the mass must be conserved, which implies [1 + δ (s) (s)]d 3 s = [1 + δ(r)]d 3 r. Thus the transformation from δ(r) to δ (s) (s) at linear order in δ and v reads, where ∇ 3 is short for ∂ ∂x3 . The two-point correlation function in Fourier space, δ (s) (k)δ (s) (k ) is then [28], where r = x − x and ∆u z = u z (x) − u z (x ). Note that we have written u instead of v to make the dependence on f explicit. In Eq. 2.7 the enhancement and damping effect of the redshift space distortions on the power spectrum are manifest. The enhancement due to RSD, also known as Kaiser effect, is produced by the the +f ∇ 3 u 3 terms in Eq. 2.2, that increase the overdensity δ(x). These terms, represent the coherent distortions by the peculiar velocities along the line-of-sight direction, and are controlled by the growth factor parameter f . On the other hand, the damping effect comes from the exponential factor in Eq. 2.7. This term is mainly due to the small scale velocity dispersion around the most clustered regions, and produces the suppression of power at small scales in the power spectrum.
In this paper we focus on two approaches: standard perturbation and resummed perturbation theory, both in Eulerian space.
Standard perturbation theory (SPT hereafter) consists of expanding the statistics of interest as a sum of infinite terms, where every term correspond to a n-loop correction. For the power spectrum the SPT prediction is written as (see appendix A for the explicit formulae of SPT terms), P SPT (k) = P (0) (k) + P (1) (k) + P (2) (k) + . . . . (2.8) The 0-loop term correction is just the linear power spectrum, P (0) (k) = P lin (k). The 1-loop term is expressed as a sum of 2 different subterms, where the subscripts i and j refer to the perturbative order of the terms δ(k) used in eq. 2.1 to compute the power spectrum P ij (k). In this case, both P 13 and P 22 requires a 2-dimensional integration (after exploiting rotational invariance). The 2-loop term is the sum of three different subterms, In this case, all these three terms, require a 5-dimensional integration (after exploiting rotational invariance). One can keep going to higher-order terms. However, the 3-loop correction term requires the computation of 8-dimensional integrals and the 4-loop correction term 11-dimensional integration. For practical and computational reasons, one does not usually go beyond the 2-loop correction terms. Thus, one has to truncate Eq. 2.8 series at some loop order. Truncating it at P (1) and at P (2) term is what is respectively called in this paper, 1L-SPT and 2L-SPT.
Renormalized perturbation theory (RPT hereafter) attempts to reorganize the perturbative series expansion of SPT and resum some of the terms into a function that can be factorized out of the series. This function is usually called the resummed propagator and we refer to it as N . In order to make the resummation possible, all the kernels of the P ( ) terms have to be expressed as a product of kernels that correspond to full-mode-coupling terms and full-propagator terms (see appendix B for details). The full-mode-coupling kernels are those kernels contained in P nn (k) terms that contain a coupling of the form, k−q 1 −. . .−q n−1 (see Eq. B.5). The full-propagator kernels are those contained in P ( ) terms of the form P 1n with no-mode coupling term (see Eq. A.4 and A.6 as examples). The resulting expression of resumming terms in this way is (see appendix B and Refs. [24,25] for a full derivation), where the term P (n−1)L nn is the part of the P nn term that describes a full-mode coupling (see Eq. B.5). In spite of the resummation, Eq. 2.11 contains an infinite series as Eq. 2.8 and has to be truncated after a certain number of loops. However, some of the infinite terms of Eq. 2.8 have now been reorganized into the N i function. We will refer as 1L-RPT-N i and 2L-RPT-N i the truncation of Eq. 2.11 at 1-and 2-loop, respectively. Finally, the form of the function N i depends on the way we approximate the kernels in the resummation process. In the case the kernels are approximated according to the Zel'dovich approximation (see Eq. B.1), they are expressed as a product of 0-loop propagators and the resulting function N 0 is, with, (2.13) When the kernels are approximated as a product of 1-loop propagator kernels (see Eq. B.17-B.19) the resulting N 1 function is, This is the expression presented in [56] (see Eq. B.31 for the correspondence). Expressing the kernels as a product of 2-loop propagators, (see Eq. B.32-B.37) yields the function N 2 , Note that the angular part of the propagator terms, P 13 , P 15 , P 17 , 3 . . . is analytically integrable for any shape of the power spectrum. Thus, the 2-dimensional integration of P 13 can be reduced to a 1-dimensional integration (see Eq. B.31). In the same way, the 5-and 8-dimensional integrations of the terms P 15 and P 17 , are reducible to 2-and 3-dimensional integration. However this is a hard task due to the symmetrized kernels, which are constructed as a sum of 5! = 120 and 7! = 5040 terms respectively (see Eq. A.14). Because of that, we stop at 2-loop. In this paper, the Nfunction is computed considering P 15 as a 5-dimensional integral. If more accuracy is needed, the 3-loop resummed propagator written in Eq. B.43 can be used. These extensions in the resummed propagator could be easily incorporated in the current public codes for the RPT [50,54,56].

Perturbation theory in redshift space
In order to describe the non-linear matter power spectrum in redshift space (Eq. 2.5) we need a model that, given the non-linear power spectrum in real space, is able to "map it" to the power spectrum in redshift space. There are several models that attempt to do this task. While in principle the same perturbation theory approach used for the dark matter could be employed to model also the velocity and the density velocity coupling yielding therefore a real-to-redshift space mapping, it has become clear in the literature that the redshift space clustering and in particular the redshift space power spectrum is not well described perturbatively. In fact, highly non-linear scales are superimposed to linear scales by the real-to-redshift space mapping, as realized in the seminal papers by [26,59]. Some of this effect is even visually apparent in the galaxy distribution as the so-called "Fingers of God" effect. In this paper we consider physically motivated, but phenomenological models, we study the Kaiser model [26], the Scoccimarro model [27] and the TNS model [28]. All these models propose a functional form of P s (k) that depends on real-space statistics. The simplest model, is the Kaiser model proposed by Nick Kaiser 25 years ago [26], where b(k) is a possibly scale-dependent biasing function, which relates the observable tracers to the underlying dark matter distribution. This expression is obtained when Eq. 2.7 is treated linearly. Because of that, in principle P δδ should refer to the matter-matter linear power spectrum. However, the prescription P lin δδ → P nl δδ has been demonstrated to work better. With this recipe, the Kaiser model is usually known as 'non-linear Kaiser'. In this paper we refer to Eq. 2.16 with a non-linear P δδ simply as Kaiser model.
The Scoccimarro model proposes that the redshift-space power spectrum is given by [27], where P δθ and P θθ are the velocity-matter and velocity-velocity power spectra in real space, respectively, and are defined by . Note that Eq. 2.17 tends to 2.16 when P δθ and P θθ tend to P δδ . This is the case for the linear regime in SPT. The TNS model [28] takes into account the cross interaction due to linear and non-linear processes. This produces two extra terms to Eq. 2.17, (2.20) where the A and B terms, arise from the interaction between the enhancement terms and the damping terms in Eq. 2.7. Also [57] have provided a complementary explanation and a possible generalization of the A and B terms. However, in this paper we use the basic formalism presented by [28], where, with k 123 ≡ k 1 + k 2 + k 3 . Since the functions A and B require an integration over the whole range of momenta q, we cannot use RPT predictions to compute them. Furthermore, A requires the crossed bispectra between δ and θ, whose computation at 1-and 2-loop requires more effort than in the power spectrum case. Since we expect A and B to be small compared to P δδ , P δθ and P θθ [28], in this paper we compute A and B assuming the leading terms for the power spectrum and bispectrum inside the integrals of Eq. 2.21 and 2.22. For the rest of the power spectrum terms of Eq. 2.16, 2.17 and 2.20 we use the perturbative approaches described in section 2.1. None of these models is able to account for non-linear effects such as Fingers of God. These effects have to be included ad hoc through a function that damps the power spectrum in redshift space at small scales.

Fingers of God
The effects of small-scale velocities are not completely included in the models presented in section 2.2. Both the Kaiser and Scoccimarro model ignore them as we have already commented. Only the TNS model takes them into account, but only as a cross term with large-scale squashing. Thus, the effect of these small-scales velocities has to be inserted as a multiplicative damping function into the models described in section 2.2 [58][59][60], P s (k, µ; z) → P s (k, µ; z)D 2 FoG (k, µ; z). (2.25) The most used prescriptions for that, are the Gaussian and the Lorentzian functions that both depend on a redshift-dependent parameter, σ(z) ≡ σ 0 D(z), Gaussian. (2.27) where σ 0 ≡ σ(z = 0). Theoretically, it has been suggested [27] that σ(z) can be computed analytically as, Such a parameter is physically motivated as it tries to enclose the effect of the velocity dispersion of dark matter particles. Similar modeling was discussed e.g. in [61] and in [27], although the two have different interpretations. In any case, the fit and theoretical numerical value for σ need not to coincide with the actual value of the velocity dispersion. The theoretical value is computed in the linear approximation and the fit value is obtained using a Gaussian approximation while the actual velocity dispersion is highly non-gaussian even on large scales. More discussion on this in [27]. In most of the cases, this modeling has shown a very poor agreement with N-body simulation results. Therefore, in this paper we always treat σ 0 as a free parameter to be fit from N-body simulations.

Simulations
The simulations used in this paper model the structure formation on very large scales within a flat ΛCDM cosmology consistent with current observational data. The adopted cosmological parameters are: Ω Λ = 0.73, Ω m = 0.27 h = 0.7, Ω b h 2 = 0.023, n s = 0.95 and σ 8 (z = 0) ≈ 0.8. Our suite of simulations consists of 160 independent runs with 768 3 particles in a box of length 2.4 h −1 Gpc. Hence, each box has a volume of V 1box = 13.8 (Gpc/h) 3 and in total, we simulate about 2, 200 (Gpc/h) 3 . The initial conditions were generated at the starting redshift z = 19 by displacing the particles according to the second-order Lagrangian PT from their initial grid points. The initial power spectrum of the density fluctuations was computed with CAMB [62]. The simulations were performed with the GADGET-2 code [63] taking only the gravitational interaction into account.
In this paper, we consider snapshots at z = 0, 0.5, 1, 1.5. In order to obtain the dark-matter field from particles we discretize each box using 512 3 grid cells. Thus the size of the Cartesian mesh is 4.68 h −1 Mpc. We assign mass to the cells using the cloud-in-cell prescription. Using a much higher resolution simulation, we checked that the power spectrum derived from the simulation data is accurate at the 1% level up to k < 0.2 hMpc −1 .
To identify the dark-matter haloes, we used the Amiga Halo Finder [64,65], which defines a halo by the bound dark-matter particles inside a spherical overdensity equal to the so-called virial overdensity. We only consider haloes which are at least resolved by 40 particle. This leads to a minimum halo mass of 10 14 M /h.
The errors associated to the statistical quantities measured from the simulations are obtained from the dispersion among the 160 runs: we report the error on the mean of the independent runs. We make the dark matter and halo power spectra publicly available and also provide the multipoles measured from the simulations used in this paper for possible comparisons 4 .

Results
In this section we compare and test the theoretical formalism presented in §2 (for both real and redshift space) with the N-body simulation data. We start with perturbation theory in real space, and we follow the theoretical predictions of different models for redshift-space power spectrum: Kaiser (Eq. 2.16), Scoccimarro (Eq. 2.17) and TNS model (Eq. 2.20). In particular, we consider here 1-and 2-loop SPT, and 2-loop RPT with N 1 and N 2 . We focus on the multipole prediction of these models but also on the accuracy on recovering the logarithmic growth rate f . We explore this for both dark matter and for massive dark matter haloes.

Performance of the perturbation theory approach: comparison to N-body simulations
In this subsection we test the formalism presented in section §2.1, both SPT and RPT, and we compare them with the outcome of the N-body simulations. In Fig. 1 -2 (left panels) we compare the power spectra obtained with different flavors of perturbation theory against the N-body simulation results at z = 0, 0.5, 1.0, and 1.5. Each power spectrum has been normalized to the linear non-wiggle power spectrum to reduce the dynamical range. In the left top subpanels the power spectrum of the N-body simulations is shown as black circles linked with a black dot-dashed line, the linear theory prediction is shown as a dotted black line, SPT in blue lines, RPT-N 0 model in red lines, RPT-N 1 model in green lines and RPT-N 2 model in orange lines; where the infinite series of Eq. 2.8 and 2.11 are truncated at 1-loop (solid lines) and 2-loop (dashed lines). In the bottom left subpanel the ratio of these models with N-body simulation data is shown with the same color and line notation. In that case, circle symbols correspond to linear prediction, square symbols to 1-loop truncation and triangle symbols to 2-loop truncation.
We see that SPT with 1-loop truncation over-predicts the N-body power spectrum at all redshifts, and can only make an accurate prediction of the power spectrum (≤ 1% deviation) at very large scales: k ≤ 0.05 h/Mpc for z = 0 and k ≤ 0.10 h/Mpc for z = 1. Going to 2-loop correction improves considerably the behavior of SPT, but one has to deal with the terms P 33 , P 24 and P 15 that require a 5-dimensional integration. In this case, 2-loop SPT makes a very good prediction of the N-body . Left bottom subpanel: ratio between the power spectrum of N-body simulation and different PT models with the same color notation. Right top subpanel: power spectrum at z = 0 up to relatively small scales (k ≤ 0.20 h/Mpc), but at z = 0.5 starts overpredicting it, and at z = 1 and z = 1.5 the over-prediction is a few percent at intermediate scales The RPT-N 0 model for both 1-and 2-loop truncation behaves accurately only at very large scales where it shows 1% deviation respect to N-body simulations both at all redshifts, but breaks down at relatively large scales: k 0.03 h/Mpc for 1-loop at z = 0 and k 0.07 h/Mpc for 2-loop at z = 0.
The RPT-N 1 model presents at large scales a very good agreement with the N-body simulation results with the advantage that the breakdown happens at smaller scales: k 0.1 h/Mpc for 1loop and k 0.15 for 2-loop at z = 0; at k 0.15 h/Mpc for 1-loop and at k 0.25 for 2-loop at z = 1. However, the for 2-loop truncation over-predicts the N-body power spectrum with a systematic ∼ 2% deviation at intermediate scales at all redshifts. This is due to the limitation in the resummed propagator N 1 , which is expressed as 1-loop expansion terms (in terms of P 13 ).
Finally, the RPT-N 2 model presents a modest improvement over previous models on large scales at z = 0, but breaks down already at k 0.07 h/Mpc and k 0.10 h/Mpc for 1-loop and 2-loop, respectively. However at z > 0, this model presents a better behavior, with < 2% accuracy, breaking down at k 0.20 h/Mpc at z = 1 and fixing the systematic ∼ 2% over-prediction observed for the RPT-N 1 model. The accuracy of this model at z ≥ 0.5 is then better than the 2-loop RPT-N 1 and 2-loop SPT.
In general, we see that the model that best describes the N-body data at z ≤ 1.5 is the RPT-N 2 model at 2-loop truncation. The 2L-RPT-N 1 model shows also good results but presents a ∼ 2%systematic over-prediction at intermediate scales. However at z = 0, RPT-N 1 is able to reach smaller scales than RPT-N 2 , which breaks down at relatively large scales. Indeed, 2L-RPT-N 2 works better than 2L-RPT-N 1 and 2L-SPT at z ≥ 0.5 but not at z = 0.
In Fig. 1 -2 (right panels) the behavior of the damping functions of RPT-N i models is shown for the same redshift range. In particular, in the right top subpanel we show the scale dependence of the damping functions of Eq. 2.11: (green dashed line) and N 2 (k) 2 (see Eq. 2.15) (orange dotted line). We see that at large scales, all N i functions converge to 1, and it is at intermediate and smaller scales (k > 0.03 h/Mpc) where the differences between these three models are significant. For the scales of interest, we see that N 0 (k) < N 2 (k) < N 1 (k). Thus, as we include higher-order propagators in the computation of N i , these functions oscillate about the 'true' damping function. This approach explains why RPT-N 0 under-predicts N-body data and why RPT-N 1 slightly over-predicts it. Every loop correction in the resummed propagator tends to the true value, but in a oscillatory way. For z ≥ 0.5, the 2-loop correction in the resummed propagator seems to be sufficient for a ≤ 1% prediction. However at z = 0, the convergence in the resummed propagator is still not reached for the 2-loop correction in N . In that case, higher-order loop corrections would be necessary to reach the ≤ 1% deviation. In the right bottom panel we show the effective σ 2 v , i.e. −2 ln[N i (k)]/k 2 , for the different orders in the resummed propagator: 0-loop (red solid line), 1-loop (blue dashed line) and 2-loop (orange dotted line). In that case, we see that at large scales the three models diverge, whereas at small scales all models seem to converge (at least for z ≥ 0.5).
In other words, it may happen that a lower-order approximation appears to work better than a higher-order approximation which, in principle, should be more accurate. This is due to a fortuitous cancellation of the truncation errors in the propagator and in the damping function. This cancellation does not hold for all redshifts and/or all cosmologies. The performance of an analytical approximation scheme must be quantified looking at different redshifts (or different cosmologies).
The formalism presented here deals with the F n kernels, that correspond to the δ-field. However, this formalism is also perfectly valid for the computation of P δθ and P θθ , only changing appropriately the F n kernels by the G n kernels as in SPT.

Dark matter multipoles
In order to obtain information about the growth rate f from the redshift-space distortions, it is convenient to work with the expansion in Legendre moments, P , defined as, where L are the Legendre polynomials of order . For the first three non-vanishing P , According to linear theory (Kaiser model with P lin δδ ) only the monopole ( = 0), the quadrupole ( = 2) and the hexadecapole ( = 4) are different from 0. In that case, for an unbiased tracer these three moments read, Hence, knowing the dark matter power spectrum in both real and redshift space, one can directly measure the growth rate f from any of these multipoles. It is interesting to note that the ratio between any of these multipoles does not depend (at large scales) on the real-space power spectrum and the ratio tends to a constant that only depends on f when k → 0. However, non-linearities produce deviations from these formulae. Depending on the ability of modeling the non-linearity in the redshift space distortions, we will be able to use information from non-linear scales to estimate f with accuracy.
In this section we focus on checking the quality of the different theoretical models in predicting the multipole power spectrum of dark matter in redshift space. We focus on the models described in section 2.2, using as real-space inputs, the PT-theory approaches described in section 2.1. Here, we assume that f is known and we only fit the FoG parameter σ 0 assuming a Lorentzian damping function (Eq. 2.26), although no significant difference is observed when a Gaussian damping function is assumed. We allow σ 0 to depend on z and on k max 5 and find the best-fit value by minimizing Here, the subscript "sims" refers to the simulations and "theo" to the theoretical models described above. σ 2 P,sims (k) is the variance of the multipoles computed from the 160 simulations and k o is the minimum k considered, which is set by the size of the simulation box but has little effect on the final results. Note that we neglect for simplicity any covariance between different k-bins, which is a good approximation for small k and broad bins.
In Fig. 3, we show the measurements of the monopole (top panels), quadrupole (middle panels) and hexadecapole (bottom panels) from N-body simulations (black empty circles) in top sub-panels. In all cases, the multipoles have been normalized to the non-wiggle linear prediction to reduce the dynamical range. In the top subpanels, different theoretical predictions are shown: the Kaiser model (dotted lines), the Scoccimarro model (dashed lines) and the TNS model (solid lines). The chosen real-space power spectrum for each of these models is: linear prediction (black dotted lines), 1L-SPT (red lines), 2L-SPT (blue lines), 2L-RPT-N 1 (green lines) and 2L-RPT-N 2 (orange lines).
In all cases, the TNS model combined with both SPT and RPT predictions is the model that describes best the N-body results. Using the TNS model it is possible to achieve < 1% accuracy for the monopole up to k ≤ 0.12 h/Mpc at z = 0 and k ≤ 0.17 h/Mpc at z = 1. The models are also very accurate for the quadrupole: at z = 0 we can describe N-body data up to scales of k = 0.12 h/Mpc and at z = 1 up to k = 0.30 h/Mpc with a deviation 2%. For the hexadecapole, the agreement is more modest: at z = 0 we can only achieve ∼ 10% accuracy up to scales of k = 0.20 h/Mpc at z = 0 and ∼ 5% at z = 1.
Both the Scoccimarro and Kaiser model provide a reasonably good approximation on large scales but both fail to give an accurate description on mildly non-linear scales where baryon acoustic oscillations (BAO) are located. The difference between the TNS model and the other models is more evident for the quadrupole and hexadecapole possibly suggesting that non-linearities become more important for higher-order multipoles.
The imprint of the BAO in the multipoles is clearly visible: note that the Scoccimarro and Kaiser models slightly over-predict the BAO amplitude, especially for the quadrupole, while the TNS model does better, although a trend towards the under-prediction is observed. All models correctly predict the BAO location. These considerations might be relevant for recovering in an unbiased way the angular and radial BAO information (separately) from forthcoming surveys.
We conclude that the TNS model with the RPT and SPT models studied here, has the ability of describing the redshift space power spectrum monopole and quadrupole at z = 0 and z = 1 within 1 − 2% for k 0.2 and the hexadecapole within about 5%.
We do not observe a crucial difference between 1-and 2-loop SPT. Also, no significant difference between using N 1 and N 2 for 2-loop RPT is detected. This indicates that on these mildly non-linear scales at redshifts 1.5 the accuracy of the modeling of redshift space distortions is more important than that of the non-linear evolution of the real-space dark matter power spectrum.
Because of that, for simplicity we focus on 1-loop SPT and 2-loop RPT with N 1 when measuring f from dark-matter multipoles in the next subsection.
In Fig. 4 we show the best-fit value for σ 0 corresponding to the fits shown in Fig. 3 using the same color notation for the different models. Additionally we show the theoretical predictions of Eq. . Multipoles corresponding to dark matter power spectrum: monopole (top panels), quadrupole (middle panels) and hexadecapole (bottom panels), for z = 0 (left panels) and z = 1 (right panels), where f is fixed to the true value, and σ0 is the only free parameter, fit to N-body data. In the upper subpanels the value of the corresponding multipole is shown, normalized to the linear, no-wiggles value to reduce the dynamical range. In the lower subpanels the ratio between the N-body simulation data and different perturbation theory predictions is shown.  2.28 using as an input the 1L-SPT prediction for P θθ (solid black line) and P lin (dot-dashed line). As indicated, top, middle and bottom panels correspond to the monopole, quadrupole and hexadecapole, whereas left panels show the result for z = 0 and right panels for z = 1.
We note that at z = 1 all the models produce a best-fit σ 0 which is close the the theoretical predictions, although an overestimate is observed for the monopole and underestimation for the hexadecapole, being the quadrupole the case which is closer to the theoretical prediction. For z = 0, the discrepancy between theory and best-fit value is larger. We will analyze again the agreement between theory and best-fit σ 0 in the next section, when will allow f also to vary.

Estimating f from dark matter multipoles
In the last section we have shown that the TNS model was able to describe well the multipoles when one free parameter was allowed to vary in order to account for the FoG effect. In this section we want to check the ability of these models to recover f from the dark-matter field. In this case, we will allow both f and σ 0 to freely vary. As before, in order to find the f and σ 0 best fit parameters, we minimize the χ 2 value.
In Fig. 5 we show the obtained values for f (top subpanels) and σ 0 (bottom subpanel) as a function of the maximum scale used in the fitting, namely k max . We show the results at different redshift: z = 0 (left panels) and z = 1 (right panels). Top, middle and bottom panels show the derived values for f and σ 0 for each of the multipoles: monopole, quadrupole and hexadecapole as indicated. As in Fig. 3, dotted lines stands for Kaiser model, dashed lines for Scoccimarro model and solid lines for TNS model. For simplicity, we only show the results corresponding to 1L-SPT (red lines) and 2L-RPT-N 1 (green lines). 2L-SPT and 2L-RPT-N 2 yield similar results. In top subpanels, horizontal solid black line shows the true value for f , whereas black dot-dashed horizontal lines show the 1% and 2% deviation, as labeled. In the bottom subpanels, the horizontal lines stands for the theoretical predictions of σ 0 according to Eq. 2.28 when P θθ (blue line) and P lin (orange line) are used as inputs. In the case of P θθ we use 1L-SPT prediction. All the error bars correspond to 1−σ errors for f and σ 0 , and have been computed from the contour in the f -σ 0 space that corresponds to ∆χ 2 = 2.3. Since our N-body sample consists of 160 realizations (∼ 13 3 [Gpc/h] 3 volume), the dispersion in the measured monopole, quadrupole and hexadecapole is small. Therefore the corresponding error bars for the recovered parameters are also small, especially in the case of the monopole and quadrupole. However, one should be aware that these statistical errors in f and σ 0 are not comparable to the  Figure 5. Estimates for f (top subpanels) and σ0 (bottom subpanels) from the dark matter multipoles of the N-body data: monopole (top panels), quadrupole (middle panels) and hexadecapole (bottom panels); for z = 0 (left panels) and z = 1 (right panels expected errors from future surveys. These errors just provide information about the uncertainties and shortcomings of the models. Note that for most of the studied models, these errors are much smaller than the expected ones from any galaxy survey. From Fig. 5, we see that the TNS model using both 1L-SPT and 2L-RPT-N 1 is the only model able to recover the value of f to 1% accuracy: for the monopole up to k 0.15 h/Mpc for z = 0 and up to k 0.20 h/Mpc for z = 1; for the quadrupole up to k 0.20 h/Mpc for z = 0 and up to k 0.25 h/Mpc for z = 1; and for the hexadecapole up to k 0.15 h/Mpc for z = 0 and up to k 0.20 h/Mpc for z = 1. However, in the case of the hexadecapole, the statistical errors are too large to be able to claim 1% accuracy of the model predictions. Also note the difference between the theoretical value of σ 0 and the best-fit value obtained from the N-body data: at z = 1 for both monopole and quadrupole, the best-fit σ 0 value for TNS model + 2L-RPT-N 1 yields a very similar result as Eq. 2.28. However at z = 0 there is a large discrepancy between these two values.

Halo biasing and stochasticity
So far, we have been able to recover the f parameter with high accuracy using the dark matter multipoles. However, galaxy redshift surveys consist of galaxies residing in dark matter haloes, which are biased and stochastic tracers of the underlying dark matter field. Furthermore, since we only consider massive isolated haloes, we do not expect any FoG effects in the halo redshift-space statistics.
In this paper, we model the biasing, i.e. the relation between the halo overdensity and the the dark matter density contrast as, where b(k) is the scale-dependent bias function and describes a stochastic field. The stochastic field stands for any physical or statistical process that produces a non-deterministic relation between the dark matter and the halo field. This includes the shot noise due to the discrete nature of haloes 6 . For a Poisson process the shot noise is inversely proportional to the mean number density, namely n h , However, the formalism used here allows for other stochastic processes. Assuming the field to be uncorrelated with δ, namely, δ = 0, the bias function can be written as, where the second equality stands for the numerator and denominator independently. Here the subscripts "h" and "m" stand for haloes and dark matter respectively. The power spectrum of the field can be computed combining Eq. 4.9 and 4.11, Finally, we define the (shot noise free) halo-halo power spectrum P hh as, Note that by this definition the following equality holds, (4.14) Hence, in this biasing scheme, the bias functions that relate the halo-matter and halo-halo power spectra to the matter-matter power spectrum are the same. This is a reasonable approximation at least on large scales [66], where the bias becomes linear. In Fig. 6 we show the scale dependence of the halo bias (left panel) and the -field power spectrum (right panel) measured from the halo catalogues at different redshifts. Remember that the minimum    halo mass is M cut = 10 14 M /h, which ensures that all the haloes have at least ∼ 40 particles. The number density of haloes with this mass cut at different redshifts is shown in Table 1. Fig. 6 shows that the bias increases with z and with k. The P increases with z and slightly decreases with k. For M cut = 10 14 M /h, the P is sub-Poissonian at low redshifts, 65% at z = 0 and 90% at = 0.5, but turns out to be very close to the Poissonian prediction at high redshifts, z = 1 and z = 1.5.
In Fig. 7 we show the signal-to-noise ratio (or P hh /P ) for the haloes studied here at different redshifts: z = 0 (red line), z = 0.5 (blue line), z = 1.0 (green line) and z = 1.5 (orange line). Horizontal dotted lines mark the values P hh /P = 1.0, 0.5 and 0.1 as a reference. In Table 2 the scale at which the signal-to-noise ratio reaches these values is written for the same z-snapshots. We see that only for z = 0 and z = 0.5 the signal-to-noise ratio is above 0.5 at scales k < 0.1h Mpc, whereas for z = 1.0 is above 0.5 only at very large scales, k < 0.03 h/Mpc and never for z = 1.5. Conservatively, in this paper we consider only the scales where P hh /P mm 0.5 to be suitable for extracting information. Hence, we do not study the halo power spectra at z = 1.5 and z = 1.0 because only at very large scales (where the behavior is linear) the signal-to-noise ratio satisfies this condition.

High-bias halo power spectrum multipoles
The multipoles for the halo power spectrum are defined by Eq. 4.1 just changing P s mm (k, µ) by P s hh (k, µ). In this case, we define the halo-power spectrum in redshift space by where we assume that P does not depend on µ. For the monopole term, the shot noise subtraction is important. However, since we are assuming that P does not depend on µ, it is irrelevant for higher-order moments. Indeed, a µ-independent offset on P s hh (k, µ) has no effect in the quadrupole and hexadecapole, only in the monopole. As in Eq. 4.5-4.7, at very large scales, the halo multipoles are written as, Note that for the case of biased tracer, f and b are degenerate in the linear regime when we treat b 2 P lin as input or when any ratio between different multipoles is used to constrain f . In this case, only the ratio among them, β(k) ≡ f /b(k) can be measured. However, as for the dark matter case, non-linearities cause deviations from these formulae. Depending on the ability of modeling the redshift space distortions, we will be able to use information from non-linear scales to estimate f with accuracy but also, in the case of biased tracers, we might be able to break the degeneracy between the bias and f . For both the Kaiser and the Scoccimarro model, f and b always appear in the β combination, only for the TNS model this degeneracy is not exact. Fig. 8 shows the halo monopole (top panels) and halo quadrupole (bottom panels) for z = 0 (left panels) and z = 0.5 (right panels). Both f and b(k) have been set to their true values. The color and line notation is the same as in Fig. 3: dashed lines stands for Scoccimarro model and solid lines for TNS model. The different real-space power spectrum inputs are: linear prediction (black dotted lines), 1L-SPT (red lines), 2L-SPT (blue lines), 2L-RPT-N 1 (green lines) and 2L-RPT-N 2 (orange lines). Since we do not expect Fingers of God for isolated haloes, the theoretical models have no free parameters in this case.  . Multipoles corresponding the halo-halo power spectrum: monopole (top panels), quadrupole (bottom panels), for z = 0 (left panels) and z = 0.5 (right panels). In the upper subpanels the value of the corresponding multipole is shown, normalized to the corresponding non-wiggle linear prediction to reduce the dynamical range. In the lower subpanels the ratio between the N-body simulation data and different PT-predictions is shown. Dashed lines corresponds to Scoccimarro model and solid lines to the TNS model. Different PT models are shown: linear prediction (black dotted lines), 1L-SPT (red lines), 2L-SPT (blue lines), 2L-RPT-N1 (green lines) and 2L-RPT-N2 (orange lines). In bottom subpanels the 2% deviation is marked with black dot-dashed horizontal lines. Vertical dot-dashed lines mark the regions where P hh /P = 1.0 and P hh /P = 0.5 as labeled.
From Fig. 8, we see considerably different results compared to the dark matter case (see Fig.  3). The accuracy of the modeling is reduced, especially for the quadrupole at z = 0. In the case of the halo-halo monopole, we see that the different PT models make very different predictions, while this was not the case for the dark matter. In particular we see that TNS + 2L-RPT-N 2 is the only model able to make sub-percent predictions at z = 0 up to k = 0.15 h/Mpc and at z = 0.5 up to k = 0.10 h/Mpc for the monopole. Any other PT theory + TNS yields worse results. We also see that Scoccimarro + 1L-SPT and Scoccimarro + 2L-RPT-N 1 provide a good description at z = 0 but not at z = 0.5. This is due to an (accidental) cancellation of two terms that go in opposite directions. The Scoccimarro model does not take into account the A and B functions of the TNS model, that add a positive contribution to P s (see Fig. 11 in the next section). On the other hand, in Fig. 1 we have seen that 1L-SPT and 2L-RPT-N 1 over-predict the true value for P δδ at z = 0 (but also for P δθ and P θθ ). For the dark matter power spectrum, both effects are small and negligible. However, for massive haloes, the high bias increases these effects. At z = 0.5, the TNS approach is clearly modeling better the monopole than the Scoccimarro model, which under-predicts the N-body data by ∼ 10%.
In the case of the halo-halo quadrupole, we see that all the PT models make similar predictions, and that the main difference is due to the RSD model chosen but the modeling breaks down at relatively large, almost linear, scales. We see that the TNS model is able to describe well the quadrupole only up to k = 0.05 h/Mpc at z = 0 and up to k = 0.10 h/Mpc at z = 0.5.
In general we see that for the monopole, all the models describe the N-body data better at z = 0 than at z = 0.5. This seems counter-intuitive, because at higher redshifts, non-linearities are less important and perturbation theory should work better. However, it could be explained by the fact that at z = 0.5 the signal-to-noise ratio is considerably less than at z = 0 (see Table 2 and Fig. 7). For instance, recall that at z = 0, P hh /P = 1.0 at k = 0.168 h/Mpc whereas at z = 0.5 this happens already at k = 0.065 h/Mpc. Also the bias of the selected haloes grows with redshift, making more apparent any systematic error in the model (of both RSD and biasing). It is reasonable to expect (but it remains to be tested) that the model performance improves for lower mass -thus less rare and less biased-haloes.

Simultaneously estimating f and b from halo multipoles
In this section we show how well the f parameter can be recovered from the halo-monopole N-body data. Since in the last section we have seen that none of the models studied here is able to reproduce the quadrupole data for haloes with sufficient accuracy at the mildly non-linear scales we are interested in, we do not try to recover f from P 2 . Instead we focus on the degeneracy between f and the bias in the monopole. According to Eq. 4.5, these two parameters are perfectly degenerate when P (k)b 2 is set from observations or when the P 2 /P 0 ratio is used to compute f and b. This is the case for the Scoccimarro and Kaiser models without Finger-of-God effects. However, non-linear terms, namely the A and B functions of the TNS model, are expected to break this degeneracy at non-linear scales even when P (k)b 2 is fixed: from Eq. 2.21-2. 22, we see that f and b do not appear always in the same combination in the A and B functions. In particular A can be expressed as b 3 A(k, µ, β) and B as b 4 B(k µ, β). Since in this paper the input is the dark matter power spectrum, the degeneracy between b and f is not perfect, and there is the possibility of recovering these parameters separately with certain accuracy even for the Kaiser and Scoccimarro models. For simplicity, we do not model or fit the scale dependence of the bias function b(k), but instead assume that we do know the scale dependence of the bias, and try to recover the growth rate f and the bias amplitude, Realistic approaches need an analytical modeling of the scale dependence of the bias, that in principle one could expand perturbatively e.g., [52,67]. Here we focus on the modeling of the dark matter and redshift space distortions and thus we measure the bias directly from the simulations. First, we assume that we also know that A b = 1 when we fit f . The fitting results are shown in Fig. 9 using the same color and line notation as in Fig. 8. As before, the error bars correspond to the interval defined by ∆χ 2 =1.0.
As we have commented for the dark matter case, these errors are much smaller than the ones which could be obtained form a real survey. They contain information of 160 realizations of a large volume, the total volume is close to that enclosed by an all-sky survey up to z = 25.
For both z = 0 and z = 0.5 the models that do best at recovering f are TNS + 2L-SPT (solid blue line) or TNS + 2L-RPT-N 2 (solid orange line), especially the latter one. In particular, at z = 0, the TNS model with 2L-RPT-N 2 implementation is able to recover the f parameter within a 5% accuracy up to scales k = 0.17 h/Mpc and up to scales k = 0.13 h/Mpc at z = 0.5. Scoccimarro models with accurate implementations of the real-space power spectra such as 2L-SPT and 2L-RPT-N 2 largely over-predict f at k > 0.05 h/Mpc.
Next we show in Fig. 10 the results of extracting both A b and f from the halo monopole at the same time. We only show the two models that performed best at extracting f when A b is known: TNS + 2L-SPT (solid blue line) and TNS + 2L-RPT-N 2 (solid orange line). Since now there are two free parameters, the 1 − σ error bars correspond to contours of ∆χ 2 = 2.3 in the f − b space. We see that in general, at z = 0 and z = 0.5, 2L-RPT-N 2 works better than 2L-SPT, as we have already observed in Fig. 9. However, 2L-RPT-N 2 tends to overestimate the bias amplitude and underestimate the f parameter. In particular, at z = 0 and z = 0.5, at scales around k 0.11 h/Mpc, this model overestimate A b by 1% and underestimate f by about 15%. However, when A b was set 1, the systematic error of f was less than 5%. Similar results were obtained by [12], when the TNS model was combined with the Closure perturbation theory of [53]. We have checked that this large underestimation in f is due to the 1% overestimation in the bias amplitude. In other words, if the bias amplitude is set by hand to a fixed value of A b = 1.01, we obtain similar results as in Fig 10. Thus, this formalism tends to underestimate f by 15% while overestimating A b by only 1%. Since the Scoccimarro model + 2-loop perturbation theory predictions were not able to predict f when A b was assumed to be 1, the TNS terms A and B are the key ingredient of the TNS model to achieve a high accuracy recovering f when the bias amplitude is assumed. However, the signal-to-noise ratio of these terms is low compared to the signal-to-noise of the whole monopole term. In Fig. 11 we show this for z = 0 (left panel) and z = 0.5 (right panel). In top subpanels, the total contribution of Scoccimarro and TNS models is shown in red and blue lines respectively. In green and orange lines, the isolated contribution of A and B of the TNS model is shown respectively. Black line shows the stochastic noise, P . In bottom subpanels, the ratio of all the signal terms with respect to P is shown. As a reference, the lines where S/N = 1, 0.5 and 0.1 are also shown as black dotted lines. We see that the signal associated to the A and B terms is much less than for the P δδ , P δθ and P θθ term. In particular, we see that for the scales of interest (0.10 h/Mpc < k < 0.15 h/Mpc) at z = 0 the signal for A and B terms is about S/N = 0.1 and even less at z = 0.5. This means that in Fig. 9, the crucial difference between red and blue lines (for a given PT model) comes from terms with low signal-to-noise ratio. In Fig. 10, when f and A b are allowed to vary, the low signal that the A and B terms have, has to be split to find one more parameter, and then, the accuracy recovering f must decrease. In order to break the degeneracy between f and b from redshift space distortions, the signal-to-noise of the non-linear A and B terms must be as high as possible and it could be optimized by selecting haloes with suitable cuts in mass. The shot noise increases with the mass cut but also does the bias. When optimizing for a measurement of the angle-averaged power spectrum, for a fixed P mm , the signal-tonoise scales as a function of mass cut like b 2 /shot noise. If the shot noise can be approximated as Poissonian then the signal-to-noise can be roughly approximated as b 2 n with n number density of tracers. When considering redshift space distortions, the signal for the A and B terms scales like b 3 and b 4 suggesting a different scaling of the signal-to-noise with mass than for the real-space power spectrum, which is favored by higher bias. These considerations might be useful when optimizing a survey selection of targets, although the bias of observable tracers might not behave as the halo bias especially when multiple galaxies occupy the same halo.

Summary and conclusions
Using a suite of 160 N-body simulations each with a volume of V box = 13.8 (Gpc/h) 3 , we have investigated the accuracy of analytic models in predicting the non-linear power spectrum of matter and dark-matter haloes in real and redshift space. The total simulated volume amounts to 2,200 (Gpc/h) 3 , much larger than the volume surveyed by any forthcoming or planned survey ensuring that statistical errors in the determination of the simulation data points is negligible. We make the dark matter and halo power spectra publicly available and also provide the multipoles measured from these simulations for possible comparisons 7 . We considered a number of theoretical schemes obtained by combining standard or resummed perturbation theory with analytical models for redshift-space distortions (based on the simplification of Eq. 2.7). To predict power spectra in real space, we have employed 1-and 2-loop standard perturbation theory, and the resummed perturbation theory proposed by [24,56] that we have generalized to account for 2-loop correction terms in the resummed propagator. For the redshift-space power spectra, we have focused on the models proposed by Kaiser [26], Scoccimarro [27] and Taruya et al. [28].
At the level of dark matter in real space, increasing the order in loop corrections for the resummed propagator improves the theoretical predictions of the power spectrum. In particular, working at 2loop correction in the resummed propagator, N 2 , an accuracy of 1% is achieved at different redshifts up to the following scales: k = 0.10h/ Mpc at z = 0; k = 0.15h/ Mpc at z = 0.5; k = 0.20h/ Mpc at z = 1.0 and k = 0.25h/ Mpc at z = 1.5 for a 2-loop truncation of the infinite series. In general, working at 2-loop correction in the resummed propagator provides a more accurate description than working at 1-loop correction (as many of the public codes do [50,54,56]).
Also, the price of working at 2-loop instead of 1-loop correction in the resummed propagator is not very high in terms of computational resources. It is true that evaluating N 2 involves the 5-dimensional integration of P 15 , but the angular part of this function can be either analytically computed, or numerically precomputed for any shape of the linear power spectrum so, in the end, one ends up with a 2-dimensional integration,which can be easily performed.
For dark matter in redshift space, our results show that the model by Taruya et al. combined with a Lorentzian damping term for the FoG effect with σ 0 as a free parameter, is able to reproduce the multipoles from N-body simulations with high accuracy. In this paper, we have fit σ 0 as a function of k max using monopole, quadrupole and hexadecapole data separately, which yields to 3 different values for σ 0 . Although is possible to fit all these 3 multipoles with the same value of σ 0 , the accuracy is expected to be reduced. When 3 different σ 0 are used, the high accuracy holds true for the monopole ( 2%) and the quadrupole ( 5%) irrespective of the flavor of perturbation theory adopted to compute the real-space power spectrum. For higher-order moments, such as the hexadecapole, the level of accuracy is more modest, ∼ 10% at z = 1. This suggests that even on mildly non-linear scales at redshifts z 1.5 the accuracy of the modeling of the redshift space distortions is more important than the modeling of the non-linear evolution of the real-space dark matter power spectrum.
The difference between the Taruya et al. model (which attempts to include the density-velocity coupling at higher orders) and the other two models become more evident as the multipole order is increased, possibly suggesting that non-linearities become more important for higher-order multipoles.
While in linear theory only the monopole, quadruple and hexadecapole are non-zero, in principle non-linearities should "excite" all higher-order multipoles, and thus cosmological signal could, in principle, be extracted from them. We find from the N-body simulations that the signal-to-noise decreases with increasing multipole, making the hexadecapole errors large even from the large suite of simulations considered here. This suggests that for cosmological applications most of the signalto-noise at these mildly non-linear scales -where analytic approaches can provide a good modelingis still enclosed in the monopole and quadrupole.
The imprint of the baryon acoustic oscillations, is also visible in the multipoles. We find that all models of redshift space distortions considered do not bias the wiggles location although the more linear models (Kaiser and Scoccimarro) over-predict their amplitude. These considerations might be relevant for recovering in an unbiased way the angular and radial BAO information (separately) from forthcoming surveys.
Overall, the accuracy of the analytic description allows measurements of the logarithmic growth rate f to percent level, when σ 0 is allowed to vary. In this case, when the Taruya et al. model is combined with the RPT prediction for the power spectrum in real space, f is recovered within 1% up to k max = 0.15 h/Mpc for the monopole and quadrupole (separately) at z = 0 and up to k max = 0.20 h/Mpc at z = 1. This indicates that the Taruya et al. model combined with the RPT real-space predictions is accurate enough to be used for precision cosmology.
Since most of the current and future redshift surveys target galaxies as tracers of the matter distribution, a more realistic way of estimating f is to use dark-matter haloes instead of the dark matter density. The limited mass resolution of our N-body simulations, allows us to consider only cluster-sized haloes M > 10 14 M /h. Dealing with isolated haloes has the advantage of eliminating the imprint of FoGs from the power spectrum. However, the effect of the scale-dependent (and possibly non-linear) bias plays an important role. In this work, we have assumed that the bias is linear and that its k-dependence is known. Under these approximations, we have been able to recover f with 5% when the amplitude of the bias is known a priori. For these massive haloes, the effect of bias is important: the degraded accuracy in recovering f indicates that, at least for these massive haloes, the modeling of biasing is crucial. In particular, given the high shot noise that the statistics of these tracers have, a modeling of its behavior both as a function of scale and halo mass is important at mildly non-linear scales.
When both the f parameter and the bias amplitude are allowed to vary, we recover the bias amplitude to 1% level in the best cases, but f is underestimated by 10% − 20% at z = 0 and slightly more at z = 0.5. Similar results are reported in the literature [12]. Most likely this is not (or not only) due to a limitation of the model for the power spectrum, but also to the poor signal-to-noise ratio of the population of haloes used to extract f . It remains to be seen whether reducing the halo mass threshold increases the signal-to-noise. The scaling with halo mass of the signal-to-noise of the RSD terms that can break the degeneracy between f and b, is different from the usual nP ∼ b 2 n used for the angle-averaged, real-space power spectrum. Our simple considerations indicate that a lower mass threshold increases the signal-to-noise, but if the number density is kept fixed, then the signal-to-noise for the A and B terms is favored by a higher bias (the signal-to-noise scales like b 3 n and b 4 n rather than like b 2 n). These considerations might be useful when optimizing a survey selection of targets, although one should keep in mind that bias of observable tracers might not behave as the bias of the host haloes.
The results found in this paper are in agreement with those found in recent works. Kwan et al. [12] found that, most of the RSD models fail at recovering f , underestimating its value even at large scales for z = 0 and z = 0.5. Considering the Taruya et al. model, relaxes the discrepancy (compared to Kaiser or Scoccimarro models), but does not completely fixes the problem. Other works, such as de la Torre & Guzzo [11], show that the Taruya et al. model (in configuration space) + numerical and phenomenological schemes to estimate the real-space spectra, are able to recover f with 5% accuracy from different low-biased galaxy population. Okumura & Jing [14] show that β can be estimated accurately from very massive haloes (M h ≥ 10 14 M /h) using linear theory (Kaiser model) at large scales (k ≤ 0.1h/ Mpc); and f can be also recovered with similar accuracy when the bias is assumed to be known, although the precision achieved is not very high.
Future surveys and missions will provide datasets about the distribution of galaxies on large scales. We envision that in order to extract useful information from these datasets, we will need more accurate theoretical models of structure formation. In this paper, we have shown that for extracting the growth of structure correctly, accurate analytic models for both real and redshift space clustering are crucial. In particular 2L-RPT-N 2 in combination with the Taruya et al. formula, seems to be able to recover f accurately from massive haloes, when the bias is assumed to be known or equivalently when one wants to recover the combination β = f /b. However, in our study we found that this approach fails when trying to recover the individual values of f and b simultaneously. This can be due to the low signal-to-noise ratio of this halo population, or a limitation in the model itself. Should the determination of f and b separately rather than in the β combination from RSD alone become a priority, studies will attempt to improve both the real and redshift-space models. In the case of real space, this can be done extending resummation theories both in Eulerian or Lagrangian spaces. In the case of redshift space, considering higher-order terms in the TNS formula should improve the model. Finally, it is also important to model correctly the stochasticity associated with the halo population and to determine correctly the scale dependence of the bias and its possible non-linearities.
According to standard perturbation theory (SPT) the power spectrum in real space can be expressed as a sum of loop corrections, where P (0) (k) = P lin (k) is the linear term. For Gaussian initial conditions, the different loop corrections read as, where, as already mentioned in the main text, the subscripts i and j refer to the perturbative order of the terms δ(k) used in eq. 2.1 to compute the power spectrum P ij (k). For the case of matter-matter power spectrum, namely P δδ , these terms are [32], × P lin (q 1 )P lin (q 2 )P lin (|k − q 1 |), In the 1-loop correction, P 22 accounts for the mode coupling between vectors with frequencies k − q and q, whereas P 13 can be interpreted as the 1-loop correction to the linear propagator. In a similar way, in the 2-loop correction term, only the second term of P 33 accounts for a full 2-loop mode coupling because is the only term that contains P lin (|k − q 1 − q 2 |). Also note that P 24 contains a term similar to a 1-loop coupling, P lin (|k − q 1 |), is similar to P 22 . P 15 and the first term of P 33 contain no coupling between k and q i and can be interpreted as a 2-loop propagators. In particular, the full n-propagator can be written as [32] P 1n (k) = n!!P lin (k) where x = (n − 1)/2. These similarities between these terms is the basis of the resummation process that is described in Appendix B. The kernels of Eqs. A.4-A.8 are expressed as, with F 1 ≡ 1 and G 1 ≡ 1. Also, k 1 ≡ q 1 + · · · + q m , k 2 ≡ q m+1 + · · · + q n , k ≡ k 1 + k 2 where the functions α and β are defined as, The symmetrization process of the kernels is given by, F s n (q 1 , . . . , q n ) = 1 n! π F n (q π(1) , . . . , q π(n) ), (A.14) where the sum is taken over all the permutations π of the set {1, . . . , n}. In particular, the expressions for F s 2 (k 1 , k 2 ) and G s 2 (k 1 , k 2 ) are, where, cos θ ≡ (k 1 · k 2 )/(k 1 k 2 ). This SPT formalism presents some shortcomings. As noted by [24], at large scales only the linear term contributes to the total power spectrum. However, at smaller scales, all loop corrections become of the same order with a significant cancellation among them. In particular, at low redshifts it can be seen that 1-loop correction overestimates the full power spectrum (from N-body simulations), whereas the 2-loop correction underestimates it. This is due to the fact that P (1) is negative on large scales while P (2) is positive, and both almost cancel out giving a remaining quantity which is close to the full power spectrum. In the same way, as we go to higher order, more cancellations come out among the different loop corrections. Thus, truncating at certain loop in SPT will naturally produce a systematic over-and under-prediction of the real-space power spectrum. A way to avoid this behavior is to resum some terms of the total SPT expansion. In [24,25] a formalism for resumming part of these terms was proposed. The resulting expansion presents a more controlled behavior because each different loop contributes only positively to the total power spectrum and acts at different scales. In Appendix B we present an alternative description (but mathematically identical) of the resummation presented by [24,25]. As an extension of current works, we write not only the 1-loop resummed propagator, but we perform our computation up to 3-loops.

B Resummation in Standard Perturbation Theory
In this section we present the derivation of Eq. 2.11, 2.14 and 2.15. We also show that the N 1 expression is identical to the one used in [56] when the propagator is perturbed at 1-loop. For completeness, we also show the result of perturbing the propagator for 3-loops. These equations come from resumming some terms in SPT under certain approximation in the kernels. In particular, under the Zel'dovich approximation, the kernels read [69], where k = k 1 + · · · + k n . As shown by [24], with this approximation the resummation process yields, P (k, z) = P lin (k, z) + P 1L 22 (k, z) + P 2L 33 (k, z) + · · · + P (n−1)L nn (k, z) + . . . N 0 (k, z) 2 , where, and σ v is a characteristic scale defined as (B.6) We will refer to Eq. B.3 as RPT-N 0 model. With this technique the behavior of PT improves, since every new loop adds a positive term that only acts on a small range of scales. Therefore, using this technique the oscillatory behavior observed in standard PT vanishes.
In this section we show that performing a slightly different approximation (not Zel'dovich) in the kernels we can end up with the same formula used in [56]. Also, depending on how we 'factorize' the kernels, we will end up with a 1-, 2-or higher-loop correction in the resummed propagator. These formulae present a notable improvement respect to Eq. B.3 for a similar computational effort. In this work we do not follow the approach of Feynman diagrams to resum the infinite terms as it is done in [24,25]. Alternatively, we present a different approach for doing this without requiring any knowledge of quantum field theory. We hope that this way of resumming is clearer for the reader who is no familiar with this kind of formalism. Furthermore, this approach allows us to easily compute the resummed propagator for higher-order loops. Our method consists of rewriting the terms of the -loop correction (where ≡ (n + m)/2 − 1), namely P nm , as a sum of subterms which can be associated to lower loop corrections as we show below 8 , P nm (k) = P 0L nm (k) + P 1L nm (k) + P 2L nm (k) + · · · + P L nm (k). (B.7) The subterm with an index 0L contains a P lin (k) and corresponds to the linear power spectrum (noloop correction), the subterm with an index 1L contains a P lin (|k − q 1 |) and therefore is similar to P 22 (k) (that corresponds to 1-loop correction). In the same way, the 2L subterm is similar to P 33 (k) (that corresponds to 2-loop correction) because contains a term P lin (|k − q 1 − q 2 |) and so on. The generic way of writing these terms is the following 9 . The 0-L subterm can be written, × F s m (k, q m 1 , −q m 1 , . . . , q m xm , −q m xm )P lin (k)P lin (q n 1 ) . . . P lin (q n xn )P lin (q m 1 ) . . . P lin (q m xm ), where x i = (i − 1)/2.
The 1-L subterm can be written, 8 The redshift dependence is understood for simplicity: it only appears through P lin . 9 Note that the terms with n odd and m even (and vice versa) vanish for Gaussian initial conditions.

B.1 1-loop factorization
If we want to end up with a resummed propagator of 1-loop correction, the prescription in factorizing the kernels is the following, 1. for 0-L subterms, 3. for 2-L subterms, F s n (q 1 , q 2 , k − q 1 − q 2 , q 3 , −q 3 , . . . , q xn , −q xn ) and similarly for higher-order loops. We factorize the n-order kernel in a product of 2-and 3-order kernels (which correspond to 1-loop correction terms) keeping the sum q 1 + q 2 + · · · = k in all the 2and 3-order kernels. Under these approximations we rewrite the subterms of Eq. B.7.  and similar for higher-order loops. As before, with this approximation we can perform an exact resummation of all terms yielding to, P (k) = P lin (k) + P 1L 22 (k) + P 2L 33 (k) + . . . N 2 (k) 2 , (B. 38) with, We will refer to Eq. B.38 as RPT-N 2 model. Note that if we perform on P 15 (k) in Eq. B.39 the approximation, F s 5 (k, q 1 , −q 1 , q 2 , −q 2 ) 1 5! 3!F s 3 (k, q 1 , −q 1 )3!F s 3 (k, q 2 , −q 2 ), (B. 41) we obtain that 2P 15 → P 13 /P lin 2 P lin and therefore N 2 → N 1 . In the same way, when we apply the Zel'dovich approximation on the kernel of P 13 , F s 3 (k, q, −q) we obtain that 2P 13 /P lin → −k 2 σ 2 v , and therefore, N 1 → N 0 . Note that for the computation of N 2 (k) one needs to compute P 15 (k) which requires the knowledge of F s 5 . This computation is a 6-dimensional integral that reduces trivially to 5-dimensional exploiting rotational invariance. In principle, one could integrate analytically the remaining 3 angles and reduce the computation of P 15 to a 2-dimensional integral in the same way P 13 is reduced to a 1-dimensional integral in Eq. B.31. However this is hard, because the symmetrized kernel F s 5 is the sum of 5! = 120 different cyclic permutations. A possible alternative, is to precompute the angular part of P 15 as a 3-dimensional integral for a wide range values of k, q 1 and q 2 , and then use this to compute P 15 as a 2-dimensional integral for any shape of P lin . Nevertheless, for practical reasons, in this paper P 15 is always computed numerically as a 5-dimensional integral. For completeness we also report the expression for the 3-loop resummed propagator, N 3 (k) function as function of the full-propagator terms P 13 , P 15 and P 17 , where the functions A,B and C are given by, We do not use this function in this paper, because it requires the computation of P 17 which is a 8-dimensional integral (after exploiting rotational invariance), which goes beyond the scope of this paper. We leave the analysis of this function for a future work.