Modelling non-linear redshift-space distortions in the galaxy clustering pattern: systematic errors on the growth rate parameter

We investigate the ability of state-of-the-art redshift-space distortions models for the galaxy anisotropic two-point correlation function \xi(r_p, \pi), to recover precise and unbiased estimates of the linear growth rate of structure f, when applied to catalogues of galaxies characterised by a realistic bias relation. To this aim, we make use of a set of simulated catalogues at z = 0.1 and z = 1 with different luminosity thresholds, obtained by populating dark-matter haloes from a large N-body simulation using halo occupation prescriptions. We examine the most recent developments in redshift-space distortions modelling, which account for non-linearities on both small and intermediate scales produced respectively by randomised motions in virialised structures and non-linear coupling between the density and velocity fields. We consider the possibility of including the linear component of galaxy bias as a free parameter and directly estimate the growth rate of structure f. Results are compared to those obtained using the standard dispersion model, over different ranges of scales.We find that the model of Taruya et al. (2010), the most sophisticated one considered in this analysis, provides in general the most unbiased estimates of the growth rate of structure, with systematic errors within 4% over a wide range of galaxy populations spanning luminosities between L>L* and L>3L*. The scale-dependence of galaxy bias plays a role on recovering unbiased estimates of f when fitting quasi non-linear scales. Its effect is particularly severe for most luminous galaxies, for which systematic effects in the modelling might be more difficult to mitigate and have to be further investigated. [...]


INTRODUCTION
The structure in the Universe grows through the competing effects of universal expansion and gravitational instability. For this reason, the large-scale spatial distribution and dynamics of galaxies, which follows in some way those of mass, provide fundamental information about the expansion history and the nature of gravity. In general, galaxies are biased tracers of the underlying mass distribution. However, they are sensitive to the same gravitational potential and their motions keep an imprint of the rate of structure growth. One manifestation of this is the observed anisotropy between the clustering of galaxies along the line-of-sight and that perpendicular to it in redshift space. These anisotropies or distortions are caused by the line-of-sight component of galaxy peculiar velocities affecting E-mail: sdlt@roe.ac.uk † Scottish Universities Physics Alliance the observed galaxy redshifts from which distances are measured. In turn, the large-scale coherent component of galaxy peculiar motions is the fingerprint of the growth rate of structure.
By mapping the large-scale structure over scales which retain this primordial information, galaxy spectroscopic surveys have become one of the most powerful probes of the cosmological model. A specific application of spectroscopic surveys involves recovering cosmological information on the expansion history H(z), by measuring the shape of the power spectrum (e.g. Tegmark et al. 2004;Cole et al. 2005;Percival et al. 2007;Reid et al. 2010) and tracking the baryonic acoustic oscillations (BAO) feature in the power spectrum or in the two-point correlation function at different redshifts (e.g Eisenstein et al. 2005;Percival et al. 2010;Kazin et al. 2010; Blake et al. 2011b). However, measurements of H(z) alone, either from BAO or Type Ia supernovae, cannot discriminate dark energy from modifications of General Relativity (e.g. Carroll et al. 2004), in order to explain the observed recent accel-eration of the expansion of the Universe (Riess et al. 1998;Perlmutter et al. 1999). This degeneracy can be lifted by measuring the growth rate at different epochs (Peacock et al. 2006;Albrecht et al. 2006). Indeed, scenarios with similar expansion history but different gravity or type of dark energy will have a different rate of structure growth resulting from different effective gravity strength in action. This makes redshift-space distortions measured from large spectroscopic surveys a very efficient probe to test cosmology, at the same level as BAO and cosmological microwave background (CMB) anisotropies. In fact, although this effect is known since the late eighties (Kaiser 1987), its usefulness as a probe of dark energy and modified gravity has been realised only recently (Guzzo et al. 2008).
Measuring the growth rate of structure from redshift-space distortions is however non trivial. The linear theory formalism for the power spectrum was first derived by Kaiser (1987) (see Hamilton 1992, for its configuration-space counterpart). Its validity is however limited to very large scales as it lacks a description of small-scale non-linear fluctuations. This model has been extended to quasi-and non-linear scales in the early nineties using the earlier ideas of the "streaming model" (Peebles 1980), in which the linear correlation function is convolved along the line-of-sight with a pairwise velocity distribution Peacock & Dodds 1994). This enables one to approximately reproduce the Fingers-of-God small-scale elongation (Jackson 1972). Fitting functions calibrated on simulations have also been proposed for this purpose (Hatton & Cole 1999;Tinker et al. 2006). Such extension of the linear model, usually referred as the "dispersion model", has been extensively used to measure the growth rate of structure f or the distortion parameter β = f /bL from redshift surveys, using measurements of both redshift-space correlation function (Peacock et al. 2001;Hawkins et al. 2003;Zehavi et al. 2005;Ross et al. 2007;Okumura et al. 2008;Guzzo et al. 2008;daÂngela et al. 2008;Cabré & Gaztañaga 2009a,b;Samushia et al. 2011b) and power spectrum (Percival et al. 2004;Tegmark et al. 2004Tegmark et al. , 2006Blake et al. 2011a). We refer the reader to Hamilton (1998) for a review of older studies. Although generally the dispersion model is found to be a good fit on linear and quasi-linear scales (Percival & White 2009;Blake et al. 2011a), it breaks down in the non-linear regime (Taruya et al. 2010;Okumura & Jing 2011). In particular, it has been shown that it introduces systematic errors of about 10 − 15% on the growth rate parameter (e.g. Taruya et al. 2010;Okumura & Jing 2011;Bianchi et al. 2012), of the order of or greater than the statistical errors expected from on-going and prospected very large spectroscopic surveys such as WiggleZ (Drinkwater et al. 2010), GAMA (Driver et al. 2011), VIPERS , BOSS , Euclid (Laureijs et al. 2011), or BigBOSS (Schlegel et al. 2011. In particular, Euclid, the recently selected ESA dark energy mission, should be able to constrain the growth rate at the percent level (Wang et al. 2010;Samushia et al. 2011a;Majerotto et al. 2012). There is therefore a strong need to go beyond the dispersion model, in order to bring systematic errors below this expected level of precision, i.e. measure the growth rate parameter in an unbiased way.
Work in this direction started since quite some time, concentrating first on describing the redshift-space clustering and dynamics of dark matter. Scoccimarro (2004) demonstrated that the dispersion model gives rise to unphysical distributions of pairwise velocities and proposed an ansatz that accounts, to some extent, for the non-linear coupling between the velocity and the density fields. This latter model has been shown to provide a better match to the observed redshift-space power spectrum in dark matter simulations (Scoccimarro 2004;Jennings et al. 2011) and has later on been refined by Taruya et al. (2010). Further approaches have also been proposed while completing this paper Reid & White 2011), but we shall not discuss them in this analysis.
All mentioned advanced redshift-space distortions models have been tested so far only on the redshift-space power spectra of dark matter and dark matter haloes, as extracted from large Nbody simulations (Kwan et al. 2011;Reid & White 2011;Nishimichi & Taruya 2011). This is quite different from a real survey, in which the most useful tracers of mass, the galaxies, are in general biased with respect to the underlying density field through a bias which is generally non-linear and scaledependent. The performance of redshift-space distortions models applied to galaxy populations with a priori unknown biases, has to be precisely investigated. This is the aim of this paper, in which we confront non-linear models of redshift-space distortions for the anisotropic two-point correlation function in the case of realistic galaxy samples. This is done in the framework of the concordant ΛCDM cosmological model. We study the ability of these models to recover the linear growth rate of structure and their range of applicability. Furthermore we investigate the effects of galaxy non-linear spatial and velocity biases and quantify how the latter affect the estimated linear growth rate for differently biased galaxy populations.
The paper is organised as follows. In Sect. 2, we present the redshift-space distortions formalism and how models can be implemented in practice. In Sect. 3, we present the comparison between the different models and study the effect of galaxy non-linear bias. In Sect. 4, we investigate the impact of neglecting galaxy velocity bias in the modelling. In Sect. 5, we summarise our results and conclude.

Fourier space
The peculiar velocity v alters objects apparent comoving position s from their true comoving position r, as where H(a) is the Hubble parameter, a is the scale factor, andê is the line-of-sight unit vector. The redshift-space density field δ s (s) can be obtained from the real-space one by requiring mass conservation, i.e.
In the following we shall work in the plane-parallel approximation and in this limit, the Jocobian of the real-to redshift-space transformation can be simply written as, where we defined u (r) = −v (r)/(f aH(a)) with f being the linear growth rate. The linear growth rate parameter is defined as the logarithmic derivative of the linear growth factor D(a) and given by f (a) = d ln D/d ln a. To a very good approximation it has a generic form (Wang & Steinhardt 1998;Linder 2005), where Ωm(a) = Ωm,0 a 3 In this parametrisation, while Ωm characterises the mass content in the Universe, the exponent γ directly relates to the theory of gravity (e.g. Linder 2004). General Relativity scenarios have γ 0.55. From Eq. 2 and Eq. 3, one can write the redshift-space density field as, One usually assumes an irrotational velocity field for which u (r) = ∂ ∆ −1 θ(r) and where θ(r) = ∇ · v(r) is the velocity divergence field and ∆ denotes the Laplacian. In that case Eq. 6 can be recast, In Fourier space, it is noticeable that ∂ 2 ∆ −1 = (k /k) 2 = µ 2 with µ being the cosine of the angle between the line-of-sight and the separation vector. Therefore, one can write the redshift-space density field (Scoccimarro et al. 1999) as, and the redshift-space power spectrum as, where in the latter equation, ∆u = u (x) − u (x ) and r = x − x . The redshift-space power spectrum given in Eq. 10 is almost exact, the only approximation which has been done is to assume that all object line-of-sight separations are parallel. This approximation is valid for samples with pairs covering angles typically lower than 10 • (Matsubara 2000). Eq. 10 captures all the different regimes of distortions. While the terms in square brackets describe the squashing effect or "Kaiser effect" which leads to an enhancement of clustering on large scales due to the coherent infall of mass towards overdensities, the exponential prefactor is responsible to some extent for the Fingers-of-God effect (FoG, Jackson 1972) which disperses objects along the line-of-sight due to random motions in virialised structures. Scoccimarro (2004) proposed a simple ansatz for the redshift-space anisotropic power spectrum by making the assumption that the exponential prefactor and the term involving the density and velocity fields can be separated in the ensemble average. In that case Eq. 10 simplifies to, where P δδ , P δθ , P θθ are respectively the non-linear mass densitydensity, density-velocity divergence, and velocity divergencevelocity divergence power spectra and σv is the pairwise velocity dispersion defined as, It is found that this model captures most of the distortion features predicted by N-body simulations (Scoccimarro 2004; Jennings et al. 2011) although it breaks down in the non-linear regime (Percival & White 2009;Taruya et al. 2010). Note that in the linear regime where P δδ = P δθ = P θθ = P and in the limit where kσv tends to zero, one recovers the original Kaiser (1987) formula, derived from linear-order calculations.
In principle, the exponential prefactor and the term involving the density and velocity fields in Eq. 10, which we will refer to as the damping and Kaiser terms in the following, cannot be treated separately. Additional terms may arise in Eq. 11 from the coupling between the exponential prefactor and the velocity divergence and density fields. Taruya et al. (2010) proposed an improved model that takes into account these couplings, adding two correction terms CA and CB to Scoccimarro (2004)'s formula such as, whose perturbative expressions are given in their appendix A. In the improved model, the exponential prefactor has been replaced by an arbitrary functional form D(kµσv) for which σv is an effective pairwise velocity dispersion parameter that can be fitted for. Taruya et al. (2010) showed that while adopting a Gaussian or a Lorentzian for the damping function and letting σv free, one improves dramatically the fit to the redshift-space power spectrum in large dark matter simulations, particularly on translinear scales. The function D(kµσv) damps the power spectra in the Kaiser term but also partially mimics the effects of the pairwise velocity distribution (PVD) in virialised systems, which translate into the FoG seen in the anisotropic power spectrum and correlation function on small scales. This is analogous to the phenomenological dispersion model proposed in the early nineties (e.g. Fisher et al. 1994;Peacock & Dodds 1994) in which the linear Kaiser model in configuration space (Hamilton 1992) is radially convolved with a PVD model to reproduce the FoG elongation on small scales, as for the early streaming model (Peebles 1980).
There is however not any simple general functional form for the PVD that matches all scales for all types of tracers. The shape of the PVD is found to depend on galaxy physical properties and halo occupation (Li et al. 2006;Tinker et al. 2006), and its associated pairwise velocity dispersion to vary with scale, in particular at small separations (e.g. Hawkins et al. 2003;Cabré & Gaztañaga 2009b). It can be shown mathematically that the PVD is in fact not a single function but rather an infinite number of PVD corresponding to different scales and angles between velocities and separation vectors (Scoccimarro 2004). In practice however, the use of an exponential distribution, a Gaussian or other forms with more degrees of freedom (e.g. Tang et al. 2011;Kwan et al. 2011) shows to be very useful to fit the residual small-scale distortions remaining once the large-scale Kaiser distortions are accounted for, unless one is interested in modelling the very small highly non-linear scales.

Configuration space
The redshift-space anisotropic two-point correlation function can be obtained by Fourier-transforming the anisotropic redshift-space power spectrum as, where ν = r /s, r ⊥ = s 2 − r 2 , and L l denote Legendre polynomials. The correlation function multipole moments ξ s l (s) are de-fined as, where j l denotes the spherical Bessel functions and In practice, it can be convenient to write the redshift-space twopoint correlation function as a convolution between the Fourier transform of the damping function D and that of the Kaiser term as, In the case of the model of Eq. 11, the three non-null correlation function multipole moments related to the Kaiser term are given by, where ξ δδ , ξ δθ , ξ θθ are the Fourier conjugate pairs of P δδ , P δθ , P θθ and (Hamilton 1992;Cole et al. 1994), The correlation function multipole moments of the Kaiser term in the case of the model of Eq. 14 are given in Appendix A. For the latter model we restrain ourselves to use only the first three nonnull multipole moments, as those of orders l = 6 and l = 8 are very poorly defined in our simulated galaxy catalogues.

From mass to galaxies
The models derived in the previous section apply in the case of perfectly unbiased tracers of mass. Real galaxies however are biased with respect to mass. Galaxy biasing is generally expected to be non-linear, scale-dependent, stochastic, and to depend on galaxy type, although it is still poorly constrained by observations. On large scales in the linear regime, one expects the bias to be a constant multiplicative factor to the mass density field as δg = bLδ. In that case, it is convenient to replace the growth rate f in the models with a "effective" growth rate (or distortion parameter) β = f /bL, which accounts for the large-scale linear bias bL of the considered galaxies. This simple model is valid on large scales where the bias asymptotes to a constant value but breaks down on small non-linear scales, where bias possibly varies with scale. Recently, Okumura & Jing (2011) showed that the scale-dependent behaviour of halo bias can strongly affect the recovery of the growth rate. While some analytical approaches have been proposed to include bias non-linearity in the model (Desjacques & Sheth 2010;Matsubara 2011), here we follow a different route and assume that the scale dependence of bias is known. In fact, the latter can be measured to some extent from the data themselves in configuration space, once the shape for the underlying non-linear mass power spectrum is assumed. General arguments may suggest that galaxy motions are also biased with respect to the mass velocity field, while observations tend to indicate that this bias is small (Tinker et al. 2006;Skibba et al. 2011). In this analysis we will neglect the galaxy velocity bias in the models but discuss and quantify its impact on the recovery of f in Section 4.

Constructing the galaxy redshift-space distortion models
We will use in this analysis different combinations of Kaiser terms, damping functions, and bias prescriptions. Although we will work in configuration space, we refer to the different models in this section as their Fourier-space counterpart for clarity. All the models we consider take the general form, where, Hereafter, we will refer as the different PK models to A, B, and C. Model A corresponds to the Kaiser (1987) model with the non-linear power spectrum instead of the linear one. It assumes a linear coupling between the density and velocity fields such that δ ∝ θ. Model B is the generalisation proposed by Scoccimarro (2004) that accounts for the non-linear coupling between the density and velocity fields, making explicitly appearing the velocity divergence auto-power spectrum and density-velocity divergence cross-power spectrum. Finally, model C is an extension of model B that contains the two additional correction terms proposed by Taruya et al. (2010) to correctly account for the coupling between the Kaiser and damping terms. Besides, we will consider two deterministic galaxy biasing prescriptions: a constant linear bias b(k) = bL and a general non-linear bias which we define as b(k) = (Pgg/P δδ ) 1/2 (k) = bLbNL(k), where Pgg is the galaxy power spectrum and bNL(k) is the scale-dependent part of the bias that tends to unity at small k. We note that the Gaussian and Lorentzian damping forms that we will consider here, correspond respectively to Gaussian and exponential functions in configuration-space.
The redshift-space distortions models necessitate P δδ , P δθ ,      and P θθ real-space power spectra as input. These can be obtained analytically using perturbation theory. Although standard perturbation theory does not describe well the shape of these power spectra on intermediate and non-linear scales, improved treatments such as Renormalised Perturbation Theory (RPT, Crocce & Scoccimarro 2006) or Closure Theory (Taruya et al. 2009) have been shown to be much more accurate (see Carlson et al. 2009, for a thorough comparison). In particular, Closure Theory predictions are found to match large N-body simulation real-space power spectra to the percent-level up to k = 0.2 for z > 0.5 (Taruya et al. 2009).
In this analysis we use the P δδ provided by CosmicEmu emulator (Lawrence et al. 2010) and the fitting functions of Jennings et al. (2011) to obtain P θθ and P δθ from P δδ . The latter fitting functions have an accuracy of 5% to k = 0.2 for both standard ΛCDM and quintessence dark energy cosmological models. In Fig. 1 and 2 we compare the P δδ , P δθ , P θθ obtained in this way with Closure Theory 2-loop analytical predictions at z = 0 and z = 1. We find that all power spectra agree very well below k 0.2 and k 0.3 respectively for the two redshifts considered, except in the case of P δθ for which they systematically differ by about 10%. While all other power spectra match on linear scales, the P δθ fitting formula from Jennings et al. (2011) stays somewhat below (dotted lines in the figures). We find that by multiplying the latter by a factor of 1.1 one obtains a good match with Closure Theory predictions on both linear and non-linear scales (solid lines in the figures). We will then adopt this correcting factor in the following when calculating the redshift-space distortions models.
It is noticeable that Closure Theory breaks down at lower k than Jennings et al. (2011) fitting functions, indicating that the latter are more suitable in practice to describe clustering on the smallest scales. It is however important to mention that the validity of these fitting functions is limited at large k (kmax 0.2, Jennings et al. 2011). One can therefore reliably describe the models involving P δθ and P θθ only down to scales of r π/kmax 10.5 h −1 Mpc. On smaller scales, P δθ and P θθ configuration-space counterparts, ξ δθ and ξ θθ , will drop rapidly and their contributions to the redshift-space distortions models will be underestimated with respect to that of ξ δδ .

Methodology
To test the redshift-space distortions models presented in section 2.4 and quantify how well the linear growth rate parameter f can be recovered from the anisotropic two-point correlation function ξ(r ⊥ , r ), we constructed a set of realistic galaxy catalogues. We populate the identified friends-of-friends haloes in the MultiDark Run 1 (MDR1) dark matter N-body simulation (Prada et al. 2011) with galaxies by specifying the Halo Occupation Distribution (HOD). MDR1 assumes a ΛCDM cosmology with (Ωm = Ω dm + Ω b , ΩΛ, Ω b , h, n, σ8) = (0.27, 0.73, 0.0469, 0.7, 0.95, 0.82) and probe a cubic volume of 1 h −3 Gpc 3 with a mass resolution of mp = 8.721 × 10 9 h −1 M . From haloes at snapshots z = 0.1 and z = 1, we built galaxy catalogues based on the current most accurate observations of the halo occupation available at these redshifts (Zheng et al. 2007;Zehavi et al. 2011;Coupon et al. 2011) and create three luminosity-threshold samples corresponding to L > L * , L > 2L * , and L > 3L * . In these catalogues, the redshift-space displacements with respect to real space were reproduced using Eq. 1 and the galaxy peculiar velocity information. The details of the procedure used to create the galaxy catalogues are given in Appendix B.
We measured the anisotropic two-point auto-correlation function ξ(r ⊥ , r ) in the redshift-space catalogues using the Landy & Szalay (1993) estimator. Because of the large number of galaxies and to keep the computational time reasonable, all pair counts have been performed using a specifically developed parallel kd-tree code following the dual-tree approach (Moore et al. 2001).
The best-fitting parameters for the different models have been determined by adopting the usual likelihood function, where Np is the number of points in the fit, ∆ is the data-model difference vector, and C is the covariance matrix of the data. The likelihood is performed on the quantity y = ln 1 + ξ(r ⊥ , r ) , rather than simply ξ(r ⊥ , r ), as to enhance the weight on large more linear scales (see Bianchi et al. 2012, for discussion). The determination of the covariance matrix is however troublesome when fitting two-dimensional correlation functions. Having only one realisation of the simulated catalogues at each redshift, we can use to this end internal estimators, such as blockwise bootstrap or jackknife resampling (e.g. Norberg et al. 2009). Using the latter method, we find that the maximum number of cubic sub-volumes that can be extracted without underestimating the variances on the scales of interest is 64. This poses a problem, since in order to have a proper estimation of the eigenvalues of the covariance matrix, this number should be at least equal to the number of degrees of freedom, which in our case ranges between 14397 and 25597 (i.e. 120 2 to 160 2 data points minus 3 free parameters) depending on the scale interval of the fit. As a result, the covariance matrix estimated in this way is degenerate and in the end not very useful. We note that this is a general problem for redshift-space distortions analysis, when one tries to fit the full anisotropic two-point correlation function or power spectrum. In principle this could be solved by using a very large number of realisations or alternatively, theoretically-motivated analytical forms for the covariance matrix. In our case, due to the limited size of the simulation, it is impracticable to define more sub-volumes than degrees of freedom, as this would result in underestimating the covariances by using too small sub-volumes. We note however that by estimating ξ(r ⊥ , r ) using bins of size 0.5h −1 Mpc in both directions in the jackknife resamplings, we are oversampling the functions so the actual number of degrees of freedom may be smaller than that associated to the number of data points (e.g. Fisher et al. 1994). We are therefore forced to perform our fits ignoring the non-diagonal elements of the covariance, i.e. use the variances only. We verified however on a test case that the best-fit values of the parameters obtained in this way do not differ significantly when using the full covariance matrix based on 64 sub-volumes.
In all cases we define f , σv, and bL as free parameters and use different scale ranges in the fit by varying r ⊥ from r min ⊥ = 1h −1 Mpc to r min ⊥ = 20h −1 Mpc and fixing r max ⊥ = r max = 80h −1 Mpc. We limit ourselves to these latter maximum scales to avoid being affected by wide-angle effects that may introduce additional distortions at larger scales (e.g. Samushia et al. 2011b). The statistical errors on the model parameters, and in particular on f , have been estimated from the 1σ dispersion of their best-fitted values among the 64 resamplings. We assume the amplitude of the input power spectra σ8 in the models to be known and fix it to that of the simulation. The fiducial values of the growth rate f f id are given by Eq. 4.

Varying the Kaiser and damping terms
Let us first study the effect of using different combinations of Kaiser term and damping function with the assumption that galaxy bias is linear and scale-independent. We use the full galaxy catalogues at z = 1 and z = 0.1 and estimate the statistical and systematic errors on the growth rate with the different models, varying

Statistical and systematic errors
The systematic errors on the growth rate at z = 1 are significant using any of the model variants, and depends on the chosen minimum scale for the fit, r min ⊥ , as shown in Fig. 3. Over the explored range of r min ⊥ , ∆f /f varies between −8% and +6% in the case of models A and B. The expected most accurate model, i.e. model C with exponential damping (C-EXP), shows a similar trend with a systematic error varying between −6% and +10% and between −21% and +9% when using a Gaussian damping (C-GAUSS). We note that model A with exponential damping (A-EXP) applied to scales r min ⊥ < 5h −1 Mpc, which is one of the most commonly used model in the literature, systematically underestimates the value of f , in agreement with recent analyses (e.g. Bianchi et al. 2012).
When looking at the z = 0.1 case in Fig. 4, model performances and associated scale-dependence of systematic errors change significantly. All models now tend to systematically underestimate the growth rate by about 10% for r min ⊥ > 5h −1 Mpc. The notable exception is model C, and more particularly model C-EXP, which performs very well and recovers the fiducial value of f with less than 3% of systematic error, even including scales as small as r min ⊥ = 5h −1 Mpc. This is consistent with the power spectrum analysis of Kwan et al. (2011), who show that for dark matter only at z = 0, C-GAUSS is the least biased model, with a systematic error of −4% when fitting up to kmax = 0.1. Our tests show however that for galaxies, model C-EXP is less biased than C-GAUSS. In fact the choice of damping function has only a significant impact on model's ability to handle small scales, with the difference diminishing with increasing r min ⊥ given the similar asymptotic behaviour of the two functional forms. Conversely, we note that the Gaussian damping produces in general slightly lower statistical errors than the exponential damping. These tend also to be about 15% smaller for models A and B than for model C.
The important result of this section is the clear difference in the systematic error on f produced by the very same models when applied to samples at different redshifts. In the next section we shall argue that this can be simply understood in terms of the evolution (weakening) of the scale-dependence of bias between z = 1 and z = 0.1 for L > L * galaxies.

Fidelity in reproducing the anisotropic two-point correlation function
We visually compare in Fig. 5 and Fig. 6 the model correlation functions at z = 1 and z = 0.1 using (a) the fiducial values of the parameters (f , σv, bL) and (b) their best-fit values obtained with the chosen model, to those measured in our simulated catalogues. For case (a), σv is set to linear theory prediction which is given by Eq. 12, replacing P θθ by the linear power spectrum. We limit this comparison to the case with exponential damping and r min ⊥ = 5h −1 Mpc. The galaxy bias is assumed to be linear and scale-independent.
As shown in the top panels, when all parameters are fixed to their fiducial value, model C (solid line) provides a very good description of the observed ξ(r ⊥ , r ) in the catalogues. Conversely, model A and B produce contours that at intermediate separations (r ⊥ < 40h −1 Mpc, r < 25 − 30h −1 Mpc) are less squashed along the line-of-sight. When (f , σv, bL) are allowed to vary (lower panels), all models are generally capable to achieve a good fit to ξ(r ⊥ , r ) above r ⊥ = 10h −1 Mpc. For model A and B, this is obtained through a lower damping than predicted by lin- Best-fit (models w/ exp. D(σ) + L-bias) r ⊥ min =5 h -1 Mpc Figure 5. Measured ξ(r ⊥ , r ) and associated models for L > L * galaxies at z = 1. In each panel the dotted, dashed, and solid curves correspond respectively to model A, B, and C with exponential damping and linear bias, while the contours correspond to the measured ξ(r ⊥ , r ) in the galaxy catalogue. The top panel shows the best-fitting model when (f ,σv,b L ) parameters are allowed to vary, while the bottom panel shows the fiducial prediction of the models. We note that in the latter case σv is fixed to its linear value. In this figure, the measured ξ(r ⊥ , r ) is smoothed using a Gaussian kernel of size 0.5h −1 Mpc. ear theory. This allows one to reproduce the significant squashing along the line-of-sight seen in the catalogues, but fails to model the FoG elongation at smaller r ⊥ . It should also be noted that at r ⊥ < 40h −1 Mpc and r > 25 − 30h −1 Mpc, model A predicts a much higher amplitude of the correlation function than models B and C. In general, low values of σv can balance the deficit of power on small scales but do not allow to recover the true value of f . This is shown in Fig. 7 and 8, where the recovered values of σv for the different models are compared to linear prediction. One notices immediately that the trends of the best-fitting values of σv as a function of r min ⊥ are correlated with those found for f in the different models: higher values of estimated f imply higher values of estimated σv. These results evidence the presence of a significant correlation among the parameters, in particular between f and σv. In this respect, model C tends to be less degenerate than A and B. Moreover, model C is able to recover more realistic values of σv, of the order of linear theory predictions. Conversely, model A provides best-fitting values strongly deviating from linear expec-  tations, leading to σv as small as 100 − 200 km · s −1 . These results are consistent with the findings of Taruya et al. (2010) and Nishimichi & Taruya (2011), who compared model A and C in the case of dark matter and dark matter halo catalogues, finding a better agreement of the recovered σv with linear theory in the case of model C.

Effect of galaxy scale-dependent bias
We now let the galaxy bias vary with scale in the models and study the effect on the recovery of the growth rate parameter. In general, the galaxy bias in configuration space can be defined as, where ξgg is the galaxy real-space auto-correlation function and bNL(r) is the non-linear scale-dependent part of the bias. It is important to stress that ξgg(r) is directly measurable from observations by deprojecting the observed projected correlation function w(r ⊥ ) (Saunders et al. 1992) as, where wi = w(ri) in the logarithmic bin centred on r ⊥,i . This procedure allows one to correctly recover the shape of ξgg(r) up to about 30h −1 Mpc (e.g. Saunders et al. 1992;Cabré & Gaztañaga 2009b). In the following we will therefore make the assumption that ξgg(r) is known and use the measured real-space ξgg(r) from the simulated catalogues to infer bNL(r) to be used in the models. In fact, it is not necessary to know the exact shape of ξgg(r) on scales larger than about 20 − 30 h −1 Mpc, where one generally finds the galaxy bias to be almost scale-independent and can thus safely assume bNL(r) = 1. A notable exception is that of more non-linear objects, for which the scale dependence may extend to larger scales (see section 3.3.1). Fig. 9 shows the non-linear scale-dependent component of galaxy bias, bNL(r), for the different galaxy populations in our simulated catalogues at the two reference redshifts considered, z = 1 and z = 0.1. In the previous section we considered only catalogues of galaxies with L > L * , while in this figure we introduce more extreme galaxy populations, which we analyse in the following section. To define bNL(r), the linear bias bL has been determined for each galaxy population by minimising the difference between ξgg and b 2 L ξ δδ on scales above r = 10h −1 Mpc. It is evident from this figure that non-linearities in the galaxy bias produce variations up to 40% in the real-space clustering on scales 1h −1 Mpc < r < 20h −1 Mpc, the strength of the effect increasing for more luminous galaxies. Let us come back to our original L > L * catalogues and repeat the analysis of the previous section now including the scale dependence of galaxy bias shown in Fig. 9. The new statistical and systematic errors on f estimated from our simulated catalogues are shown in Figs. 10 and 11. In general, one sees that including the bias scale-dependence information improves the estimation of f , in particular when using the smallest scales in the fit (r min ⊥ < 5h −1 Mpc), as expected. This is particularly true at z = 1 although at z = 0.1 the effect is almost negligible. This is clearly related to the scale dependence of the bias being stronger at earlier epochs, as visible in Fig. 9. This is a very important result, as it implies that neglecting bias scale-dependence in analysis of deep redshift surveys, would produce estimates of f which are differently affected by systematic errors at different redshifts. Accounting for this dependence on scale tends also to reduce the dependence of the systematic error on the minimum fitted scale: the retrieved value is almost constant down to r min ⊥ = 1h −1 Mpc for all considered models. Finally, it is important to note that model C-EXP provides the best performance at both considered redshifts, recovering the fiducial value of the growth rate with less than 3% systematics while using clustering information down to r ⊥ = 5h −1 Mpc. Remarkably, including the bias scale-dependence also reduces the statistical error on f by about 15% for all models.

The case of highly-biased galaxies
One might naively expect highly biased objects to be more suitable to constrain f from observations, by virtue of their higher clustering signal that should in principle reduce statistical uncertainty (see Bianchi et al. 2012). However, highly biased objects, which are more likely to reside in the most massive haloes, have undergone a stronger non-linear evolution, this explaining their stronger bias scale-dependence seen in Fig. 9. As such, the inclusion of bNL in the models may become even more critical for these objects, if one's goal is to accurately measure the growth rate parameter (e.g. Cabré & Gaztañaga 2009a). In this section, we explicitly test this hypothesis, extending the model comparison to the regime of highly biased galaxy populations. Let us first define higher-luminosity galaxies from our simulated catalogues. We shall consider two sub-samples, defined as having galaxies with L > 2L * and L > 3L * . We complement these with a catalogue of Luminous Red Galaxies (LRG) drawn from the "LasDamas" suite of simulations, meant to accurately reproduce the galaxy clustering in the SDSS-DR7 release (McBride et al., in preparation). More precisely, we shall use 100 mock realisations of what are defined as "faint LRG" (Mr < −21.2).
We compare in Fig. 12 the relative systematic error on the growth rate obtained for these highly biased galaxies. We use here only model C-EXP and include the scale-dependence of bias as described in Section 3.3. We find that the trends are similar to those observed for L > L * galaxies, with a dependence on r min ⊥ that varies with the considered redshift. Overall, the systematic error remains within ±5% for r min ⊥ = 5 − 10h −1 Mpc for all considered galaxy populations, but the LRG sample. In this case the growth rate is overestimated by on average about 10%, suggesting additional non-linear effects that are not accounted for in our models. A cautionary remark in interpreting this result is that the simulated LRG sample is the only one which has not been built by us on purpose for this work. Although this sample has been already used for other investigations, we have no way to verify the details of the HOD implementation and its impact on our results. Assuming the effect is indeed real, one may hope however to reduce the systematic error by including in the fit scales greater than our current upper limit of r ⊥ = 80h −1 Mpc (Samushia et al. 2011a), at the expense of a worse precision in the measurement.
These results suggest that for highly biased galaxies, additional systematic effects arise which could be due to an incorrect inclusion of scale-dependence of bias or a too simplistic approach for the FoG. From a theoretical point-of-view, bias non-linearity changes the correction terms CA and CB in model C, which might need to be modified to properly include galaxy bias non-linearities and scale-dependence (see Tang et al. 2011;Nishimichi & Taruya 2011). We plan to investigate in more details these aspects in a future paper. Finally, we note that FoG modelling could be potentially improved by adding more freedom in the damping function (e.g. Kwan et al. 2011) or including scale dependence in the pairwise velocity dispersion (e.g. Hawkins et al. 2003), but at the price of increasing the statistical error on the growth rate estimate.

EFFECT OF GALAXY VELOCITY BIAS
In the framework of understanding the impact of non-linear effects on the accuracy of growth rate estimates from redshift-space distortions, we investigate in this section the possible impact of velocity bias, an effect which is is usually neglected. The galaxy catalogues that we used so far, assume that the radial distribution of satellite galaxies inside dark-matter halos follows that of mass, as described by a Navarro et al. (1996) (NFW) radial density profile. Moreover, central galaxies have been defined as being at rest at the centre of their dark matter halo, inheriting its mean velocity. These assumptions make galaxy velocities unbiased with respect to the mass velocity field. However, there are some observational evidences that galaxies does not exactly follow the same radial distribution as dark matter and exhibit some velocity bias. In particular, recent smallscale clustering measurements in the SDSS suggest that relatively luminous galaxies have a steeper radial density profile than NFW, with inner slopes close to −2 and lower concentration parameters (Watson et al. 2011). In addition, galaxy groups and galaxy clusters analysis tend to indicate that central galaxies might in general not be exactly at rest at the centre of the potential well of dark matter haloes (e.g. van den Bosch et al. 2005;Skibba et al. 2011). Differences between the radial distribution of galaxies and that of mass implies the presence of additional spatial and velocity biases. This has a direct impact on the description of the observed redshift-space distortions. In this section we provide a first quantitative assessment of the systematic uncertainty on f that it can introduce if not accounted for in the models.
To this end, we include in the catalogues some amount of velocity bias coming from either central or satellite galaxies. For central galaxies we follow van den Bosch et al. (2005) and assign halocentric positions assuming a radial density distribution of the form, where m is the halo mass, rv is the halo virial radius, and fr is a free parameter that controls the amount of spatial and velocity bias introduced. This effectively offsets central galaxies from their halo centre of mass. By solving Jeans equation one obtains the associated one-dimensional velocity dispersion, where ψ(r) is gravitational potential. One can thus define the velocity bias of central galaxies as, where the halo-averaged velocity dispersions are given by, σcen|m = 4π m rv (m) 0 ρcen(r|m)σcen(r|m)r 2 dr (31) (32) and the NFW radial density profile is defined as, where c dm is the dark matter concentration parameter (see Appendix B for its definition). Galaxy halo-centric velocities are drawn from a Gaussian distribution in each dimension with velocity dispersion given by Eq. 29. We choose the value of fr in order to have b cen v (m) = 0.2. The latter is an average value motivated by recent observations (e.g. Coziol et al. 2009;Skibba et al. 2011).
In the case of satellite galaxies we followed a similar methodology. Based on the recent results of Watson et al. (2011) in the SDSS for Mr < −20.5 galaxies, we radially distribute satellite galaxies as to reproduce a generalised radial density profile of the form, with csat(m) = c dm (m) 2 and γ = 2. As for central galaxies, we assign satellite galaxy velocities from the associated one-dimensional velocity dispersion, This leads to a velocity bias, where σsat|m = 4π m rv (m) 0 ρsat(r|m)σsat(r|m)r 2 dr.
In that case, b sat v (m) slowly increases from b sat v = 1 at m = 10 9 h −1 M to b sat v 1.2 at m = 10 16 h −1 M . Assuming that the haloes are spherical and isotropic, all previous integrals can be solved analytically (see e.g. Appendix B and van den Bosch et al. 2005). Fig. 13 and 14 show the percent variation on the estimated f , when a velocity bias is included in the simulated catalogues according to the previously described procedure. The curves show the systematic effect of velocity bias coming from either central or satellite galaxies, when f is estimated from models B-EXP and C-EXP assuming a scale-dependent spatial bias. We find that the impact of central galaxy velocity bias is smaller than that of satellite galaxies: at z = 0.1 a velocity bias of b cen v = 0.2 introduces a negative systematic error of about 1%, independently of the model. At z = 1 the effect is larger, with a systematic error that can reach about −3% over the range 1h −1 Mpc < r min ⊥ < 20h −1 Mpc. We should note that the test at z = 1 may be pessimistic, as the amount of velocity bias introduced in the catalogues corresponds to observational constraints from the local Universe and it is plausible that at z = 1 the actual velocity bias is lower. In the case of satellite galaxies, the introduction of a velocity bias has a stronger impact: at z = 1 the systematic error remains within 1 − 2%, while at z = 0.1 it is of about 2 − 3%, depending on the model used. In the latter case, we note that model C tend to be less affected.   Figure 13. Percent systematic variation of the measured growth rate (f vb − f novb )/f novb induced by including some amount of velocity bias from central galaxies in the simulated catalogues. Galaxies with L > L * at the two reference redshifts are considered, applying models B-EXP and C-EXP with scale-dependent bias.  Figure 14. Same as Fig. 13 but when including velocity bias coming from satellite galaxies in the simulated catalogues.
These admittedly simple tests suggest that velocity bias, if not accounted for in the models, can introduce additional systematic errors of the order of 1 − 3%, i.e. of the same order of statistical errors expected to be reachable by future large surveys. Although not dramatic, this additional source of systematic error needs to be kept in mind and possibly accounted for in future models. In principle, galaxy velocity bias can be included in the models by adding an effective velocity bias factor in front of the terms involving the velocity divergence field. At first approximation, one could assume this velocity bias factor to be constant and set it as a free parameter while fitting redshift-space distortions. Introducing more degrees of freedom in the models would however inevitably increase statistical errors.

SUMMARY AND CONCLUSIONS
Measurements of the growth rate of structure from redshift-space distortions have come to be considered as one of the most promising probes for future massive redshift surveys that aim at solving the dilemma of cosmic acceleration. A notable example is the survey planned for the approved ESA Euclid space mission (Laureijs et al. 2011). In this paper, we have investigated in some detail how well non-linear effects can be accounted for by current models when performing such measurements on galaxy catalogues. The question is how to optimally use observations and models, in order to minimise statistical and systematic errors alike. The actual extraction of the redshift-space distortions signal from real data is subject to two competing requirements. On one hand, one would like to use the simple linear description of redshift-space distortions induced by large-scale coherent motions, thus limiting the measurements and modelling to very large scales. On such scales fluctuations are close to be linear and systematic effects appear to be reduced (e.g. Samushia et al. 2011a). On the other hand, the clustering signal on those scales is so weak that statistical errors on the measured growth rate remain large, even when using samples probing large volumes. One may therefore prefer to extend the modelling to smaller scales, enlarging the range of analysed scales and thus reducing the statistical error.
The most effective compromise to exploit the size and statistics of future surveys seems thereby that of including intermediate quasi-nonlinear scales, at the expense of a more complicated modelling effort. As discussed in the introduction, most recent developments in this direction have so far concentrated on improving the description of non-linear effects in the redshift-space clustering of matter fluctuations in Fourier space. Here we have investigated how these models perform when applied to catalogues of realistic galaxies in configuration space, in both the local Universe at z = 0.1 and the more distant Universe at z = 1. For this purpose, we have reformulated in terms of the anisotropic two-point correlation function of galaxies ξ(r ⊥ , r ), the analytical model for the redshift-space anisotropic power spectrum by Scoccimarro (2004) as well as the recent improvements proposed by Taruya et al. (2010), in addition to the standard dispersion model. In these models, we have included the possibility of using a realistic non-linear scale-dependent galaxy bias, the latter being in principle measurable from the observations. At variance with the usual habit of fitting for the distortion parameter β = f /bL, we have considered the possibility of including the linear component of galaxy bias as a free parameter and directly estimate the growth rate of structure f . Our key findings and results can be summarised as follows: (i) When applied to the galaxy anisotropic correlation function, Taruya et al. (2010)'s model, the most sophisticated model considered in this analysis, generally provide the most unbiased estimates of the growth rate of structure f .
(ii) The inclusion of the scale-dependence of bias in the models is an important ingredient for minimising the systematic error on f . This is particularly significant in the case of galaxy populations with non-negligible non-linearity of the bias, as seen in the z = 1 samples of the present analysis.
(iii) Systematic errors vary with the degree of non-linearity in the bias of the considered galaxy population, which in turn is a function of redshift. Accounting for it is therefore particularly relevant when "slicing" deep surveys to measure f (z) over a significant redshift range.
(iv) Taruya et al. (2010)'s model with scale-dependent bias allows to reach the 4% accuracy on the recovery of the growth rate of structure independently of how biased galaxies are. The only exception is that of very luminous and biased galaxies such as LRGs, for which f can be over-estimated by more than 10%. This shows the need to further improving the way bias is included in redshiftspace distortions models.
(v) Galaxy velocity bias could represent, if not accounted for, an additional source of systematic error. By implementing realistic prescriptions for galaxy velocity bias coming from either central or satellite galaxies in our simulated samples, we estimate that neglecting it yields to an underestimate of the recovered f by 1 − 3%.
Overall, these results emphasize the need for careful modelling non-linear effects, if redshift-space distortions have to be used as a precision cosmology probe. They nevertheless indicate promising venues along which to develop further methods to overcome systematic biases and support the hope that redshift-space distortions in future surveys will provide us with a unique test of the cosmological model. currence relations and integral forms of spherical Bessel functions (see Toyoda & Ozaki 2010, for details).
Although the model predicts non-null multipole moments of orders l = 6 and l = 8, we neglected these terms in this analysis, including only correlation function multipole moments ξ s 0 (s), ξ s 2 (s), and ξ s 4 (s) (see main text).